+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: HP Prime (/forum-5.html)
+--- Thread: Continuous fractions in CAS (/thread-15622.html)
Continuous fractions in CAS - pinkman - 09-23-2020 11:22 AM
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!
RE: Continuous fractions in CAS - roadrunner - 09-23-2020 11:33 AM
This works:
a:=[2]
[1,2,a]
-road
RE: Continuous fractions in CAS - Albert Chan - 09-23-2020 02:04 PM
(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]
RE: Continuous fractions in CAS - pinkman - 09-23-2020 09:28 PM
(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.
RE: Continuous fractions in CAS - pinkman - 09-23-2020 10:05 PM
(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.
RE: Continuous fractions in CAS - Joe Horn - 09-24-2020 01:12 AM
Omit the outer brackets. dfc2f(1,2,[2]) --> sqrt(2) (after simplifying).
RE: Continuous fractions in CAS - pinkman - 09-24-2020 07:14 AM
(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]] ?
RE: Continuous fractions in CAS - Albert Chan - 09-24-2020 06:19 PM
(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
RE: Continuous fractions in CAS - Albert Chan - 03-04-2021 05:06 PM
(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 !
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.
RE: Continuous fractions in CAS - Han - 03-05-2021 01:38 AM
(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.