Divisors And Totatives
November 26, 2010
We have examined functions to compute the factors of a number in several previous exercises. In today’s exercise we examine functions to compute divisors and totients, which are concepts of number theory closely related to factors.
The divisors of a number are those numbers that divide it evenly; for example, the divisors of 60 are 1, 2, 3, 4, 5, 6, 10, 12, 15, 20, 30, and 60. The sum of the divisors of 60 is 168, and the number of divisors of 60 is 12.
The totatives of a number are those numbers less than the given number and coprime to it; two numbers are coprime if they have no common factors other than 1. The number of totatives of a given number is called its totient. For example, the totatives of 30 are 1, 7, 11, 13, 17, 19, 23, and 29, and the totient of 30 is 8. The totient function was discovered by Leonhard Euler, and it pops up in all kinds of strange places in number theory; for instance, the totatives of 30 are the spokes on a 2,3,5 factorization wheel.
Your task is to write a small library of five functions that compute the divisors of a number, the sum and number of its divisors, the totatives of a number, and its totient. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
My Haskell solution (see http://bonsaicode.wordpress.com/2010/11/26/programming-praxis-divisors-and-totatives/ for a version with comments):
I just learned something about Scheme while reading Remco’s solution in Haskell. The result of
subsequences . primeFactorsis the same in Haskell as the expression(subsets (factors n))in Scheme, and includes the null set as one of the subsets. I thought that mapping the product function over that list would produce an error when applied to the null set, but apparently it doesn’t in Haskell. So I checked Scheme, and found out that(apply * '())doesn’t produce an error, as I thought it did, but instead returns 1. That means I can simplify thedivisors1function so it is the same as Remco’s version:(define (divisors1 n)(define (prod xs) (apply * xs))
(unique = (sort < (map prod (subsets (factors n))))))
Now I’m curious: how does Scheme know the value for the base case? In Haskell, product is defined as foldl (*) 1. The equivalent of (apply * xs), which would appear to be foldl1 (*) xs, produces an error when given an empty list. I looked at your Scheme Prelude, but couldn’t find a definition of apply.
I didn’t know this either until I looked it up. I thought (*) was an error.
The magic is in the * operator, not the apply function. The product of an empty list is 1, by definition. See Section 6.2.5 of R5RS. You can give any number of arguments to the * function, where any number includes zero. If you say (* 2 3 4) the answer is 24. If you say (* 2 3) the answer is 6. If you say (* 2) the answer is 2. And if you say (*) the answer is 1. Likewise, addition with no arguments (+) is 0.
Apply is defined in Section 6.4 of R5RS, not the Standard Prelude. An operational definition of apply, using denotational semantics, is in Section 7.2.4 of R5RS; I won’t try to re-type all the funny characters and Greek letters, so you’ll have to go look at it. Apply is the often-forgotten half of the eval-apply yin-yang.
It is amazing how often Scheme just “does the right thing” with your program. It was truly designed and refined by experts.
I think the Haskell equivalent of Scheme’s
(apply * xs)must befoldl (*) 1 xs, notfoldl1 (*) xs. (Or whatever order the operator and base case appear in.) Of coursefoldl1 (*) xsproduces an error becausefoldl1is not defined on the null list — it has nothing to do with the(*)function.I wrote a set of one-liners equivalent to Remco’s version. You can run the program at http://programmingpraxis.codepad.org/poR5avz2. Follow the link to see all the prelude functions that are required; here is the meat of the library:
(define (divisors n)(unique = (sort < (map product (subsets (factors n))))))
(define (sumdiv n) (sum (divisors n)))
(define (numdiv n)
(apply * (map add1 (map cdr (uniq-c = (factors n))))))
(define (totatives n)
(filter (lambda (x) (= (gcd n x) 1)) (range 1 n)))
(define (totient n)
(product (cons n (map (lambda (x) (- 1 (/ x))) (unique = (factors n))))))
Pretty busy this weekend, so I went for the obvious (and slow) answers in Python, mostly written with clarity in mind.
def gcd(a, b): while b: a, b = b, a % b return a def divides(k, n): if n % k == 0: return True else: return False def coprime(k, n): if gcd(k, n) == 1: return True else: return False def divisors(n): return [k for k in xrange(1, n + 1) if divides(k, n)] def sum_divs(n): return sum(divisors(n)) def num_divs(n): return len(divisors(n)) def totatives(n): return [k for k in xrange(1, n) if coprime(k, n)] def totient(n): return len(totatives(n))In F#..
let divisors n = List.filter (fun x -> n % x = 0) [1..n] let sum_of_divisors n = divisors n |> List.sum let number_of_divisors n = divisors n |> List.length let totatives n = let coprime a b = Set.intersect (divisors a |> Set.ofList) (divisors b |> Set.ofList) |> Set.count = 1 List.filter (fun x -> coprime n x ) [1..n] let totient n = totatives n |> List.lengthIf you want to factor numbers a little bit larger than trial division can comfortably handle, here is a version of John Pollard’s rho algorithm that is short, simple, and quick; it can find factors up to ten or twelve digits without inconvenience:
(define (factors n)(let fact ((n n) (c 1) (fs '()))
(define (rho x c) (modulo (+ (* x x) c) n))
(cond ((prime? n) (sort < (cons n fs)))
((even? n) (fact (/ n 2) 1 (cons 2 fs)))
(else (let loop ((t 2) (h 2) (f 1))
(cond ((= f 1)
(let ((t (rho t c)) (h (rho (rho h c) c)))
(loop t h (gcd (- t h) n))))
((= f n) (fact n (+ c 1) fs))
((prime? f) (fact (/ n f) c (cons f fs)))
(else (append fs (fact f c '())
(fact (/ n f) c '()))))))))))
I was wrong to say that there is no simpler way to compute the sum of the divisors of the number than to enumerate the divisors and compute their sum. This function is due to Berndt (see Eq 14 at http://mathworld.wolfram.com/DivisorFunction.html) and requires only the factors of n, not the list of divisors:
(define (sumdiv n)(define (f x) (/ (- (expt (car x) (+ (cdr x) 1)) 1) (- (car x) 1)))
(product (map f (uniq-c = (factors n)))))
Here is my take in APPLE 1 BASIC
#!/usr/bin/apple1basic
100 PRINT “ENTER THE NUMBER”
200 INPUT N
300 REM * HOLDS SUM OF DIVISORS
310 LET S = 0
400 REM * HOLDS NUMBER OF DIVISORS
410 LET D = 0
505 FOR I = 1 TO N
510 REM IF NOT EXACT DIVISIBLE NEXT NUMBER
600 IF (N/I)*I N THEN NEXT I
610 PRINT I
620 S = S + I: D = D + 1
700 NEXT I
800 PRINT “SUM = “;S,”NUMBER OF DIVISORS = “;D
900 END
For SmallTalk, I created an intermediate object (a Factorization object), the idea being that one may wish to work call several methods (e.g. totient, divisors, etc.) without having a need to constantly refactor, just save the Factorization object if you need to resuse it… For the divisors method, I take subsets by noting that there are 2^n subsets in a set of n objects, and they can be enumerated just by counting from 0 .. 2^n-1. For each number I do bit twiddling to select all subsets. (e.g. get the lowbit, which is first element to count, mask out the lowbit, keep going until zero.)
Object subclass: Factorization [ |factors n| factorize: aNumber [ |p m exponent| factors := Bag new. n := aNumber. m := n abs. p := 2. [p squared > m] whileFalse: [ (m rem: p) = 0 ifTrue: [ m := m / p. factors add: p. ] ifFalse: [ p := p nextPrime. ]. ]. factors add: m. ] asInteger [ ^n. ] printOn: aStream [ <category: 'printing'> aStream nextPutAll: self class storeString , ' ('. n < 0 ifTrue: [ aStream print: -1; space ]. factors asSet sorted do: [:fac| |exp| exp := factors occurrencesOf: fac. exp = 1 ifTrue: [ aStream print: fac ] ifFalse: [ aStream print: fac; nextPut: $^; print: exp]. aStream space. ]. aStream nextPut: $). ] divisors [ |mults divs val| mults := factors asArray. divs := Set new. 0 to: (2 raisedTo: mults size) - 1 do: [:m| |n| val := 1. n := m. [ n ~= 0 ] whileTrue: [ val := val * (mults at: n lowBit). n := n bitAnd: n-1. ]. divs add: val. ]. ^divs sorted. ] divisorCount [ ^ factors asSet inject: 1 into: [:a :b| a * ((factors occurrencesOf: b) + 1)]. ] divisorSigma [ ^ factors asSet inject: 1 into: [:x :p| |e| e := factors occurrencesOf: p. ((p raisedTo: e+1) - 1) / (p - 1) * x]. ] totient [ n <= 0 ifTrue: [^0]. n = 1 ifTrue: [^1]. ^ n * (factors asSet inject: 1 into: [:x :p| x * (1 - (1/p))]). ] totatives [ ^(1 to: n abs - 1) select: [:m| (n gcd: m) = 1]. ] ] Factorization class extend [ new: n [ ^ (super new) factorize: n. ] ] Integer extend [ factored [ ^ Factorization new: self ] ]Examples of usage:
[…] first solution to this problem, or you can factor n and compute the divisors, as we have done in a previous exercise. But if you have to find the divisors of a bunch of numbers, in sequence, you can sieve for them; […]