Re: request for a few minutes CPU-time
- To: mathgroup at smc.vnet.net
- Subject: [mg79521] Re: [mg79477] request for a few minutes CPU-time
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sat, 28 Jul 2007 05:37:43 -0400 (EDT)
- References: <200707271001.GAA03340@smc.vnet.net>
Peter Pein wrote:
> Dear group,
>
> I've written code, which looks for numbers for which a smaller natural
> number exists which has the same sum of cubes of it's divisors:
>
> Block[{spa, nmax = 6*106, expo = 3},
> Reap[For[n = 1, n <= nmax, n++,
> (If[Head[#1] === spa, #1 = n, Sow[{n, #1}]] & )[
> spa[DivisorSigma[expo, n]]]]][[2,1]]]
>
> the name spa is an artefact; I tried this with SparseArrays, but the
> allowed range of indices has not been sufficient. I now use spa as an
> initially undefined function (Block[{spa..}]). For each n to test I look
> wether spa[sigma(r,n)] has been defined. If not, the Head is still spa
> and I set spa[sigma(r,n)] to n; else the remembered value together with
> n will go to the result via the Sow-Reap mechanism.
>
> If you've got RAM (4GB or so) than I do (1.5 GB),
> could you please run this code with, say nmax=10 or 20 million? On my
> machine it swapped heavily with nmax=6 million and I had to kill
> MathKernel as I tried nmax=10^7. The lines above took ~181 seconds to
> evaluate (nmax=10^7 has been stopped by me after 15 minutes). I do not
> expect any runtimes of more than ~7-10 minutes. Would this be possible,
> please?
>
> Alternatively any hints how to calculate these sequences more efficient
> would be highly appreciated (AFAIK there exists no kind of "inverse
> function" to sigma(r,n) w.r.t. n which could be calculated without this
> brute-force method).
>
> Thank you for your attention and in advance for CPU-time,
>
> Peter
Here is what happens for me at 10^6.
len = 10^6;
Timing[divisorsumcubes = Table[{j, DivisorSigma[3, j]}, {j, len}];]
Out[12]= {16.473, Null}
Your method, more or less:
In[13]:= Timing[l1 = Reverse[Last[Last[Module[{j, seen},
Reap[Do[j = seen[divisorsumcubes[[n, 2]]];
If[Head[j] === seen,
seen[divisorsumcubes[[n, 2]]] = divisorsumcubes[[n, 1]],
Sow[j]],
{n, len, 1, -1}]]]]]]]
Out[13]= {12.6128, {194315, 295301, 590602}}
A different approach is to sort on the second values, then split based
on that value and keep only the runs of length larger than one.
In[14]:= Timing[
l2 = Map[#[[-1, -1]] &,
Select[Split[
Sort[Map[Reverse, divisorsumcubes]], #1[[1]] === #2[[1]] &],
Length[#] > 1 &]]]
Out[14]= {9.5966, {194315, 295301, 590602}}
At 10^7 the first approach brought my machine to its knees due to memory
use. I do not know if this is something bad in our implementation of
pattern free down values (which uses hashing) or if there really is a
huge memory requirement at play.
For the second approach, well, the common first step of computing all
the DivisorSigma[3,xxx]'s takes around 200 seconds. Step two, the
split-sort thing, is as below.
Out[17]= {105.73, {194315, 295301, 590602, 1181204, 1476505,
1886920, 2067107, 2362408, 2526095, 2953010, 3248311,
3691985, 3838913, 4134214, 4469245, 4724816, 5020117,
5610719, 5635135, 6023765, 5906020, 6496622, 6791923,
7382525, 7677826, 7966915, 8355545, 8563729, 8268428,
9132805, 9449632}}
Some time in the not-too-far-distant future I expect Mathematica will
support 64-bit integers in packed arrays, and approach #2 to this
computation will become less stressful in terms of memory requirement on
64-bit machines.
Daniel Lichtblau
Wolfram Research
- References:
- request for a few minutes CPU-time
- From: Peter Pein <petsie@dordos.net>
- request for a few minutes CPU-time