Continuous fractions in CAS

[pinkman]

I wanted to improve my knowledge about continuous fractions, to be able to compete with Gerson and Albert.
I noticed that the CAS of the Prime has the two functions dfc and dfc2f, really fun to use to simplify the comprehension of continuous fractions.

Well indeed it works well, ie. dfc(π) returns [3,7,15,1,292,1,1,1,2], and dfc2f([3,7,15,1,292,1,1,1,2]) returns 833719/265381, which is a good approximation for π (3.14159265358).

Let’s continue with sqrt(2)
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]]

Here are my questions!
1- What does this last result [1,2,[2]] mean? I guess it means “repeat 2 until the end of times”, but I just discover this syntax. And why isn’t it also the answer for dfc(sqrt(2)) (I let the default epsilon value in CAS params) ?
2- How can I create such an array? I tried to edit it manually: (shift) (5[]) (1) (right) (2) (shift) (5[]) but it doesn’t open new brackets. When I use an external editor, as simple as the Notes app of the Prime, I can textually create [1,2,[2]], and with the copy-paste function I can paste it in the CAS environnement.

Thanks!

[roadrunner]

This works:

a:=[2]
[1,2,a]

-road

[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]]

sqrt(n) has repeating CF coefs is because radical can be flipped to the bottom:

√n = p + (√(n) - p) = p + q / (√(n) + p) = p + 1 / ((√(n) + p)/q)

where p = floor(√n), q = n - p²

If q = 0, n is perfect square, and we are done.
If q = 1, FP((√(n) + p)/q) = FP(√n), thus will repeat itself. (here, FP(x) meant x - floor(x))

Example:

XCas> dfc(sqrt(7),1)       → [2, (sqrt(7)+2)/3]
XCas> dfc(sqrt(7),4)       → [2,1,1,1, sqrt(7)+2]          // q=1, next one will repeat
XCas> dfc(sqrt(7),5)       → [2,1,1,1,4, (sqrt(7)+2)/3] ≡ [2, [1,1,1,4]]

Code:

def CFsqrt(n, repeat=False):    # Generate CF coefs of sqrt(n)
    p = sqrtn = int(n**0.5)
    yield sqrtn                 # 1st continued fraction coef
    q = n - p * p
    if not q: return            # perfect square, nothing to do!
    while q != 1 or repeat:     # q==1 signal recurring about to start
        coef = (sqrtn + p) // q # continued fraction coef
        yield coef
        p = coef * q - p        # update next p and q
        q = (n - p*p) // q
    yield (sqrtn + p) // q      # last repeating coef

>>> list(CFsqrt(7))         # \(\sqrt{7} = [2; \overline{1, 1, 1, 4}]\)
[2, 1, 1, 1, 4]
>>> list(CFsqrt(77))       # \(\sqrt{77} = [8; \overline{1, 3, 2, 3, 1, 16}]\)
[8, 1, 3, 2, 3, 1, 16]

[pinkman]

(09-23-2020 11:33 AM)roadrunner Wrote:  This works:

a:=[2]
[1,2,a]

-road

Yes good trick, but is also another roundabout, but not really straightforward.

[pinkman]

(09-23-2020 02:04 PM)Albert Chan Wrote:  sqrt(n) has repeating CF coefs is because radical can be flipped to the bottom:

√n = p + (√(n) - p) = p + q / (√(n) + p) = p + 1 / ((√(n) + p)/q)

where p = floor(√n), q = n - p²

If q = 0, n is perfect square, and we are done.
If q = 1, FP((√(n) + p)/q) = FP(√n), thus will repeat itself. (here, FP(x) meant x - floor(x))

Thank you for this detailed explanation, really enjoying it.

Example:

(09-23-2020 02:04 PM)Albert Chan Wrote:  XCas> dfc(sqrt(7),5)       → [2,1,1,1,4, (sqrt(7)+2)/3] ≡ [2, [1,1,1,4]]

Prime’s CAS gives:

dfc(sqrt(7),5) -> [2,1,1,1,4,1.54858377036]
No inner brackets.

(09-23-2020 02:04 PM)Albert Chan Wrote:  >>> list(CFsqrt(7))         # \(\sqrt{7} = [2; \overline{1, 1, 1, 4}]\)
[2, 1, 1, 1, 4]

This gave me the idea of inserting a list in the array, and I thought it was a success:
[1,2,{2}] Enter gives:
[1 2 [2]] in the history (hooray!)
[1 2 table()] in the result field
And it is unusable.

If someone has another clue, in addition to road’s one, I’ll take it.

[Joe Horn]

Omit the outer brackets. dfc2f(1,2,[2]) --> sqrt(2) (after simplifying).

[pinkman]

(09-24-2020 01:12 AM)Joe Horn Wrote:  Omit the outer brackets. dfc2f(1,2,[2]) --> sqrt(2) (after simplifying).

Perfect, thanks! (Even if incoherent)

Now that this is solved, I also don’t understand why dfc(sqrt(2),5) returns [1,2,[2]] and dfc(sqrt(2)) returns [1,2,2,2,2,2,2,2,2,2,2,2,2,2,2].
Why both don’t return [1,2,[2]] ?

[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

[Albert Chan]

(09-24-2020 06:19 PM)Albert Chan Wrote:  From giac source, I learned a nice trick to estimate continued fraction convergent error.
Without knowing the convergent ! (see prog.cc, float2continued_frac())

Let y[i] = [a[i]; a[i+1], a[i+2], ..., a[n] + x]. This is the effect of dropping x of y[0]:

\( \displaystyle ε = |y_0 - [a_0;a_1,a_2,\;...,\;a_n] |
≈ \left| {y_1 - [a_1;a_2,a_3,\;...,\;a_n] \over y_1^2 } \right|
≈\;...\;≈ {x \over y_1^2\;y_2^2 \;...\;y_n^2}\)

Flip RHS denominator to LHS, we have rough estimate of when to stop generating CF coefficients.

[Han]

(09-24-2020 07:14 AM)pinkman Wrote:  

(09-24-2020 01:12 AM)Joe Horn Wrote:  Omit the outer brackets. dfc2f(1,2,[2]) --> sqrt(2) (after simplifying).

Perfect, thanks! (Even if incoherent)

The CAS on the HP Prime is based off of Giac/Xcas, which -- as far as I can recall -- has always treated comma-delimited input as a list/vector even if typed without the [] or {} delimiters.