Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

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

2 Answers 2

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

share|improve this answer

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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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