The Next Prime
March 26, 2010
Save-primes loops through the list of pre-computed primes, writing a byte every time k changes; bits accumulates the 1-bits representing the primes:
(define (save-primes n file-name)
(with-output-to-file file-name
(lambda ()
(let loop ((ps (primes n)) (k 0) (bits 0))
(cond ((null? ps) (display (integer->char bits)))
((< k (quotient (car ps) 30))
(display (integer->char bits))
(do ((k (+ k 1) (+ k 1)))
((= k (quotient (car ps) 30)))
(display (integer->char 0)))
(loop ps (quotient (car ps) 30) 0))
(else (case (modulo (car ps) 30)
((1) (loop (cdr ps) k (+ bits 1)))
((7) (loop (cdr ps) k (+ bits 2)))
((11) (loop (cdr ps) k (+ bits 4)))
((13) (loop (cdr ps) k (+ bits 8)))
((17) (loop (cdr ps) k (+ bits 16)))
((19) (loop (cdr ps) k (+ bits 32)))
((23) (loop (cdr ps) k (+ bits 64)))
((29) (loop (cdr ps) k (+ bits 128)))
(else (loop (cdr ps) k bits)))))))))
Load-primes reads the compressed primes into a global vector prime-bits; note that n can be smaller than the number of primes stored on disk, if you don’t need them all for your application:
(define prime-bits #f)
(define (load-primes n file-name)
(with-input-from-file file-name
(lambda ()
(let ((k-max (+ (quotient n 30) (if (zero? (modulo n 30)) 0 1))))
(set! prime-bits (make-vector k-max))
(do ((k 0 (+ k 1))) ((= k k-max))
(vector-set! prime-bits k (char->integer (read-char))))))))
Next-prime looks more complicated than it really is; Scheme’s poor support for bit-twiddling hurts us. We use a 2,3,5-wheel, as in the previous exercise, for values larger than the maximum prime, which means that we test only eight of each thirty values greater than n rather than the fifteen we would test if we examined each odd number. Next-bit returns the index and offset for the next 30k±1,7,11,13 greater than n. Bit-value converts a bit position within a byte to the corresponding 30k offset. Get-wheel positions the wheel for the beginning of iteration when n exceeds the maximum prime; it uses last-pair and cycle from the wheel factorization exercise. Here’s next-prime:
(define (next-prime n)
(define (next-bit n)
(let ((index (quotient n 30))
(offset (modulo n 30)))
(case offset
((0) (values index 1))
((1 2 3 4 5 6) (values index 2))
((7 8 9 10) (values index 4))
((11 12) (values index 8))
((13 14 15 16) (values index 16))
((17 18) (values index 32))
((19 20 21 22) (values index 64))
((23 24 25 26 27 28) (values index 128))
((29) (values (+ index 1) 1)))))
(define (bit-value offset)
(case offset
((1) 1) ((2) 7) ((4) 11) ((8) 13)
((16) 17) ((32) 19) ((64) 23) ((128) 29)))
(define (last-pair xs)
(if (null? (cdr xs)) xs
(last-pair (cdr xs))))
(define (cycle . xs)
(set-cdr! (last-pair xs) xs) xs)
(define (get-wheel n)
(let ((base (* (quotient n 30) 30))
(offset (modulo n 30)))
(case offset
((0) (values (+ base 1) (cycle 6 4 2 4 2 4 6 2)))
((1 2 3 4 5 6) (values (+ base 7) (cycle 4 2 4 2 4 6 2 6)))
((7 8 9 10) (values (+ base 11) (cycle 2 4 2 4 6 2 6 4)))
((11 12) (values (+ base 13) (cycle 4 2 4 6 2 6 4 2)))
((13 14 15 16) (values (+ base 17) (cycle 2 4 6 2 6 4 2 4)))
((17 18) (values (+ base 19) (cycle 4 6 2 6 4 2 4 2)))
((19 20 21 22) (values (+ base 23) (cycle 6 2 6 4 2 4 2 4)))
((23 24 25 26 27 28) (values (+ base 29) (cycle 2 6 4 2 4 2 4 6)))
((29) (values (+ base 31) (cycle 6 4 2 4 2 4 6 2))))))
(cond ((< n 2) 2) ((< n 3) 3) ((< n 5) 5)
((< n max-prime)
(let-values (((index offset) (next-bit n)))
(let loop ((index index) (offset offset))
(cond ((= offset 256) (loop (+ index 1) 1))
((zero? (logand (vector-ref prime-bits index) offset))
(loop index (* offset 2)))
(else (+ (* index 30) (bit-value offset)))))))
(else (let-values (((k wheel) (get-wheel n)))
(let loop ((k k) (wheel wheel))
(if (prime? k) k (loop (+ k (car wheel)) (cdr wheel))))))))
When saving and loading primes, it is best to make n divisible by 30, to fill out the last byte. To load the primes to a million, say:
> (save-primes 1000020 "prime.bits")
> (load-primes 1000020 "prime.bits")
> (define max-prime 1000003)
Timing tests show that it pays to pre-compute and store primes. On my machine, generating the primes to a million using a sieve takes 47 milliseconds, and cdr’ing through them takes 15 milliseconds. On the other hand, next-prime takes 0 milliseconds to load the pre-computed primes and 46 milliseconds to iterate through them; thus, unless your program needs the list of primes for something else, next-prime beats sieving. And next-prime saves space; the compressed primes occupy 33,334 bytes; the list of primes generated by the sieve takes four bytes for each of the 78,498 primes less than a million, plus a four-byte pointer from each prime to the next, plus a four-byte pointer to nil, for a total of 635,988 bytes (assuming a 32-bit machine), which is over nineteen times larger than the compressed data structure.
As an example of the use of next-prime, here is the goldbach function of a previous exercise, rewritten using next-prime:
(define (goldbach n)
(let loop ((p 2))
(if (prime? (- n p))
(list p (- n p))
(loop (next-prime p)))))
For example:
> (goldbach 986332)
(353 985979)
We used primes and prime? from previous exercises, and logand from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/KwCJJfYH.
Programming Praxis 
Another illustration of
next-primeis this alternate implementation of Pollard’s p-1 factorization algorithm, which we have studied in two previous exercises:(define (pollard n b1 b2)(let stage1 ((a 2) (p 2))
(if (< p b1)
(stage1 (expm a (expt p (ilog pb1)) n) (next-prime p))
(let ((d (gcd (- a 1) n)))
(if (< 1 d n) (list 'stage1 d)
(let stage2 ((p (next-prime b1)))
(if (< b2 p) #f
(let ((d (gcd (- (expm a p n) 1) n)))
(if (< 1 d n) (list 'stage2 d)
(stage2 (next-prime p)))))))))))
In addition to
next-prime, we useexpmandilogfrom the Standard Prelude. For example:> (pollard 15770708441 150 180)(stage2 135979)
A more serious effort is this factorization of the sixty-ninth repunit R69 = (10^69 – 1) / 9. The small factors 3, 37 and 277 are quickly found by trial division. Then Pollard’s p-1 method with bounds B1 = 30000 and B2 = 600000 find the fifteen digit factor 203864078068831. Applying the p-1 method again with the same bounds finds the twenty-three digit factor 11111111111111111111111, another repunit. The remaining twenty-eight digit factor 1595352086329224644348978893 is prime.
If you are interested in the factorization of repunits, see Makoto Kamada’s web page.
A million primes will not compress down to 33334 bytes. However, all primes up to one million compress down to this file size.
[…] : well, you might find this link of interest …it speaks to timing farther down the narrative https://programmingpraxis.com/2010/03/26/the-next-prime/Re: original ideas, the only things that come to mind is publishing via arxiv at some point if […]
A clojure solution, includes the Miller Rabin primality test, modified to use known deterministic tests when possible. To access the table, we use the fact that the expression x & -x masks out all bits in x except for the lower bit, so this is a constant time operation to get the smallest bit. We then index into a table mod 11 to convert the lowbit (which is a power of two) to the appropriate offset (1, 7, 11, 13, 17, 19, 23, or 29) I use an array of masks (based on the input mod 30,) to mask out the bits we will not consider prior to calculating x & -x. For primes > 1,000,000, I pretty much follow the reference solution.
(load "miller-rabin") (defn load-primes "Load a file of compressed primes" [filename] (let [input (java.io.DataInputStream. (java.io.BufferedInputStream. (java.io.FileInputStream. filename))), next-byte (fn [] (try (.readByte input) (catch java.io.EOFException e nil)))] (loop [data (vector-of :byte)] (let [x (next-byte)] (if (nil? x) data (recur (conj data x))))))) (defn lowbit "Get low order 1 bit in 1 byte value in constant time and convert to the prime mod 30" [x] (let [mod_30 (vector-of :int -1 1 7 -1 11 17 -1 29 13 23 19)] (mod_30 (mod (bit-and x (- x)) 11)))) (def masks (reduce into [ (vector-of :int) (repeat 1 2r11111111) ; 0 (repeat 6 2r11111110) ; 1 - 6 (repeat 4 2r11111100) ; 7 - 10 (repeat 2 2r11111000) ; 11 - 12 (repeat 4 2r11110000) ; 13 - 16 (repeat 2 2r11100000) ; 17 - 18 (repeat 4 2r11000000) ; 19 - 22 (repeat 6 2r10000000) ; 23 - 28 (repeat 1 2r00000000)])) ; 29 (def wheel235 (cycle [6 4 2 4 2 4 6 2])) (def start-wheel ; [offset of next possible prime & # of times to spin the wheel] (reduce into [ [] (repeat 1 [1 0]) ; 0 (repeat 6 [7 1]) ; 1 - 6 (repeat 4 [11 2]) ; 7 - 10 (repeat 2 [13 3]) ; 11 - 12 (repeat 4 [17 4]) ; 13 - 16 (repeat 2 [19 5]) ; 17 - 18 (repeat 4 [23 6]) ; 19 - 22 (repeat 6 [29 7]) ; 23 - 28 (repeat 1 [31 0])])) ; 29 (def prime-table (load-primes "prime1e6.dat")) (defn next-prime "Given n, return next prime > n For large n (> 3.4e14) tests are probabilistic" [n] (let [q (quot n 30), r (rem n 30)] (cond (< n 2) 2 (< n 3) 3 (< n 5) 5 (<= n 1000000) (loop [index q, mask (masks r)] (let [offset (lowbit (bit-and (prime-table index) mask))] (if (> offset 0) (+ (* index 30) offset) (recur (inc index) 16rFF)))) :otherwise (let [[offset, spins] (start-wheel r), base (* 30 q)] (loop [n (+ base offset), wheel (drop spins wheel235)] (if (prime? n) n (recur (+ (first wheel) n) (rest wheel)))))))) (def primes (iterate next-prime 2)) (defn primes-after [n] (drop 1 (iterate next-prime n)))To generate the file, there is another Clojure program that uses the (slow) O’Neil sieve to generate primes up to 1,000,030, and create a binary file.
(load "lazy-sieve") (defn bit_pos "Return byte # and bit position of a prime number" [n] (let [bit_offset {1 1, 7 2, 11 4, 13 8, 17 16, 19 32, 23 64, 29 128}] [(quot n 30), (bit_offset (mod n 30))])) (defn adjust_max [n] ; table size should be multiple of 30 (let [r (mod n 30)] (if (> r 0) (+ n (- 30 r)) n))) (defn to_signed_byte [b] ; java work-around, byte must be in range -128..127 (if (< b 128) b (- b 256))) (defn make-table "Create a compressed table of prime numbers with given size. (excluding 2,3,5)" [size] (let [max_prime (adjust_max size) table (into (vector-of :byte) (repeat (quot max_prime 30) 0)) primes (take-while #(< %1 max_prime) (drop 3 lazy-primes))] (loop [l primes, t table] (if (empty? l) t (let [[index, offset] (bit_pos (first l))] (recur (rest l) (assoc t index (to_signed_byte (bit-or (t index) offset))))))))) (def primes (make-table (int 1e6))) (def output (java.io.DataOutputStream. (java.io.BufferedOutputStream. (java.io.FileOutputStream. "prime1e6.dat")))) (doseq [p primes] (.writeByte output (int p))) (.close output)And some testing…