Add a new clifford.taylor_expansions module (#371)

This contains series expansions of trigonometric, exponential, and hyperbolic functions.

This also tweaks the implementation of the hitzer inverse to behave slightly better on degenerate algebras, and fixes a crash when computing an inverse in `Cl(0, 0, 0)`.

`cf.general_exp(mv)` has been deprecated in favor of `cf.taylor_expansions.exp(mv)`, or the more convenient spelling `mv.exp()`. The new function is jitted, unlike the old one.

Finally, this overloads pow in numba to allow it to be used in jitted code.

clifford/__init__.py

@@ -93,9 +93,9 @@
 
 from ._version import __version__  # noqa: F401
 from . import _numba_utils
-from . import _settings
 
 from ._settings import pretty, ugly, eps, print_precision  # noqa: F401
+import clifford.taylor_expansions as taylor_expansions
 
 # For backwards-compatibility. New code should import directly from `clifford.operator`
 from .operator import gp, op, ip  # noqa: F401
@@ -108,6 +108,11 @@
     NUMBA_PARALLEL = not bool(NUMBA_DISABLE_PARALLEL)
 
 
+def general_exp(x, **kwargs):
+    warnings.warn("cf.general_exp is deprecated. Use `mv.exp()` or `np.exp(mv)` on multivectors, or `cf.taylor_expansions.exp(x)` on arbitrary objects", DeprecationWarning, stacklevel=2)
+    return taylor_expansions.exp(x, **kwargs)
+
+
 def linear_operator_as_matrix(func, input_blades, output_blades):
     """
     Return a matrix that performs the operation of the provided linear
@@ -232,43 +237,6 @@ def grade_obj_func(objin_val, gradeList, threshold):
     return np.argmax(modal_value_count)
 
 
-def general_exp(x, max_order=15):
-    """
-    This implements the series expansion of e**mv where mv is a multivector
-    The parameter order is the maximum order of the taylor series to use
-    """
-
-    result = 1.0
-    if max_order == 0:
-        return result
-
-    # scale by power of 2 so that its norm is < 1
-    max_val = int(np.max(np.abs(x.value)))
-    scale = 1
-    if max_val > 1:
-        max_val <<= 1
-    while max_val:
-        max_val >>= 1
-        scale <<= 1
-
-    scaled = x * (1.0 / scale)
-
-    # taylor approximation
-    tmp = 1.0 + 0.0*x
-    for i in range(1, max_order):
-        if np.any(np.abs(tmp.value) > _settings._eps):
-            tmp = tmp*scaled * (1.0 / i)
-            result += tmp
-        else:
-            break
-
-    # undo scaling
-    while scale > 1:
-        result *= result
-        scale >>= 1
-    return result
-
-
 def grade_obj(objin, threshold=0.0000001):
     '''
     Returns the modal grade of a multivector

clifford/_layout.py

@@ -577,11 +577,11 @@ def _hitzer_inverse(self):
         """
         Performs the inversion operation as described in the paper :cite:`Hitzer_Sangwine_2017`
         """
-        tot = np.sum(self.sig != 0)
+        tot = len(self.sig)
         @_numba_utils.njit
         def hitzer_inverse(operand):
             if tot == 0:
-                numerator = operand.layout.scalar
+                numerator = 1 + 0*operand
             elif tot == 1:
                 # Equation 4.3
                 mv_invol = operand.gradeInvol()

clifford/_multivector.py

@@ -6,7 +6,7 @@
 import numpy as np
 
 import clifford as cf
-from . import general_exp
+import clifford.taylor_expansions as taylor_expansions
 from . import _settings
 from ._layout_helpers import layout_short_name
 
@@ -99,7 +99,25 @@ def _newMV(self, newValue=None, *, dtype: np.dtype = None) -> 'MultiVector':
     # binary
 
     def exp(self) -> 'MultiVector':
-        return general_exp(self)
+        return taylor_expansions.exp(self)
+
+    def cos(self) -> 'MultiVector':
+        return taylor_expansions.cos(self)
+
+    def sin(self) -> 'MultiVector':
+        return taylor_expansions.sin(self)
+
+    def tan(self) -> 'MultiVector':
+        return taylor_expansions.tan(self)
+
+    def sinh(self) -> 'MultiVector':
+        return taylor_expansions.sinh(self)
+
+    def cosh(self) -> 'MultiVector':
+        return taylor_expansions.cosh(self)
+
+    def tanh(self) -> 'MultiVector':
+        return taylor_expansions.tanh(self)
 
     def vee(self, other) -> 'MultiVector':
         r"""
@@ -301,7 +319,7 @@ def __rpow__(self, other) -> 'MultiVector':
         # else.
 
         # pow(x, y) == exp(y * log(x))
-        newMV = general_exp(math.log(other) * self)
+        newMV = taylor_expansions.exp(math.log(other) * self)
 
         return newMV
 

clifford/numba/_multivector.py

@@ -250,6 +250,22 @@ def impl(a, b):
         return impl
 
 
+@numba.extending.overload(operator.pow)
+def ga_pow(a, b):
+    if isinstance(a, MultiVectorType) and isinstance(b, types.Integer):
+        gmt_func = a.layout_type.obj.gmt_func
+        def impl(a, b):
+            if b < 0:
+                raise NotImplementedError('Negative powers are currently not implemented')
+            if b == 0:
+                return 1 + 0*a
+            op = a.value
+            for i in range(1, b):
+                op = gmt_func(op, a.value)
+            return a.layout.MultiVector(op)
+        return impl
+
+
 @numba.extending.overload(operator.truediv)
 def ga_truediv(a, b):
     if isinstance(a, MultiVectorType) and isinstance(b, types.abstract.Number):

clifford/taylor_expansions.py

@@ -0,0 +1,164 @@
+"""
+.. currentmodule:: clifford.taylor_expansions
+
+=====================================================
+taylor_expansions (:mod:`clifford.taylor_expansions`)
+=====================================================
+
+.. versionadded:: 1.4.0
+
+This file implements various Taylor expansions for useful functions of multivectors.
+For some algebra signatures there may exist closed forms of these functions which would likely be faster
+and more accurate. Nonetheless, having pre-written taylor expansions for the general case is useful.
+
+.. note::
+    Many of these functions are also exposed as :class:`~clifford.MultiVector` methods,
+    such as :meth:`clifford.MultiVector.sin`. This means that ``mv.sin()`` or even ``np.sin(mv)`` can be used
+    as a convenient interface to functions in this module, without having to import it directly.
+
+    For example::
+
+    >>> from clifford.g3 import *
+    >>> import numpy as np
+    >>> np.sin(np.pi*e12/4)
+    (0.86867^e12)
+
+Implemented functions
+----------------
+
+.. autofunction:: exp
+.. autofunction:: sin
+.. autofunction:: cos
+.. autofunction:: tan
+.. autofunction:: sinh
+.. autofunction:: cosh
+.. autofunction:: tanh
+
+"""
+import math
+import numpy as np
+
+from . import _numba_utils
+from . import _settings
+
+
+@_numba_utils.njit
+def exp(x, max_order=15):
+    """
+    This implements the series expansion of :math:`\exp x` where :math:`x` is a multivector
+    The parameter `max_order` is the maximum order of the taylor series to use
+    """
+
+    result = 1.0 + 0.0*x
+    if max_order == 0:
+        return result
+
+    # scale by power of 2 so that its norm is < 1
+    max_val = int(np.max(np.abs(x.value)))
+    scale = 1
+    if max_val > 1:
+        max_val <<= 1
+    while max_val:
+        max_val >>= 1
+        scale <<= 1
+
+    scaled = x * (1.0 / scale)
+
+    # taylor approximation
+    tmp = 1.0 + 0.0*x
+    for i in range(1, max_order):
+        if np.any(np.abs(tmp.value) > _settings._eps):
+            tmp = tmp*scaled * (1.0 / i)
+            result = result + tmp
+        else:
+            break
+
+    # undo scaling
+    while scale > 1:
+        result = result*result
+        scale >>= 1
+    return result
+
+
+@_numba_utils.njit
+def sin(X, max_order=30):
+    """
+    A taylor series expansion for sin
+    The parameter `max_order` is the maximum order of the taylor series to use
+    """
+    op = +X
+    X2 = X*X
+    X2np1 = X
+    for n in range(1, max_order):
+        X2np1 = X2np1 * X2
+        op = op + ((-1) ** (n) / math.gamma(2 * n + 2)) * X2np1
+    return op
+
+
+@_numba_utils.njit
+def cos(X, max_order=30):
+    """
+    A taylor series expansion for cos
+    The parameter `max_order` is the maximum order of the taylor series to use
+    """
+    op = 1 + 0*X
+    X2 = X * X
+    X2n = 1 + 0*X
+    for n in range(1, max_order):
+        X2n = X2n*X2
+        op = op + ((-1) ** (n) / math.gamma(2 * n + 1)) * X2n
+    return op
+
+
+def tan(X, max_order=30):
+    """
+    The tan function as the ratio of sin and cos
+    The parameter `max_order` is the maximum order of the taylor series to use
+
+    .. note::
+        It would probably be better to implement this as its own taylor series. This function
+        is not JITed as currently we do not overload the truediv operator for multivectors.
+    """
+    return sin(X, max_order) / cos(X, max_order)
+
+
+@_numba_utils.njit
+def sinh(X, max_order=30):
+    """
+    A taylor series expansion for sinh
+    The parameter `max_order` is the maximum order of the taylor series to use
+    """
+    op = +X
+    X2 = X * X
+    X2np1 = X
+    for n in range(1, max_order):
+        X2np1 = X2np1 * X2
+        op = op + (1 / math.gamma(2 * n + 2)) * X2np1
+    return op
+
+
+@_numba_utils.njit
+def cosh(X, max_order=30):
+    """
+    A taylor series expansion for cosh
+    The parameter `max_order` is the maximum order of the taylor series to use
+    """
+    op = 1 + 0 * X
+    X2 = X * X
+    X2n = 1 + 0 * X
+    for n in range(1, max_order):
+        X2n = X2n * X2
+        op = op + (1 / math.gamma(2 * n + 1)) * X2n
+    return op
+
+
+def tanh(X, max_order=30):
+    """
+    The tanh function as the ratio of sinh and cosh
+    The parameter `max_order` is the maximum order of the taylor series to use
+
+    .. note::
+        It would probably be better to implement this as its own taylor series. This function
+        is not JITed as currently we do not overload the truediv operator for multivectors.
+    """
+    return sinh(X, max_order) / cosh(X, max_order)