Re: mpz_powm crashes with modulus>2^2100000

[Paul Zimmermann]

   From: "Jens Kruse Andersen" <address@hidden>
   Date: Wed, 12 Mar 2003 06:48:32 +0100

   I am a new GMP user and new to MinGW and gcc. I have built GMP 4.1.2 under
   Windows XP with MinGW and MSYS using "./configure --enable-fft" and "make" on
   my Athlon XP with 256 MB ram.

   My program crashed (Windows reported an error had occured and ended the
   program) after 4 days computation involving the largest known non-Mersenne
   prime 3*2^2145353+1. The crash happened just as mpz_powm was called for the
   first time. The Windows error report starts:
   Exception information
   Code: 0xc00000fd       Flags: 0x00000000

It is my fault. I'm the original author of the current mpz_powm code.
I wrote it to optimize the computer time, without taking into account
the memory used. For your example (exponent and modulus of 2145355 bits)
the mpz_powm function will first create an auxiliary table containing
x^(2i+1) for i < 2^15, i.e. 2^14 values each of 268Ko, i.e. a total of
about 4.4Go!!!

The workaround (which I think will be in the next gmp release) is to have
a hard limit for the parameter 'k' in mpz_powm. For example kmax=7 seems
reasonable, which gives a table of 2^6 elements at most, i.e. about 17Mo
in your case. Since the cost is proportional to 1+1/(k+1), with k=7 you
loose only about 6% wrt the "optimal" k=15.

Below is a patch that limits k to 7, and in addition prints some information
to show the progress in the computation. For computing 3^p mod p for
p=3*2^2145353+1 on an old PIII-500, I estimate it will take about 120 days.

Regards,
Paul Zimmermann

leibniz% diff -u powm.c powm.c.new
--- powm.c      Wed Mar 12 14:04:11 2003
+++ powm.c.new  Wed Mar 12 13:52:38 2003
@@ -242,11 +242,12 @@
   enb = en * GMP_NUMB_BITS - e_zero_cnt; /* number of bits of exponent */
   k = 1;
   K = 2;
-  while (2 * enb > K * (2 + k * (3 + k)))
+  while (2 * enb > K * (2 + k * (3 + k)) && k < 7)
     {
       k++;
       K *= 2;
     }
+  printf ("k=%u\n", k);
 
   tp = TMP_ALLOC_LIMBS (2 * mn + 1);
   qp = TMP_ALLOC_LIMBS (mn + 1);
@@ -295,6 +296,7 @@
 
   /* Compute xx^i for odd g < 2^i.  */
 
+  printf ("computing x^2\n");
   xp = TMP_ALLOC_LIMBS (mn);
   mpn_sqr_n (tp, gp, mn);
   if (use_redc)
@@ -302,8 +304,10 @@
   else
     mpn_tdiv_qr (qp, xp, 0L, tp, 2 * mn, mp, mn);
   this_gp = gp;
+  printf ("start computing xx^i for odd g < 2^%u\n", k);
   for (i = 1; i < K / 2; i++)
     {
+      printf ("i=%u\n", i);
       mpn_mul_n (tp, this_gp, xp, mn);
       this_gp += mn;
       if (use_redc)
@@ -312,6 +316,8 @@
        mpn_tdiv_qr (qp, this_gp, 0L, tp, 2 * mn, mp, mn);
     }
 
+  printf ("finished computing xx^i for odd g < 2^%u\n", k);
+
   /* Start the real stuff.  */
   ep = PTR (e);
   i = en - 1;                          /* current index */
@@ -337,6 +343,7 @@
   MPN_COPY (xp, gp + mn * (c >> 1), mn);
   while (--j >= 0)
     {
+      printf ("j=%u\n", j);
       mpn_sqr_n (tp, xp, mn);
       if (use_redc)
        redc (xp, mp, mn, invm, tp);
@@ -354,6 +361,7 @@
          if (i > 0)
            {
              i--;
+              printf ("i=%u\n", i);
              c <<= (-sh);
              sh += GMP_NUMB_BITS;
              c |= ep[i] >> sh;
@@ -384,6 +392,7 @@
          else
            {
              i--;
+              printf ("i=%u\n", i);
              sh = GMP_NUMB_BITS - 1;
              c = (c << 1) + (ep[i] >> sh);
            }