[Bug-gmp] Re: mpz_powm bug!

[Paul Zimmermann]

Dear Kengo and Torbjo"rn,

   From: Torbjorn Granlund <address@hidden>
   Date: 20 Jul 2000 18:06:26 +0200

   Kengo Nakajima <address@hidden> writes:

     See this code :

     #include <stdio.h>
     #include <gmp.h>
     int main(){
           mpz_t a,b,c,res;
           mpz_init_set_str(a,"3",10);
           mpz_init_set_str(b,"2",10);
           mpz_init_set_str(c,"9",10);
           mpz_powm(res,a,b,c);
           printf( mpz_get_str( NULL,10,res));
     }

   I can reproduce the problem if I add an mpz_init(res) before mpz_powm.

Sorry for this horrible bug for which I'm responsible!

   Here is a patch.  Paul, please check that this patch is sufficient!

Maybe it is not, however I could not prove it. The only thing I can prove is
that the output c of mpz_redc(c, a, b, m) satisfies:

(1) c - a*b/R^n = 0 (mod m)
(2) 0 <= c < R^n

where R = 2^BITS_PER_MP_LIMB and n is the number of limbs of the modulus m.
Thus it may be that several copies of m have to be subtracted. However the
worst case I could find using several random trials is c=m, and it appears
only when c is odd (of course otherwise mpz_redc is not used), and non 
square-free. So a safe patch would be to add the line:

   mpz_mod(xx, xx, mod);

instead of the '+' lines below. If it is worth investigating further, I could
try to prove that c < 2*m, or perhaps even c <= m.

Sorry again for that stupid bug!

Paul Zimmermann

   2000-07-20  Torbjorn Granlund  <address@hidden>

           * mpz/powm.c (mpz_powm): After final mpz_redc call, subtract `mod'
           from result if it is greater than `mod'.

   Index: mpz/powm.c
   ===================================================================
   RCS file: /home/cvsfiles/gmp/mpz/powm.c,v
   retrieving revision 1.11
   diff -c -r1.11 powm.c
   *** powm.c   2000/06/07 17:46:36     1.11
   --- powm.c   2000/07/20 15:39:28
   ***************
   *** 352,357 ****
   --- 352,359 ----
         {
           mpz_set_ui (g[0], 1);
           mpz_redc (xx, xx, g[0], mod, invm);
   +       if (mpz_cmp (xx, mod) >= 0)
   +    mpz_sub (xx, xx, mod);
         }
       mpz_set (res, xx);


   -- 
   Torbjörn