nvec ← {tolerance} ##.cfract numb ⍝ Continued fraction approximation of real ⍵.
The result of [cfract] is an integer vector, whose continued fraction reduction
is within comparison tolerance of the original argument. For example, the sequ-
ence: 3 7 15 1 292 1 1 1 2 1 3 1, represents a continued fraction approxim-
ation to pi in that:
pi ~= 3 + 1
─────
7 + 1
──────
15 + 1
─────
1 + 1
───────
292 + 1
─────
1 + 1
─────
1 + 1
─────
1 + 1
─────
2 + 1
─────
1 + 1
─────
3 + 1
─
1 + ...
Argument [numb] is easily reconstituted from the resulting vector [nvec]:
cfract 1.234 ⍝ continued fraction representation of 1.234
1 4 3 1 1 1 10
+∘÷/ cfract 1.234 ⍝ round-trip.
1.234
Technical notes:
The coding of the function is as follows:
cfract←{ ⍝ Continued fraction approximation of real ⍵.
⍺←⎕CT ⋄ ⎕CT←⍺ ⍝ default comparison tolerance.
There is a misleading symmetry about this line. ⍺←··· is a special syntax that
supplies a default value for a missing left argument. This line means: If no
left argument is given, default it to the current value of ⎕CT (comparison tol-
erance), 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.
real←⍵ ⍝ real number.
As the subject value does not change during the computation, it is named here
for reference from within the inner function. An alternative would be to bind
the value as an operand to an inner _operator_, where it would be referenced as
⍺⍺. This ploy would be hampered slightly by having to separate the operand value
from the inner function's left argument as in Dyalog, array-array binding is
stronger than that of operand-operator. As the inner function occupies several
lines, the opening and closing parentheses around (operand operator) would wind
up on separate lines, an ugliness that currently (Sept 2001), offends the funct-
ion editor's syntax-colouring mechanism.
⍬{ ⍝ starting with null sequence,
The left argument of the tail-recursive inner function is used as an accumulator
for the resulting vector of continued fraction terms.
whole part←0 1⊤⍵ ⍝ whole and fractional part of number.
separates the whole and fractional part of ⍵. For example 123.456 → 123, 0.456.
seq←⍺,whole ⍝ sequence extended with next term.
real=+∘÷/seq:seq ⍝ continued fraction approximates real: done.
An irrational number, such as (2*÷2), will have an infinite continued fraction,
the reduced value (+∘÷/seq) of a sufficient number of terms of which, will even-
tually lie within comparison tolerance of the original.
Finally, in the recursive call:
seq ∇÷part ⍝ otherwise: continued accumulation of terms.
the extended sequence is passed in a recursive tail call. Notice that the recip-
rocal of the fractional part must have a value greater than 1.
}⍵ ⍝ 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://www.dyalog.com/dfnsdws/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 3 1
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 1 1
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 1 1
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 1 1
14 : 3 1 2 1 6 1 2 1 6 1 2 1 6 1 2 1 6 1 2 1
15 : 3 1 6 1 6 1 6 1 6 1 6 1 6 1 6 1
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 17 1
6 : 1 1 4 2 7 3 508 1 5 5 1
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 1 1
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 4272943 5419351
1 7 106 113 33102 33215 66317 99532 265381 364913 1360120 1725033
⍝ CFs, rational numbers, fibonacci sequences
⍝ (and, oh yeah, Euclid's GCD algorithm)
⍝ seem to be deeply intertwingled:
fib←{⍬⍴ rational +∘÷/ 0,⍵/1} ⍝ ⍵th Fibonacci number.
fib¨ 0 to 10
0 1 1 2 3 5 8 13 21 34 55
See also: rational factors sieve gcd root fibonacci
Back to: contents
Back to: Dyalog APL
Trouble seeing APL font?