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 = 22 · 31 · 51; 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 = 22 · 31 · 51. This time we look at the distinct factors 2, 3, and 5. Each is turned into a fraction using the formula 1 − 1/f, where f is the factor. Then the totient is n times each of the fractions, so φ(60) = 60 · 1/2 · 2/3 · 4/5 = 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.

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 616 other followers

%d bloggers like this: