The Divisor Function Sigma

January 10, 2014

The divisor function first factors n, by any convenient method, then sorts the factors and applies the uniq-c function of the Standard Prelude, which collects the factors and their multiplicities; for instance, (uniq-c = (sort < (factors 60))) returns the list ((2 . 2) (3 . 1) (5 . 1)). Then the appropriate definition of the divisor function is applied, using map instead of explicit loops:

(define (divisor x n)
  (define (prod xs) (apply * xs))
  (if (= n 1) 1
    (let ((fs (uniq-c = (sort < (factors n)))))
      (if (zero? x)
          (prod (map add1 (map cdr fs)))
          (prod (map (lambda (p a)
                       (/ (- (expt p (* (+ a 1) x)) 1)
                          (- (expt p x) 1)))
                     (map car fs) (map cdr fs)))))))

The summatory function calculates the number of times each divisor appears and applies the appropriate powering operation:

(define (summatory-divisor x n)
  (do ((d 1 (+ d 1))
       (s 0 (+ s (* (expt d x) (quotient n d)))))
      ((< n d) s)))

It may help to see visually what is happening. Consider what happens when n = 10; the divisors of the numbers from 1 to 10 are:

1: 1
2: 1 2
3: 1 3
4: 1 2 4
5: 1 5
6: 1 2 3 6
7: 1 7
8: 1 2 4 8
9: 1 3 9
10: 1 2 5 10

If you rearrange the divisors, you can see how many times each divisor appears:

1:  (10) 1 2 3 4 5 6 7 8 9 10
2:  (5)  2 4 6 8 10
3:  (3)  3 6 9
4:  (2)  4 8
5:  (2)  5 10
6:  (1)  6
7:  (1)  7
8:  (1)  8
9:  (1)  9
10: (1) 10

Thus, each divisor d appears (quotient n d) times, which the function then raises to the appropriate power and sums over all divisors. Here are some examples:

> (divisor 0 60)
12
> (divisor 1 60)
168
> (divisor 2 60)
5460
> (summatory-divisor 0 60)
261
> (summatory-divisor 1 60)
3014
> (summatory-divisor 2 60)
89286

You can run the program at http://programmingpraxis.codepad.org/gpTcbMQr.

Today’s exercise was inspired by Project Euler 401, which you are now equipped to solve. But if you want to solve it within the one-minute rule, you will need to improve the summatory-divisor function, because even though it operates in linear time, doing a quadrillion of anything will take a while. Out of respect for Project Euler, I won’t give the code, but the idea is to combine adjacent divisors with identical counts and compute the sum of their squares using the formula for square pyramidal numbers.

About these ads

Pages: 1 2

One Response to “The Divisor Function Sigma”

  1. Paul said

    In Python.

    from collections import Counter
    from ma.primee import td_factors
    
    def divisor_sigma(n, x):
        res = 1
        if x  == 0 :
            for pi, ai in Counter(td_factors(n)).iteritems():
                if ai:
                    res *= (ai + 1)
        else:
            for pi, ai in Counter(td_factors(n)).iteritems():
                if ai:
                    res *= (pi ** ((ai + 1) * x) - 1) // (pi ** x - 1)
        return res
        
    def sum_divisor_sigma(n, x):
        if x == 0:
            return sum(n // d for d in xrange(1, n + 1))
        else:
            return sum(n // d * d ** x for d in xrange(1, n + 1))
    

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

%d bloggers like this: