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