-USING: accessors arrays classes.tuple io kernel locals math math.functions
- math.ranges prettyprint project-euler.common sequences ;
+USING: accessors arrays classes.tuple io kernel locals math
+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
dup tuple>array rest cont-frac slots>tuple ;
: create-cont-frac ( n -- n cont-frac )
- dup sqrt >fixnum
- [let :> root
- root
- root
- 1
- ] <cont-frac> ;
+ dup sqrt >fixnum dup 1 <cont-frac> ;
: step ( n cont-frac -- n cont-frac )
swap dup
drop new-whole new-num-const new-denom <cont-frac>
] ;
-: loop ( c l n cont-frac -- c l n cont-frac )
- [let :> cf :> n :> l :> c
- n cf step
- :> new-cf drop
- c 1 + l n new-cf
- l new-cf = [ ] [ loop ] if
- ] ;
+:: loop ( c l n cf -- c l n cf )
+ n cf step :> new-cf drop
+ c 1 + l n new-cf
+ l new-cf = [ loop ] unless ;
: find-period ( n -- period )
0 swap
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
: 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 a>> neg
- fr n>> fr a>> sq - fr b>> /
- <cfrac>
- ;
+ fr >cfrac< :> ( n a b )
+ n
+ a neg
+ n a sq - b /
+ <cfrac> ;
:: split ( fr -- k fr' )
- fr n>> sqrt fr a>> + fr b>> / >integer
- dup fr n>> swap
- fr b>> * fr a>> swap -
- fr b>>
- <cfrac>
- ;
+ fr >cfrac< :> ( n a b )
+ n sqrt a + b / >integer
+ dup n swap
+ b * a swap -
+ b
+ <cfrac> ;
: pure ( n -- fr )
- 0 1 <cfrac>
- ;
+ 0 1 <cfrac> ;
: next ( fr -- fr' )
- reciprocal split nip
- ;
-
-:: period ( n -- per )
- n pure split nip :> start
- n sqrt >integer sq n =
- [ 0 ]
- [ 1 start next
- [ dup start = not ]
- [ next [ 1 + ] dip ]
- while
- drop
- ] if
- ;
+ reciprocal split nip ;
+
+:: period ( n -- period )
+ n sqrt >integer sq n = [ 0 ] [
+ n pure split nip :> start
+ 1 start next
+ [ dup start = not ]
+ [ next [ 1 + ] dip ]
+ while drop
+ ] if ;
PRIVATE>
: euler064b ( -- ct )
- 1 10000 [a,b]
- [ period odd? ] count
- ;
+ 10000 [1,b] [ period odd? ] count ;
+
+SOLUTION: euler064b