extra/project-euler/050/050.factor

   1 ! Copyright (c) 2008 Aaron Schaefer.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: arrays kernel locals math math.primes sequences ;
   4 IN: project-euler.050
   5
   6 ! http://projecteuler.net/index.php?section=problems&id=50
   7
   8 ! DESCRIPTION
   9 ! -----------
  10
  11 ! The prime 41, can be written as the sum of six consecutive primes:
  12
  13 !     41 = 2 + 3 + 5 + 7 + 11 + 13
  14
  15 ! This is the longest sum of consecutive primes that adds to a prime below
  16 ! one-hundred.
  17
  18 ! The longest sum of consecutive primes below one-thousand that adds to a
  19 ! prime, contains 21 terms, and is equal to 953.
  20
  21 ! Which prime, below one-million, can be written as the sum of the most
  22 ! consecutive primes?
  23
  24
  25 ! SOLUTION
  26 ! --------
  27
  28 ! 1) Create an sequence of all primes under 1000000.
  29 ! 2) Start summing elements in the sequence until the next number would put you
  30 !    over 1000000.
  31 ! 3) Check if that sum is prime, if not, subtract the last number added.
  32 ! 4) Repeat step 3 until you get a prime number, and store it along with the
  33 !    how many consecutive numbers from the original sequence it took to get there.
  34 ! 5) Drop the first number from the sequence of primes, and do steps 2-4 again
  35 ! 6) Compare the longest chain from the first run with the second run, and store
  36 !    the longer of the two.
  37 ! 7) If the sequence of primes is still longer than the longest chain, then
  38 !    repeat steps 5-7...otherwise, you've found the longest sum of consecutive
  39 !    primes!
  40
  41 <PRIVATE
  42
  43 :: sum-upto ( seq limit -- length sum )
  44     0 seq [ + dup limit > ] find
  45     [ swapd - ] [ drop seq length swap ] if* ;
  46
  47 : pop-until-prime ( seq sum -- seq prime )
  48     over length 0 > [
  49         [ unclip-last-slice ] dip swap -
  50         dup prime? [ pop-until-prime ] unless
  51     ] [
  52         2drop { } 0
  53     ] if ;
  54
  55 ! a pair is { length of chain, prime the chain sums to }
  56
  57 : longest-prime ( seq limit -- pair )
  58     dupd sum-upto dup prime? [
  59         2array nip
  60     ] [
  61         [ head-slice ] dip pop-until-prime
  62         [ length ] dip 2array
  63     ] if ;
  64
  65 : longest ( pair pair -- longest )
  66     2dup [ first ] bi@ > [ drop ] [ nip ] if ;
  67
  68 : continue? ( pair seq -- ? )
  69     [ first ] [ length 1- ] bi* < ;
  70
  71 : (find-longest) ( best seq limit -- best )
  72     [ longest-prime longest ] 2keep 2over continue? [
  73         [ rest-slice ] dip (find-longest)
  74     ] [ 2drop ] if ;
  75
  76 : find-longest ( seq limit -- best )
  77     { 1 2 } -rot (find-longest) ;
  78
  79 : solve ( n -- answer )
  80     [ primes-upto ] keep find-longest second ;
  81
  82 PRIVATE>
  83
  84 : euler050 ( -- answer )
  85     1000000 solve ;
  86
  87 ! [ euler050 ] 100 ave-time
  88 ! 291 ms run / 20.6 ms GC ave time - 100 trials
  89
  90 MAIN: euler050