Add general trig functions via taylor series

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,142 @@
+"""
+.. currentmodule:: clifford.taylor_expansions
+
+=================================================
+taylor_expansions (:mod:`clifford.taylor_expansions`)
+=================================================
+
+.. versionadded:: 1.3.2
+
+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.
+
+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 e**mv where mv is a multivector
+    The parameter 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
+    """
+    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
+    """
+    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
+    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
+    """
+    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
+    """
+    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
+    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)

clifford/test/test_trig_functions.py

@@ -0,0 +1,86 @@
+
+from .. import Cl
+import numpy as np
+import pytest
+
+
+class TestScalarProperties:
+    @pytest.fixture()
+    def element(self):
+        alg, blades = Cl(0, 0, 0)
+        return alg.scalar
+
+    @pytest.mark.parametrize('np_func', [np.sin,
+                                         np.cos,
+                                         np.tan,
+                                         np.sinh,
+                                         np.cosh,
+                                         np.tanh])
+    def test_trig(self, element, np_func):
+        for x in np.linspace(0, 2*np.pi, 100):
+            assert abs(np_func(x*element).value[0] - np_func(x)) < 1E-10
+
+
+class TestDualNumberProperties:
+    @pytest.fixture()
+    def element(self):
+        alg, blades = Cl(0, 0, 1)
+        return alg.scalar, blades['e1']
+
+    @pytest.mark.parametrize('func, deriv_func', [(np.sin, np.cos),
+                                                  (np.cos, lambda x: -np.sin(x)),
+                                                  (np.tan, lambda x: (1/np.cos(x)**2)),
+                                                  (np.sinh, np.cosh),
+                                                  (np.cosh, np.sinh),
+                                                  (np.tanh, lambda x: (1 - np.tanh(x)**2))])
+    def test_derivatives(self, element, func, deriv_func):
+        for x in np.linspace(0, 2*np.pi, 10):
+            result = func(x * element[0] + element[1])
+            assert abs(result.value[0] - func(x)) < 1E-10
+            assert abs(result.value[1] - deriv_func(x)) < 1E-10
+
+
+class TestComplexNumberProperties:
+    @pytest.fixture()
+    def Cl010element(self):
+        alg, blades = Cl(0, 1, 0)
+        return alg.scalar, blades['e1']
+
+    @pytest.fixture()
+    def Cl000element(self):
+        alg, blades = Cl(0, 0, 0)
+        return alg.MultiVector(0*1j + alg.scalar.value), alg.MultiVector(1j*alg.scalar.value)
+
+    @pytest.mark.parametrize('np_func', [np.sin,
+                                         np.cos,
+                                         np.tan,
+                                         np.sinh,
+                                         np.cosh,
+                                         np.tanh])
+    def test_trig_Cl010(self, Cl010element, np_func):
+        """
+        This tests the a clifford algebra isomorphic to the complex numbers
+        """
+        for x in np.linspace(0, 2 * np.pi, 10):
+            for y in np.linspace(0, 2 * np.pi, 10):
+                complex_mv = x * Cl010element[0] + y * Cl010element[1]
+                complex_value = x + 1j * y
+                result = np_func(complex_mv)
+                assert abs(result.value[0] + 1j * result.value[1] - np_func(complex_value)) < 1E-10
+
+    @pytest.mark.parametrize('np_func', [np.sin,
+                                         np.cos,
+                                         np.tan,
+                                         np.sinh,
+                                         np.cosh,
+                                         np.tanh])
+    def test_trig_CxCl000(self, Cl000element, np_func):
+        """
+        This tests the complexified clifford algebra of only the scalars
+        """
+        for x in np.linspace(0, 2 * np.pi, 10):
+            for y in np.linspace(0, 2 * np.pi, 10):
+                complex_mv = x * Cl000element[0] + y * Cl000element[1]
+                complex_value = x + 1j * y
+                result = np_func(complex_mv)
+                assert abs(result.value[0] - np_func(complex_value)) < 1E-10

docs/api/index.rst

@@ -10,3 +10,4 @@ API
     tools
     operator
     transformations
+    taylor_expansions

docs/api/taylor_expansions.rst

@@ -0,0 +1 @@
+.. automodule:: clifford.taylor_expansions