project-euler.064: adding description and SOLUTION:.

extra/project-euler/064/064.factor

@@ -1,7 +1,59 @@
-USING: accessors arrays classes.tuple io kernel locals math math.functions
+USING: accessors arrays classes.tuple io kernel locals math
-    math.ranges prettyprint project-euler.common sequences ;
+math.functions math.ranges prettyprint project-euler.common
+sequences ;
 IN: project-euler.064
 
+! http://projecteuler.net/index.php?section=problems&id=64
+
+! DESCRIPTION
+! -----------
+
+! All square roots are periodic when written as continued
+! fractions and can be written in the form:
+
+! √N=a0+1/(a1+1/(a2+1/a3+...))
+
+! For example, let us consider √23:
+
+! √23=4+√(23)−4=4+1/(1/(√23−4)=4+1/(1+((√23−3)/7)
+
+! If we continue we would get the following expansion:
+
+! √23=4+1/(1+1/(3+1/(1+1/(8+...))))
+
+! The process can be summarised as follows:
+
+! a0=4, 1/(√23−4) = (√23+4)/7 = 1+(√23−3)/7
+! a1=1, 7/(√23−3) = 7*(√23+3)/14 = 3+(√23−3)/2
+! a2=3, 2/(√23−3) = 2*(√23+3)/14 = 1+(√23−4)/7
+! a3=1, 7/(√23−4) = 7*(√23+4)/7 = 8+√23−4
+! a4=8, 1/(√23−4) = (√23+4)/7 = 1+(√23−3)/7
+! a5=1, 7/(√23−3) = 7*(√23+3)/14 = 3+(√23−3)/2
+! a6=3, 2/(√23−3) = 2*(√23+3)/14 = 1+(√23−4)/7
+! a7=1, 7/(√23−4) = 7*(√23+4)/7 = 8+√23−4
+
+! It can be seen that the sequence is repeating. For
+! conciseness, we use the notation √23=[4;(1,3,1,8)], to
+! indicate that the block (1,3,1,8) repeats indefinitely.
+
+! The first ten continued fraction representations of
+! (irrational) square roots are:
+
+! √2=[1;(2)] , period=1
+! √3=[1;(1,2)], period=2
+! √5=[2;(4)], period=1
+! √6=[2;(2,4)], period=2
+! √7=[2;(1,1,1,4)], period=4
+! √8=[2;(1,4)], period=2
+! √10=[3;(6)], period=1
+! √11=[3;(3,6)], period=2
+! √12=[3;(2,6)], period=2
+! √13=[3;(1,1,1,1,6)], period=5
+
+! Exactly four continued fractions, for N <= 13, have an odd period.
+
+! How many continued fractions for N <= 10000 have an odd period?
+
 <PRIVATE
 
 TUPLE: cont-frac
@@ -15,12 +67,7 @@ C: <cont-frac> cont-frac
     dup tuple>array rest cont-frac slots>tuple ;
 
 : create-cont-frac ( n -- n cont-frac )
-    dup sqrt >fixnum
+    dup sqrt >fixnum dup 1 <cont-frac> ;
-    [let :> root
-        root
-        root
-        1
-    ] <cont-frac> ;
 
 : step ( n cont-frac -- n cont-frac )
     swap dup
@@ -54,13 +101,10 @@ C: <cont-frac> cont-frac
         drop new-whole new-num-const new-denom <cont-frac>
     ] ;
 
-: loop ( c l n cont-frac -- c l n cont-frac )
+:: loop ( c l n cf -- c l n cf )
-    [let :> cf :> n :> l :> c
+    n cf step :> new-cf drop
-        n cf step
-        :> new-cf drop
     c 1 + l n new-cf
-        l new-cf = [ ] [ loop ] if
+    l new-cf = [ loop ] unless ;
-    ] ;
 
 : find-period ( n -- period )
     0 swap
@@ -70,7 +114,8 @@ C: <cont-frac> cont-frac
     loop
     drop drop drop ;
 
-: try-all ( -- n ) 2 10000 [a,b]
+: try-all ( -- n )
+    2 10000 [a,b]
     [ perfect-square? not ] filter
     [ find-period ] map
     [ odd? ] filter
@@ -81,52 +126,51 @@ PRIVATE>
 : euler064a ( -- n ) try-all ;
 
 <PRIVATE
+
 ! (√n + a)/b
 TUPLE: cfrac n a b ;
 
 C: <cfrac> cfrac
 
+: >cfrac< ( fr -- n a b )
+    [ n>> ] [ a>> ] [ b>> ] tri ;
+
 ! (√n + a) / b = 1 / (k + (√n + a') / b')
 !
 ! b / (√n + a) = b (√n - a) / (n - a^2) = (√n - a) / ((n - a^2) / b)
 :: reciprocal ( fr -- fr' )
-    fr n>>
+    fr >cfrac< :> ( n a b )
-    fr a>> neg
+    n
-    fr n>> fr a>> sq - fr b>> /
+    a neg
-    <cfrac>
+    n a sq - b /
-    ;
+    <cfrac> ;
 
 :: split ( fr -- k fr' )
-    fr n>> sqrt fr a>> + fr b>> / >integer
+    fr >cfrac< :> ( n a b )
-    dup fr n>> swap
+    n sqrt a + b / >integer
-    fr b>> * fr a>> swap -
+    dup n swap
-    fr b>>
+    b * a swap -
-    <cfrac>
+    b
-    ;
+    <cfrac> ;
 
 : pure ( n -- fr )
-    0 1 <cfrac>
+    0 1 <cfrac> ;
-    ;
 
 : next ( fr -- fr' )
-    reciprocal split nip
+    reciprocal split nip ;
-    ;
 
-:: period ( n -- per )
+:: period ( n -- period )
+    n sqrt >integer sq n = [ 0 ] [
         n pure split nip :> start
-    n sqrt >integer sq n =
+        1 start next
-    [ 0 ]
-    [ 1 start next
         [ dup start = not ]
         [ next [ 1 + ] dip ]
-      while
+        while drop
-      drop
+    ] if ;
-    ] if
-    ;
 
 PRIVATE>
 
 : euler064b ( -- ct )
-    1 10000 [a,b]
+    10000 [1,b] [ period odd? ] count ;
-    [ period odd? ] count
+
-    ;
+SOLUTION: euler064b