extra/project-euler/027/027.factor

   1 ! Copyright (c) 2008 Aaron Schaefer.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: kernel math math.primes math.ranges project-euler.common sequences ;
   4 IN: project-euler.027
   5
   6 ! http://projecteuler.net/index.php?section=problems&id=27
   7
   8 ! DESCRIPTION
   9 ! -----------
  10
  11 ! Euler published the remarkable quadratic formula:
  12
  13 !     n² + n + 41
  14
  15 ! It turns out that the formula will produce 40 primes for the consecutive
  16 ! values n = 0 to 39. However, when n = 40, 402 + 40 + 41 = 40(40 + 1) + 41 is
  17 ! divisible by 41, and certainly when n = 41, 41² + 41 + 41 is clearly
  18 ! divisible by 41.
  19
  20 ! Using computers, the incredible formula n² - 79n + 1601 was discovered, which
  21 ! produces 80 primes for the consecutive values n = 0 to 79. The product of the
  22 ! coefficients, -79 and 1601, is -126479.
  23
  24 ! Considering quadratics of the form:
  25
  26 !     n² + an + b, where |a| < 1000 and |b| < 1000
  27
  28 !     where |n| is the modulus/absolute value of n
  29 !     e.g. |11| = 11 and |-4| = 4
  30
  31 ! Find the product of the coefficients, a and b, for the quadratic expression
  32 ! that produces the maximum number of primes for consecutive values of n,
  33 ! starting with n = 0.
  34
  35
  36 ! SOLUTION
  37 ! --------
  38
  39 ! b must be prime since n = 0 must return a prime
  40 ! a + b + 1 must be prime since n = 1 must return a prime
  41 ! 1 - a + b must be prime as well, hence >= 2. Therefore:
  42 !    1 - a + b >= 2
  43 !        b - a >= 1
  44 !            a < b
  45
  46 <PRIVATE
  47
  48 : source-027 ( -- seq )
  49     1000 [0,b) [ prime? ] filter [ dup [ neg ] map append ] keep
  50     cartesian-product [ first2 < ] filter ;
  51
  52 : quadratic ( b a n -- m )
  53     dup sq -rot * + + ;
  54
  55 : (consecutive-primes) ( b a n -- m )
  56     3dup quadratic prime? [ 1 + (consecutive-primes) ] [ 2nip ] if ;
  57
  58 : consecutive-primes ( a b -- m )
  59     swap 0 (consecutive-primes) ;
  60
  61 : max-consecutive ( seq -- elt n )
  62     dup [ first2 consecutive-primes ] map dup supremum
  63     over index [ swap nth ] curry bi@ ;
  64
  65 PRIVATE>
  66
  67 : euler027 ( -- answer )
  68     source-027 max-consecutive drop product ;
  69
  70 ! [ euler027 ] 100 ave-time
  71 ! 111 ms ave run time - 6.07 SD (100 trials)
  72
  73 ! TODO: generalize max-consecutive/max-product (from #26) into a new word
  74
  75 SOLUTION: euler027