totient and aliquot rewrite + corrected tests for n=1 (issue ashinn#751…

… cont.)

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))
@@ -87,6 +88,7 @@
       (test '(2 2 2 3 3) (factor 72))
       (test '(3 3 3 5 7) (factor 945))
 
+      (test 0 (aliquot 1))
       (test 975 (aliquot 945))
 
       (do ((i 3 (+ i 2)))

lib/chibi/math/prime.scm

@@ -189,6 +189,48 @@
          ((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 arguments \var{put2} and \var{put-1}, if provided,
+;; are the specialized versions of \var{put} that get called
+;; for \var{p}=2 (as \scheme{(\var{put2} \var{r} \var{k})})
+;; and \var{p}=-1 (as \scheme{(\var{put-1} \var{r})}), respectively.
+(define make-factorizer
+  (case-lambda
+   ((r1 put)
+    (make-factorizer
+     r1 put
+     (lambda (r k) (put r  2 k))
+     (lambda (r)   (put r -1 1))))
+   ((r1 put put2 put-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 monotonically
 ;;> increasing list of primes.
 (define (factor n)
@@ -210,30 +252,38 @@
 	     (+ ii (* (+ i 1) 4))
 	     n
 	     res)))))))
+;; this version is slightly slower
+;;(define factor
+;;  (let ((rfactor (make-factorizer '()
+;;                  (lambda (l p k) (cons (make-list k p) l)))))
+;;    (lambda (n) (apply append! (rfactor n)))))
+
+
+;;> The Euler totient φ(\var{n}) is the number of positive
+;;> integers less than or equal to \var{n} that are
+;;> relatively prime to \var{n}.
+(define totient
+  (make-factorizer 1
+   (lambda (tot p k)
+     (* tot (- p 1) (expt p (- k 1))))
+   (lambda (tot k)
+     (arithmetic-shift tot (- k 1)))
+   (lambda (_)
+     (error "totient of negative number?"))))
 
-;;> Returns the Euler totient function, the number of positive
-;;> integers less than \var{n} that are relatively prime to \var{n}.
-(define (totient n)
-  (let ((limit (exact (ceiling (sqrt n)))))
-    (let lp ((i 2) (count 1))
-      (cond ((> i limit)
-             (if (= count (- i 1))
-                 (- n 1)                ; shortcut for prime
-                 (let lp ((i i) (count count))
-                   (cond ((>= i n) count)
-                         ((= 1 (gcd n i)) (lp (+ i 1) (+ count 1)))
-                         (else (lp (+ i 1) count))))))
-            ((= 1 (gcd n i)) (lp (+ i 1) (+ count 1)))
-            (else (lp (+ i 1) count))))))
+;;> The aliquot sum s(\var{n}) is
+;;> 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,6 +1,6 @@
 
 (define-library (chibi math prime)
-  (import (scheme base) (scheme inexact) (srfi 27))
+  (import (scheme base) (scheme inexact) (scheme case-lambda) (srfi 27))
   (cond-expand
    ((library (srfi 151)) (import (srfi 151)))
    ((library (srfi 33)) (import (srfi 33)))