20200605.float-approximations-to-pi.sage

## Search for floating point numbers which are close to a multiple of pi.

## Double precision
mantissa_size = 52
exponent_size = 11
## Single precision
# mantissa_size = 23
# exponent_size = 8
## 80-bit extended precision
# mantissa_size = 63
# exponent_size = 15

mantissa_bound = 2^(mantissa_size+1)

best_i = None
best_k = None
best_n = None
best_eps = 1

def testappr(i, k, n):
    # Check if n*2^i is closer to k*pi than the best we've had so far.
    global best_i, best_k, best_n, best_eps
    eps = abs(N(k * pi - n * 2^i, 20000))
    if eps < best_eps:
        best_i = i
        best_k = k
        best_n = n
        best_eps = eps
        print "best so far: ", i, n * 2^i, k, N(eps,100)

def transform_continued_fraction(a,b,c,d, infrac):
    # Compute continued fraction expansion of (a+b*x)/(c+d*x) in function of
    # that of x (infrac).
    intab = list(infrac)
    outtab = []
    while True:
        # print "transforming by",a,b,c,d
        if c != 0 and d != 0 and floor(a/c) == floor(b/d):
            q = floor(a/c)
            outtab.append(q)
            (a,b,c,d) = (c, d, a-c*q, b-d*q)
        else:
            if len(intab)==0:
                if d!=0:
                    outtab.extend(list(continued_fraction(b/d)))
                return outtab
            p = intab.pop(0)
            (a,b,c,d) = (b, a+b*p, d, c+d*p)

def testexp(i, f):
    # Try to find n and k so that n*2^i is close to k*pi, by trying to get n/k
    # close to pi/2^i.  Here f is the list of coefficients of the continued
    # fraction expansion of pi/2^i.
    # First, compute continued fraction approximations to pi/2^i.
    j = 0
    (prev_k, prev_n) = (1, 0)
    (cur_k, cur_n) = (0, 1)
    while True:
        # Compute convergents until numerator exceeds mantissa_bound
        next_k = cur_k * f[j] + prev_k
        next_n = cur_n * f[j] + prev_n
        if next_n >= mantissa_bound:
            break
        j += 1
        (prev_k, prev_n) = (cur_k, cur_n)
        (cur_k, cur_n) = (next_k, next_n)
    if j<=1:
        return
    testappr(i, cur_k, cur_n)
    # Now try best rational approximations which aren't convergents
    for m in range(floor(f[j]/2), f[j]):
        test_k = cur_k * m + prev_k
        test_n = cur_n * m + prev_n
        if test_n >= mantissa_bound:
            break
        testappr(i, test_k, test_n)

# Now try all exponents sucessively
lst = list(continued_fraction(N(pi*2^mantissa_size, 2^(exponent_size+1) + 2*mantissa_size + 100)))
for ei in range(2^(exponent_size-1) + 1):
    if ei%100==0:
        print "ei=", ei
    testexp(ei - mantissa_size, lst)
    lst = transform_continued_fraction(0,1,2,0, lst)