feat: Implement MultiVector.log() principal branch inverse of exp()

[sunkmechie]

Closes #131

Summary

This PR adds MultiVector.log() as the principal inverse of exp() for normalized simple rotors.

It covers the three standard cases for a rotor r = c + B:

• circular case when B^2 < 0
• hyperbolic case when B^2 > 0
• null/translation case when B^2 = 0

What changed

• added MultiVector.log() in kingdon/multivector.py
• added tests for:
  • PGA translations
  • PGA rotations
  • hyperbolic boosts
  • symbolic inputs
  • array-valued inputs
  • principal-branch behavior
  • invalid inputs such as scaled rotors and -1
• added docs for log() usage

Notes

log() is implemented for normalized simple rotors over real scalar/array inputs and symbolic expressions.

Complex-valued inputs are intentionally not supported in this PR. The current implementation relies on the real-valued split between B^2 < 0, B^2 > 0, and B^2 = 0, together with real atan2 / atanh behavior. A correct complex implementation would be difficult and is best done in a separate PR.

[tBuLi]

Thank you for this great PR, this is definitely something we want to have.

Overall I think this is a great PR and already looks very promising. I do however have some comments:

• I agree that we only want to support simple rotors for now, but I do not think the normalized condition should be imposed. Since arctan functions are ratios anyway, the overall scale does not matter.
• The idea behind the sqrt, cosh, and sinhc functions for exp was that once these have been provided, the logic of taking the exp can be written purely in terms of these functions, without reference being made to math/sympy/numpy. For the log it seems to me that the minimal set of functions is sqrt and arctanh2 (where arctanh2 is a function $f(y, x)$ that takes two scalars as input and produces the angle), such that the log of a rotor $R = c(B) + s(B)$ is always

$$ \text{log}(R) = \text{arctanh2}(\sqrt{s^2(B)}, c(B)) \frac{s(B)}{\sqrt{s^2(B)}} $$

For a rotation you would then set arctanh2 = arctan2 and sqrt = lambda x: (-x) ** 0.5, etc. The reason I went for hyperbolic functions as the primary choice btw is because I think of the exp as being the exponential of a bivector $B$, and so the series expansion naturally suggests hyperbolic functions:

$$ e^B = \sum_{n=0}^\infty \frac{B^n}{n!} = \cosh(\sqrt{B^2}) + B \frac{\sinh(\sqrt{B^2})}{\sqrt{B^2}} $$

• Add a test in 3DPGA or possibly 3DCGA as well to make sure the filter steps work correctly, because only from a 4D algebra onwards can we have non-simple bivectors, so those checks are not being tested properly. We should also add a test to see that a non-simple bivector does raise the right exception, and likewise for the other exceptions.
• I love the fact that you also thought about the array case, which is indeed a bit of an outlier which is why I didn't do it yet for the exp.
  • I would check for this scenario by asking if len(self) > 0, since this is non-zero when the coefficients are arrays.
  • If we are in this scenario, make the sqrt and arctanh2 functions accept the masks as input and deal with it accordingly.
  • Replace e.g. np.any(mask_neg) by mask_neg.any() so that it works on other array types like torch tensors as well.
  • Instead of explicitly using np.zeros_like, make zeros_like a keyword-only argument to the function so that users can dependency inject their own preferred zeros. Alternativelly, I suppose you could also do res = 0 * c_val so we do not even need to know how to make a zero array and don't have to worry about dtype etc.

Looking forward to your reply.

[sunkmechie]

Thank you for the detailed feedback and the kind words.

I will refactor the PR to:

• Remove the normalization check and make it general.

• Implement the unified formula replacing the if/else statements.

• Update the array handling using len(self)>0 and .any() for better Torch/Tensor support.

• Add the 3DCGA/3DPGA tests for non-simple bivectors and exception handling.

I love the res = 0 * c_val suggestion for the zero check, makes it GPU ready.

For the non-simple test case, we can perhaps use two skew lines?

I'll push these updates shortly!

[sunkmechie]

Thanks again for the detailed review, Martin! I’ve pushed an update addressing the points you raised.

Key Changes:

• MultiVector.log() now accepts non-normalized simple rotors and correctly recovers the bivector generator.
• Switched to the minimal helper set of sqrt and arctanh2.
• len(self) > 0 for array detection, .any() for masks, and the 0 * value trick for zero-initialization
• Added a positive 4D round-trip case and an explicit non-simple bivector rejection.

I used $e_{03} + e_{12}$. for the non-simple case as shown below:

>>> from kingdon import Algebra
>>> pga3d = Algebra.fromname('3DPGA')
>>> B = pga3d.multivector(e03=1, e12=1)
>>> 
>>> print((B * B).filter())      # scalar + grade-4 part
-1 + 2 e0123
>>> print((B ^ B).filter())      # nonzero grade-4 part
2 e0123
>>> (1 + B).log()                # Raises NotImplementedError as expected
Traceback (most recent call last):
  ...
