Mv class return 1 for all negative integer power

[mammalwong]

The __ pow __ method claims supporting integer power but it is returning 1 for any negative integer power, it also didn't raise exception if it do not support negative integer power, the __ pow __ method should check if n < 0 and return self.inv() to the power abs(n)

[utensil]

@mammalwong Hi, thanks for reporting these issues! Sorry that I just noticed them.

Which version of GAlgebra are you using? Are these bugs recent regression or they existed for a long time?

P.S. I wonder what are you using GAlgebra for? Physics, computer graphics, learning GA etc.?

[mammalwong]

@utensil the version is 0.5.1. I was running the code on google colab, and it has to execute pip install galgebra whenever a new cloud machine is assigned to run my code. I have only start to use galgebra 2 days before I report this issue so I have no idea whether this issue is a recent regression or is a long existing issue.

My background is computer science and I had personal studies of traditional vector algebra based computer graphics and physics (which are mostly physics for games which are discretely approximated to just "convince" the player), I had also self studied basic geometric algebra (no geometric calculus and only GA with a diagonal metric tensor with 1/-1/0 elements, and primarily Cl(3,0,1) which is 3d projective GA for games).

For what I am using galgebra for, I have implemented the methods in section 15 and 20 of https://www.cphysics.org/article/86175.pdf to generate matrix representation of non-degenerate GA, and extended it to also generate matrix rep of degenerate GA, I am using galgebra to crosscheck and validate whether the operations on matrix presentations are correct or not, and my main focus is to see how computational effective are GA compare to traditional vector algebra (which have lots of existing hardware optimization). Besides it, geometric algebra neural layer are also one of the raising topic in the deep-learning/AI field which relate to my profession.

I wonder if you will implement matrix representation of GA? as I found it is pretty easy to calculate fractional power and exponential of a multivector using its matrix rep with matrix exponential, and it is quite easy to inverse map a matrix rep back to its multivector form.

[utensil]

Good to learn about Matrix Representations of Clifford Algebras, it would be nice to have matrix representations implemented but not a priority for me at present.

Would you be interested in discussing GA in a discord with many GA enthusiasts ? https://discord.gg/2tQZVYkD

[mammalwong]

Thanks for the discord invitation, yeah I will join it, is it the discord of bivector.net?

By the way, do you interest in how am I extending the non-degenerate GA only matrix representation generation to include also degenerate GA?

[utensil]

Yes, it is.

Yes, I'm also interested in degenerate GAs, and looking forwards to tests or new features in this direction.

[mammalwong]

I extended the non-degenerate GA matrix rep generation to include also degenerate GA like this: for Cl(p,q,r), I generate the matrix rep of Cl( p+r, q+r, 0) instead, then for each positive negative extra r basis pair a and b, the degenerate basis can be created as either $R_+ = (a+b)/2$ or $R_- = (a-b)/2$, choose only one of them and the other is discarded.

I think you should know them already, $R_+ R_-$ and $R_- R_+$ defines a projector pair and the anticommutator {R_+,R_-} = 1 as stated in this video https://youtu.be/bfuxhLnDzlQ?si=lxXQnY_KjIOMbUZJ

[utensil]

@mammalwong Can you also post minimal working example for this issue? Thanks!

[mammalwong]

here is an example to reproduce this issue:

ga = Ga('e', g=[1,1,1], coords=S.symbols(f"0:{3}", real=True), wedge=False)
ex,ey,ez = ga.mv()
ex**-1, ga.mv(1)/ex

returns:

1, e0

expected correct output:

e0, e0

[utensil]

Thanks, this can be confirmed as a bug. The current implementation is naive and didn't consider negative cases. Will address it later.

[utensil]

@mammalwong This should also be fixed in the branch now.

[mammalwong]

@utensil I have tested this issue with the branch again, I can confirm this issue is fixed.

[utensil]

Thanks, it's also added to the test cases as examples/issue-511.ipynb to prevent future regression.

[mammalwong]

@utensil, looks like the current implementation of __pow__ is having $O(n)$ time complexity when raise a multivector to power n, I suggest to adopt this implementation which takes only $O(log_2 n)$:

def pow(x, y):
  if y < 0:
    x, y = 1/x, -y

  rval = 1
  y_int = int(y)
  y_fractional = y - y_int
  # y_fractional is always between 0 and 1
  if y_fractional > 0:
    rval *= x**y_fractional

  y = y_int
  while True:
    if y & 1 != 0:
      rval *= x
    y >>= 1
    if y == 0:
      break
    x = x * x

  return rval

# test
X,Y = 2,5.5
pow(X, Y), X**Y

[utensil]

Would the implementation also support non-integer powers of multivectors? It doesn't look correct to me in that part.

As a symbolic library, my personal take is to prefer simple readable implementation rather than over-optimized ones. I would also prefer to delegate such optimizations to low-level implementation in SymPy as much as possible and keep GAlgebra minimal.

Although in this case it's simple enough, but to adopt such an implementation, much more unit tests upon more GAs and more blades are required. Your test only verified numeric x, y so far.

[mammalwong]

@utensil that implementation was copied from my previous code for other purpose, let me rewrite the code to compare with the current implementation of galgebra:

def __pow__(self, n):  # Integer power operator
        if not isinstance(n, int):
            raise ValueError('!!!!Multivector power can only be to integer power!!!!')

        x, y = self, n
        if y < 0:
            x, y = x.inv(), -y

        # compute result = x**y
        result = Mv(S.One, 'scalar', ga=self.Ga)
        while True:
            if y & 1 != 0:
                result *= x
            y >>= 1
            if y == 0:
                break
            x = x * x

        return result

This should be still simple compare to the current O(n) implementation in galgebra

[utensil]

I'll be happy to take a PR to address this, with proper test coverage, and some investigation about the possibility to delegate such details to SymPy.

GAlgebra is built as a reference symbolic implementation and my usage of it relies on its correctness heavily. So I'll need to be very sure about the correctness (i.e. proper test coverage, or delegating to maintained dependencies e.g. SymPy which also breaks sometimes).

The implementation you are suggesting is simple and it's particularly good for GA, besides it's $O(\log_2 n)$, x*x would be simplified a lot due to the quadratic form of GA.

But unfortunately I need to note that GAlgebra is in maintenance mode, as its maintainers are preoccupied from time to time, like I was in the last few years (when GAlgebra was almost idle, except for the changes being updated in #510 ). Lately I have more time to invest into it as I need it again for cross-library validations, so it will be actively maintained for a while.

Being in maintenance mode means that maintainers would accept PRs, fix bugs with their best effort, but will not be actively developing new features unless they need them themselves too.

In the case of the changes you're proposing, to ensure its correctness, I need 10x work to add the tests etc. for different GAs and blades, and it's not going to be a priority in weeks at least. Sorry.

[mammalwong]

Understandable and thanks for the explanation of the situation 👍

[utensil]

Thanks, and I've opened separate issues for them. This issue will be closed when the fixes of the original issues are merged into master.

[utensil]

Thank you, @mammalwong , the fix should be in master now as #510 is now merged.

[mammalwong]

@utensil thanks for the acknowledgement!