# 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:

sum(n=1,64*10^6,issquare(sigma(n,2))*n)


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

• 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
Feb 15, 2013 at 4:44
• @kale that is much faster but the output doesn't appear to be correct. Feb 15, 2013 at 5:37
• @Mr.Wizard, Is that a machine precision thing or a mistake in my formula?
– kale
Feb 15, 2013 at 14:22
• 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. Feb 15, 2013 at 17:23
• @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. Feb 15, 2013 at 18:44

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
First@a/First@b

(*True*)
19.9746841

• You ran the full thing on your machine and the results agree? I guess I need to try this again. Feb 16, 2013 at 2:46
• @Mr.Wizard, Yessir.
– kale
Feb 16, 2013 at 2:53
• 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. Feb 16, 2013 at 3:06
• @Mr.Wizard I've verified the correct answer with Project Euler. Bizarre. I'm on v9, BTW.
– kale
Feb 16, 2013 at 4:06
• @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]]] Mar 7, 2013 at 5:13

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.

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


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