extra/project-euler/064/064.factor

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