cleanup, add test cases

janet-lang · Jul 27, 2023 · e74b793 · e74b793
1 parent 1b3a3c3
commit e74b793
Show file tree

Hide file tree

Showing 2 changed files with 194 additions and 18 deletions.
diff --git a/spork/math.janet b/spork/math.janet
@@ -1212,8 +1212,8 @@
   [2 325 9375 28178 450775 9780504 1795265022])
 
 (def- prime-prod
-  # (* ;(take 13 small-primes))
-  304250263527210)
+  #(* ;(array/slice small-primes 1 14))
+  6541380665835015)
 
 (defn miller-rabin-prp?
   ``Performs a Miller-Rabin probable prime test on `n` against all given bases.
@@ -1241,33 +1241,32 @@
     true))
 
 (defn prime?
-  ``A primality test, deterministic for all `n` less than 2^63.``
+  "A primality test, deterministic for all `n` less than 2^63."
   [n]
   (cond
     (compare<= n 211) (do
                         (def m (binary-search n small-primes compare<))
                         (compare= n (in small-primes m)))
+    (zero? (mod n 2)) false
     (= 0 (jacobi n prime-prod)) false
     (miller-rabin-prp? n)))
 
 (defn next-prime
   "Returns the next prime number strictly larger than `n`."
   [n]
-  (label result
-    (def x (+ n 1))
-    (if (compare<= x 2) (return result 2))
-    (when (compare<= x 211)
-      (def m (binary-search x small-primes compare<))
-      (return result (in small-primes m)))
-    (def j (mod x 210))
-    (def m (binary-search j soe-indices compare<))
-    (var i (+ x (- (in soe-indices m) j)))
-    (def offs [;(slice soe-offsets m) ;(slice soe-offsets 0 m)])
-    (forever
-      (each o offs
-        (if (and (not= 0 (jacobi i prime-prod)) (miller-rabin-prp? i))
-          (return result i))
-        (+= i o)))))
+  (def x (+ n 1))
+  (if (compare<= x 211)
+    (in small-primes (binary-search x small-primes compare<))
+    (label result
+      (def j (mod x 210))
+      (def m (binary-search j soe-indices compare<))
+      (var i (+ x (- (in soe-indices m) j)))
+      (def offs [;(slice soe-offsets m) ;(slice soe-offsets 0 m)])
+      (forever
+        (each o offs
+          (if (and (not= 0 (jacobi i prime-prod)) (miller-rabin-prp? i))
+            (return result i))
+          (+= i o))))))
 
 (defn primes
   "A boundless prime number generator."

diff --git a/test/suite-math.janet b/test/suite-math.janet
@@ -421,4 +421,181 @@
                       @[-3 3 -1]]))
         "permanent")
 
+(def test-primes
+  [100123456789
+   107928278317
+   113131311401
+   125411328001
+   137438691329
+   150614187107
+   168409140869
+   191198888863
+   232911191713
+   260389232731
+   313473008141
+   344980016453
+   411379717319
+   526858348381
+   581485876661
+   657835997711
+   701234567897
+   777775777777
+   977973373171
+   1000008000001
+   1110110110111
+   1146890986411
+   1344409044431
+   1534139560403
+   1686910196861
+   1918960000861
+   2748779069441
+   3484606209571
+   4615392897979
+   5599297703999
+   6987191424553
+   7777775552353
+   9977770001777
+   10861196119801
+   12345678987431
+   15154262241479
+   18826507658281
+   24738041398529
+   29202114663949
+   35557777777987
+   50000010000001
+   60818091990661
+   74596893730427
+   88929267773197
+   100033000330001
+   106111115118119
+   123123454321321
+   139239739439839
+   176860696068671
+   233444000011111
+   303160419086407
+   333555577577777
+   485398038695407
+   727777887889889
+   799333555511111
+   973369606963379
+   1003229774283941
+   1235711131175321
+   1510553619999637
+   1889080110806881
+   2357353373727757
+   3391382115599173
+   5111111111111119
+   7005264275346131
+   9007199254740847])
+
+(def pseudoprimes
+  [1373653
+   1530787
+   1987021
+   2284453
+   3116107
+   5173601
+   6787327
+   11541307
+   13694761
+   15978007
+   16070429
+   16879501
+   25326001
+   27509653
+   27664033
+   28527049
+   54029741
+   61832377
+   66096253
+   74927161
+   80375707
+   101649241
+   161304001
+   960946321
+   1157839381
+   3215031751
+   3697278427
+   5764643587
+   6770862367
+   14386156093
+   15579919981
+   18459366157
+   19887974881
+   21276028621
+   27716349961
+   29118033181
+   37131467521
+   41752650241
+   42550716781
+   43536545821
+   118670087467
+   307768373641
+   315962312077
+   354864744877
+   457453568161
+   528929554561
+   546348519181
+   602248359169
+   1362242655901
+   1871186716981
+   2152302898747
+   2273312197621
+   2366338900801
+   3343433905957
+   3461715915661
+   3474749660383
+   3477707481751
+   4341937413061
+   4777422165601
+   5537838510751])
+
+(defn trial-division
+  [n]
+  (cond
+    (< n 2) false
+    (= n 2) true
+    (= 0 (mod n 2)) false
+    (label result
+      (var p 3)
+      (while (<= (* p p) n)
+        (if (= 0 (mod n p)) (return result false))
+        (+= p 2))
+      true)))
+
+(def pg (take 10000 (primes)))
+(assert (all prime? pg) "small primes")
+(assert (= 104729 (last pg)) "10000th prime")
+
+(assert (all = (map next-prime pg) (drop 1 pg)) "next-prime")
+
+(assert
+  (all |(= (trial-division $) (prime? $)) (range 1 100000))
+  "prime? vs trial division")
+
+(assert
+  (all |(= (trial-division $) (not= 0 (jacobi $ 6541380665835015))) (range 49 2209 2))
+  "jacobi vs trial division")
+
+(assert (all prime? test-primes) "test primes")
+(when-let [s64 int/s64
+           u64 int/u64]
+  (assert (all prime? (map s64 test-primes)) "test primes s64")
+  (assert (all prime? (map u64 test-primes)) "test primes u64"))
+
+(assert (all |(miller-rabin-prp? $ 2 3) pseudoprimes) "pseudoprimes base 2 3")
+(assert (not (some miller-rabin-prp? pseudoprimes)) "pseudoprimes miller-rabin")
+
+(assert
+  (all |(= 1 (mulmod $ (invmod $ 100123456789) 100123456789)) (range 1 1000))
+  "invmod")
+
+(assert
+  (all |(= 1 (mulmod
+               (powmod $ $ 100123456789)
+               (powmod $ (- $) 100123456789)
+               100123456789))
+       (range 1 1000))
+  "powmod")
+
 (end-suite)