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 all

;; 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 all

> > (time primes 1000000)
937.407

Ok within a second so how much primes are there?

Code: Select all

> (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 all

(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 all

> (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 all

(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 all

> (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 all

(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 all

> (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 all

(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 all

> (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 all

(define (comp x y) 
    (>= (last x) (last y)))

And use it:

Code: Select all

(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, hope you enjoyed this little exercise!