Fun with primes (and newLISP)
Posted: Tue May 14, 2019 10:36 pm
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:
Great a fast factor function! So how long this take for all the primes up to 1,000,000 ?
Ok within a second so how much primes are there?
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:
I really like my '(unless (rest (factor(.. bit here, so what about the time?
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:
But is it faster then using the built in factor function?
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:
So a generic function which we can use to calculate the differences between consecutive list members:
Ok great no we want to count all the differences, we'll use this function:
It is so nice to have all this functions in newLISP like count!
So applying this we get:
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:
And use it:
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!
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)
Code: Select all
> > (time primes 1000000)
937.407
Code: Select all
> (length (primes 1000000))
78498
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
)
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
)
Code: Select all
> (time (sieve 1000000))
259.257
Code: Select all
(define (funlist lst func rev)
(if rev
(map func (chop lst) (rest lst))
(map func (rest lst) (chop lst))))
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)
Code: Select all
(define (freq lst)
(let (ulist (unique (sort lst)))
(map list ulist (count ulist lst))))
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))
Code: Select all
(define (comp x y)
(>= (last x) (last y)))
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)