| Ruud H.G. van Tol on Tue, 06 Oct 2026 17:46:55 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| logint: estimate first, one power instead of binary splitting |
As a test, I asked Anthropic Claude Fable 5.1 to make a faster logint. This is what it came up with: - - - - - - - - - - - - - - - - - - - Hello pari-dev,logintall() (src/basemath/ispower.c) finds e = floor(log_y B) for multiword B by binary splitting: the squares y^(2^i), then one trial product per bit of e.
The attached patch (against dfb0ec9, 17 lines in, 54 out) replaces that and the naive loop for small e by an estimate:
t = dbllog2(B) / dbllog2(y), e = floor(t)1. Without the third argument, if t is further than (t+1) * 1e-12 from an integer, e is returned at once. dbllog2 is good to about 2^-52 relative, so the margin is wide; no multi-precision work is done.
2. Otherwise q = powiu(y, e) and one comparison with B; e is corrected by steps of 1 as long as needed, so this path does not depend on the accuracy of t. The product q * y is skipped when its bit length already exceeds that of B.
Timings, ns per call, GMP kernel, x86-64 (old -> new): B random B B = y^k + 1 same, with &q 10^25 312 -> 114 219 -> 209 298 -> 283 10^100 459 -> 122 369 -> 277 569 -> 369 (2^1000+1)^6 5383 -> 216 4306 -> 1603 4446 -> 1973 3^1000 715 -> 115 735 -> 435 825 -> 905 3^10000 29550 -> 100 27700 -> 5550 28250 -> 5800 3^100000 1054000 -> <100 1107000 -> 172000 1169000 -> 177000 10^100000 2752000 -> <100 3026000 -> 530000 3006000 -> 572000 3^1000000 25825000 -> <100 25300000 -> 5125000 22675000 -> 5050000Checks: src/test/in/logint gives the same output; 60000 random cases around powers (bases up to 200 bits, with and without &q) agree with the definition and with the current code (logint_check.gp).
If relying on the accuracy of dbllog2 in case 1 is unwanted, dropping the early return keeps the gains of the middle column for every call.
Regards, Anthropic Claude Fable 5.1. - - - - - - - - - - - - - - - - - - - See atttachments: logint_01.patch logint_bench_01.gp logint_check_01.gp -- Ruud
diff --git a/src/basemath/ispower.c b/src/basemath/ispower.c
index 8609f4e..a145dae 100644
--- a/src/basemath/ispower.c
+++ b/src/basemath/ispower.c
@@ -61,8 +61,8 @@ long
logintall(GEN B, GEN y, GEN *ptq)
{
pari_sp av;
- long ey, e, emax, i, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
- GEN q, pow2;
+ long ey, e, emax, fl, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
+ GEN q;
if (lgefint(B) == 3)
{
@@ -83,60 +83,23 @@ logintall(GEN B, GEN y, GEN *ptq)
}
av = avma;
ey = expi(y);
- /* eB/(ey+1) - 1 < e <= eB/ey */
- emax = eB/ey;
- if (emax <= 13) /* e small, be naive */
- {
- GEN r = y, r2 = gen_1;
- for (e=1;; e++)
- { /* here, r = y^e, r2 = y^(e-1) */
- long fl = cmpii(r, B);
- if (fl >= 0)
- {
- if (fl) { e--; cgiv(r); r = r2; }
- if (ptq) *ptq = gc_INT(av, r); else set_avma(av);
- return e;
- }
- r2 = r; r = mulii(r,y);
- }
- }
- /* e >= 13 ey / (ey+1) >= 6.5 */
-
- /* binary splitting: compute bits of e one by one */
- /* compute pow2[i] = y^(2^i) [i < crude upper bound for log_2 log_y(B)] */
- pow2 = new_chunk((long)log2(eB)+2);
- gel(pow2,0) = y;
- for (i=0, q=y;; )
- {
- GEN r = gel(pow2,i); /* r = y^2^i */
- long fl = cmpii(r,B);
- if (!fl)
- {
- e = 1L<<i;
- if (ptq) *ptq = gc_INT(av, r); else set_avma(av);
- return e;
- }
- if (fl > 0) { i--; break; }
- q = r;
- if (1L<<(i+1) > emax) break;
- gel(pow2,++i) = sqri(q);
+ { /* e = floor(t) unless t is close to an integer; t has relative error
+ * < 2^-50, see dbllog2 */
+ double t = dbllog2(B) / dbllog2(y), s = (t + 1) * 1e-12;
+ e = (long)floor(t);
+ if (!ptq && t - e > s && e + 1 - t > s) return e;
}
-
- for (e = 1L<<i;;)
- { /* y^e = q < B < r = q * y^(2^i) */
- pari_sp av2 = avma;
- long fl;
- GEN r;
- if (--i < 0) break;
- r = mulii(q, gel(pow2,i));
+ emax = eB/ey; /* eB/(ey+1) - 1 < e <= eB/ey */
+ if (e > emax) e = emax;
+ q = powiu(y, e);
+ fl = cmpii(q, B);
+ while (fl > 0) { q = diviiexact(q, y); e--; fl = cmpii(q, B); }
+ while (fl < 0 && expi(q) + ey <= eB)
+ { /* y^e = q < B and q * y < 2^(eB+2) */
+ GEN r = mulii(q, y);
fl = cmpii(r, B);
- if (fl > 0) set_avma(av2);
- else
- {
- e += (1L<<i);
- q = r;
- if (!fl) break; /* B = r */
- }
+ if (fl > 0) break;
+ q = r; e++;
}
if (ptq) *ptq = gc_INT(av, q); else set_avma(av);
return e;
\\ ns per call: [far, near, with power]; far = B random, near = B = y^k + 1
tm(f, n) = my(t = getabstime()); for(i = 1, n, f()); (getabstime() - t) * 10^6 \ n;
row(name, y, k, n) = {
my(Bn = y^k + 1, Bf = y^k * 3 \ 2 + 12345, q);
if(y == 3, Bf = y^k * 5 \ 3);
printf("%-22s %10d %10d %10d\n", name,
tm(() -> logint(Bf, y), n), tm(() -> logint(Bn, y), n), tm(() -> logint(Bn, y, &q), n));
}
row("10^25", 10, 25, 2*10^6);
row("3^50", 3, 50, 2*10^6);
row("3^100", 3, 100, 10^6);
row("10^100", 10, 100, 10^6);
row("(2^64+13)^5", 2^64+13, 5, 10^6);
row("(2^1000+1)^6", 2^1000+1, 6, 3*10^5);
row("3^1000", 3, 1000, 2*10^5);
row("10^1000", 10, 1000, 2*10^5);
row("3^10000", 3, 10000, 2*10^4);
row("3^100000", 3, 100000, 1000);
row("10^100000", 10, 100000, 500);
row("(2^64+13)^10000", 2^64+13, 10000, 200);
row("3^1000000", 3, 10^6, 40);
quit
\\ results of logint with and without the power, as one checksum
chk(N) = {
my(s = 0, q);
setrand(1);
for(i = 1, N,
my(y = if(i%3, random(40) + 2, random(2^(1 + random(200))) + 2),
k = random(400) + 1,
B = y^k + (random(7) - 3) * if(i%2, 1, random(y^(k\2 + 1)) + 1));
if(B < 1, B = 1);
my(e = logint(B, y), f = logint(B, y, &q));
if(e != f || q != y^e || y^e > B || y^(e+1) <= B, error([B, y, e, f]));
s += e);
s;
}
print(chk(60000));
\\ real arguments and edge cases
print([logint(Pi*10^50, 10), logint(2^64, 2^64), logint(2^64 - 1, 2^64), logint(2^200, 2^100 + 1), logint(3^1000, 3), logint(3^1000 - 1, 3), logint(3^1000 + 1, 3), logint(10^100000 - 1, 10), logint(2^(2^20), 2^(2^19) + 1)]);
quit