Help stabilising mini-gmp

[Niels Möller]

tg at gmplib.org (Torbjörn Granlund) writes:

> nisse at lysator.liu.se (Niels Möller) writes:
>
>   I've tracked this one down to incorrect results from mpn_invert_3by2.
>   I'll investigate, and in the process try to document how it works, and
>   maybe rewrite it (it's a bit too hairy for its own good).

Bug found. The details: Problem was the computation of the 2/1 inverse
done by mpn_invert_3by2 (later adjusted to be a 3/2 inverse).

This was implemented as follows:

1. Compute the (approximate) high half limb of the inverse, using a plain division.

2. Adjust it to become a 3/2 inverse on half limbs, i.e, with 64 bit
   limbs, that would be floor((2^96-1) / u1) - 2^32.

3. Use this for a 3/2 division (on half-limbs) to compute the low half
   of the inverse. Similar logic used on full limbs for gmp's
   invert_4by2 (div_qr_2.c).

In the 3/2 division algorithm, candidate remainder is computed as

  r1       <-- (u1 - q1 d1) (mod B)
  {t1, t0} <-- d0 q1
  {r1, r0} <-- ({r1, u0} - {t1, t0} - {d1, d0}) (mod B^2)
  q1       <-- (q1 + 1) (mod B)

Multiplication and the extra subtraction of {d1, d0} must be done,
before incrementing q1, because incrementing might overflow. That's no
problem for the rest of the algorithm, but it would break the important
high half of the product d0 q1.

Now when working with half limbs, q1 fills only half a word so this can
be simplified by incrementing q1 first. The problem was that the
variable wasn't typed as a full limb (mp_limb_t) but as an unsigned,
which on 64-bit archs typically is exactly half a limb. And bug
triggered when low half of the inverse is 2^32 - 1.

> I am very worried about this, as well as past and lingering arithmetic
> bugs in mini-gmp.  Important lessons need to be learned, and full
> disclosure to our user needs to be made.

I think a bug like this ought to have been caught by asserts in the
division code. So that using an incorrect inverse results in a crash
rather than silent miscomputation.

> I have not looked at the mini-gmp sources, except very briefly.  I
> understood the project as a plain, safe implementation of core functions
> of GMP.

It's not always obvious where to draw the line. The objective is to be
at most 10 times slower than real gmp, and then I think division
corresponding to gmp:s _pi1 schoolbook division is reasonable. And the
way to do it is pretty well understood. The "novelty" is to compute the
inverse using a single-limb divide instruction for the high half, and a
3/2 for the low half.

> If we still believe in the idea of mini-gmp as part of a known robust
> library as GMP, I suggest that we do the following: Do not write any new
> code such as mpn_invert_3by2.  Only grab code from GMP, extracting the
> most fundamental algorithms.  Sometimes this will require editing to
> simplify the code; this work needs to be done carefully and
> thoughtfully.  Perhaps this is boring, but I am sure that our users
> prefer the boring property of correct computation results over never so
> surprising miscalculations!

To this, I'd like to add that we should use asserts whenever practical.
E.g., mpn_div_qr_pi1 only uses assert to check input requirements, but
it should be pretty easy to verify that each quotient used is sane
(correct or one off in the way that is expected). In mini-gmp, I think
it is reasonable to spend a % or two of the cpu resources on sanity
checks.

>   It could also use some units test of its own; result is easy to
>   validate.
>   
> Perhaps better tests are needed.

There already were unit tests for mpn_invert_3by2
(mini-gmp/tests/t-invert). But we didn't get sufficient coverage in
inputs until testing was enabled in gmp night builds.

And the t-div and t-powm failures I reproduced yesterday seem to also
have been caused by the broken inverse.

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.