Merge pull request #752 from wrog/math_prime_fixes

(chibi math prime) fix miller-rabin-composite?, factor, etc (issue #751), add factor-alist

lib/chibi/math/prime-test.sld

@@ -47,6 +47,7 @@
       (test 5 (prime-below 7))
       (test 797 (prime-below 808))
 
+      (test 1 (totient 1))
       (test 1 (totient 2))
       (test 2 (totient 3))
       (test 2 (totient 4))
@@ -56,6 +57,7 @@
       (test 4 (totient 8))
       (test 6 (totient 9))
       (test 4 (totient 10))
+      (test-error (totient 0))
 
       (test #f (perfect? 1))
       (test #f (perfect? 2))
@@ -71,7 +73,7 @@
       (test #t (perfect? 496))
       (test #t (perfect? 8128))
 
-      (test '(1) (factor 1))
+      (test '() (factor 1))
       (test '(2) (factor 2))
       (test '(3) (factor 3))
       (test '(2 2) (factor 4))
@@ -86,8 +88,16 @@
       (test '(2 3 3) (factor 18))
       (test '(2 2 2 3 3) (factor 72))
       (test '(3 3 3 5 7) (factor 945))
+      (test-error (factor 0))
 
+      (test '() (factor-alist 1))
+      (test '((2 . 3) (3 . 2)) (factor-alist 72))
+      (test '((3 . 3) (5 . 1) (7 . 1)) (factor-alist 945))
+      (test-error (factor-alist 0))
+
+      (test 0 (aliquot 1))
       (test 975 (aliquot 945))
+      (test-error (aliquot 0))
 
       (do ((i 3 (+ i 2)))
           ((>= i 101))
@@ -107,4 +117,7 @@
                         5772301760555853353
                         (* 2936546443 3213384203)))
 
+      (test "Miller-Rabin vs. Carmichael prime"
+            #t (miller-rabin-composite? 118901521))
+
       (test-end))))

lib/chibi/math/prime.scm

@@ -4,12 +4,13 @@
 
 ;;> Prime and number theoretic utilities.
 
-;;> Returns a pair whose car is the power of 2 in the factorization of
+;; Given \var{n} and a continuation \var{return},
-;;> n, and whose cdr is the product of all remaining primes.
+;; returns (\var{return} \var{k2} \var{n2}) where
-(define (factor-twos n)
+;; \var{k2} is the power of 2 in the factorization of \var{n}, and
-  (do ((p 0 (+ p 1))
+;; \var{n2} is product of all other prime powers dividing \var{n}
-       (r n (arithmetic-shift r -1)))
+(define (factor-twos n return)
-      ((odd? r) (cons p r))))
+  (let ((b (first-set-bit n)))
+    (return b (arithmetic-shift n (- b)))))
 
 ;;> Returns the multiplicative inverse of \var{a} modulo \var{b}.
 (define (modular-inverse a b)
@@ -73,22 +74,36 @@
 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 ;; Probable primes.
 
-(define (modular-root-of-one? twos odd a n neg1)
+;; Given \var{n}, return a predicate that tests whether
-  ;; Returns true iff any (modular-expt a odd*2^i n) for i=0..twos-1
+;; its argument \var{a} is a witness for \var{n} not being prime,
-  ;; returns 1 modulo n.
+;; either (1) because \var{a}^(\var{n}-1)≠1 mod \var{n}
+;; \em{or} (2) because \var{a}'s powers include
+;; a third square root of 1 beyond {1, -1}
+(define (miller-rabin-witnesser n)
+  (let ((neg1 (- n 1)))
+    (factor-twos neg1
+     (lambda (twos odd)
+       (lambda (a)
          (let ((b (modular-expt a odd n)))
            (let lp ((i 0) (b b))
-      (cond ((or (= b 1) (= b neg1))) ; in (= b 1) case we could factor
+             (cond ((= b neg1)
-            ((>= i twos) #f)
+                    ;; found -1 (expected sqrt(1))
-            (else (lp (+ i 1) (remainder (* b b) n)))))))
+                    #f)
+                   ((= b 1)
+                    ;; !! (previous b)^2=1 and was not 1 or -1
+                    (not (zero? i)))
+                   ((>= i twos)
+                    ;; !! a^(n-1)!=1 mod n
+                    )
+                   (else
+                    (lp (+ i 1) (remainder (* b b) n)))))))))))
 
-;;> Returns true if we can show \var{n} to be composite by finding an
+;;> Returns true if we can show \var{n} to be composite
-;;> exception to the Miller Rabin lemma.
+;;> using the Miller-Rabin test (i.e., finding a witness \var{a}
+;;> where \var{a}^(\var{n}-1)≠1 mod \var{n} or \var{a} reveals
+;;> the existence of a 3rd square root of 1 in \b{Z}/(n))
 (define (miller-rabin-composite? n)
-  (let* ((neg1 (- n 1))
+  (let* ((witness? (miller-rabin-witnesser n))
-         (factors (factor-twos neg1))
-         (twos (car factors))
-         (odd (cdr factors))
          ;; Each iteration of Miller Rabin reduces the odds by 1/4, so
          ;; this is a 1 in 2^40 probability of false positive,
          ;; assuming good randomness from SRFI 27 and no bugs, further
@@ -97,11 +112,10 @@
          (rand-limit (if (< n 341550071728321) fixed-limit 20)))
     (let try ((i 0))
       (and (< i rand-limit)
-           (let ((a (if (< i fixed-limit)
+           (or (witness? (if (< i fixed-limit)
                              (vector-ref prime-table i)
-                        (+ (random-integer (- n 3)) 2))))
+                             (+ (random-integer (- n 3)) 2)))
-             (or (not (modular-root-of-one? twos odd a n neg1))
+               (try (+ i 1)))))))
-                 (try (+ i 1))))))))
 
 ;;> Returns true if \var{n} has a very high probability (enough that
 ;;> you can assume a false positive will never occur in your lifetime)
@@ -175,57 +189,84 @@
          ((prime? n) n)
          (else (lp (+ n 2)))))))))
 
+;; Given an initial value \var{r1} representing the (empty)
+;; factorization of 1 and a procedure \var{put}
+;; (called as \scheme{(\var{put} \var{r} \var{p} \var{k})})
+;; that, given prior representation \var{r},
+;; adds a prime factor \var{p} of multiplicity \var{k},
+;; returns a factorization function which returns the factorization
+;; of its non-zero integer argument \var{n} in this representation.
+;; The optional 3rd and 4th arguments, if provided, specialize \var{put}
+;; for particular primes:
+;; \var{put2}  for \var{p}=2,  called as \scheme{(\var{put2} \var{r} \var{k})})
+;; \var{put-1} for \var{p}=-1, called as \scheme{(\var{put-1} \var{r})}).
+(define (make-factorizer r1 put . o)
+  (let-optionals o ((put2  (lambda (r k) (put r  2 k)))
+                    (put-1 (lambda (r)   (put r -1 1))))
+    (lambda (n)
+      (when (zero? n)
+        (error "cannot factor 0"))
+      (factor-twos
+       n
+       (lambda (k2 n)
+         (let lp ((i 3) (ii 9)
+                  (n (abs n))
+                  (res (let ((res (if (negative? n) (put-1 r1) r1)))
+                         (if (zero? k2) res (put2 res k2)))))
+           (let next-i ((i i) (ii ii))
+             (cond ((> ii n)
+                    (if (= n 1) res (put res n 1)))
+                   ((not (zero? (remainder n i)))
+                    (next-i (+ i 2) (+ ii (* (+ i 1) 4))))
+                   (else
+                    (let rest ((n (quotient n i))
+                               (k 1))
+                      (if (zero? (remainder n i))
+                          (rest (quotient n i) (+ k 1))
+                          (lp (+ i 2) (+ ii (* (+ i 1) 4))
+                              n (put res i k)))))))))))))
+
+;;> Returns the factorization of \var{n} as a list of
+;;> elements of the form \scheme{(\var{p} . \var{k})},
+;;> where \var{p} is a prime factor
+;;> and \var{k} is its multiplicity.
+(define factor-alist
+  (let ((rfactor (make-factorizer '()
+                  (lambda (l p k) (cons (cons p k) l)))))
+    (lambda (n) (reverse (rfactor n)))))
+
 ;;> Returns the factorization of \var{n} as a monotonically
 ;;> increasing list of primes.
-(define (factor n)
+(define factor
-  (cond
+  (let ((rfactor (make-factorizer '()
-   ((negative? n)
+                  (lambda (l p k) (cons (make-list k p) l)))))
-    (cons -1 (factor (- n))))
+    (lambda (n) (concatenate! (reverse (rfactor n))))))
-   ((<= n 2)
-    (list n))
-   (else
-    (let lp ((n n)
-             (res (list)))
-      (cond
-       ((even? n)
-        (lp (quotient n 2) (cons 2 res)))
-       ((= n 1)
-        (reverse res))
-       (else
-        (let lp ((i 3) (n n) (limit (exact (ceiling (sqrt n)))) (res res))
-          (cond
-           ((= n 1)
-            (reverse res))
-           ((> i limit)
-            (reverse (cons n res)))
-           ((zero? (remainder n i))
-            (lp i (quotient n i) limit (cons i res)))
-           (else
-            (lp (+ i 2) n limit res))))))))))
 
-;;> Returns the Euler totient function, the number of positive
+;;> The Euler totient φ(\var{n}) is the number of positive
-;;> integers less than \var{n} that are relatively prime to \var{n}.
+;;> integers less than or equal to \var{n} that are
-(define (totient n)
+;;> relatively prime to \var{n}.
-  (let ((limit (exact (ceiling (sqrt n)))))
+(define totient
-    (let lp ((i 2) (count 1))
+  (make-factorizer 1
-      (cond ((> i limit)
+   (lambda (tot p k)
-             (if (= count (- i 1))
+     (* tot (- p 1) (expt p (- k 1))))
-                 (- n 1)                ; shortcut for prime
+   (lambda (tot k)
-                 (let lp ((i i) (count count))
+     (arithmetic-shift tot (- k 1)))
-                   (cond ((>= i n) count)
+   (lambda (_)
-                         ((= 1 (gcd n i)) (lp (+ i 1) (+ count 1)))
+     (error "totient of negative number?"))))
-                         (else (lp (+ i 1) count))))))
+
-            ((= 1 (gcd n i)) (lp (+ i 1) (+ count 1)))
+;;> The aliquot sum s(\var{n}) is
-            (else (lp (+ i 1) count))))))
+;;> the sum of proper divisors of a positive integer \var{n}.
+(define aliquot
+  (let ((aliquot+n
+         (make-factorizer 1
+          (lambda (aliq p k)
+            (* aliq (quotient (- (expt p (+ k 1)) 1) (- p 1))))
+          (lambda (aliq k)
+            (- (arithmetic-shift aliq (+ k 1)) aliq))
+          (lambda (_)
+            (error "aliquot of negative number?")))))
+    (lambda (n) (- (aliquot+n n) n))))
 
-;;> The aliquot sum s(n), equal to the sum of proper divisors of an
-;;> integer n.
-(define (aliquot n)
-  (let ((limit (+ 1 (quotient n 2))))
-    (let lp ((i 2) (sum 1))
-      (cond ((> i limit) sum)
-            ((zero? (remainder n i)) (lp (+ i 1) (+ sum i)))
-            (else (lp (+ i 1) sum))))))
 
 ;;> Returns true iff \var{n} is a perfect number, i.e. the sum of its
 ;;> divisors other than itself equals itself.

lib/chibi/math/prime.sld

@@ -1,11 +1,12 @@
 
 (define-library (chibi math prime)
-  (import (scheme base) (scheme inexact) (srfi 27))
+  (import (scheme base) (scheme inexact) (chibi optional) (srfi 1) (srfi 27))
   (cond-expand
    ((library (srfi 151)) (import (srfi 151)))
    ((library (srfi 33)) (import (srfi 33)))
    (else (import (srfi 60))))
-  (export prime? nth-prime prime-above prime-below factor perfect?
+  (export prime? nth-prime prime-above prime-below
+          factor factor-alist perfect?
           totient aliquot
           provable-prime? probable-prime?
           random-prime random-prime-distinct-from