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)