Fun with primes (and newLISP)

[fdb]

I think newLISP is a great tool for tinkering with some ideas in your head, AKA 'exploratory computing'. If you like this as well maybe you'll also like the little write up below.



Recently I had to explain 'what is a prime' to my eldest daughter, while explaining 'there are an infinite number of primes' I thought, what about the distance between primes.., maybe there is a certain distribution?



I resisted the easiest way (google it) and went for the best (explore it in newLisp).



Ok , so we need to the primes, how to generate them ? Look in the documentation and I got this:

Code Select Expand


;; primes.lsp - return all primes in a list, up to n

(define (primes n , p)
  (dotimes (e n)
    (if (= (length (factor e)) 1)
      (push e p -1))) p)

Great a fast factor function! So how long this take for all the primes up to 1,000,000 ?

Code Select Expand

> > (time primes 1000000)
937.407
Ok within a second so how much primes are there?

Code Select Expand

> (length (primes 1000000))
78498
So about 80000 primes /second, can we do better? Well an obvious improvement is to step through the odd numbers only and start with 2:

Code Select Expand

(define (primes-to n , (p-list '(2)))
(for (x 3 n 2)
(unless (rest (factor x))
(push x p-list -1)))
p-list
)
I really like my '(unless (rest (factor(.. bit here, so what about the time?

Code Select Expand

> (time (primes-to 1000000))
569.302

Aha! as expected almost twice as fast, now one last improvement before we go to the main topic, of course to generate primes has been known since antiquity so how does the  'Sieve of Eratosthenes' look in newLisp? Here is my version:

Code Select Expand

(define (sieve n)
(set 'arr (array (+ n 1)) 'lst '(2))
(for (x 3 n 2)
(when (not (arr x))
(push x lst -1)
(for (y (* x x) n (* 2 x) (> y n))
(setf (arr y) true))))
lst
)
But is it faster then using the built in factor function?

Code Select Expand

> (time (sieve 1000000))
259.257
Not bad for antiquity, we'll be using this function from now on.So what we want to do now is calculate the distance between primes. After some tinkering I came up with this:

Code Select Expand

(define (funlist lst func rev)
(if rev
(map func (chop lst) (rest lst))
(map func (rest lst) (chop lst))))
So a generic function which we can use to calculate the differences between consecutive list members:

Code Select Expand

> (funlist (sieve 1000) -)
(1 2 2 4 2 4 2 4 6 2 6 4 2 4 6 6 2 6 4 2 6 4 6 8 4 2 4 2 4 14 4 6 2 10 2 6 6 4 6
 6 2 10 2 4 2 12 12 4 2 4 6 2 10 6 6 6 2 6 4 2 10 14 4 2 4 14 6 10 2 4 6 8 6 6 4
 6 8 4 8 10 2 10 2 6 4 6 8 4 2 4 12 8 4 8 4 6 12 2 18 6 10 6 6 2 6 10 6 6 2 6 6 4
 2 12 10 2 4 6 6 2 12 4 6 8 10 8 10 8 6 6 4 8 6 4 8 4 14 10 12 2 10 2 4 2 10 14 4
 2 4 14 4 2 4 20 4 8 10 8 4 6 6 14 4 6 6 8 6)
Ok great no we want to count all the differences, we'll use this function:

Code Select Expand

(define (freq lst)
(let (ulist (unique (sort lst)))
(map list ulist (count ulist lst))))
It is so nice to have all this functions in newLISP like count!

So applying this we get:

Code Select Expand

> (freq (funlist (sieve 1000) -))
((1 1) (2 35) (4 40) (6 44) (8 15) (10 16) (12 7) (14 7) (18 1) (20 1))
Hmm for all the primes below 1000, 6 is the most frequent difference,i didn't expect that, how about all the primes below 1,000,000 ? We probably need to sort the list to get the highest first, so lets make a function to use in the sort:

Code Select Expand

(define (comp x y)
    (>= (last x) (last y)))
And use it:

Code Select Expand

(sort(freq(funlist (sieve 1000000) -)) comp)
((6 13549) (2 8169) (4 8143) (12 8005) (10 7079) (8 5569) (18 4909) (14 4233) (16
  2881)
 (24 2682)
 (20 2401)
 (22 2172)
 (30 1914)
 (28 1234)
 (26 1175)
 (36 767)
 (34 557)
 (32 550)

Definitely looks like number 6, does this take forever? Have a look at https://en.wikipedia.org/wiki/Prime_gap">//https://en.wikipedia.org/wiki/Prime_gap, hope you enjoyed this little  exercise!