libtom / libtommath

LibTomMath is a free open source portable number theoretic multiple-precision integer library written entirely in C.
https://www.libtom.net
Other
650 stars 194 forks source link

s_mp_invmod_odd returns wrong result for negative numbers #534

Closed friedrichsenm closed 1 year ago

friedrichsenm commented 2 years ago

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:

https://github.com/libtom/libtommath/blob/4b47368501321c795d5b54d87a5bab35a21a7940/s_mp_invmod_odd.c#L31-L32 The comment is assuming that mp_mod gives the absolute value, but it doesn't. It will make the value positive though.

Then, in https://github.com/libtom/libtommath/blob/4b47368501321c795d5b54d87a5bab35a21a7940/s_mp_invmod_odd.c#L98-L109 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 commented 2 years ago

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 commented 2 years ago

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.