a property of multisets: unit fractions with a twist

[Maximilian Hasler]

Hugo,
thanks for your code.
I made some simplifications, and got a minor speedup
(might take us from 300 years down to 200 years...),
achieved in writing the body of the "big" fordiv() loop using several
... & next type statements.
(Note: loop variables are automatically local in PARI, so if you
declare them outside the loop with local(...) they are in fact
declared twice.)

To make it self contained, I completed for the case n=1 and n=2,
this does not take up many microseconds.

It can be seen that the function often calculates the same thing (for
a(6), the n=3 piece is done 344 times, but there are only 113
different results, and in 344-228 cases the result is 0 or 1). Thus,
it could be worth while implementing an (intelligent) "remember"
feature and/or some additional check.

Maximilian

A118085(n,p=2,q=1,cmin=0) = {
 local(cmax = floor(1 / (sqrtn(p / q, n) - 1)),
       count = 0, p2, q2, d, fmin);
 cmin = max( cmin, q \ (p - q) + 1);
 if (n==3,/* p*cur > q*(cur+1) <=> (p-q)cur>q <=> cur>=[ q/(p-q) ]+1 */
    for(cur = cmin, cmax, /* this ensures d > 0 */
       if( 1 < d = gcd( p2 = p*cur, q2 = q*(cur + 1)), p2\=d; q2\=d);
       d = p2 - q2;
       p2 *= q2;
       fmin = cur * d - q2;
       if (d > 1, /* else d==1 since we know d>0 */
         fordiv (p2, f,
           f < fmin & next; f * f > p2 & break;
           (f + q2) % d & next; (p2 / f + q2) % d | count++
         )
       , /* here d==1 */
         count += numdiv(p2) \ 2;
         if (fmin > 1,
           fordiv (p2, f, f >= fmin & break; count-- )
         );
       );
    ); count
 ,
    if( n<3, 1 /* result for n=1 and n=2 */
    ,  sum( cur = cmin, cmax, A118085(n-1, p*cur, q*(cur+1), cur))
    );
 )}

values of "count" calculated for a(6):
5413,3047,2046,1449,1535,1265,666,761,576,781,602,799,345,311,365,
385,160,204,125,246,84,43,57,99,38,32,29,18,21,34,12,6,11,6,10,3,1,1,5,
3,0,1,0,0,1,0,2332,1244,781,667,502,445,422,213,272,232,124,124,49,44,
89,19,15,37,14,2,7,0,2,2,0,0,1,0,1535,799,563,331,240,149,256,258,58,50,
37,28,18,25,8,3,4,2,0,1,1,0,934,502,398,271,206,175,148,64,58,36,12,20,
6,13,1,1,0,1,0,180,256,282,127,70,36,25,47,3,2,8,0,0,1,0,0,180,143,125,54,
31,28,25,10,2,2,0,1,0,57,103,23,9,13,6,0,1,1,0,0,17,24,10,7,5,1,0,1,0,7,25,
4,3,0,0,0,2,7,0,0,0,0,2,4,0,0,1,1,0,0,0,1809,993,595,407,347,398,229,217,
186,83,58,116,34,24,12,23,14,9,2,5,0,1,0,1,0,799,398,256,180,143,125,54,
31,28,25,10,2,2,0,1,0,445,258,183,130,108,33,23,17,6,6,1,1,0,142,108,87,
21,16,8,2,1,2,0,44,54,9,7,1,3,0,0,22,21,6,1,1,0,2,5,0,0,1,1,0,0,0,576,256,86,
91,29,16,16,6,2,1,1,1,130,87,69,22,21,6,1,1,0,38,43,7,1,2,0,0,10,7,1,1,0,1,
1,0,0,1,0,33,43,11,6,8,0,0,18,20,4,1,0,1,3,0,0,1,1,0,3,11,0,0,4,2,0,0,0,0,0,0,0

On Dec 18, 2007 6:45 AM,  <hv at crypt.org> wrote:
> hv at crypt.org wrote:
> :Using this approach, I calculate a(7) = 13005235 with the code below
> :(which takes 860 minutes on my computer - for comparison, a(6) takes
> :just under 12s).
> :
> :Given my track record so far on this, I would welcome confirmation. :)
> :
> :Further optimisations are possible [...]
>
> After trying a variety of them, and then converting the whole thing
> into PARI code (and then finding that under PARI the recursive version
> was after all faster than the stack-based version), the timings are now:
>
> ? print(A118085(6));
> 49513
> time = 1,711 ms.
> ? print(A118085(7));
> 13005235
> time = 2h, 19mn, 38,180 ms.
>
...
>
> Based on these timings and the dominant factor, I'd expect a(8) to take
> somewhere between 300 and 85000 CPU years on my machine.
>
> Hugo



"Maximilian Hasler" <maximilian.hasler at gmail.com> wrote:
:Hugo,
:thanks for your code.
:I made some simplifications, and got a minor speedup

Excellent stuff. I have a longer reply to your comments, but I feel we
must be getting a bit boring for the rest of seqfan, so I'll send that
offlist.

Hugo