extra/rosetta-code/continued-fraction/continued-fraction.factor

   1 ! Copyright (c) 2012 Anonymous
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: arrays combinators io kernel locals math math.functions
   4 math.ranges prettyprint sequences ;
   5 IN: rosetta-code.continued-fraction
   6
   7 ! http://rosettacode.org/wiki/Continued_fraction
   8
   9 ! A number may be represented as a continued fraction (see
  10 ! Mathworld for more information) as follows:
  11
  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:
  16
  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.
  20
  21 ! Every continued fraction must implement these two words.
  22 GENERIC: cfrac-a ( n cfrac -- a )
  23 GENERIC: cfrac-b ( n cfrac -- b )
  24
  25 ! square root of 2
  26 SINGLETON: sqrt2
  27 M: sqrt2 cfrac-a
  28     ! If n is 1, then a_n is 1, else a_n is 2.
  29     drop { { 1 [ 1 ] } [ drop 2 ] } case ;
  30 M: sqrt2 cfrac-b
  31     ! Always b_n is 1.
  32     2drop 1 ;
  33
  34 ! Napier's constant
  35 SINGLETON: napier
  36 M: napier cfrac-a
  37     ! If n is 1, then a_n is 2, else a_n is n - 1.
  38     drop { { 1 [ 2 ] } [ 1 - ] } case ;
  39 M: napier cfrac-b
  40     ! If n is 1, then b_n is 1, else b_n is n - 1.
  41     drop { { 1 [ 1 ] } [ 1 - ] } case ;
  42
  43 SINGLETON: pi
  44 M: pi cfrac-a
  45     ! If n is 1, then a_n is 3, else a_n is 6.
  46     drop { { 1 [ 3 ] } [ drop 6 ] } case ;
  47 M: pi cfrac-b
  48     ! Always b_n is (n * 2 - 1)^2.
  49     drop 2 * 1 - 2 ^ ;
  50
  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
  56     ] each ;
  57
  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
  62     :> fraction
  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 ;
  68
  69 <PRIVATE
  70 : main ( -- )
  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
  75     "               Pi: " write
  76     pi 950 cfrac-estimate 10 decimalize print ;
  77 PRIVATE>
  78
  79 MAIN: main