extra/project-euler/087/087.factor

   1 ! Copyright (C) 2009 Kye W. Shi.
   2 ! See https://factorcode.org/license.txt for BSD license.
   3 USING: math math.functions math.primes project-euler.common
   4 sequences sets ;
   5 IN: project-euler.087
   6
   7 ! https://projecteuler.net/problem=87
   8
   9 ! DESCRIPTION
  10 ! -----------
  11
  12 ! The smallest number expressible as the sum of a prime square,
  13 ! prime cube, and prime fourth power is 28. In fact, there are
  14 ! exactly four numbers below fifty that can be expressed in such
  15 ! a way:
  16
  17 ! 28 = 2^2 + 2^3 + 2^4
  18 ! 33 = 3^2 + 2^3 + 2^4
  19 ! 49 = 5^2 + 2^3 + 2^4
  20 ! 47 = 2^2 + 3^3 + 2^4
  21
  22 ! How many numbers below fifty million can be expressed as the
  23 ! sum of a prime square, prime cube, and prime fourth power?
  24
  25 <PRIVATE
  26
  27 :: prime-powers-less-than ( primes pow n -- prime-powers )
  28     primes [ pow ^ ] map [ n <= ] filter ;
  29
  30 ! You may think to make a set of all possible sums of a prime
  31 ! square and cube and then subtract prime fourths from numbers
  32 ! ranging from 1 to 'n' to find this. As n grows large, this is
  33 ! actually more inefficient!
  34 !
  35 ! Prime numbers grow ~ n / log n
  36 !
  37 ! Thus there are (n / log n)^(1/2) prime squares <= n,
  38 !                (n / log n)^(1/3) prime cubes   <= n,
  39 !            and (n / log n)^(1/4) prime fourths <= n.
  40 !
  41 ! If we compute the cartesian product of these, this takes
  42 ! O((n / log n)^(13/12)).
  43 !
  44 ! If we instead precompute sums of squares and cubes, and
  45 ! iterate up to n, checking each fourth power against it, this
  46 ! takes:
  47 !
  48 ! O(n * (n / log n)^(1/4)) = O(n^(5/4)/(log n)^(1/4)) >> O((n / log n)^(13/12))
  49 !
  50 ! When n = 50,000,000, the first equation is approximately 10
  51 ! million and the second is approximately 2 billion.
  52
  53 :: prime-triples ( n -- answer )
  54     n sqrt primes-upto                :> primes
  55     primes 2 n prime-powers-less-than :> primes^2
  56     primes 3 n prime-powers-less-than :> primes^3
  57     primes 4 n prime-powers-less-than :> primes^4
  58     primes^2 primes^3 [ + ] cartesian-map concat
  59              primes^4 [ + ] cartesian-map concat
  60     [ n <= ] filter members length ;
  61
  62 PRIVATE>
  63
  64 :: euler087 ( -- answer )
  65     50,000,000 prime-triples ;
  66
  67 SOLUTION: euler087