Improve the accuracy of sine_taylor_dekker

Author: pearu

functional_algorithms/expr.py

@@ -622,8 +622,11 @@ def __str__(self):
         return Printer().tostring(self)
 
     def __repr__(self):
-        operands = tuple(f"{o.kind}:{o.intkey}" for o in self.operands)
-        return f"{type(self).__name__}({self.kind}, {operands}, {self.props})"
+        if self.kind in {"symbol"}:
+            operands = tuple(repr(o) for o in self.operands)
+        else:
+            operands = tuple(f"{o.kind}:{o.intkey}" for o in self.operands)
+        return f"{type(self).__name__}({self.kind!r}, {operands}, {self.props})"
 
     def __abs__(self):
         return self.context.absolute(self)

functional_algorithms/floating_point_algorithms.py

@@ -318,18 +318,6 @@ def mul_dekker(ctx, x, y, C=None):
     return mul_dw(ctx, x, y, xh, xl, yh, yl)
 
 
-def _mul_series_scalar(ctx, x, y):
-    assert type(y) is not tuple
-    terms = tuple(x_ * y for x_ in x[1:])
-    return ctx._series(terms, dict(unit_index=x[0][0], scaling_exp=x[0][1]))
-
-
-def _mul_scalar_series(ctx, x, y):
-    assert type(x) is not tuple
-    terms = tuple(x * y_ for y_ in y[1:])
-    return ctx._series(terms, dict(unit_index=y[0][0], scaling_exp=y[0][1]))
-
-
 def _div_series_scalar(ctx, x, y):
     assert type(y) is not tuple
     terms = tuple(x_ / y for x_ in x[1:])
@@ -343,16 +331,13 @@ def _mul_series_series(ctx, x, y):
 
     terms = []
     for n in range(len(terms1) + len(terms2) - 1):
-        xy = None
         for i, x in enumerate(terms1):
             for j, y in enumerate(terms2):
                 if i + j == n:
-                    if xy is None:
-                        xy = x * y
+                    if ctx.parameters.get("series_uses_dekker"):
+                        _terms_add(ctx, terms, n, *mul_dekker(ctx, x, y))
                     else:
-                        xy += x * y
-        assert xy is not None
-        terms.append(xy)
+                        _terms_add(ctx, terms, n, x * y)
 
     return ctx._series(tuple(terms), dict(unit_index=index1 + index2, scaling_exp=sexp1))
 
@@ -364,10 +349,10 @@ def mul_series(ctx, x, y):
     if type(x) is tuple:
         if type(y) is tuple:
             return _mul_series_series(ctx, x, y)
-        return _mul_series_scalar(ctx, x, y)
+        return _mul_series_series(ctx, x, ((0, 0), y))
     elif type(y) is tuple:
-        return _mul_scalar_series(ctx, x, y)
-    return x * y
+        return _mul_series_series(ctx, ((0, 0), x), y)
+    return _mul_series_series(ctx, ((0, 0), x), ((0, 0), y))
 
 
 def div_series(ctx, x, y):
@@ -1415,8 +1400,14 @@ def sine_taylor_dekker(ctx, x, order=7):
     C, f = [x], 1
     for i in range(3, order + 1, 2):
         f *= -i * (i - 1)
-        f1 = ctx.constant(1 / f, x)
-        C.append(mul_series(ctx, x, f1))
+        # See sine_taylor:
+        if i >= 5:
+            f1 = ctx.constant(1 / f, x)
+            C.append(mul_series(ctx, x, f1))
+        else:
+            fh, fl = split_veltkamp(ctx, ctx.constant(f, x))
+            a = fl / fh
+            C.append(ctx.series(x / fh, -a * x / fh))
     xx = mul_series(ctx, x, x)
     # Horner's scheme is most accurate
     return fast_polynomial_dekker(

functional_algorithms/tests/test_floating_point_algorithms.py

@@ -825,6 +825,7 @@ def test_sine_taylor(dtype, func, fma):
                     paths=[fpa],
                     dtype=dtype,
                     debug=(1.5 if size <= 10 else 0),
+                    parameters=dict(series_uses_dekker=True, series_uses_2sum=True),
                 )
                 def sin_dekker_func(ctx, x):
                     return fpa.sine_taylor_dekker(ctx, x, order=order)
@@ -835,7 +836,8 @@ def sin_dekker_func(ctx, x):
                     paths=[fpa],
                     dtype=dtype,
                     debug=(1.5 if size <= 10 else 0),
-                    rewrite_parameters=dict(optimize_cast=False, fma_backend=fma, series_uses_2sum=True),
+                    parameters=dict(series_uses_2sum=True),
+                    rewrite_parameters=dict(optimize_cast=False, fma_backend=fma),
                 )
                 def sin_func(ctx, x):
                     return fpa.sine_taylor(ctx, x, order=order, split=False)