Add cosine_taylor_dekker

Author: pearu

functional_algorithms/expr.py

@@ -624,6 +624,8 @@ def __str__(self):
     def __repr__(self):
         if self.kind in {"symbol"}:
             operands = tuple(repr(o) for o in self.operands)
+        elif self.kind in {"series", "constant"}:
+            operands = (repr(self.operands[0]),) + tuple(f"{o.kind}:{o.intkey}" for o in self.operands[1:])
         else:
             operands = tuple(f"{o.kind}:{o.intkey}" for o in self.operands)
         return f"{type(self).__name__}({self.kind!r}, {operands}, {self.props})"

functional_algorithms/floating_point_algorithms.py

@@ -305,8 +305,8 @@ def mul_dw(ctx, xh, xl, yh, yl):
 def div_dw(ctx, xh, xl, yh, yl):
     """Dekker's division"""
     q = xh / yh
-    th, tl = mul_dw(q, yh)
-    l = (xh - th - tl + xl - q * yl) / yh
+    th, tl = mul_dekker(ctx, q, yh)
+    l = ((((xh - th) - tl) + xl) - q * yl) / yh
     rh = q + l
     rl = q - rh + l
     return rh, rl
@@ -1447,26 +1447,11 @@ def sine_taylor(ctx, x, order=7, split=False):
     return p0, p1 * xxl
 
 
