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.

About these ads

Pages: 1 2

12 Responses to “Divisors And Totatives”

  1. 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]
    
  2. programmingpraxis said

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

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

  4. programmingpraxis said

    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.

  5. programmingpraxis said

    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.

  6. programmingpraxis said

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

  7. Graham said

    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))
    
  8. Khanh Nguyen said

    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
    
  9. programmingpraxis said

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

  10. programmingpraxis said

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

  11. Shankar Chakkere said

    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

  12. David said

    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
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 615 other followers

%d bloggers like this: