1 ! Copyright (c) 2012 Anonymous
2 ! See https://factorcode.org/license.txt for BSD license.
3 USING: arrays combinators io kernel math math.functions
4 ranges prettyprint sequences ;
5 IN: rosetta-code.continued-fraction
7 ! https://rosettacode.org/wiki/Continued_fraction
9 ! A number may be represented as a continued fraction (see
10 ! Mathworld for more information) as follows:
12 ! The task is to write a program which generates such a number
13 ! and prints a real representation of it. The code should be
14 ! tested by calculating and printing the square root of 2,
15 ! Napier's Constant, and Pi, using the following coefficients:
17 ! For the square root of 2, use a0 = 1 then aN = 2. bN is always 1.
18 ! For Napier's Constant, use a0 = 2, then aN = N. b1 = 1 then bN = N − 1.
19 ! For Pi, use a0 = 3 then aN = 6. bN = (2N − 1)2.
21 ! Every continued fraction must implement these two words.
22 GENERIC: cfrac-a ( n cfrac -- a )
23 GENERIC: cfrac-b ( n cfrac -- b )
28 ! If n is 1, then a_n is 1, else a_n is 2.
29 drop { { 1 [ 1 ] } [ drop 2 ] } case ;
37 ! If n is 1, then a_n is 2, else a_n is n - 1.
38 drop { { 1 [ 2 ] } [ 1 - ] } case ;
40 ! If n is 1, then b_n is 1, else b_n is n - 1.
41 drop { { 1 [ 1 ] } [ 1 - ] } case ;
45 ! If n is 1, then a_n is 3, else a_n is 6.
46 drop { { 1 [ 3 ] } [ drop 6 ] } case ;
48 ! Always b_n is (n * 2 - 1)^2.
51 :: cfrac-estimate ( cfrac terms -- number )
52 terms cfrac cfrac-a ! top = last a_n
53 terms 1 - 1 [a..b] [ :> n
54 n cfrac cfrac-b swap / ! top = b_n / top
55 n cfrac cfrac-a + ! top = top + a_n
58 :: decimalize ( rational prec -- string )
59 rational 1 /mod ! split whole, fractional parts
60 prec 10^ * ! multiply fraction by 10 ^ prec
61 [ >integer unparse ] bi@ ! convert digits to strings
63 "." ! push decimal point
64 prec fraction length -
65 dup 0 < [ drop 0 ] when
66 "0" <repetition> concat ! push padding zeros
67 fraction 4array concat ;
71 " Square root of 2: " write
72 sqrt2 50 cfrac-estimate 30 decimalize print
73 "Napier's constant: " write
74 napier 50 cfrac-estimate 30 decimalize print
76 pi 950 cfrac-estimate 10 decimalize print ;