Project Euler problem 46

[dougo]

On Friday, aarondf posted a call for help with problem 46 of Project Euler. It turned out to be a pretty simple problem, but aarondf had found a misleading definition of "odd composite" via Google which led him to miss the answer. I've become slightly obsessed with this problem, and have come up with a number of ways to write solution programs. You should try to write one yourself before reading further!

My first thought when it came to dealing with primes was to use infinite streams as described in Structure and Interpretation of Computer Programs (SICP). Fortunately, there is a SRFI for infinite streams, so I only needed to copy a few definitions from SICP. (Note that cons-stream in SICP is called stream-cons in SRFI 40.)

(require srfi/40)

(define (integers-starting-from n)
  (stream-cons n (integers-starting-from (+ n 1))))

(define integers (integers-starting-from 1))

(define (divisible? x y) (= (remainder x y) 0))

(define (square x) (* x x))

(define primes
  (stream-cons
   2
   (stream-filter prime? (integers-starting-from 3))))

(define (prime? n)
  (define (iter ps)
    (cond ((> (square (stream-car ps)) n) true)
          ((divisible? n (stream-car ps)) false)
          (else (iter (stream-cdr ps)))))
  (iter primes))

Yes, the definitions for primes and prime? are mutually recursive! Go back and read SICP if this is too mind-boggling.

There is also a nice stream utility package for PLT Scheme that adds a bunch of convenience functions for manipulating streams. With these, it's fairly straightforward to build up an infinite stream of counterexamples to Goldbach's conjecture. (Note that this is not the Goldbach conjecture.)

(require (planet "stream.ss" ("dherman" "stream.plt" 1 1)))

(define doubled-squares
  (stream-map (lambda (x) (* 2 x x)) integers))

(define odd-composites
  (stream-remove prime? (stream-filter odd? integers)))

;; Is there a doubled square m < n such that n = p + m for some prime p?
(define (goldbach? n)
  (stream-any (lambda (m) (prime? (- n m)))
              (stream-take-while (lambda (m) (< m n)) doubled-squares)))

(define goldbach-counterexamples
  (stream-remove goldbach? odd-composites))

You can then just get the answer with (stream-first goldbach-counterexamples). It turns out that there are only two counterexamples, or at least only two that are reasonably small—(stream-third goldbach-counterexamples) runs for longer than I am willing to wait.

tinhorn2 posted a very brief solution for Mathematica that displays all (i.e., both) solutions up to 100,000. I translated it into idiomatic Scheme:

(for ((n (in-range 4 50000)))
  (do ((i 0 (+ i 1)))
      ((not (and (not (prime? (- (+ (* 2 n) 1) (* 2 i i)))) (< (* i i) n)))
       (when (> (* i i) (- n 2))
         (display (+ (* 2 n) 1)) (newline)))))

(Note that the for construct for iteration over sequences is new in PLT Scheme v3.99.) It's comparatively ugly when compared to the original due to the lack of infix notation, but it's still pretty compact. There is a subtle shortcut here: instead of searching through odd composites and doubled squares greater than 1, it searches through all the odd naturals while considering 0 to be a doubled square. This way, any prime can be expressed as the sum of a prime and a doubled square, i.e. itself and zero, so no prime will be found as a counterexample, so we don't need to filter out primes ahead of time. This led me to a shorter solution using infinite streams:

(stream-find
 (lambda (n)
   (not (stream-any (lambda (i) (prime? (- n (* 2 i i))))
                    (stream-iota (ceiling (sqrt (/ n 2)))))))
 (stream-filter odd? (stream-iota)))

This takes another shortcut: instead of making a stream of doubled squares and stopping the stream at n, it makes a stream of naturals stopping at the square root of n/2. It also uses stream-find as an abbreviation for stream-first of stream-filter.

