extra/project-euler/064/064.factor

   1 USING: accessors arrays classes.tuple io kernel locals math
   2 math.functions ranges prettyprint project-euler.common
   3 sequences ;
   4 IN: project-euler.064
   5
   6 ! http://projecteuler.net/index.php?section=problems&id=64
   7
   8 ! DESCRIPTION
   9 ! -----------
  10
  11 ! All square roots are periodic when written as continued
  12 ! fractions and can be written in the form:
  13
  14 ! √N=a0+1/(a1+1/(a2+1/a3+...))
  15
  16 ! For example, let us consider √23:
  17
  18 ! √23=4+√(23)−4=4+1/(1/(√23−4)=4+1/(1+((√23−3)/7)
  19
  20 ! If we continue we would get the following expansion:
  21
  22 ! √23=4+1/(1+1/(3+1/(1+1/(8+...))))
  23
  24 ! The process can be summarised as follows:
  25
  26 ! a0=4, 1/(√23−4) = (√23+4)/7 = 1+(√23−3)/7
  27 ! a1=1, 7/(√23−3) = 7*(√23+3)/14 = 3+(√23−3)/2
  28 ! a2=3, 2/(√23−3) = 2*(√23+3)/14 = 1+(√23−4)/7
  29 ! a3=1, 7/(√23−4) = 7*(√23+4)/7 = 8+√23−4
  30 ! a4=8, 1/(√23−4) = (√23+4)/7 = 1+(√23−3)/7
  31 ! a5=1, 7/(√23−3) = 7*(√23+3)/14 = 3+(√23−3)/2
  32 ! a6=3, 2/(√23−3) = 2*(√23+3)/14 = 1+(√23−4)/7
  33 ! a7=1, 7/(√23−4) = 7*(√23+4)/7 = 8+√23−4
  34
  35 ! It can be seen that the sequence is repeating. For
  36 ! conciseness, we use the notation √23=[4;(1,3,1,8)], to
  37 ! indicate that the block (1,3,1,8) repeats indefinitely.
  38
  39 ! The first ten continued fraction representations of
  40 ! (irrational) square roots are:
  41
  42 ! √2=[1;(2)] , period=1
  43 ! √3=[1;(1,2)], period=2
  44 ! √5=[2;(4)], period=1
  45 ! √6=[2;(2,4)], period=2
  46 ! √7=[2;(1,1,1,4)], period=4
  47 ! √8=[2;(1,4)], period=2
  48 ! √10=[3;(6)], period=1
  49 ! √11=[3;(3,6)], period=2
  50 ! √12=[3;(2,6)], period=2
  51 ! √13=[3;(1,1,1,1,6)], period=5
  52
  53 ! Exactly four continued fractions, for N <= 13, have an odd period.
  54
  55 ! How many continued fractions for N <= 10000 have an odd period?
  56
  57 <PRIVATE
  58
  59 TUPLE: cont-frac
  60     { whole integer }
  61     { num-const integer }
  62     { denom integer } ;
  63
  64 C: <cont-frac> cont-frac
  65
  66 : deep-copy ( cont-frac -- cont-frac cont-frac )
  67     dup tuple>array rest cont-frac slots>tuple ;
  68
  69 : create-cont-frac ( n -- n cont-frac )
  70     dup sqrt >fixnum dup 1 <cont-frac> ;
  71
  72 : step ( n cont-frac -- n cont-frac )
  73     swap dup
  74     ! Store n
  75     [let :> n
  76         ! Extract the constant
  77         swap dup num-const>>
  78         :> num-const
  79
  80         ! Find the new denominator
  81         num-const 2 ^ n swap -
  82         :> exp-denom
  83
  84         ! Find the fraction in lowest terms
  85         dup denom>>
  86         exp-denom simple-gcd
  87         exp-denom swap /
  88         :> new-denom
  89
  90         ! Find the new whole number
  91         num-const n sqrt + new-denom / >fixnum
  92         :> new-whole
  93
  94         ! Find the new num-const
  95         num-const new-denom /
  96         new-whole swap -
  97         new-denom *
  98         :> new-num-const
  99
 100         ! Finally, update the continuing fraction
 101         drop new-whole new-num-const new-denom <cont-frac>
 102     ] ;
 103
 104 :: loop ( c l n cf -- c l n cf )
 105     n cf step :> new-cf drop
 106     c 1 + l n new-cf
 107     l new-cf = [ loop ] unless ;
 108
 109 : find-period ( n -- period )
 110     0 swap
 111     create-cont-frac
 112     step
 113     deep-copy -rot
 114     loop
 115     drop drop drop ;
 116
 117 : try-all ( -- n )
 118     2 10000 [a..b]
 119     [ perfect-square? ] reject
 120     [ find-period ] map
 121     [ odd? ] filter
 122     length ;
 123
 124 PRIVATE>
 125
 126 : euler064a ( -- n ) try-all ;
 127
 128 <PRIVATE
 129
 130 ! (√n + a)/b
 131 TUPLE: cfrac n a b ;
 132
 133 C: <cfrac> cfrac
 134
 135 : >cfrac< ( fr -- n a b )
 136     [ n>> ] [ a>> ] [ b>> ] tri ;
 137
 138 ! (√n + a) / b = 1 / (k + (√n + a') / b')
 139 !
 140 ! b / (√n + a) = b (√n - a) / (n - a^2) = (√n - a) / ((n - a^2) / b)
 141 :: reciprocal ( fr -- fr' )
 142     fr >cfrac< :> ( n a b )
 143     n
 144     a neg
 145     n a sq - b /
 146     <cfrac> ;
 147
 148 :: split ( fr -- k fr' )
 149     fr >cfrac< :> ( n a b )
 150     n sqrt a + b /i
 151     dup n swap
 152     b * a swap -
 153     b
 154     <cfrac> ;
 155
 156 : pure ( n -- fr )
 157     0 1 <cfrac> ;
 158
 159 : next ( fr -- fr' )
 160     reciprocal split nip ;
 161
 162 :: period ( n -- period )
 163     n sqrt >integer sq n = [ 0 ] [
 164         n pure split nip :> start
 165         1 start next
 166         [ dup start = not ]
 167         [ next [ 1 + ] dip ]
 168         while drop
 169     ] if ;
 170
 171 PRIVATE>
 172
 173 : euler064b ( -- ct )
 174     10000 [1..b] [ period odd? ] count ;
 175
 176 SOLUTION: euler064b