5
$\begingroup$

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

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

2 Answers 2

Reset to default
4
$\begingroup$

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

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.