NotImplementedError: Currently only rotors with a simple bivector part can be logarithmized.

I’ve also added tests for the other exception cases and reran the full test suite.

To keep log() itself lean, I had to introduce new helper functions.

Looking forward to your thoughts!

[tBuLi]

Thanks for your changes! Just some first comments, I'll do a more detailed review and reply later.

[tBuLi]

I would put scale outside as a single multiply to prevent mistakes. R_arr = scale * alg.multivector(...)

[sunkmechie]

Thanks for the feedback, I've refactored to have the scale as a common multiplier for better clarity.

[tBuLi]

This has a non-simple bivector part, but it is not a valid rotor, as that should be the exp of bivector and this is not. You need to have a test input of the form $e^{b_1}e^{b_2}$ where $b_1$ and $b_2$ are simple bivectors.

[sunkmechie]

Great point on the validity of the rotor manifold. I've updated the test case to use a valid screw motion constructed from the product of two skew exponentials ($e^{b_1} e^{b_2}$).

All the tests are passing locally, looking forward to your review.

[tBuLi]

Also check for a translation

[sunkmechie]

Good point. I'll include in in the tests in next update.

[tBuLi]

I actually think complex values should be supported, but I am not sure how to interpret this test. Would this error also be raised if you take the log of a more interesting complex simple rotor? In our paper we have an example in R2,2 where complex numbers are needed in the invariant decomposition of certain rotors, so this is something we absolutely need to support.
https://www.researchgate.net/publication/370750268_Graded_Symmetry_Groups_Plane_and_Simple

[sunkmechie]

Thanks for pointing this out and for linking the paper.

The current implementation assumes real coefficients (might actually have been a regression, seeing how exp() supported them). I'll revisit this restriction while refactoring and check how it should behave for more general complex simple rotors.

[tBuLi]

I do not like all these new staticmethods on MultiVector that do not have anything to do with multivector per se. I would propose that you make a new file, log.py and place everything there, and then just add a simple call to the log entry point on MultiVector itself. That way we keep MultiVector as clean as possible.

[sunkmechie]

That makes sense. I introduced the helpers to keep the main log() implementation close in style to exp(), but I agree it ended up cluttering MultiVector.

I'll move the implementation into a separate log.py module and keep MultiVector.log() as a thin entry point.

[tBuLi]

If we need a helper like this to call a function then something is going very wrong :P. Your first PR looked man-made but now I am beginning to doubt it, tell Claude to take a chill-pill :P.

[sunkmechie]

Fair point. I will remove the helper and simplify in the next update.

[tBuLi]

See my other comment about complex.

[sunkmechie]

I started restructuring the implementation as suggested. For now I moved the structure into log.py and left placeholders while I refactor the logic. I'll fill in the implementation in the next commits.

[sunkmechie]

Thank you for your comprehensive feedback and the time you spent reviewing. Following your detailed review @tBuLi , I refactored the log() func.

Breakdown of the changes:

1. Complex Rotors & $R_{2,2}$ Invariant Decompositions

• I removed the restrictive checks. get_default_log_helpers() now screens inputs through _is_complex_value() and passes complex pathways to cmath and numpy.sqrt/arctanh.
• The $R_{2,2}$ Test: I added an integration test targeting $R_{2,2}$ using a complex simple bivector $B = (1+1j)e_{12} + (2-0.5j)e_{13}$. The principal log cleanly recovers the original generator from exp(B) down to $10^{-15}$ .

2. Arrays

• Arrays are detected using array_valued = len(rotor) > 0 and filtered via .any() (no np.any()).
• Iimplemented your coefficient = 0 * sample trick.
• Multi-Branch Matrix Test: We explicitly test a unified numpy batch mapping simultaneously to a Circular branch ($B^2 < 0$), a Hyperbolic branch ($B^2 > 0$), and a Null/Translation branch ($B^2 = 0$).

3. Misc

• Used sympy.Ne to make sure both SymPy rotations and translations work (using a direct x != 0 check causes Python to drop the translation logic entirely.)

• log() uses eps = 1e-12 tolerance (unlike exp() which uses exact equality (x == 0), dodging possible $0/0$ explosions of ($\frac{angle}{root}$).

Question

Regarding the's Riesz example (Example 6.4) from your paper:

Right now, Kingdon's exp() throws an error if we try to pass it a non-simple bivector. Should the next PR try to implement exp() and log() for those non-simple cases too?

Let me know if there's anything else you'd like me to modify.

[tBuLi]

Hi @sunkmechie, thanks for the updated PR! I will review it later, but I just wanted to answer your question already: non-simple bivectors/rotors are absolutely on the roadmap for future PR's.