Calculation of big powers of big numbers in a modular ring

We've seen that powers in a modular ring are important when using primality tests.  In fact, powers in modular rings will also have a lot of cryptographic applications.  Of course, a trivial way to calculate a power like aN would consist in multiplying N times a with itself.  However, this is not feasible when we are dealing with very large numbers and very large powers.  The number of iterations needed to do so would make the calculation last for ever.

However, there's a much faster way to calculate such a power.  We can calculate a2 = a * a ; a4 = a2 * a2 ; a8 = a4 * a4 ; a16 = a8 * a8 ; a32 = a16 * a16 etc.   As we see, to calculate a32, we do not need 31 multiplications, as was the case with the trivial power algorithm, but rather 5 multiplications, as spelled out.   If we need to calculate, say, a38 we can write this as a32 * a4 * a2 ; so with 2 extra multiplications we can calculate a38.

In fact, it turns out that this scheme is universal.  To calculate a power N of a number a, we calculate the powers of a which have exponents which are all the successive powers of 2, and we multiply together certain of those factors.  In fact, the factors that are to be taken, are given by the binary representation of N: there were a bit "1" occurs in the binary representation of N, we have to take the corresponding pre-calculated factor.

No matter in what base we have written N, we can find its binary representation by taking the successive Euclidean divisions of N by 2.  The successive rest values are the bits (the binary digits) of N (starting at the lsb).  For instance, if we consider (in the decimal system), the number 23, then:

  • 23 divided by 2 gives 11, rest 1.  1 is hence the lsb of 23.
  • 11 divided by 2 gives 5, rest 1.  1 is hence the second lsb of 23.
  • 5 divided by 2 gives 2, rest 1.  1 is hence the third lsb of 23
  • 2 divided by 2 gives 1, rest 0.  0 is hence the  fourth lsb of 23.
  • 1 divided by 2 gives 0, rest 1.  1 is hence the msb of 23.

The binary representation of 23 is: 10111

In order to do these calculations, we introduce an auxiliary function which is "Euclidean division by 2".

We would like to recall: Code examples

void divide2(number* x, number* q, int* r)
/* Euclidean division by two */
  {
  int i;
  int rest;
  rest = 0;
  int divi;
  for (i = nofdigits - 1 ; i >= 0 ; i--)
    {
    divi = rest * basis + x->digit[i];
    q->digit[i] = divi / 2;
    rest = divi % 2;
    }
  *r = rest;
  }  // end divide2

With the help of this function, we can finally introduce a function which calculates a number to the power another number in a modular ring:

void power(number* x, number* y, number* xy)
/* number x to the power y */
  {
  int i;
  int j;
  number p;
  number s;
  number zero;
  number half;
  number full;
  int uneven;
  int bitsequence[basisbits * nofdigits];
  int bitcounter;
  for (i = 0 ; i < nofdigits ; i++)
    {
    xy->digit[i] = 0;
    zero.digit[i] = 0;
    full.digit[i] = y->digit[i];
    }
  xy->digit[0] = 1;
  bitcounter = 0;
  /* we find the binary representation of the power */
  do
    {
    divide2(&full,&half,&uneven);
    bitsequence[bitcounter] = uneven;
    for (i = 0 ; i < nofdigits ; i++)
      {
      full.digit[i] = half.digit[i];
      }
    bitcounter++;
    }
  while (isbigger(&half,&zero) > 0);
  /* the square-and-multiply algorithm */
  for (i = bitcounter - 1 ; i >= 0 ; i--)
    {
    timesnum(xy,xy,&s);  // we square
    timesnum(&s,x,&p);   // we multiply
    if (bitsequence[i] == 1)
      {
      for (j = 0 ; j < nofdigits ; j++)
        {
        xy->digit[j] = p.digit[j];
        }
      }
    else
      {
      for (j = 0 ; j < nofdigits ; j++)
        {
        xy->digit[j] = s.digit[j];
        }
      }
    }
  }  // end power

In this function, we first determine the binary representation of the exponent.  Next, we apply the square-and-multiply system of the base, using this binary representation to decide whether we "multiply and square", or simply square.

Fermat primality test

Now that we can take powers in modular rings, the implementation of the Fermat primality test is not difficult:

