include proof program for _arb_atan_taylor_rs

arb/atan_taylor_rs.c

@@ -25,6 +25,112 @@
 
 #include "arb.h"
 
+/*
+Python code for generating the coefficients and verifying that the
+rectangular splitting code is correct. Assumptions: N < 256, xn >= 1, x <= 1/16,
+FLINT_BITS = 32 or 64. Executing the verification code proves that:
+
+* The final 2 ulp error bound is correct
+* The fixed-point arithmetic never overflows (i.e. the value of the single
+  integral limb never exceeds 2^32-1 or 2^64-1)
+* After temporarily working with twos complemented negative values, we
+  always get nonnegative values back in the end
+
+
+from gmpy import mpq, lcm, denom, numer
+
+def coefficients(NN, bits):
+    ps = []
+    qs = []
+    temp = []
+    Q = 1
+    for k in range(2*NN+50):
+        p = 1
+        q = 2*k+1
+        if lcm(Q, q) < 2**bits:
+            temp.append(mpq(p,q))
+            Q = lcm(Q, q)
+        else:
+            for a in temp:
+                ps.append(int(a * Q))
+                qs.append(int(Q))
+            Q = q
+            temp = [mpq(p,q)]
+    return ps[:NN], qs[:NN]
+
+# verifies error bound in ulps
+def verify_rs(N, PS, QS, bits):
+    # by assumption, |x| <= xbound
+    xbound = mpq(1,16)
+    # step length
+    m = 2
+    while m * m < N:
+        m += 2
+    # bound for sum (ignoring error)
+    sbound = 0
+    # error bound
+    error = 0
+    power = (N - 1) % m
+    for k in range(N-1, -1, -1):
+        c = PS[k]
+        new_denom = QS[k]
+        old_denom = QS[k+1]
+        if new_denom != old_denom and k < N-1:
+            # if alternating, adding old_denom must give a nonnegative number
+            assert sbound + error/mpq(2**bits) <= old_denom
+            # adding old_denom must not overflow
+            assert sbound + old_denom + error/mpq(2**bits) < 2**bits - 1
+            # change denominator
+            sbound = sbound * new_denom / mpq(old_denom)
+            # new error from the change of denominator
+            error = error * (new_denom / mpq(old_denom)) + 1
+            # must not overflow
+            assert sbound + new_denom + error/mpq(2**bits) < 2**bits - 1
+        if power == 0:
+            # if alternating, adding c must give a nonnegative number
+            assert sbound + error/mpq(2**bits) <= c
+            sbound += c
+            # check no overflow
+            assert sbound + error/mpq(2**bits) < 2**bits - 1
+            if k != 0:
+                # multiply by x^m
+                xmerror = 2            # error of x^m is <= 2 ulps
+                xmbound = xbound ** m  # magnitude
+                # (s + err(s))*(x^m + err(x^m)) - s * x^m = 
+                # s*err(x^m) + err(s)*x^m + err(s)*err(x^m)
+                nerror = 0
+                nerror += error * xmbound       # ulps
+                nerror += xmerror * sbound      # ulps
+                nerror += error * xmerror / mpq(2**bits) # ulps
+                nerror += 1     # 1 ulp rounding
+                error = nerror
+                sbound = sbound * xmbound
+                assert sbound + error/mpq(2**bits) < 2**bits - 1
+            power = m - 1
+        else:
+            sbound += c * xbound**power
+            # error of x^power is <= 0, 1 or 2 ulps
+            if power == 2:
+                error += c
+            elif power > 2:
+                error += 2 * c
+            # check no overflow
+            assert sbound + error/mpq(2**bits) < 2**bits - 1
+            power -= 1
+    # final division and multiplication by x
+    error = error / mpq(QS[0]) + 1
+    assert sbound + error/mpq(2**bits) < 2**bits - 1
+    error = error * xbound + 1
+    assert sbound + error/mpq(2**bits) < 2**bits - 1
+    assert error <= 2
+
+for bits in [32, 64]:
+    PS, QS = coefficients(256, bits)
+    for N in range(256):
+        verify_rs(N, PS, QS, bits)
+
+*/
+
 int _arf_get_integer_mpn(mp_ptr y, mp_srcptr x, mp_size_t xn, long exp);
 
 #define TMP_ALLOC_LIMBS(size) TMP_ALLOC((size) * sizeof(mp_limb_t))