nvec ← {tolerance←⎕ct} ##.cfract numb ⍝ Continued fraction approx of real ⍵.
Result [nvec] is a simple numeric vector, whose continued fraction reduction is
within comparison tolerance of scalar argument [numb].
Any real number ⍵ may be respresented by a continued fraction (CF) of the form:
⍵ = n + 1
───── where n is an integer and
i + 1 where i j k ... are positive integers.
─────
j + 1
─────
k + 1
─────
...
or:
⍵ = n + 1 ÷ i + 1 ÷ j + 1 ÷ k + 1 ÷ ...
or:
⍵ = n + ÷ i + ÷ j + ÷ k + ÷ ...
or:
⍵ = n + ∘ ÷ i + ∘ ÷ j + ∘ ÷ k + ∘ ÷ ...
or:
⍵ = +∘÷/ n i j k ...
Given a (possibly infinite) CF, we can reconstitute a rational number arbitrar-
ily close to the real number it represents by reducing successively more terms.
Rational numbers have finite CFs and irrational numbers have infinite though
(after ignoring the first term) possibly regular ones. For example:
cfract 123÷345 ⍝ rational number: finite CF.
0 2 1 4 8
cfract *1 ⍝ irrational number: infinite, regular CF.
2 1 2 1 1 4 1 1 6 1 1 8 1 1 10 1 1 12 ...
cfract ○1 ⍝ irrational number: infinite, irregular CF.
3 7 15 1 292 1 1 1 2 1 3 1 14 ...
Reducing successively longer leading sequences of the CF, produces increasingly
closer approximations. Here are reductions of the first 10 terms of the CF for
e = 2.718281828459...:
+∘÷/¨ ,\ 10↑cfract *1
2 3 2.6666667 2.75 2.7142857 2.71875 2.7179487 2.7183099 2.7182796 2.7182836
Here are the (base-10) log differences for successive terms of e, Pi sqrt(2) and
Phi. The last rows in each column show that reducing all terms in the CF
produces a difference that is within comparison tolerance; in this case with a
default value of 10*¯14:
seq ← {+∘÷/¨,\cfract ⍵} ⍝ successive reductions.
err ← {10⍟|(⍵-⍺⍺ ⍵)÷⍵} ⍝ base 10 log error factor
col ← {6 2⍕⍪⍵} ⍝ formatted column.
vals ← (*1) (○1) (2*÷2) (0.5×1+5*÷2) ⍝ exp, Pi, sqrt 2, phi.
col∘(seq err)¨ vals ⍝ 10∘⍟ errors in successive approximations.
¯0.58 ¯1.35 ¯0.53 ¯0.42
¯0.98 ¯3.40 ¯1.22 ¯0.63
¯1.72 ¯4.58 ¯2.00 ¯1.14
¯1.93 ¯7.07 ¯2.76 ¯1.52
¯2.83 ¯9.74 ¯3.53 ¯1.95
¯3.76 ¯9.98 ¯4.29 ¯2.37
¯3.91 ¯10.41 ¯5.06 ¯2.79
¯4.99 ¯11.03 ¯5.82 ¯3.20
¯6.08 ¯11.56 ¯6.59 ¯3.62
¯6.19 ¯12.29 ¯7.35 ¯4.04
¯7.39 ¯14.15 ¯8.12 ¯4.46
¯8.61 ¯8.89 ¯4.88
¯8.69 ¯9.65 ¯5.29
¯9.99 ¯10.42 ¯5.71
¯11.30 ¯11.18 ¯6.13
¯11.37 ¯11.95 ¯6.55
¯12.75 ¯12.71 ¯6.97
¯14.13 ¯13.48 ¯7.38
¯14.24 ¯7.80
¯8.22
¯8.64
¯9.05
¯9.47
¯9.89
¯10.31
¯10.73
¯11.14
¯11.56
¯11.98
¯12.40
¯12.82
¯13.24
¯14.07
Technical notes:
This recurrence relation defines the CF of a rational number P÷Q. The '.' intro-
duces definitions, which are local to an exdented antecedent (see max.dws), and
% is integer quotient {⌊⍺÷⍵}.
cf p÷1 = p ⍝ whole number.
cf p÷q = n, cf q÷r ⍝ q>1
. n = p%q ⍝ integer quotient,
. r = q|p ⍝ and remainder.
which translates directly into code:
cfract←{ ⍝ Continued fraction approximation of real ⍵.
⍺←⎕CT ⋄ ⎕CT←⍺ ⍝ default comparison tolerance.
,↑{ ⍝ cf from rational ⍺÷⍵:
⍵=1:⍺ ⍝ whole number: finished.
n r←0 ⍵⊤⍺ ⍝ next term and remainder.
n,⍵ ∇ r ⍝ next term and cf of remainder.
}/⌊⍵ 1÷1∨⍵ ⍝ whole number ratio. See →rational←
}
The first line:
⍺←⎕CT ⋄ ⎕CT←⍺ ⍝ default comparison tolerance.
has a misleading symmetry. "⍺←···" is a special syntax, which supplies a default
value for a missing left argument. The line means: If no left argument is given,
default it to the current value of ⎕CT (comparison tolerance), and then make a
local system variable ⎕CT with the explicit or assumed value. In other words:
if a left argument is given, use its value for comparison tolerance; otherwise,
use the current value of ⎕CT.
With finite-precision representation of real numbers, such as IEEE 64-bit float-
ing point, [cfract] necessarily returns a finite vector (though see the "muse"
section below). Therefore, the result always represents a _rational_ number ⍺-
tolerably close to its "real" argument.
The left and right arguments of the inner function, which is applied by reduct-
ion:
,↑{ ⍝ cf from rational ⍺÷⍵:
...
}/⌊⍵ 1÷1∨⍵ ⍝ whole number ratio. See →rational←
are integers whose quotient is a rational number ⎕CT-tolerably close to ⍵. This
inner function could be made tail-recursive by passing the rational pair as
right argument, thus freeing the left argument to be an accumulator:
cfract←{ ⍝ Continued fraction approximation of real ⍵.
⍺←⎕CT ⋄ ⎕CT←⍺ ⍝ default comparison tolerance.
⍬{ ⍝ sequence accumulator.
p q←⍵ ⍝ rational argument pair.
q=1:⍺,p ⍝ whole number: finished.
n r←0 q⊤p ⍝ next term and remainder.
(⍺,n)∇ q r ⍝ extended accumulator and new ratio.
}⌊⍵ 1÷1∨⍵ ⍝ whole number ratio. See →rational←
}
but in this case, as the number of recursive calls is relatively small, the
saving does not pay for the additional local name assignment and so the function
as a whole, is marginally slower.
Here is an alternative, slower coding, which generates the CF directly, compar-
ing its value against the subject number at the generation of each term:
cfract←{ ⍝ Continued fraction approximation of real ⍵.
⍺←⎕CT ⋄ ⎕CT←⍺ ⍝ default comparison tolerance.
real←⍵ ⍝ "real" number.
⍬{ ⍝ starting with null sequence.
whole part←0 1⊤⍵ ⍝ whole and fractional part of number.
seq←⍺,whole ⍝ sequence extended with next term.
real=+∘÷/seq:seq ⍝ continued fraction approximates real: done.
seq ∇÷part ⍝ otherwise: continued accumulation of terms.
}⍵ ⍝ starting with real number.
}
Ref: http://en.wikipedia.org/wiki/Continued_fraction
(muse:
Continued fractions have many properties that make them an attractive alter-
native to regular decimals (3.1415...) and rationals (355÷113).
Bill Gosper's treatment of continued fraction arithmetic suggests that val-
ues be represented as objects, each of which, when referenced, returns the
next term in the CF sequence.
Implementing "closures" in Dyalog would give us an elegant way to do this.
For example, a CF representation of "e" (*1) is the infinite regular sequ-
ence:
1 0 1 1 2 1 1 4 1 1 6 1 1 8 1 1 10 1 1 12 ...
A closure-based approach to providing an object that returns the next ⍵
terms in this sequence might be:
cfe∆←{ ⍝ continued fraction stream for *1.
S←⍵ ⍝ initial stream.
{ ⍝ closure:
⍵<⍴S:(S↓⍨←⍵)⊢⍵↑S ⍝ draw ⍵ items from stream.
S,←1 1,2+¯1↑S ⍝ afix more items and
∇ ⍵ ⍝ draw ⍵ items from extended stream.
}
}
then:
e_nxt← cfe∆ 1 0 ⍝ closure with a state for "stream".
e_nxt 10 ⍝ first 10 terms in series.
1 0 1 1 2 1 1 4 1 1
e_nxt 10 ⍝ following 10 terms.
6 1 1 8 1 1 10 1 1 12
e_nxt 10 ⍝ ... and so on.
1 1 14 1 1 16 1 1 18 1
(*1)=+∘÷/(cfe∆ 1 0) 20 ⍝ first 20 terms sufficient for default ⎕ct.
1
Similarly, continued fractions for (irrational) square-roots have an initial
value followed by a repeated sequence. We could write a more general clos-
ure-generator for all square-root values.
cfseq∆←{ ⍝ return sequence-generating closure.
seq←⍵ ⍝ regular sequence.
S←⍺,⍵ ⍝ initial stream.
{ ⍝ closure:
⍵<⍴S:(S↓⍨←⍵)⊢⍵↑S ⍝ draw ⍵ items from stream.
S,←seq ⍝ afix more items and try again.
∇ ⍵
}
}
Then:
r2 ← 1 cfseq∆ 2 ⍝ Seqence for sqrt 2
r3 ← 1 cfseq∆ 1 2 ⍝ .. .. 3
...
r2 20 ⍝ first 20 terms of CF for sqrt 2.
1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
(1 seq∆ 1 2)20 ⍝ first 20 terms of CF for sqrt 19.
1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1
Finally, we could generalise our sequence generator to return a closure for
_any_ regular sequence. This operator takes, as right argument, the initial
terms of the sequence and as left operand, a function that generates further
terms on demand:
seq∆←{ ⍝ generate regular sequence closure.
S←,⍵ ⍝ initial sequence.
s←⍴S ⍝ minimum size.
⍺⍺{ ⍝ closure:
⍵<(⍴S)-s:(S↓⍨←⍵)⊢⍵↑S ⍝ ⍵ items from sequence.
S∘←⍺⍺ S ⍝ extended stream,
∇ ⍵ ⍝ and try again.
}
}
Then:
fibs ← {⍵,+/¯2↑⍵} seq∆ 0 1 ⍝ Fibonacci numbers.
fibs 20 ⍝ first 20 Fibonacci numbers.
0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181
cfe ← {⍵,1 1,2+¯1↑⍵} seq∆ 1 0 ⍝ continued fraction sequence for e.
cfe 20 ⍝ first 20 CF terms for e.
1 0 1 1 2 1 1 4 1 1 6 1 1 8 1 1 10 1 1 12
Rational numbers, whose continued fraction representation is finite, pose a
small problem for this technique. On exhausting the finite part of the sequ-
ence, the generator must continue with terms that do not affect the result;
in other words, terms whose reciprocal is _effectively_ zero. Happily, (⌊/⍬)
does the trick:
cfract 93÷113 ⍝ finite continued fraction for rational number.
0 1 4 1 1 1 5 1
r93_113 ← {⍵,⌊/⍬} seq∆ cfract 93÷113 ⍝ generator for CF 93÷113
r93_113 20 ⍝ first 20 terms of "infinite" CF for 93÷113
0 1 4 1 1 1 5 1 1.797693135E308 1.797693135E308 1.797693135E308
1.797693135E308 1.797693135E308 1.797693135E308 1.797693135E308
1.797693135E308 1.797693135E308 1.797693135E308 1.797693135E308
1.797693135E308
rational +∘÷/ r93_113 20 ⍝ check round-trip of 93÷113
93 113
rational +∘÷/ ({⍵,⌊/⍬} seq∆ cfract 93÷113) 20 ⍝ check round-trip.
93 113
In fact, we can abstract a CF-closure-generator for rational numbers:
cf_rat ← {{⍵,⌊/⍬}seq∆ cfract ⍺÷⍵} ⍝ CF generator for ⍺÷⍵
(11 cfrat 15) 8 ⍝ first 8 terms for 11÷15
0 1 2 1 2 1 1.797693135E308 1.797693135E308
Ref: http://mathworld.wolfram.com/ContinuedFraction.html
http://www.tweedledum.com/rwg/cfup.htm
Details of an experimental version of Dyalog, which implements closures may
be found at: http://dfns.dyalog.com/downloads/fre.pdf
)
Examples:
cfract 5÷8 ⍝ rational numbers have a finite CF.
0 1 1 1 2
cfract *1 ⍝ e has an infinite but regular CF.
2 1 2 1 1 4 1 1 6 1 1 8 1 1 10 1 1 12
+∘÷/ cfract *1 ⍝ reduction reconstitutes e.
2.718281828
+∘÷/⎕←cfract ○1 ⍝ irregular infinite CF for pi.
3 7 15 1 292 1 1 1 2 1 4
3.141592654
+∘÷/⎕←cfract root 2 ⍝ CF for sqrt(2).
1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
1.414213562
+∘÷/⎕←cfract 0.5×1+root 5 ⍝ CF for golden mean (phi ⌽)
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
1.618033989
cfract 3○1 ⍝ regular CF for tan(1)
1 1 1 3 1 5 1 7 1 9 1 11 1 13 1 15
{⎕←⍵,':',cfract root ⍵}¨⍳20 ⍝ (regular) CFs for first 20 sqrts.
1 : 1
2 : 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
3 : 1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2
4 : 2
5 : 2 4 4 4 4 4 4 4 4 4 4 4
6 : 2 2 4 2 4 2 4 2 4 2 4 2 4 2 4
7 : 2 1 1 1 4 1 1 1 4 1 1 1 4 1 1 1 4 1 1 1 4 1 2
8 : 2 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4
9 : 3
10 : 3 6 6 6 6 6 6 6 6 6
11 : 3 3 6 3 6 3 6 3 6 3 6 3
12 : 3 2 6 2 6 2 6 2 6 2 6 2 6
13 : 3 1 1 1 1 6 1 1 1 1 6 1 1 1 1 6 1 1 1 1 6 1 2
14 : 3 1 2 1 6 1 2 1 6 1 2 1 6 1 2 1 6 1 3
15 : 3 1 6 1 6 1 6 1 6 1 6 1 6 1 7
16 : 4
17 : 4 8 8 8 8 8 8 8
18 : 4 4 8 4 8 4 8 4 8 4
19 : 4 2 1 3 1 2 8 2 1 3 1 2 8 2 1 3 1 2
20 : 4 2 8 2 8 2 8 2 8 2 8 2
{⎕←⍵,':',cfract 3 root ⍵}¨⍳20 ⍝ irregular CFs for first 20 cube roots.
1 : 1
2 : 1 3 1 5 1 1 4 1 1 8 1 14 1 10 2 1 4
3 : 1 2 3 1 4 1 5 1 1 6 2 5 8 3 3
4 : 1 1 1 2 2 1 3 2 3 1 3 1 30 1 4 1 2 9
5 : 1 1 2 2 4 3 3 1 5 1 1 4 10 18
6 : 1 1 4 2 7 3 508 1 5 6
7 : 1 1 10 2 16 2 1 4 2 1 21 1 3 5
8 : 2
9 : 2 12 2 18 1 1 1 1 4 1 1 24 1 9
10 : 2 6 2 9 1 1 2 4 1 12 1 1 1 1 57
11 : 2 4 2 6 1 1 2 1 2 9 88 2 1 2
12 : 2 3 2 5 15 7 3 1 1 3 1 1 96
13 : 2 2 1 5 1 1 43 3 2 1 1 3 10 7
14 : 2 2 2 3 1 1 5 5 9 6 21 2
15 : 2 2 6 1 8 1 10 8 12 1 721
16 : 2 1 1 12 10 18 1 6 1 21 1 2 2
17 : 2 1 1 3 138 1 1 3 2 3 1 1 207
18 : 2 1 1 1 1 1 3 22 1 2 2 2 24 64
19 : 2 1 2 63 1 2 2 2 1 95 2 1 1 2
20 : 2 1 2 1 1 154 6 1 1 1 6 232
rational +∘÷/¨ ,\ cfract ○1 ⍝ successive rational approximations to Pi.
3 22 333 355 103993 104348 208341 312689 833719 1146408 5419351
1 7 106 113 33102 33215 66317 99532 265381 364913 1725033
See also: rational factors sieve gcd root fibonacci
Back to: contents
Back to: Workspaces