project-euler: add solutions to 064, 087

extra/project-euler/064/064.factor

@@ -0,0 +1,132 @@
+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
+    ;

extra/project-euler/087/087.factor

@@ -0,0 +1,40 @@
+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 ;