From e2e88c21fe021f645f012b073fd7c56243a6ce68 Mon Sep 17 00:00:00 2001
From: Hugo Sansaqua <privet.kitty99@gmail.com>
Date: Wed, 15 Nov 2023 07:59:43 +0900
Subject: [PATCH] Improve prime-p

---
 module/primality.lisp      | 44 +++++++++++++++++++-------------------
 module/test/primality.lisp |  3 +++
 2 files changed, 25 insertions(+), 22 deletions(-)

diff --git a/module/primality.lisp b/module/primality.lisp
index 427568a..b16cd67 100644
--- a/module/primality.lisp
+++ b/module/primality.lisp
@@ -8,31 +8,32 @@ This module is tuned for SBCL on x86-64, i.e., (INTEGER 0
 #.MOST-POSITIVE-FIXNUM) is assumed to be (UNSIGNED-BYTE 62)"))
 (in-package :cp/primality)
 
+;; NOTE: not sufficiently tested
+
 (deftype uint () '(unsigned-byte 62))
 
 (defun %strong-probable-prime-p (n base)
   (declare (optimize (speed 3))
            (uint n base))
-  (or (= n base)
-      (labels ((mod-power (base power)
-                 (declare (uint base power))
-                 (loop with res of-type uint = 1
-                       while (> power 0)
-                       when (oddp power)
-                       do (setq res (mod (* res base) n))
-                       do (setq base (mod (* base base) n)
-                                power (ash power -1))
-                       finally (return res))))
-        (let* ((d (floor (- n 1) (logand (- n 1) (- 1 n))))
-               (y (mod-power base d)))
-          (declare (uint y))
-          (or (= y 1)
-              (= y (- n 1))
-              (let ((s (tzcount (- n 1))))
-                (loop repeat (- s 1)
-                      do (setq y (mod (* y y) n))
-                         (when (<= y 1) (return nil))
-                         (when (= y (- n 1)) (return t)))))))))
+  (labels ((mod-power (base power)
+             (declare (uint base power))
+             (loop with res of-type uint = 1
+                   while (> power 0)
+                   when (oddp power)
+                   do (setq res (mod (* res base) n))
+                   do (setq base (mod (* base base) n)
+                            power (ash power -1))
+                   finally (return res))))
+    (let* ((d (floor (- n 1) (logand (- n 1) (- 1 n))))
+           (y (mod-power base d)))
+      (declare (uint y))
+      (or (<= y 1)
+          (= y (- n 1))
+          (let ((s (tzcount (- n 1))))
+            (loop repeat (- s 1)
+                  do (setq y (mod (* y y) n))
+                     (when (<= y 1) (return nil))
+                     (when (= y (- n 1)) (return t))))))))
 
 ;; https://primes.utm.edu/prove/prove2_3.html
 ;; https://miller-rabin.appspot.com/
@@ -42,8 +43,7 @@ This module is tuned for SBCL on x86-64, i.e., (INTEGER 0
   (cond ((<= n 1) nil)
         ((evenp n) (= n 2))
         ((< n 49141)
-         (or (= n 331)
-             (%strong-probable-prime-p n 921211727)))
+         (%strong-probable-prime-p n 921211727))
         ((< n 4759123141)
          (loop for base in '(2 7 61)
                always (%strong-probable-prime-p n base)))
diff --git a/module/test/primality.lisp b/module/test/primality.lisp
index 6ec3011..0172c91 100644
--- a/module/test/primality.lisp
+++ b/module/test/primality.lisp
@@ -13,4 +13,7 @@
         (assert (eq (if (prime-p x) 1 0) (aref table x)))))
     (loop repeat 5000
           for x = (+ 4759123141 (* 2 (random 10000000 state)))
+          do (is (eq (prime-p x) (sb-int:positive-primep x))))
+    (loop repeat 5000
+          for x = (+ 2152302898747 (* 2 (random 10000000 state)))
           do (is (eq (prime-p x) (sb-int:positive-primep x))))))