Then I got to thinking about how to write this using just for instead of streams. With the help of escape continuations, it came out fairly elegant. (You can think of calling an escape continuation as being analogous to raising an exception—it's a form of non-local return.)

(let/ec return
  (for ((n (in-naturals)) #:when (odd? n))
    (let/ec continue
      (for ((i (in-naturals)))
        (let ((p (- n (* 2 i i))))
          (cond ((negative? p) (return n))
                ((prime? p) (continue))))))))

This version avoids having to compute square roots, yet also avoids doing the doubled-square computation in both the end test and the prime test. However, going back to the square root version and using some alternate forms of for led to the most compact version I could come up with:

(for/or ((n (in-naturals)) #:when (odd? n))
  (for/and ((i (in-range 0 (ceiling (sqrt (/ n 2))))))
    (if (prime? (- n (* 2 i i))) false n)))

After all that, though, I think my favorite answer just sticks with explicit recursion using good old named-let:

(let loop ((n 1))
  (or (let loop ((i 1) (p n))
        (cond ((negative? p) n)
              ((prime? p) false)
              (else (loop (+ i 1) (- n (* 2 i i))))))
      (loop (+ n 2))))

It's probably not a coincidence that this is also the fastest solution, taking about 28ms on my machine.

Meanwhile, ahkond posted his solution in T-SQL, with an explanation (at my request). It's fascinating to see dynamic programming made concrete by constructing an actual table in a database.

Edit: Ha! I just noticed that prime? has a bug, and if you fix it, all the above solutions break (except the one translated from tinhorn2). Can you spot it?

Comments

[gemini6ice.livejournal.com]

stream-car reads the next item in the stream, and stream-cdr pulls off the next item in the stream? Oh, i barely remember this stuff.

[gemini6ice.livejournal.com]

Okay, the bug is that prime? 1 evaluates to true. I don't understand how fixing it would break your solutions though.

Can't you put a

((= n 1) false)

into your cond?

[dougo.livejournal.com]

The problem is that all the programs start their search at 1 (except tinhorn2's, which starts at 9). If 1 is composite, then it's a counterexample, because it cannot be expressed as a prime plus twice a square. (In fact, 1 is neither prime nor composite.)

Fortunately, they're all easy to fix. The stream-based programs can use (integers-starting-with 3) instead of integers or (stream-iota); the sequence-based programs can use (in-naturals 3); and the loop program can use (let ((n 3)) ...).

[dougo.livejournal.com]

That last one should be (let loop ((n 3)) ...).

[gemini6ice.livejournal.com]

just to be difficult, 1 = 5 + (2i)^2

;)

Another solution would be to simply change prime? to composite? and reverse your true and false statements:

(define primes
(stream-cons
2
(stream-filter not(composite? (integers-starting-from 3)))))

(define (composite? n)
(define (iter ps)
(cond ((> (square (stream-car ps)) n) false)
((divisible? n (stream-car ps)) true)
(else (iter (stream-cdr ps)))))
(iter primes))

[gemini6ice.livejournal.com]

in fact, i think you could do:

(define primes
(stream-filter not(composite? (integers-starting-from 2))))


huzzzah!

[gemini6ice.livejournal.com]

or maybe not, whoops. the mutually recursive definitions would cause an infinite loop here, maybe, if the first prime isn't explicitly defined.

[dougo.livejournal.com]

stream-car returns the first element of the stream, and stream-cdr returns the rest of the stream. The trick is that the elements of the stream aren't actually computed until stream-car is evaluated.

[darius.livejournal.com]

FWIW, here's my brute-force solution. I screwed up initially with (- s n) instead of (- n s).

(define (problem-46)
(let loop ((k 3))
(if (or (prime? k)
(sum-of-a-prime-and-twice-a-square? k))
(loop (+ k 2))
k)))

(define (sum-of-a-prime-and-twice-a-square? n)
(any prime? (map (lambda (s) (- n s))
(twice-squares-less-than n))))

(define (twice-squares-less-than n)
(let loop ((k 1) (results '()))
(let ((s (* 2 k k)))
(if (< s n)
(loop (+ k 1) (cons s results))
results))))

(define (prime? n)
(or (= n 2)
(and (not (even? n))
(let testing ((k 3))
(or (< n (* k k))
(and (not (= 0 (remainder n k)))
(testing (+ k 2))))))))

Yet another solution

[(Anonymous)]

It is nice to see the different angles of attack for this problem.
Here is a variation, where the inner loop iterates over primes.
If next-prime uses a precomputed prime table, then potentially
this could save some calls to prime?. But whether this kicks
in this example where n is relatively small, is another question.

/soegaard

The problem:
  Find first odd, composite n where 
         n = p + 2 x^2 , p prime, 
  doesn't have a solution.

Note:  If p<n is prime , then:
                    n = p + 2 x^2 has a solution
               <=>  (n-p)/2 is a square


(require (planet "math.ss" ("soegaard" "math.plt")))

(define (perfect-square? n)
  (let ([sqrt-n (integer-sqrt n)])
    (zero? (- n (* sqrt-n sqrt-n)))))

(let n-loop ([n 3])
  (if (prime? n)
      (n-loop (+ n 2))
      (let p-loop ([p 2])
        (cond [(> p n) 
               n]
              [(perfect-square? (quotient (- n p) 2))
               (n-loop (+ n 2))]
              [else
               (p-loop (next-prime p))]))))

Re: Yet another solution

[dougo.livejournal.com]

Yeah, I did it this way at first, looping through the primes and testing for square-ness. But then it seemed better to do it the other way, though it probably doesn't make much difference.

By the way, I tried to use math.plt, but it wouldn't compile for v3.99:

setup-plt: number-theory.ss:1578:21: compile: unbound variable in module in: reverse!

Re: Yet another solution

[(Anonymous)]

Yeah, use 372 for now.

But most likely you can get it
working on 399 with a search and
replace from reverse! to reverse.

/soegaard

[aarondf.livejournal.com]

Glad you're having fun with this.

you are evil

[stoneself.livejournal.com]

i just start project euler.