5

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 – 2013-02-15T04:44:29.307@kale that

ismuch 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

1

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

– whuber – 2013-02-15T17:23:45.250realgains, you need to think of ways to rule out the majority of those testsa 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 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

aminterested in seeing better algorithms and perhapsMathematicacannot be made faster here otherwise. – Mr.Wizard – 2013-02-15T18:46:43.863Relevant 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.613I 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