--- /dev/null
+USING: accessors arrays classes.tuple io kernel locals math math.functions
+ math.ranges prettyprint project-euler.common sequences ;
+IN: project-euler.064
+
+<PRIVATE
+
+TUPLE: cont-frac
+ { whole integer }
+ { num-const integer }
+ { denom integer } ;
+
+C: <cont-frac> cont-frac
+
+: deep-copy ( cont-frac -- cont-frac 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> ;
+
+: step ( n cont-frac -- n cont-frac )
+ swap dup
+ ! Store n
+ [let :> n
+ ! Extract the constant
+ swap dup num-const>>
+ :> num-const
+
+ ! Find the new denominator
+ num-const 2 ^ n swap -
+ :> exp-denom
+
+ ! Find the fraction in lowest terms
+ dup denom>>
+ exp-denom simple-gcd
+ exp-denom swap /
+ :> new-denom
+
+ ! Find the new whole number
+ num-const n sqrt + new-denom / >fixnum
+ :> new-whole
+
+ ! Find the new num-const
+ num-const new-denom /
+ new-whole swap -
+ new-denom *
+ :> new-num-const
+
+ ! Finally, update the continuing fraction
+ 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
+ ] ;
+
+: find-period ( n -- period )
+ 0 swap
+ create-cont-frac
+ step
+ deep-copy -rot
+ loop
+ drop drop drop ;
+
+: try-all ( -- n ) 2 10000 [a,b]
+ [ perfect-square? not ] filter
+ [ find-period ] map
+ [ odd? ] filter
+ length ;
+
+PRIVATE>
+
+: euler064a ( -- n ) try-all ;
+
+<PRIVATE
+! (√n + a)/b
+TUPLE: cfrac n a b ;
+
+C: <cfrac> cfrac
+
+! (√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>
+ ;
+
+:: split ( fr -- k fr' )
+ fr n>> sqrt fr a>> + fr b>> / >integer
+ dup fr n>> swap
+ fr b>> * fr a>> swap -
+ fr b>>
+ <cfrac>
+ ;
+
+: pure ( n -- fr )
+ 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
+ ;
+
+PRIVATE>
+
+: euler064b ( -- ct )
+ 1 10000 [a,b]
+ [ period odd? ] count
+ ;
--- /dev/null
+USING: locals math math.primes sequences math.functions sets kernel ;
+IN: project-euler.087
+
+<PRIVATE
+
+: remove-duplicates ( seq -- seq' )
+ dup intersect ;
+
+:: prime-powers-less-than ( primes pow n -- prime-powers )
+ primes [ pow ^ ] map [ n <= ] filter ;
+
+! You may think to make a set of all possible sums of a prime square and cube
+! and then subtract prime fourths from numbers ranging from 1 to 'n' to find
+! this. As n grows large, this is actually more inefficient!
+! Prime numbers grow ~ n / log n
+! Thus there are (n / log n)^(1/2) prime squares <= n,
+! (n / log n)^(1/3) prime cubes <= n,
+! and (n / log n)^(1/4) prime fourths <= n.
+! If we compute the cartesian product of these, this takes
+! O((n / log n)^(13/12)).
+! If we instead precompute sums of squares and cubes, and iterate up to n,
+! checking each fourth power against it, this takes
+! O(n * (n / log n)^(1/4)) = O(n^(5/4)/(log n)^(1/4)) >> O((n / log n)^(13/12))
+!
+! When n = 50000000, the first equation is approximately 10 million and
+! the second is approximately 2 billion.
+
+:: prime-triples ( n -- answer )
+ n sqrt primes-upto :> primes
+ primes 2 n prime-powers-less-than :> primes^2
+ primes 3 n prime-powers-less-than :> primes^3
+ primes 4 n prime-powers-less-than :> primes^4
+ primes^2 primes^3 [ + ] cartesian-map concat
+ primes^4 [ + ] cartesian-map concat
+ [ n <= ] filter remove-duplicates length ;
+
+PRIVATE>
+
+:: euler087 ( -- answer )
+ 50000000 prime-triples ;