add math/factor

spork/math.janet

@@ -1295,3 +1295,61 @@
             (put sieve (+ q (* p (in soe-table (mod p 105)))) p)
             (set p (resume pg))
             (set q (* p p))))))))
+
+(defn- pollard-rho
+  "pollard's rho algorithm, with brent's cycle detection"
+  [n &opt c]
+  (default c 1)
+  (def m 16)
+  (var [x q r] [2 1 m])
+  (var xs nil)
+  (var y nil)
+  (while (not= 0 (jacobi q n))
+    (set y x)
+    (repeat r
+      (set x (+ (mulmod x x n) c)))
+    (loop [k :range [0 r m]
+           :while (not= 0 (jacobi q n))
+           :before (set xs x)
+           :repeat m]
+      (set x (+ (mulmod x x n) c))
+      (set q (mulmod q (- x y) n)))
+    (*= r 2))
+  (if (zero? q)
+    (while true
+      (set xs (+ (mulmod xs xs n) c))
+      (set q (mod (- xs y) n))
+      (if (= 0 (jacobi q n)) (break))))
+  (if (zero? q)
+    (pollard-rho n (+ c 1))
+    (do
+      (set x n)
+      (while (not (zero? q))
+        (set* [x q] [q (mod x q)]))
+      x)))
+
+(defn- factor-pollard
+  [n]
+  (if (miller-rabin-prp? n)
+    [n]
+    (do
+      (def p (pollard-rho n))
+      (merge-sorted
+        (factor-pollard p)
+        (factor-pollard (/ n p))))))
+
+(defn factor
+  "Returns an array containing the prime factors of `n`."
+  [n]
+  (def res @[])
+  (def cmp-fn (get n :compare cmp))
+  (def ident (- n n))
+  (when (>= (cmp-fn n 2) 0)
+    (var x n)
+    (each p small-primes
+      (while (= (cmp-fn (mod x p) 0) 0)
+        (array/push res (+ p ident))
+        (/= x p)))
+    (if (> (cmp-fn x 1) 0)
+      (array/concat res (factor-pollard x))))
+  res)

test/suite-math.janet

@@ -597,4 +597,15 @@
     (all |(= 1 (mulmod (powmod $ $ p) (powmod $ (- $) p) p)) (range 1 1000))
     (string "powmod " p)))
 
+(defn check-factor [n]
+  (def res (factor n))
+  (and 
+    (all prime? res)
+    (all = res (sorted res))
+    (= n (* ;res))))
+
+(assert (all check-factor (range 1 10000)) "factor small integers")
+(assert (all check-factor test-primes) "factor primes")
+(assert (all check-factor pseudoprimes) "factor pseudoprimes")
+
 (end-suite)