Continuous fractions in CAS

[Albert Chan]

(09-23-2020 11:22 AM)pinkman Wrote:  dfc(sqrt(2)) returns [1,2,2,2,2,2,2,2,2,2,2,2,2,2,2] (good!) and dfc(sqrt(2),5) returns... [1,2,[2]]
... why isn’t it also the answer for dfc(sqrt(2))

Timings suggested dfc(x) evaluate CF using floating point:

XCas> time(dfc(pi))                   → 2.18e-05
XCas> time(dfc(approx(pi)))       → 2.18e-05

This may generate catastrophic cancellation errors. Thus, if ε is too small, it does nothing.

XCas> dfc(pi, 1e-12)       → [pi]

It is too bad that XCas use the same name for dfc(pi, n), it is really different function.
For dfc(pi, n), last number is the remainder, not necessarily integer.

XCas> dfc(pi, 1e-11)       → \([3,7,15,1,292,1,1,1,2]\)
XCas> dfc(pi, 8)             → \([3,7,15,1,292,1,1,1,\large\frac{-66317\cdot \pi +208341}{99532\cdot \pi -312689}]\)
XCas> simplify(dfc2f(ans()))       → \(\large\pi\)

dfc2f(dfc(x,n)) = x, does not matter what n is

---

From giac source, I learned a nice trick to estimate continued fraction convergent error.
Without knowing the convergent ! (see prog.cc, float2continued_frac())

Code:

def cf(x, eps=1e-11, max_coefs=100):
    print 'coef\t\tfrac\t\teps'
    for i in range(max_coefs):
        c = int(x)
        x = x-c
        print 'c(%d) = %d\t%8g\t%g' % (i, c, x, eps)
        if x < eps: return
        x = 1/x
        eps *= x*x

>>> import math
>>> cf(math.pi)
coef           frac           eps
c(0) = 3       0.141593       1e-11
c(1) = 7       0.0625133      4.98791e-10
c(2) = 15      0.996594       1.27636e-07
c(3) = 1       0.00341723     1.2851e-07
c(4) = 292     0.634591       0.0110049
c(5) = 1       0.575818       0.0273275
c(6) = 1       0.736659       0.0824194
c(7) = 1       0.357481       0.151879
c(8) = 2       0.797351       1.18848

XCas> float(pi - dfc2f([3,7,15,1,292,1,1,1,2]))       → 8.71525074331e-12

This estimate trick is not fool-proof, but is amazingly good.
Another example where it slightly under-estimated required CF coefficients.

XCas> lst := dfc(e, 1e-11)      → [2,1,2,1,1,4,1,1,6,1,1,8,1,1,10,1]
XCas> float(e - dfc2f(lst))       → -1.15378817611e-11
XCas> lst := append(lst, 1)
XCas> float(e - dfc2f(lst))       → 4.82280881897e-13