Re; [PATCH] mini-gmp: avoid a mpz_sizeinbase call in mpz_get_str

[David Sparks]

> So we have a rather narrow range
>
>   1 + floor ((k-1) / log_2 b) <= n <= 1 + floor (k / log_2 b)
>
> So it seems natural to use an approximation of 1 / log_2 b. E.g., say we
> put w = ceil (2^{32} / log_2 b), and
>
>   n' = 1 + floor (k w / 2^{32})
>
> Then n' >= n, but can n' be too large? The accuracy of w is limited, so
> for large enough k we'll surely get n' > n+1, but how large must k be
> for that to become a problem?

The rest of the derivation is straightforward.
(TL;DR: The first failure occurs when b=5 and k=2655232169.)

Let's consider the two multiplicative factors:
   m = 1 / log_2 b = log_b 2   (exact)
   m' = w / 2^{32} = m + \eps  (approximation)
where 0 <= \eps < 2^{-32}.

If we are looking for n' >= n + 2, that implies that
   floor(k m') >= floor((k-1) m) + 2,
which implies that
   k m' > (k-1) m + 1.
(The worst case occurs if (k-1) m' is minutely smaller than
an integer.)
Subtracting (k-1) m from both sides, this means
   k \eps + m > 1
   k > (1-m)/\eps
If we use a fixed \eps < 2^{-32}, this is worst when m is
large, i.e. b = 3.  However, we can compute the actual \eps for
small b and see where it tapers off.  Given a bound for a
specific b, we can search until m is so small that, even if
\eps were 2^32, it wouldn't limit us.

Making the search easy, the worst case happens very early, in base 5.
2^32 * log_5 2 = 1849741732.079068, rounding up by \eps = 0.920932
to w = 1849741733.  From there, it's easy to see that we don't need
to consider bases greater than 6.144 (although I checked up to 1000
in my first draft code).

In base 5, we might fail for k > 2655163919.987920.
Brute-force searching starting there, the first actual
failure occurs at k=2655232169, when we have:

2^(k-1) = 5^1143546250.999993333...
n = 1143546251
n' = k*w >> 32 = 1143546253

There's a close call at k = 2653058431, where
2^k-1 = 5^1142610072.9999999917837433812576137004
n = 5^1142610073
n' = 1 + floor(k * w * 2^-32)
   = 1 + floor(1142610073.9995485458057373762130737305)
   = 1142610074
With this k below the bound, the error introduced (the difference
between the pre-floored values) by the overestimate in w
is less than 1, so the floor() can't make it worse.

(I did the bulk searching in IEEE double, with anything within
1e-5 of an integer double-checked manually with Pari/GP.)

Just to be sure, I re-checked brute-force searching from k=1
(it's amazing how fast a naive, scalar, single-threaded loop
can do 2e9 iterations on a mediocre computer) and that is
indeed the first failure.

FWIW, here are the first 10 failures in base-5:
Gotcha!  k=2655232169 n=1143546251 n'=1143546253
Gotcha!  k=2655281640 n=1143567557 n'=1143567559
Gotcha!  k=2655311024 n=1143580212 n'=1143580214
Gotcha!  k=2655340408 n=1143592867 n'=1143592869
Gotcha!  k=2655360495 n=1143601518 n'=1143601520
Gotcha!  k=2655389879 n=1143614173 n'=1143614175
Gotcha!  k=2655409966 n=1143622824 n'=1143622826
Gotcha!  k=2655419263 n=1143626828 n'=1143626830
Gotcha!  k=2655439350 n=1143635479 n'=1143635481
Gotcha!  k=2655459437 n=1143644130 n'=1143644132

2.6 billion bits is probably safely outside mini-GMP's
domain of usefulness.