Divisors And Totatives

November 26, 2010

All of these functions rely on knowing the factors of the given number. You can use any method to calculate the factors; we will use trial division, knowing that it can be replaced by a more sophisticated method if needed (we assume below that factors returns a list in increasing order):

(define (factors n) (if (even? n) (cons 2 (factors (/ n 2))) (let loop ((n n) (f 3) (fs '())) (cond ((< n (* f f)) (reverse (cons n fs))) ((zero? (modulo n f)) (loop (/ n f) f (cons f fs))) (else (loop n (+ f 2) fs))))))

One way to find the divisors of a number is to take the product of all the subsets of its factors; the product of the null subset is 1. We use subsets from a previous exercise:

(define (divisors1 n) (define (prod xs) (apply * xs)) (define (cons1 xs) (cons 1 xs)) (unique = (sort < (map prod (map cons1 (subsets (factors n)))))))

Another way to find the divisors of a number is to iterate over its factors. Each factor can either be included or excluded, so append both possibilities to the list of divisors, recursively calculating the divisors of the remaining factors:

(define (divisors2 n) (define (times x) (lambda (y) (* x y))) (let divs ((fs (factors n))) (unique = (sort < (if (null? fs) '(1) (let ((ds (divs (cdr fs)))) (append ds (map (times (car fs)) ds))))))))

There is no better way to compute the sum of the divisors of a number than to compute its divisors and calculate their sum; we arbitrarily choose divisors1:

(define divisors divisors1)

(define (sumdiv n) (apply + (divisors n)))

The number of divisors can be computed as the length of the list of divisors, but it is possible to calculate the count without generating the divisors by looking at the exponents in the factorization of the number. For instance, 60 = 2 · 2 · 3 · 5 = 2² · 3¹ · 5¹; the exponents in the factorization are 2, 1 and 1. Add 1 to each exponent, giving 3, 2, and 2, and take their product: 3 · 2 · 2 = 12, which is the number of divisors of 60. We make the same calculation somewhat differently below:

(define (numdiv n) (let ((fs (factors n))) (let loop ((prev (car fs)) (fs (cdr fs)) (f 2) (d 1)) (cond ((null? fs) (* d f)) ((= (car fs) prev) (loop prev (cdr fs) (+ f 1) d)) (else (loop (car fs) (cdr fs) 2 (* d f)))))))

There is no other way to compute the totatives of a number than by iterating through the positive integers less than the number and testing each for coprimality:

(define (totatives n) (let loop ((t n) (ts '())) (cond ((= t 1) (cons 1 ts)) ((= (gcd t n) 1) (loop (- t 1) (cons t ts))) (else (loop (- t 1) ts)))))

The totient of a number is just the length of the list of totatives, but a better way to calculate the totient is similar to the method for calculating the number of divisors. We illustrate with 60 = 2 · 2 · 3 · 5 = 2² · 3¹ · 5¹. This time we look at the distinct factors 2, 3, and 5. Each is turned into a fraction using the formula 1 − ¹/_f, where f is the factor. Then the totient is n times each of the fractions, so φ(60) = 60 · ¹/₂ · ²/₃ · ⁴/₅ = 16. Our function performs this calculation:

(define (totient n) (let loop ((fs (unique = (factors n))) (t n)) (if (null? fs) t (loop (cdr fs) (* t (- 1 (/ (car fs))))))))

We used sort and unique from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/SmASe1eC.

Posted by programmingpraxis

Filed in Exercises

13 Comments »

13 Responses to “Divisors And Totatives”

Remco Niemeijer said

November 26, 2010 at 10:17 AM

My Haskell solution (see http://bonsaicode.wordpress.com/2010/11/26/programming-praxis-divisors-and-totatives/ for a version with comments):

import Data.List
import Data.Numbers.Primes
import Data.Ratio

divisors :: Integral a => a -> [a]
divisors = nub . sort . map product . subsequences . primeFactors

divisorSum :: Integral a => a -> a
divisorSum = sum . divisors

divisorCount :: Integral a => a -> Int
divisorCount = product . map (succ . length) . group . primeFactors

totatives :: Integral a => a -> [a]
totatives n = filter ((== 1) . gcd n) [1..n]

totient :: (Integral a, Integral b) => a -> b
totient n = floor . product $ n % 1 : [1 - 1 % f | f <- nub $ primeFactors n]

programmingpraxis said
November 26, 2010 at 3:04 PM
I just learned something about Scheme while reading Remco’s solution in Haskell. The result of subsequences . primeFactors is 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 the divisors1 function so it is the same as Remco’s version:

(define (divisors1 n) (define (prod xs) (apply * xs)) (unique = (sort < (map prod (subsets (factors n))))))
Remco Niemeijer said
November 26, 2010 at 7:33 PM
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.
programmingpraxis said
November 26, 2010 at 7:58 PM
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.
programmingpraxis said
November 26, 2010 at 8:11 PM
I think the Haskell equivalent of Scheme’s (apply * xs) must be foldl (*) 1 xs, not foldl1 (*) xs. (Or whatever order the operator and base case appear in.) Of course foldl1 (*) xs produces an error because foldl1 is not defined on the null list — it has nothing to do with the (*) function.
programmingpraxis said
November 27, 2010 at 12:10 AM
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))))))

Graham said

November 29, 2010 at 2:01 AM

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))

Khanh Nguyen said

November 29, 2010 at 9:01 PM

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.length

programmingpraxis said
December 6, 2010 at 4:28 PM
If 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 '()))))))))))
programmingpraxis said
July 14, 2011 at 8:52 PM
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)))))
Shankar Chakkere said
February 22, 2012 at 4:10 AM
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

David said

September 20, 2013 at 7:59 PM

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:

st> a := 7007 factored
Factorization (7^2 11 13 )
st> a totient
5040
st> a divisors
(1 7 11 13 49 77 91 143 539 637 1001 7007 )
st> a divisorCount
12
st> a divisorSigma
9576
st> 12 factored totatives
(1 5 7 11 )
st> 12 factored totient
4

Highly Composite Numbers, A Sieving Approach | Programming Praxis said
July 19, 2016 at 9:01 AM
[…] 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; […]

Programming Praxis