How can this DivisorSigma code be made fast?


Since Project Euler problems are now fair game for questions I have a question of my own.
A certain problem* states:

For a positive integer n, let σ2(n) be the sum of the squares of its divisors. For example,

σ2(10) = 1 + 4 + 25 + 100 = 130.

Find the sum of all n, 0 < n < 64,000,000 such that σ2(n) is a perfect square.

This Mathematica code takes something like an hour to run on a modern machine:

Sum[If[IntegerQ @ Sqrt @ DivisorSigma[2, i], i, 0], {i, 64*^6 - 1}] ~Monitor~ i // Timing

The similarly naive PARI/GP code takes a minute or two:


Is there some way to make the Mathematica code fast, or otherwise solve the problem quickly in Mathematica?

Using a faster perfect square test helps quite a bit but it is still far from the PARI/GP performance.

Compilation does not seem possible as numbers exceed the maximum machine-size integer.

(*To foil search engines please do not mention the number of the Project Euler problem related to this question. Thanks.)


Posted 2013-02-15T02:17:21.900

Reputation: 259 163

I know this probably isn't in the spirit of the problem (since they're technically not integers) but I'm having success with ParallelSum[If[Mod[DivisorSigma[2., i]^0.5, 1] == 0, i, 0.],{i, 64000000.}]. On my 2010 Core 2 Macbook (and on battery power) it takes ~7 min. – kale – 2013-02-15T04:44:29.307

@kale that is much faster but the output doesn't appear to be correct. – Mr.Wizard – 2013-02-15T05:37:13.637

@Mr.Wizard, Is that a machine precision thing or a mistake in my formula? – kale – 2013-02-15T14:22:50.917


Testing all $n=2^{26}$ values is at least a $O(n)$ algorithm, so any such approach is just trying to improve an implicit constant. (It reminds me of Churchill's take-down of a lady: "we already know what you are; we're merely haggling about the price.") If you want to achieve real gains, you need to think of ways to rule out the majority of those tests a priori. That suggests exploiting the multiplicative structure of $\sigma_2$ and building a solution based on its values at prime powers: that holds out the hope of a $O(\log(n))$ solution.

– whuber – 2013-02-15T17:23:45.250

@whuber I should remove the "without involved mathematical reasoning" line from my question because I would like to see you demonstrate such an approach, even if I have trouble following it. That is, if you have the time and interest of course. – Mr.Wizard – 2013-02-15T18:44:42.970

@MichaelE2 I removed that line from the question because I am interested in seeing better algorithms and perhaps Mathematica cannot be made faster here otherwise. – Mr.Wizard – 2013-02-15T18:46:43.863

Relevant oeis. The formula sigma_k(n) = Product_p ((p^((e(p)+1)k))-1)/(p^k-1)* seems like a good place to start reducing the the computational cost since you can reuse large parts – ssch – 2013-02-15T19:40:42.613

I have found an intriguing approach based on partially factoring $\sigma_2(q)$ for all primes and prime powers up to a suitable size. The idea is that for $n=p_1^{e_1}\cdots p_k^{e_k}$, $\sigma_2(n)$ = $\sigma_2(p_1^{e_1})\cdots\sigma_2(p_k^{e_k})$. Writing $\sigma_2(p_i^e)$ = $q_1^{f_{i1}}\cdots q_j^{f_{ij}}$, we find that $\sigma_2(n)$ is a square iff the sums over $i$ of the $f_{ij}$ are zero modulo $2$: the solutions can be found, e.g., with NullSpace using sparse matrices. It's still a delicate matter, though, to improve on $O(n)$ performance--perhaps $O(n/\log(n)^2)$ is possible. – whuber – 2013-02-16T02:37:56.173

@Mr.Wizard Deleted my comment; will delete this one shortly. – Michael E2 – 2013-02-16T15:32:08.107



Here is my attempt. A 20X speedup on my computer...

a = Sum[If[IntegerQ@Sqrt@DivisorSigma[2, i], i, 0], {i, 64*^6 - 1}] //Timing (*Mr.Wizard*)
b = Rationalize[ParallelSum[If[Mod[DivisorSigma[2., i]^0.5, 1] == 0, i, 0.],
  {i, 64*^6 - 1}]] // AbsoluteTiming (*Proposed method*)
Last@a == Last@b



Posted 2013-02-15T02:17:21.900

Reputation: 10 290

Today I confirmed that this yields the correct result in version 10.1.0. I guess machine precision changes make the difference. A belated +1 and Accept. – Mr.Wizard – 2015-06-18T06:07:05.780

@Mr.Wizard Firing up the "Way-Back machine"! – kale – 2015-06-19T01:19:47.907

You ran the full thing on your machine and the results agree? I guess I need to try this again. – Mr.Wizard – 2013-02-16T02:46:05.653

@Mr.Wizard, Yessir. – kale – 2013-02-16T02:53:41.953

I really want this to work because it's a wonderful speed-up but I get a result of 167020264273 which doesn't even resemble the correct answer. It seems to match for small values but apparently on my system something goes wrong with larger numbers. I'll continue to explore this. – Mr.Wizard – 2013-02-16T03:06:11.383

@Mr.Wizard I've verified the correct answer with Project Euler. Bizarre. I'm on v9, BTW. – kale – 2013-02-16T04:06:27.777

@Mr.Wizard Building on the fast solution by @kale, the following gives another 40% speed-up Total[ParallelCombine[Pick[#, Mod[DivisorSigma[2., #]^0.5, 1], 0.\] &, Range[64*10^6-1]]]` – KennyColnago – 2013-03-07T05:13:00.077


I found testing for an odd number of divisors can be more efficient sometimes. Try

Sum[If[OddQ[DivisorSigma[0,DivisorSigma[2, i]]],i,0],{i,10^5-1}]

The following is 2.5 times faster than your method.


If you have more processors, there is a way of using ParallelCombine as in the following.

Total[ParallelCombine[Pick[#,OddQ[DivisorSigma[0,DivisorSigma[2, #]]]]&,Range[10^5-1]]]]

Of course with "involved" mathematical reasoning our results would improve...


Posted 2013-02-15T02:17:21.900

Reputation: 14 269