From b89bd9f8896626de1af4a998ce9d08bd18faa54c Mon Sep 17 00:00:00 2001
From: Roger Crew <crew@cs.stanford.edu>
Date: Sun, 27 Jun 2021 03:58:58 -0700
Subject: [PATCH] faster factor, (factor 1) = () (issue #751 cont.)

no need to go up to sqrt(n), Instead track i^2 and quit when that gets
larger than the (remaining) n (i.e., not the original n)
---
 lib/chibi/math/prime-test.sld |  3 ++-
 lib/chibi/math/prime.scm      | 42 +++++++++++++++--------------------
 2 files changed, 20 insertions(+), 25 deletions(-)

diff --git a/lib/chibi/math/prime-test.sld b/lib/chibi/math/prime-test.sld
index b006439a..c4a3ae0a 100644
--- a/lib/chibi/math/prime-test.sld
+++ b/lib/chibi/math/prime-test.sld
@@ -71,7 +71,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,6 +86,7 @@
       (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 975 (aliquot 945))
 
diff --git a/lib/chibi/math/prime.scm b/lib/chibi/math/prime.scm
index 54184e39..4ec5a061 100644
--- a/lib/chibi/math/prime.scm
+++ b/lib/chibi/math/prime.scm
@@ -192,30 +192,24 @@
 ;;> Returns the factorization of \var{n} as a monotonically
 ;;> increasing list of primes.
 (define (factor n)
-  (cond
-   ((negative? n)
-    (cons -1 (factor (- 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))))))))))
+  (when (zero? n)
+    (error "cannot factor 0"))
+  (factor-twos
+   n
+   (lambda (b n)
+     (let lp ((i 3) (ii 9) (n (abs n))
+              (res (append (make-list b 2)
+                           (if (negative? n) (list -1) '()))))
+       (cond
+        ((> ii n)
+         (reverse (if (= n 1) res (cons n res))))
+        ((zero? (remainder n i))
+         (lp i ii (quotient n i) (cons i res)))
+        (else
+         (lp (+ i 2)
+             (+ ii (* (+ i 1) 4))
+             n
+             res)))))))
 
 ;;> Returns the Euler totient function, the number of positive
 ;;> integers less than \var{n} that are relatively prime to \var{n}.