Calculate the sum of divisors

Let's say you have an arbitrary number \(n\in \mathbb{N}\). If you want to calculate the sum of it's divisors, you would come up with a loop which is \(O(n)\):

var sigma(n) {
  var s = 0;
  for (var i = 1; i <= n; i++) {
    if (n % i==0)
      s+= i;
    }
    return s;
}

Here is how you can make it much faster.Let \(\delta(n)\) be the number of divisors and \(\sigma(n)\) the sum of the divisors of number \(n\). For any prime \(p\), it's obvious that \(\sigma(p)=p+1\), as it's only divisors can be \(1\) and \(p\). If you go ahead, it's also obvious that \(\delta(p^a)=a+1\), and:

\[\begin{array}{rl} \sigma(p^a) &= \sum\limits_{i=0}^{\delta(p^a)}p^i \\ &= 1 + p + p^2 + ... + p^a \\ p\sigma(p^a) &= p + p^2 + p^3 + ... + p^{a+1} \\ p\sigma(p^a) - \sigma(p^a) &= p + p^2 + p^3 + ... + p^{a+1} - (1 + p + p^2 + ... + p^a) \\ (p - 1) \sigma(p^a) &= p^{a+1}-1 \\ \sigma(p^a) &= \frac{p^{a+1}-1}{p-1} \end{array}\]

We now use the multiplicativity property of \(\sigma\), which means that \(\sigma(x\cdot y\cdot ...)=\sigma(x)\cdot\sigma(y)\cdot...\) with \(x, y, ...\) relatively prime and can say that

\[\begin{array}{rl} \sigma(n) &= \sigma(p_1^{a_1} \cdot p_2^{a_2} \cdot ...) \\ &= \sigma(p_1^{a_1}) \cdot \sigma(p_2^{a_2}) \cdot ...\\ &= \frac{p_1^{a_1+1}-1}{p_1-1} \cdot \frac{p_2^{a_2+1}-1}{p_2-1} \cdot... \\ &= \prod\limits_{p_i \in \mathbb{P}} \frac{p_i^{a_i+1}-1}{p_i-1} \end{array}\]

This loop goes over all primes, which isn't needed though. The largest prime to be checked actually is only the square root of the number, which ultimately results in the following algorithm, if a set of primes is already known:

function sigma(n) {
  var sum = 1;
  var primes = [2, 3, 5, 7, ...];

  if (n < 4) {
    return 1;
  }

  for (var i = 0; i < primes.length; i++) {

    var p = primes[i];

    if (0 === n % p) {

      var t = p * p;
      n /= p;
      while (0 === n % p) {
        t *= p;
        n /= p;
      }
      sum = sum * (t - 1) / (p - 1);
    }
    if (p * p > n) {
      break;
    }
  }

  if (n > 1) {
    sum*= n + 1;
  }
  return sum;
}

This is practically useful only with a prime factor table available. A primal search can be combined with the summation of the factors, which leads to the following derivation:

\[\begin{array}{rl} \sigma(n) &= n+\sum\limits_{i=1}^n \begin{cases} i & \text{if } i | n\\ 0 & \text{else} \end{cases}\\ &= 1+n+\sum\limits_{i=2}^{\sqrt{n}} \begin{cases} i + n / i & \text{if } i | n\\ 0 & \text{else} \end{cases}\\ \end{array} \]

For a general case, the square root of n must be subtracted from the result if \(n\) is a perfect square. Implemented in JavaScript this can then look like:

function sigma(n) {

  var sum = 1 + n;
  var sqrt = Math.sqrt(n);

  for (var i = 2; i <= sqrt; i++) {

    if (0 === n % i) {
      sum+= i + n / i;
    }
  }

  if (0 === sqrt % 1) {
    sum-= sqrt;
  }
  return sum;
}

This algorithm isn't as good as the previous algorithm. So generating a prime factor table with a sieve beforehand would be the best option.

Calculate sum of divisors from 1 to n

If we go even further and wonder what \(\Delta(n) = \sum\limits_{i=1}^n \sigma(i)\) is, we can do much better, without factorization. Let's say that \(\phi(i,j) = \begin{cases} 1, & \text{if}\ i|j \\ 0, & \text{else} \end{cases}\), then

\[\begin{array}{rl} \Delta(n) &= \sum\limits_{i=1}^n \sigma(i)\\ &= \sum\limits_{j=1}^n \sum\limits_{i=1}^n i \cdot \phi(i,j)\\ &= \sum\limits_{i=1}^n \sum\limits_{j=1}^n i \cdot \phi(i,j)\\ &= \sum\limits_{i=1}^n i \cdot \sum\limits_{j=1}^n \phi(i,j)\\ &= \sum\limits_{i=1}^n i \lfloor\frac{n}{i} \rfloor \\ &= n^2 - \sum\limits_{i=1}^n mod(n, i) \end{array}\]