int isprimefermat(number* x, int k)
/* applies the Fermat primality test to the number x with k trials randomly selected in the range 2..x-1    */
/* the function returns 1 when the number x is most probably a prime, and 0 when it is for sure not a prime */
  {
  number oldmodulo;
  number zero;
  number one;
  number testaraw;
  number q;
  number testa;
  number xminus1;
  number testresult;
  int kcnt;
  int isstillprime;
  int i;

  srand((unsigned int)time(NULL));

  /* we put temporarily the modulus equal to x */
  for (i = 0 ; i < nofdigits ; i++)
    {
    oldmodulo.digit[i] = modulo.digit[i];
    modulo.digit[i] = x->digit[i];
    zero.digit[i] = 0;
    one.digit[i] = 0;
    }
  one.digit[0] = 1;
  difnum(x,&one,&xminus1);
  isstillprime = 0;
  for (kcnt = 0; kcnt < k ; kcnt++)
    {
    /* we generate a new random test number that is not zero */
    do
      {
      for (i = 0 ; i < nofdigits ; i++)
        {
        testaraw.digit[i] = rand() % basis ; 
        }
      dividenum(&testaraw,&modulo,&q,&testa);
      }
    while (isbigger(&testa,&zero) == 0);

    /* we calculate a^(x-1)*/
    power(&testa,&xminus1,&testresult);
    /* we do the Fermat test */
    if (isbigger(&testresult,&one) != 0)
      {
      isstillprime = 1;
      }
    } // end of loop over number of random tests
    
  /* we put the original modulo parameter back */
  for (i = 0 ; i < nofdigits ; i++)
    {
    modulo.digit[i] = oldmodulo.digit[i];
    }
  return 1 - isstillprime;
  }  // end isprimefermat

There is a first trick to replace whatever is the active modulus, into the candidate prime number, because all calculations in the test are to be in the modular ring of order the candidate prime number.  At the end, the original modulus is put back.

A random number is generated.  The quality of this random number is not perfect: first of all we use the C-library rand() function, which is known not to be the best random number generator.  Next, we simply use the C modulus function to restrict an arbitrary random integer to the base.  Finally, we "fold" this number in the modular ring.  These simple methods can harm somewhat the uniform distribution over the desired interval, but in fact, we don't care much about a perfect uniformity.  The random test number shouldn't be 0 of course ; if it is, we re-draw a random test number for the Fermat test.

The Fermat test itself is in fact simply a calculation of the power, and an equality test to 1.

Miller-Rabin primality test

The Miller-Rabin primality test is somewhat more complicated, but doesn't need any more specific tools than we already have.

 
int isprimerabin(number* x, int k)
/* applies the Miller-Rabin primality test to the number x with k trials randomly selected in the range 2..x-1 */
/* the function returns 1 when the number x is most probably a prime, and 0 when it is for sure not a prime    */
  {
  number oldmodulo;
  number zero;
  number one;
  number testaraw;
  number q;
  number testa;
  number testa2;
  number xminus1;
  number rabin;
  number testresult;
  int kcnt;
  int isstillprime;
  int i;
  int u;
  int uu;
  int r;

  srand((unsigned int)time(NULL));

  /* we put temporarily the modulus equal to x */
  for (i = 0 ; i < nofdigits ; i++)
    {
    oldmodulo.digit[i] = modulo.digit[i];
    modulo.digit[i] = x->digit[i];
    zero.digit[i] = 0;
    one.digit[i] = 0;
    }
  one.digit[0] = 1;
  /* we obtain the composition of x - 1 = 2^u * r */
  difnum(x,&one,&rabin);
  difnum(x,&one,&xminus1);
  sumnum(&rabin,&zero,&q);
  u = 0;
  do
    {
    sumnum(&q,&zero,&rabin);
    divide2(&rabin,&q,&r);
    u++;
    }
  while (r == 0);
  u--;
  /* at this point, the rabin exponent r is in rabin, and u has the power of 2 */
  isstillprime = 0;
  for (kcnt = 0; kcnt < k ; kcnt++)
    {
    /* we generate a new random test number that is not zero */
    do
      {
      for (i = 0 ; i < nofdigits ; i++)
        {
        testaraw.digit[i] = rand() % basis ; 
        }
      dividenum(&testaraw,&modulo,&q,&testa);
      }
    while (isbigger(&testa,&zero) == 0 || isbigger(&testa,&one) == 0 || isbigger(&testa,&xminus1) == 0);
    /* we start the Miller-Rabin test */
    power(&testa,&rabin,&testresult);
    /* the first part of the test verifies whether this test results in a 1 or a x-1 */
    if (isbigger(&testresult,&one) != 0 && isbigger(&testresult,&xminus1) != 0)
      {
      uu = 0;
      do
        {
        uu++;
        // we take the square of the a^m and put it back in testa
        timesnum(&testresult,&testresult,&testa2);
        sumnum(&testa2,&zero,&testresult);
        if (isbigger(&testresult,&one) == 0)
          {
          isstillprime = 1;
          }
        }  // end of loop over squaring
      while (uu < u && isbigger(&testresult,&xminus1) != 0);
      // final test on last entity
      if (isbigger(&testresult,&xminus1) != 0)
        {
        isstillprime = 1;
        }       
      }  // end need for a loop test
    }  // next random number
    
  /* we put the original modulo parameter back */
  for (i = 0 ; i < nofdigits ; i++)
    {
    modulo.digit[i] = oldmodulo.digit[i];
    }
  return 1 - isstillprime;
  }  // end isprimefermat

The first part of the code finds the decomposition of p - 1 as a power of 2 times an uneven number, as the Miller-Rabin test prescribes.  Next, we implement the Miller-Rabin algorithm, using a single power, and a number of squarings, after we have implemented a random drawing of a test number, just as in the case of the Fermat test.  However, we also don't want the number to be 1 or p - 1 in this case.