-def sine_taylor_dekker(ctx, x, order=7):
-    """Return sine of x using Taylor series approximation and Dekker's product.
-
-    See also sine_taylor.
-    """
-    C, f = [x], 1
-    for i in range(3, order + 1, 2):
-        f *= -i * (i - 1)
-        C.append(div_series(ctx, x, ctx.constant(f, x)))
-    xx = mul_series(ctx, x, x)
-    # Horner's scheme is most accurate
-    return fast_polynomial_dekker(
-        ctx, xx, C, reverse=False, scheme=[None, horner_scheme, estrin_dac_scheme, canonical_scheme][1]
-    )
-
-
 def cosine_taylor(ctx, x, order=6, split=False, drop_leading_term=False):
     """Return sine of x using Taylor series approximation:
 
     C(x) = 1 - x ** 2 / 2 + x ** 4 / 24 - ... + O(x ** (order + 2))
+         = 1 + x**2 * P(x**2, C=[-1/2!, 1/4!, -1/6!, ...])
 
     If split is True, return `ch, cl` such that
 
@@ -1556,6 +1541,51 @@ def cosine_taylor(ctx, x, order=6, split=False, drop_leading_term=False):
     return ctx.constant(1, x) + p0, p1 * xxl
 
 
+def sine_taylor_dekker(ctx, x, order=7):
+    """Return sine of x using Taylor series approximation and Dekker's product.
+
+    See also sine_taylor.
+    """
+    C, f = [x], 1
+    for i in range(3, order + 1, 2):
+        f *= -i * (i - 1)
+        C.append(div_series(ctx, x, ctx.constant(f, x)))
+    xx = mul_series(ctx, x, x)
+    # Horner's scheme is most accurate
+    return fast_polynomial_dekker(
+        ctx, xx, C, reverse=False, scheme=[None, horner_scheme, estrin_dac_scheme, canonical_scheme][1]
+    )
+
+
+def cosine_taylor_dekker(ctx, x, order=7, drop_leading_term=False):
+    """Return cosine of x using Taylor series approximation and Dekker's product.
+
+    See also cosine_taylor.
+    """
+    one = ctx.constant(1, x)
+    xx = mul_series(ctx, x, x)
+    if not drop_leading_term:
+        # P(x**2, Ce=[1, -1/2!, 1/4!, -1/6!, ...])
+        C, f = [one], 1
+        for i in range(2, order + 1, 2):
+            f *= -i * (i - 1)
+            C.append(div_series(ctx, one, ctx.constant(f, x)))
+            # Horner's scheme is most accurate
+        return fast_polynomial_dekker(
+            ctx, xx, C, reverse=False, scheme=[None, horner_scheme, estrin_dac_scheme, canonical_scheme][1]
+        )
+    else:
+        # x**2 * P(x**2, C=[-1/2!, 1/4!, -1/6!, ...])
+        C, f = [], 1
+        for i in range(2, order + 1 - 2, 2):
+            f *= -i * (i - 1)
+            C.append(div_series(ctx, xx, ctx.constant(f, x)))
+            # Horner's scheme is most accurate
+        return fast_polynomial_dekker(
+            ctx, xx, C, reverse=False, scheme=[None, horner_scheme, estrin_dac_scheme, canonical_scheme][1]
+        )
+
+
 def sine_pade(ctx, x, variant=None):
     # See https://math.stackexchange.com/questions/2196371/how-to-approximate-sinx-using-pad%C3%A9-approximation
 

functional_algorithms/rewrite.py

@@ -2,7 +2,7 @@
 import numpy
 from collections import defaultdict
 from . import expr as _expr
-from .utils import number_types, value_types, float_types, complex_types, boolean_types
+from .utils import number_types, value_types, float_types, complex_types, boolean_types, warn_once
 
 
 class Printer:
@@ -478,6 +478,8 @@ def _binary_op(self, expr, op):
             yvalue, ylike = y.operands
             if isinstance(xvalue, number_types) and isinstance(yvalue, number_types):
                 r = op(xvalue, yvalue)
+                if numpy.isfinite(xvalue) and numpy.isfinite(yvalue) and not numpy.isfinite(r):
+                    warn_once(f"{expr} evaluation resulted a non-finite value `{r}`")
                 return expr.context.constant(r, xlike)
 
     def add(self, expr):

functional_algorithms/tests/test_floating_point_algorithms.py

@@ -894,66 +894,117 @@ def test_sine_taylor_dekker(dtype):
 
 
 @pytest.mark.parametrize(
-    "func,fma", [("cos", "upcast"), ("cos", "mul_add"), ("cos", "native"), ("cosm1", "upcast"), ("numpy.cos", None)]
+    "func,fma",
+    [
+        ("cos", "upcast"),
+        ("cos", "mul_add"),
+        ("cos", "native"),
+        ("cos_dekker", "native"),
+        ("cos_numpy", None),
+        ("cosm1_dekker", "native"),
+        ("cosm1", "upcast"),
+        ("cosm1_sin", "upcast"),
+        ("cosm1_sin_numpy", None),
+        ("cosm1_numpy", None),
+    ],
 )
 def test_cosine_taylor(dtype, fma, func):
     import mpmath
     from collections import defaultdict
 
     t_prec = utils.get_precision(dtype)
     working_prec = {11: 50 * 4, 24: 50 * 4, 53: 74 * 16}[t_prec]
-    optimal_order = {11: 9, 24: 11, 53: 17}[t_prec]
+    optimal_order = {11: 9, 24: 13, 53: 19}[t_prec]
     size = 1000
     samples = list(utils.real_samples(size, dtype=dtype, min_value=dtype(0), max_value=dtype(numpy.pi / 4)))
     size = len(samples)
     with mpmath.mp.workprec(working_prec):
         mpctx = mpmath.mp
         for order in [optimal_order, 1, 3, 5, 7, 9, 11, 13, 17, 19][:1]:
 
-            @fa.targets.numpy.jit(
-                paths=[fpa],
-                dtype=dtype,
-                debug=(1.5 if size <= 10 else 0),
-                rewrite_parameters=dict(optimize_cast=False, fma_backend=fma),
-            )
-            def cos_func(ctx, x):
-                return fpa.cosine_taylor(ctx, x, order=order, split=False)
-
-            @fa.targets.numpy.jit(
-                paths=[fpa],
-                dtype=dtype,
-                debug=(1.5 if size <= 10 else 0),
-                rewrite_parameters=dict(optimize_cast=False, fma_backend=fma),
-            )
-            def cosm1_func(ctx, x):
-                return fpa.cosine_taylor(ctx, x, order=order, split=False, drop_leading_term=True)
+            if func.startswith("cosm1"):
+
+                def f_expected(x):
+                    return mpctx.cos(x) - 1
+
+            else:
+
+                def f_expected(x):
+                    return mpctx.cos(x)
+
+            if func in {"cos_dekker", "cosm1_dekker"}:
+
+                @fa.targets.numpy.jit(
+                    paths=[fpa],
+                    dtype=dtype,
+                    debug=(1.5 if size <= 10 else 0),
+                    parameters=dict(series_uses_dekker=True, series_uses_2sum=True),
+                )
+                def f(ctx, x):
+                    return fpa.cosine_taylor_dekker(ctx, x, order=order, drop_leading_term=func.startswith("cosm1"))
+
+            elif func == "cos":
+
+                @fa.targets.numpy.jit(
+                    paths=[fpa],
+                    dtype=dtype,
+                    debug=(1.5 if size <= 10 else 0),
+                    rewrite_parameters=dict(optimize_cast=False, fma_backend=fma),
+                )
+                def f(ctx, x):
+                    return fpa.cosine_taylor(ctx, x, order=order, split=False)
+
+            elif func == "cosm1":
+
+                @fa.targets.numpy.jit(
+                    paths=[fpa],
+                    dtype=dtype,
+                    debug=(1.5 if size <= 10 else 0),
+                    rewrite_parameters=dict(optimize_cast=False, fma_backend=fma),
+                )
+                def f(ctx, x):
+                    return fpa.cosine_taylor(ctx, x, order=order, split=False, drop_leading_term=True)
+
+            elif func == "cosm1_sin":
+
+                @fa.targets.numpy.jit(
+                    paths=[fpa],
+                    dtype=dtype,
+                    debug=(1.5 if size <= 10 else 0),
+                    rewrite_parameters=dict(optimize_cast=False, fma_backend=fma),
+                )
+                def f(ctx, x):
+                    two = ctx.constant(2, x)
+                    sn = fpa.sine_taylor(ctx, x / two, order=order, split=False)
+                    return -two * sn * sn
+
+            elif func == "cos_numpy":
+
+                f = numpy.cos
+
+            elif func == "cosm1_numpy":
+
+                def f(x):
+                    return numpy.cos(x) - dtype(1)
+
+            elif func == "cosm1_sin_numpy":
+
+                def f(x):
+                    two = dtype(2)
+                    sn = numpy.sin(x / two)
+                    return -two * sn * sn
+
+            else:
+                assert 0, func  # not impl
 
             ulp = defaultdict(int)
             for x in samples:
-                expected_cs = utils.mpf2float(dtype, mpctx.cos(utils.float2mpf(mpctx, x)))
-                if func == "numpy.cos":
-                    cs = numpy.cos(x)
-                    u = utils.diff_ulp(cs, expected_cs)
-                elif func == "cos":
-                    cs = cos_func(x)
-                    """
-                        csh, csl = fpa.cosine_taylor(ctx, x, order=order, split=True)
-                        cs2 = utils.mpf2float(dtype, utils.float2mpf(mpctx, csh) + utils.float2mpf(mpctx, csl))
-                        assert cs == cs2, (cs, cs2, expected_cs)
-                    """
-                    u = utils.diff_ulp(cs, expected_cs)
-                elif func == "cosm1":
-                    expected_csm1 = utils.mpf2float(dtype, mpctx.cos(utils.float2mpf(mpctx, x)) - 1)
-                    cs = cosm1_func(x)
-                    u = utils.diff_ulp(cs, expected_csm1, flush_subnormals=True)
-                else:
-                    assert 0, func  # unreachable
+                expected = utils.mpf2float(dtype, f_expected(utils.float2mpf(mpctx, x)))
+                cs = f(x)
+                u = utils.diff_ulp(cs, expected)
                 ulp[u] += 1
 
-                if 0:
-                    c = "." if u == 0 else ("v" if cs < expected_cs else "^")
-                    print(c, end="", flush=True)
-
+            print()
             show_ulp(ulp)