Saturday, September 08, 2012

A Small Experiment

One proof of the infinitude of prime numbers goes as follows:

Assume there were k primes, and we had them all in a list. Now consider the number one more than the product of all these primes. None of the primes in our list divides this number. Either it is prime, or it has a prime divisor not in our list. So there is no such k, that is, the primes are without number.

I started to think, which primes would we discover if this was our only method? I fiddled around for a few minutes with a piece of lisp code to generate these "euclidean" primes (Euclid included this proof long ago in the Elements).


(defun factor (n)
  "return a list of factors of n"
  (let ((test 3)
(limit (+ 1 (floor (sqrt n)))))
    (cond
      ((eq n 1) nil)
      ((evenp n) (cons 2 (factor (/ n 2))))
      (t (do()
 ((or (zerop (mod n test)) (> test limit)))
  (incf test 2))
(if (> test limit)
    (list n)
    (cons test  (factor (/ n test))))))))



(defun euclidean-primes (n)
  "generate the n primes generated by euclids proof"
  (let ((primes (list 1)))
    (do ((count 0 (1+ count)))
        ((= count n) (remove 1 primes))
      (let ((eprimorial (1+ (reduce #'* primes))))
        (setf primes (union primes (factor eprimorial)))))))

So I ran this, an up to (euclidean-primes 8) there is no issue, with a nice result:

(31051 1033 607 547 139 13 7 2 3 43 3263443 29881 67003 9119521 6212157481)

But running this for 9 proved time consuming, as the value to be factored was an unwieldy 
(FACTOR 12864938683278671740537145998360961546653259485195807)
which I found in the debugger when I was wondering why my CPU had pegged at 100% for several minutes. The overhead of bignum-truncate and bignum-normalize were swamping the process.

Now, I bet that my factorization is unnecessarily complicated, and that it's not the most efficient program for doing this, but it was an interesting thought experiment. Also, it appears union does some funny thing to element ordering, while I was expecting these to be ascending or descending, it looks like the smallest values tend to propagate to the center.

Just for a full picture:

(loop for i from 1 upto 8 do (print (euclidean-primes i)))
(2) 
(2 3) 
(3 2 7) 
(7 2 3 43) 
(43 3 2 7 13 139) 
(139 13 7 2 3 43 3263443) 
(3263443 43 3 2 7 13 139 547 607 1033 31051) 
(31051 1033 607 547 139 13 7 2 3 43 3263443 29881 67003 9119521 6212157481) 

No comments: