# HG changeset patch # User Steve Losh # Date 1509067119 14400 # Node ID 813d3e887794ed287f3c2380737c54bc6583d4ec # Parent cb34c7bc34fc57ef8b79a298fbac533cf74cd2f5 Clean up the divisor situation diff -r cb34c7bc34fc -r 813d3e887794 src/primes.lisp --- a/src/primes.lisp Fri Oct 13 01:32:55 2017 -0400 +++ b/src/primes.lisp Thu Oct 26 21:18:39 2017 -0400 @@ -76,10 +76,13 @@ (push n result)) (nreverse result))) -(defun prime-factorization-pairs (n) +(defun prime-factorization-pairs (n &key include-zeros) "Return the prime factorization of `n` as a list of prime/exponent conses. - The result will be a list of `(prime . exponent)` conses. For example: + The result will be a list of `(prime . exponent)` conses. If `include-zeros` + is given, even primes whose exponent turns out to be zero will be included. + + For example: (prime-factorization-pairs 60) ; => @@ -88,6 +91,8 @@ " (let ((result (list))) (iterate (while (evenp n)) ; handle 2, the only even prime factor + (else (when include-zeros + (push (cons 2 0) result))) (if-first-time (push (cons 2 1) result) (incf (cdar result))) @@ -95,6 +100,8 @@ (iterate (for i :from 3 :to (sqrt n) :by 2) ; handle odd (prime) divisors (iterate (while (dividesp n i)) + (else (when (and include-zeros (primep i)) + (push (cons i 0) result))) (if-first-time (push (cons i 1) result) (incf (cdar result))) diff -r cb34c7bc34fc -r 813d3e887794 src/problems.lisp --- a/src/problems.lisp Fri Oct 13 01:32:55 2017 -0400 +++ b/src/problems.lisp Thu Oct 26 21:18:39 2017 -0400 @@ -452,12 +452,10 @@ ;; 71 and 142; so d(284) = 220. ;; ;; Evaluate the sum of all the amicable numbers under 10000. - (labels ((sum-of-divisors (n) - (sum (proper-divisors n))) - (amicablep (n) - (let ((other (sum-of-divisors n))) - (and (not= n other) - (= n (sum-of-divisors other)))))) + (flet ((amicablep (n) + (let ((other (sum-of-divisors n :proper t))) + (and (not= n other) + (= n (sum-of-divisors other :proper t)))))) (sum (remove-if-not #'amicablep (range 1 10000))))) (defun problem-22 () diff -r cb34c7bc34fc -r 813d3e887794 src/utils.lisp --- a/src/utils.lisp Fri Oct 13 01:32:55 2017 -0400 +++ b/src/utils.lisp Thu Oct 26 21:18:39 2017 -0400 @@ -269,27 +269,25 @@ (sort sequence #'<)) -(defun unsorted-divisors (n) - (iterate (for i :from 1 :to (sqrt n)) - (for (values divisor remainder) = (truncate n i)) - (when (zerop remainder) - (collect i) - ;; don't collect the square root twice - (unless (= i divisor) - (collect divisor))))) - (defun-inlineable divisors-up-to-square-root (n) (loop :for i :from 1 :to (floor (sqrt n)) :when (zerop (rem n i)) :collect i)) -(defun divisors (n) - (sort< (unsorted-divisors n))) -(defun proper-divisors (n) - (remove n (divisors n))) +(defun unsorted-divisors (n &key proper) + (iterate (for i :from 1 :to (sqrt n)) + (for (values divisor remainder) = (truncate n i)) + (when (zerop remainder) + (collect i) + (unless (or (= i divisor) ;; don't collect the square root twice + (and proper (= i 1))) ;; don't collect n if proper + (collect divisor))))) -(defun count-divisors (n) +(defun divisors (n &key proper) + (sort< (unsorted-divisors n :proper proper))) + +(defun count-divisors (n &key proper) ;; From *Recreations in the Theory of Numbers: The Queen of Mathematics ;; Entertains* by Albert Beiler. ;; @@ -302,8 +300,29 @@ ;; factorization. This is because there are k₁+1 options for how many of the ;; first prime to include (because 0 is an option), k₂+1 options for the ;; second prime, etc. - (iterate (for (nil . exponent) :in (prime-factorization-pairs n)) - (multiplying (1+ exponent)))) + (- (iterate (for (nil . exponent) :in (prime-factorization-pairs n)) + (multiplying (1+ exponent))) + (if proper 1 0))) + +(defun sum-of-divisors (n &key proper) + "Return the sum of the divisors of `n`. + + If `proper` is given, only include the proper divisors (i.e. not `n` itself). + + " + (sum (divisors n :proper proper))) + +(defun product-of-divisors (n &key proper) + ;; From *Recreations in the Theory of Numbers: The Queen of Mathematics + ;; Entertains* by Albert Beiler. + ;; + ;; If `N` is a positive integer with `d` divisors, the product of all the + ;; divisors of `N` is `N^(d/2)`. + ;; + ;; If we only want the product of the proper divisors, just divide `N` back + ;; out at the end. + (/ (expt n (/ (count-divisors n) 2)) + (if proper n 1))) (defmacro-driver (FOR var IN-COLLATZ n) @@ -370,13 +389,13 @@ (defun perfectp (n) - (= n (sum (proper-divisors n)))) + (= n (sum-of-divisors n :proper t))) (defun abundantp (n) - (< n (sum (proper-divisors n)))) + (< n (sum-of-divisors n :proper t))) (defun deficientp (n) - (> n (sum (proper-divisors n)))) + (> n (sum-of-divisors n :proper t))) (defun multiplicative-order (integer modulus)