s_mp_invmod_odd returns wrong result for negative numbers

[friedrichsenm]

Using mp_invod can return incorrect results for negative numbers. Using it to find the inverse of -1 mod 7 yields -6 instead of 6. The problem looks like it could be because of a couple of places:

libtommath/s_mp_invmod_odd.c

Lines 31 to 32 in 4b47368

	/* we need y = |a| */
	if ((err = mp_mod(a, b, &y)) != MP_OKAY) goto LBL_ERR;

The comment is assuming that mp_mod gives the absolute value, but it doesn't. It will make the value positive though.

Then, in

libtommath/s_mp_invmod_odd.c

Lines 98 to 109 in 4b47368

	sign = a->sign;
	while (mp_isneg(&D)) {
	if ((err = mp_add(&D, b, &D)) != MP_OKAY) goto LBL_ERR;
	}
	/* too big */
	while (mp_cmp_mag(&D, b) != MP_LT) {
	if ((err = mp_sub(&D, b, &D)) != MP_OKAY) goto LBL_ERR;
	}
	mp_exch(&D, c);
	c->sign = sign;

the sign of a is used to flip the sign of the inverse, but the sign should already be correct since mp_mod is using a positive member of the equivalence class for a. Removing this sign change should fix it.

[czurnieden]

That output is correct.

LibTomMath computes two different remainders. If we take (-1)/7 as the example, we have two solutions:, mp_div computes the quotient q = 0 and the remainder r = -1; mp_mod computes the remainder r = 6, the quotient is not computed but would be q = -1, which concludes the second solution. The second solution is also the one used in standard mathematics, so PARI/GP gives [-1,6] for the input divrem(-1,7).

The comment about the absolute value /* we need y = |a| */ was valid last in version 0.33, a couple of years ago.

mp_invmod computes a*c equiv 1 ( mod b ) for the input c = mp_invmod(a,b) according to the documentation.

I think Pari/GP can be trusted to be correct so let it check s_mp_invmod_odd's result c = -6 for the input a = -1, b = 7:

? (-1 * -6) == Mod(1,7)
%1 = 0

That is correct, 6 notequiv 1 (mod 7) but we expected that.

Pari/GP gives 6 as the result:

? Mod(1/(-1),7)
%4 = Mod(6, 7)
? -1 == Mod(6,7)
%5 = 1
? 6 == Mod(-1,7)
%6 = 1
? -6 == Mod(1,7)
%7 = 1

That is correct, (-1 * 6) equiv 1 (mod 7) but we expected that.

The (slower but general) function s_mp_invmod gives the correct result.

So what happened and why did it not get caught?
The latter is answered quickly: there is only a single test for mp_invmod and that is a check for a bugfix for an edgecase and nothing more.
I think that just bit back.

I cannot see a reason for that sign to get copied over and would concur with @friedrichsenm and remove that part completely, including the mentioned comment. Additionally: add some tests to test.c.

Might take an hour or two.

(Sorry, @ssinai1 bugfixes first, but you are not forgotten! )

[czurnieden]

While writing the tests and making a small table with Pari/GP I have been reminded of the fact that mp_invmod does not make use of the symmetry (-n * c) mod -m = (-n * c) mod m for c, n,m in N\{0}

The generated table:

m = -10, -9, -8, -7, -6, -5, -4, -3, -2, -1,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10 
   { -1,  8, -1,  2, -1, -1, -1,  2, -1,  0, -1,  0, -1,  2, -1, -1, -1,  2, -1,  8, -1 }  n =  -10
   {  1, -1,  7,  3, -1,  1,  3, -1,  1,  0, -1,  0,  1, -1,  3,  1, -1,  3,  7, -1,  1 } -9
   { -1,  1, -1,  6, -1,  3, -1,  1, -1,  0, -1,  0, -1,  1, -1,  3, -1,  6, -1,  1, -1 } -8
   {  7,  5,  1, -1,  5,  2,  1,  2,  1,  0, -1,  0,  1,  2,  1,  2,  5, -1,  1,  5,  7 } -7
   { -1, -1, -1,  1, -1,  4, -1, -1, -1,  0, -1,  0, -1, -1, -1,  4, -1,  1, -1, -1, -1 } -6
   { -1,  7,  3,  4,  1, -1,  3,  1,  1,  0, -1,  0,  1,  1,  3, -1,  1,  4,  3,  7, -1 } -5
   { -1,  2, -1,  5, -1,  1, -1,  2, -1,  0, -1,  0, -1,  2, -1,  1, -1,  5, -1,  2, -1 } -4
   {  3, -1,  5,  2, -1,  3,  1, -1,  1,  0, -1,  0,  1, -1,  1,  3, -1,  2,  5, -1,  3 } -3
   { -1,  4, -1,  3, -1,  2, -1,  1, -1,  0, -1,  0, -1,  1, -1,  2, -1,  3, -1,  4, -1 } -2
   {  9,  8,  7,  6,  5,  4,  3,  2,  1,  0, -1,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9 } -1
   { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,  0, -1, -1, -1, -1, -1, -1, -1, -1, -1 }  0
   {  1,  1,  1,  1,  1,  1,  1,  1,  1,  0, -1,  0,  1,  1,  1,  1,  1,  1,  1,  1,  1 }  1
   { -1,  5, -1,  4, -1,  3, -1,  2, -1,  0, -1,  0, -1,  2, -1,  3, -1,  4, -1,  5, -1 }  2
   {  7, -1,  3,  5, -1,  2,  3, -1,  1,  0, -1,  0,  1, -1,  3,  2, -1,  5,  3, -1,  7 }  3
   { -1,  7, -1,  2, -1,  4, -1,  1, -1,  0, -1,  0, -1,  1, -1,  4, -1,  2, -1,  7, -1 }  4
   { -1,  2,  5,  3,  5, -1,  1,  2,  1,  0, -1,  0,  1,  2,  1, -1,  5,  3,  5,  2, -1 }  5
   { -1, -1, -1,  6, -1,  1, -1, -1, -1,  0, -1,  0, -1, -1, -1,  1, -1,  6, -1, -1, -1 }  6
   {  3,  4,  7, -1,  1,  3,  3,  1,  1,  0, -1,  0,  1,  1,  3,  3,  1, -1,  7,  4,  3 }  7
   { -1,  8, -1,  1, -1,  2, -1,  2, -1,  0, -1,  0, -1,  2, -1,  2, -1,  1, -1,  8, -1 }  8
   {  9, -1,  1,  4, -1,  4,  1, -1,  1,  0, -1,  0,  1, -1,  1,  4, -1,  4,  1, -1,  9 }  9
   { -1,  1, -1,  5, -1, -1, -1,  1, -1,  0, -1,  0, -1,  1, -1, -1, -1,  5, -1,  1, -1 } 10

I faintly remember there was an argument against it but forgot and couldn't find it.