Continuous fractions in CAS
09-23-2020, 11:22 AM
Post: #1
 pinkman Senior Member Posts: 432 Joined: Mar 2018
Continuous fractions in CAS
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!

Thibault - not collector but in love with the few HP models I own - Also musician : http://walruspark.co
09-23-2020, 11:33 AM
Post: #2
 roadrunner Senior Member Posts: 446 Joined: Jun 2015
RE: Continuous fractions in CAS
This works:

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

-road
09-23-2020, 02:04 PM (This post was last modified: 09-23-2020 02:30 PM by Albert Chan.)
Post: #3
 Albert Chan Senior Member Posts: 2,650 Joined: Jul 2018
RE: Continuous fractions in CAS
(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]
09-23-2020, 09:28 PM
Post: #4
 pinkman Senior Member Posts: 432 Joined: Mar 2018
RE: Continuous fractions in CAS
(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.

Thibault - not collector but in love with the few HP models I own - Also musician : http://walruspark.co
09-23-2020, 10:05 PM
Post: #5
 pinkman Senior Member Posts: 432 Joined: Mar 2018
RE: Continuous fractions in CAS
(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.

Thibault - not collector but in love with the few HP models I own - Also musician : http://walruspark.co
09-24-2020, 01:12 AM
Post: #6
 Joe Horn Senior Member Posts: 2,002 Joined: Dec 2013
RE: Continuous fractions in CAS
Omit the outer brackets. dfc2f(1,2,[2]) --> sqrt(2) (after simplifying).

<0|ɸ|0>
-Joe-
09-24-2020, 07:14 AM
Post: #7
 pinkman Senior Member Posts: 432 Joined: Mar 2018
RE: Continuous fractions in CAS
(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]] ?

Thibault - not collector but in love with the few HP models I own - Also musician : http://walruspark.co
09-24-2020, 06:19 PM
Post: #8
 Albert Chan Senior Member Posts: 2,650 Joined: Jul 2018
RE: Continuous fractions in CAS
(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
03-04-2021, 05:06 PM (This post was last modified: 03-04-2021 05:13 PM by Albert Chan.)
Post: #9
 Albert Chan Senior Member Posts: 2,650 Joined: Jul 2018
RE: Continuous fractions in CAS
(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.
03-05-2021, 01:38 AM (This post was last modified: 03-05-2021 01:39 AM by Han.)
Post: #10
 Han Senior Member Posts: 1,882 Joined: Dec 2013
RE: Continuous fractions in CAS
(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.

Graph 3D | QPI | SolveSys
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)