Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Re: Re: How do little quickest

Daniel Lichtblau wrote:
> Artur wrote:
>> Now the quickest procedure is:
>>
>> aa = {}; Do[Print[x]; rmin = 10^10; k = 2^x; w = Floor[(k - 1)/2];
>>  Do[m = FactorInteger[k n (k - n)]; rad = 1;
>>     Do[rad = rad m[[s]][[1]], {s, 1, Length[m]}];
>>    If[rad < rmin, rmin = rad], {n, 1, w, 2}];
>>  AppendTo[aa, rmin], {x, 2, 30}]; aa
>>
>> because all possible numbers was odd <{n, 1, w, 2} inspite 
>> {n, 1, w}> and GCD checking was delated because GCD[2^x-(2n-1),(2n-1)]=1 
>>
>> What to do yet because time is still expotential and for x=27 my computer need whole day to count.
>>
>> Best wishes
>> Artur
>>
>> Mathematical problem is following:
>>
>> I need help with finding rule on follwing problem
>>
>> Sequences:
>>
>> B={1, 3, 7, 15, 27, 63, 125, 243, 343, 999, 1805, 3721, 8181, 16335,
>>
>> 32761, 65533, 112847, 190269, 519375, 1046875, 1953125}
>>
>> C={2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384,
>>
>> 32768, 65536, 131072, 262144, 524288, 1048576, 2097152}
>>
>> A={1, 1, 1, 1, 5, 1, 3, 13, 169, 25, 243, 375, 11, 49, 7, 3, 18225,
>>
>> 71875, 4913, 1701, 144027}
>>
>> C(n) = 2^n
>>
>> A(n) = C(n) - B(n)
>>
>> A(n) < B(n) < C(n)
>>
>> is obtained by algorhitm such
>>
>> that GCD[A(n),B(n),C(n)] = 1 and product of distinct prime divisors of
>>
>> A(n)*B(n)*C(n) have minimal value.
>>
>> I'm looking for formula for B(n).
>>
>> Who can help?
> 
> [...]
> 
> The code below is around 3.5 times faster than what I last sent, for the 
> size range in question. It avoids direct factorization, instead using a 
> standard sieving approach. I opted to work with sums of (approximated) 
> logs rather than products of factors, just to be certain I could use 
> packed arrays to sufficient size. The disadvantage is that this is 
> memory intensive and will require a 64 bit machine in order to have any 
> hope for handling the full range you have in mind (2^30).
> [...]

I missed a few optimizations, both in speed and memory.

For speed, it turns out that NextPrime computations were the bottleneck. 
A bit surprising, since that's not in an inner loop, but there it is. 
Easy to repair though, because we get those just by searching the sieve 
for the next position that still has a zero.

For memory I did two things. One is to combine everything into one 
Compile. This is because Compile has, in effect, call-by-value 
semantics. Thus the second invocation would start by making a copy of 
the large sieve, in effect doubling memory use. The second change was to 
not use list operations that involve Take/Drop. This unfortunately makes 
things a bit slower than otherwise, but also avoids making copies of 
large arrays.

Running on a 64 bit machine, I can now handle at least n=29.

minprods2 = Compile[{{n,_Integer}}, Module[
   {fx, len=2^(n-1), pr=3, start=2, logpr,
     n2, parts, min=100.},
   fx = Table[0.,{len}];
   While[start<=len,
     logpr = Log[N[pr]];
     Do [fx[[j]] += logpr, {j,start,len,pr}];
     While[start<=len && fx[[start]]!=0., start++];
     pr = 2*start-1;
     ];
   Round[2*Exp[Table[
     n2 = 2^(j-1);
     min = 100.;
     Do[min = Min[min,fx[[k]]+fx[[2*n2-k+1]]], {k,n2}];
     min
     , {j,n-1}]]]
   ]]

In[2]:= Timing[minprods2[28]]

Out[2]= {194.669, {6, 14, 30, 30, 42, 30, 78, 182, 1110,
  570, 1830, 6666, 2310, 2534, 5538, 9870, 20010, 141270,
  14070, 480090, 155490, 334110, 1794858, 2463270,
  2132130, 2349390, 3238746}}

In[3]:= Timing[minprods2[29]]

Out[3]= {396.514, {6, 14, 30, 30, 42, 30, 78, 182, 1110, 570,
   1830, 6666, 2310, 2534, 5538, 9870, 20010, 141270, 14070,
   480090, 155490, 334110, 1794858, 2463270, 2132130,
   2349390, 3238746, 4206930}}

The speed improvement was another factor of 3.5 or so (I'm on a slightly 
slower machine than in the previous post), so I handle n=28 size almost 
as fast as the n=26 case from the last attempt.

Daniel Lichtblau
Wolfram Research
  • Prev by Date: Re: Storing and Loading Definitions / Emulating Associative Arrays
  • Next by Date: Re: Re: Re: How do little quickest ?
  • Previous by thread: Re: Re: How do little quickest ?
  • Next by thread: Re: How do little quickest