From e5016951d06b87104317d4c0b499ddb769d24e20 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jens=20Axel=20S=C3=B8gaard?= Date: Tue, 11 Dec 2012 19:43:31 +0100 Subject: [PATCH] Improved performance of prime? for small numbers --- .../private/number-theory/number-theory.rkt | 25 ++--- .../private/number-theory/small-primes.rkt | 93 +++++++++++++++++++ 2 files changed, 107 insertions(+), 11 deletions(-) create mode 100644 collects/math/private/number-theory/small-primes.rkt diff --git a/collects/math/private/number-theory/number-theory.rkt b/collects/math/private/number-theory/number-theory.rkt index 604de881c3..f019e2942c 100644 --- a/collects/math/private/number-theory/number-theory.rkt +++ b/collects/math/private/number-theory/number-theory.rkt @@ -3,7 +3,8 @@ (require "../base/base-random.rkt" "divisibility.rkt" "modular-arithmetic.rkt" - "types.rkt") + "types.rkt" + "small-primes.rkt") (require/typed typed/racket [integer-sqrt/remainder (Natural -> (Values Natural Natural))]) @@ -55,10 +56,9 @@ (define prime-strong-pseudo-trials (integer-length (assert (/ 1 prime-strong-pseudo-certainty) integer?))) -(define *SMALL-PRIME-LIMIT* 10000) -; (define *SMALL-PRIME-LIMIT* 1000) ; use 1000 for coverage testing -; Determines the size of the pre-built table of small primes -(define *SMALL-FACORIZATION-LIMIT* *SMALL-PRIME-LIMIT*) +(define *VERY-SMALL-PRIME-LIMIT* 1000) +; Determines the size of the pre-built table of very small primes +(define *SMALL-FACORIZATION-LIMIT* *VERY-SMALL-PRIME-LIMIT*) ; Determines whether to use naive factorization or Pollards rho method. @@ -193,7 +193,7 @@ (define prime? (let () ; TODO: Only store odd integers in this table - (define N *SMALL-PRIME-LIMIT*) + (define N *VERY-SMALL-PRIME-LIMIT*) (define ps (make-vector (+ N 1) #t)) (define ! vector-set!) (! ps 0 #f) @@ -204,10 +204,13 @@ (! ps m #f)))) (lambda (n) (let ([n (abs n)]) - (if (< n N) - (vector-ref ps n) - (prime-strong-pseudo? n)))))) - + (cond + [(< n N) + (vector-ref ps n)] + [(< n *SMALL-PRIME-LIMIT*) + (small-prime? n)] + [else + (prime-strong-pseudo? n)]))))) (: next-prime : (case-> (Natural -> Natural) (Integer -> Integer))) @@ -304,7 +307,7 @@ (: factorize : Natural -> (Listof (List Natural Natural))) (define (factorize n) - (if (< n *SMALL-PRIME-LIMIT*) ; NOTE: Do measurement of best cut + (if (< n *SMALL-FACORIZATION-LIMIT*) ; NOTE: Do measurement of best cut (factorize-small n) (factorize-large n))) diff --git a/collects/math/private/number-theory/small-primes.rkt b/collects/math/private/number-theory/small-primes.rkt new file mode 100644 index 0000000000..b95b5f82ce --- /dev/null +++ b/collects/math/private/number-theory/small-primes.rkt @@ -0,0 +1,93 @@ +#lang typed/racket +(require/typed + data/bit-vector + [#:opaque BitVector bit-vector?] + [make-bit-vector (Integer Boolean -> BitVector)] + [bit-vector-set! (BitVector Integer Boolean -> Void)] + [bit-vector-ref (BitVector Integer -> Boolean)]) + +(provide small-prime? + *SMALL-PRIME-LIMIT*) + +; The moduli mod 60 that 2, 3 and 5 do not divide are: +(define non-235 '(1 7 11 13 17 19 23 29 31 37 41 43 47 49 53 59)) +; The differences of these numbers are: +(define deltas '( 6 4 2 4 2 4 6 2 6 4 2 4 2 4 6 2)) +; Note that there are exactly 16 of these moduli, so they fit in a u16. +; That is, a single u16 can represent a block of 60 numbers. + +(define mod60->bits (make-vector 60 (cast #f (U #f Integer)))) +(for ([x (in-list non-235)] + [b (in-naturals)]) + (vector-set! mod60->bits x b)) + +(define-syntax-rule (mod60->bit m) (vector-ref mod60->bits m)) + +(define *number-of-groups* 17000) ; each group contain 16 numbers +(define *SMALL-PRIME-LIMIT* (- (* 60 *number-of-groups*) 1)) + +; primes holds (* 60 *number-of-groups*) bits each +; representing a number not congruent to 2, 3, 5 +(define primes (make-bit-vector (* 60 *number-of-groups*) #t)) + +(define: (set-bit! [x : Integer]) : Void + (define-values (q r) (quotient/remainder x 60)) + (define b (mod60->bit r)) + (when b (bit-vector-set! primes (+ (* q 16) b) #t))) + +(define: (clear-bit! [x : Integer]) : Void + (define-values (q r) (quotient/remainder x 60)) + (define b (mod60->bit r)) + (when b (bit-vector-set! primes (+ (* q 16) b) #f))) + +(define: (bit [x : Integer]) : Boolean + (define-values (q r) (quotient/remainder x 60)) + (define b (mod60->bit r)) + (if b + (bit-vector-ref primes (+ (* q 16) b)) + #f)) + +(clear-bit! 1) ; 1 is not prime + +(define: (mark-composites [x : Integer]) : Void + ; x is prime => mark 2*x, 3*x, 4*x, 5*x, 6*x, 7*x, ... as prime + ; Well 2*x, 3*x, 4*x, 5*x, 6*x are not in our table, + ; so the first number to mark is 7*x . + ; Use the deltas to figure out which to mark. + (define y x) + (define delta*x 0) + (let loop ([ds deltas]) + ; (for ([delta (in-cycle deltas)] ... + (when (empty? ds) (set! ds deltas)) + (let ([delta (car ds)]) + (set! delta*x (* delta x)) + (cond + [(> y (- *SMALL-PRIME-LIMIT* delta*x)) + (void)] + [else + (set! y (+ y delta*x)) + (clear-bit! y) + (loop (cdr ds))])))) + +(define: (sieve) : Void + (define x 1) + (let loop ([ds deltas]) + ; (for ([delta (in-cycle deltas)] ... + (when (empty? ds) (set! ds deltas)) + (let ([delta (car ds)]) + (cond + [(> (* x x) (- *SMALL-PRIME-LIMIT* delta)) + (void)] + [else + ; x runs through all numbers incongruent to 2, 3 and 5 + (set! x (+ x delta)) + (when (bit x) ; x is prime + (mark-composites x)) + (loop (cdr ds))])))) + +(sieve) + +(define: (small-prime? [x : Integer]) : Boolean + (or (= x 2) (= x 3) (= x 5) + (and (mod60->bit (modulo x 60)) + (bit x))))