Karim Belabas on Sun, 08 Jan 2023 23:48:04 +0100


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: lookup-optimization


* Ruud H.G. van Tol [2023-01-08 21:48]:
> 
> In the below code I build up a set 's',
> against which I first check each new case
> for "congruence".
> 
> If no such case was found,
> then the new case gets added to the set.
> 
> Each element of 's' is a tuple
> of a power-of-2 and an offset.
> 
> (1) I wonder if Mod() is better for that. Or store the sum?
> (2) Would vecsearch() be good to use for this?
> (3) Any ideas on doing things differently, are welcome.

Nothing really new, I just streamlined the code to understand what it
was doing. I also 
- removed the costly (useless) 'Set'
- avoided the computation of 2^p  (check whether valuation at 2 >= p instead).
- used bit operations when possible

All this only gains a small constant factor, though. (Less than 2.)

If I were to sort the lookup set, I would sort wrt the ѕecond entry to
maximize the likelihood of an early abort in the 'valuation' loop.
(And use a dynamic data structure so as to directly insert the new
element in the already sorted list. Didn't bother.)

A243115(N = 575) =
{ my(r= [], s= []);
  forstep(v2 = 3, N, 4,
    my(found = 0);
    for (i = 1, #s,
      if (valuation(s[i][2] - v2, 2) >= s[i][1], found = 1; break));
    if (found, next);

    my(p2 = 0, v3 = v2);
    until (v3 <= v2,
      if (bitand(v3,1), v3 += v3 >> 1 + 1;
                      , v3 >>= 1);
      p2++);
    r = concat(r, v2);
    s = concat(s, [ [p2, v2] ]));
  r;
}

Before:
? A243115(10^5);
time = 3,824 ms.

? A243115(10^6);
time = 3min, 23,506 ms.


After:
? A243115(10^5);
time = 2,056 ms.

? A243115(10^6);
time = 2min, 689 ms.

Cheers,

    K.B.

P.S. A slightly faster variant using an installed function:

install(vali,lG)

A243115(N = 575) =
{ my(r= [], s= []);
  forstep(v2 = 3, N, 4,
    my(found = 0);
    for (i = 1, #s,
      if (vali(s[i][2] - v2) >= s[i][1], found = 1; break));
    if (found, next);
 
    my(p2 = 0, v3 = v2);
    until (v3 <= v2,
      if (bitand(v3,1), v3 += v3 >> 1 + 1;
                      , v3 >>= 1);
      p2++);
    r = concat(r, v2);
    s = concat(s, [ [p2, v2] ]));
  r;
}

? A243115(10^5);
time = 1,856 ms.

? A243115(10^6);
time = 1min, 49,184 ms.

--
Pr Karim Belabas, U. Bordeaux, Vice-président en charge du Numérique
Institut de Mathématiques de Bordeaux UMR 5251 - (+33) 05 40 00 29 77
http://www.math.u-bordeaux.fr/~kbelabas/
`