6.2 M-ary Exponentiation

Through a generalization of the binary algorithm on page 80 the number of modular multiplications for exponentiation can be reduced even further. The idea is to represent the exponent in a base greater than 2 and to replace multiplication by a in step 3 by multiplication by powers of a. Thus let the exponent e be given by e = (e_n−1 e_n−2...e₀)_M, to a base M yet to be determined. The following algorithm calculates the powers a^e mod m.

M-ary algorithm for exponentiation a^e mod m

1. Calculate and store a² mod m, a³ mod m,..., a^M−1 mod m as a table.

2. Set p ← mod m and i ← n − 2.

3. Set p ← p^M mod m.

4. If e_i ≠ 0, set p ← mod m.

5. Set i ← i − 1; if i ≥ 0, go to step 3.

6. Output p.

The number of necessary multiplications evidently depends on the number of digits of the exponent e and thus on the choice of M. Therefore, we would like to determine M such that the exponentiation in step 3 can be computed to the greatest extent possible by means of squaring, as in the example above for 2¹⁶, and such that the number of multiplications by the precomputed powers of a be minimized to a justifiable cost of storage space for the table.

The first condition suggests that we choose M as a power of 2: M = 2^k. In view of the second condition we consider the number of modular multiplications as a function of M:

We require

(6.1)

squares in step 3 and on average

(6.2)

modular multiplications in step 4, where

is the probability that a digit e_i of e is nonzero. If we include the M − 2 multiplications for the computation of the table, then the M-ary algorithm requires on average

(6.3)

(6.4)

modular squarings and multiplications.

For exponents e and moduli m of, say, 512 binary places and M = 2^k we obtain the numbers of modular multiplications for the calculation of a^e mod m as shown in Table 6.1. The table shows as well the memory requirement for the precomputed powers of a mod m, which result from the product (2k − 2)CLINTMAXSHORT · sizeof(USHORT).

We see from the table that the average number of multiplications reaches a minimum of 640 at k = 5, where the required memory for each larger k grows by approximately a factor of 2. But what are the time requirements for other orders of magnitude of the exponents?

Table 6.2 gives information about this. It displays the requirements for modular multiplication for exponentiation with exponents with various numbers of binary digits and for various values of M = 2^k. The exponent length of 768 digits was included because it is a frequently used key length for the RSA cryptosystem (see Chapter 16). The favorable numbers of multiplications appear in boldface type.

In consideration of the ranges of numbers for which the FLINT/C package was developed, it appears that with k = 5 we have found the universal base M = 2^k, with which, however, there is a rather high memory requirement of 15 kilobytes for the powers a², a³,..., a³¹ to base a that are to be precomputed. The M-ary algorithm can be improved, however, according to [Cohe], Section 1.2, to the extent that we can employ not M − 2, but only M/2, premultiplications and thus require only half the memory. The task now additionally consists in the calculation of the power a^e mod m, where e = (e_n−1e_n−2...e₀)_M is the representation of the exponent to the base M = 2^k.

M-ary Algorithm for exponentiation with reduced number of premultiplications

1. Compute and store a³ mod m, a⁵ mod m, a⁷ mod m,..., mod m.

2. If e_n−1 = 0, set p ← 1.

   If e_n−1 ≠ 0, factor e_n−1 = 2^tu with odd u. Set p ← a^u mod m.

   In each case set i ← n − 2.

3. If e_i = 0, set p ← mod m by calculating mod m (k-fold squaring modulo m).

   If e_i ≠ 0, factor e_i = 2^tu with odd u; set p ← mod m and then p ← pa^u mod m; now set p ← mod m.

4. Set i ← i − 1; if i ≥ 0, go to step 3.

5. Output p.

The trick of this algorithm consists in dividing up the squarings required in step 3 in a clever way, such that the exponentiation of a is taken care of together with the even part 2^t of e_i. Within the squaring process the exponentiation of a by the odd part u of e_i remains. The balance between multiplication and squaring is shifted to the more favorable squaring, and only the powers of a with odd exponent need to be precomputed and stored.

For this splitting one requires the uniquely determined representation e_i = 2^tu, u odd, of the exponent digit e_i. For rapid access to t and u a table is used, which, for example, for k = 5 is displayed in Table 6.3.

To calculate these values we can use the auxiliary function twofact_l(), which will be introduced in Section 10.4.1. Before we can program the improved M-ary algorithm there remains one problem to be solved: How, beginning with the binary representation of the exponent or the representation to base B = 2¹⁶, do we efficiently obtain its representation to base M = 2^k for a variable k > 0? It will be of use here to do a bit of juggling with the various indices, and we can "mask out" the required digits e_i to base M from the representation of e to base B. For this we set the following: Let (ε_r−1ε_r−2...ε₀)₂ be the representation of the exponent e to base 2 (we need this on account of the number r of binary digits). Let (e_u−1e_u−2...e₀)_B be the representation of e as a CLINT type to base B = 2¹⁶, and let be the representation of e to the base M = 2^k, k ≤ 16 (M should not be greater than our base B). The representation of e in memory as a CLINT object e_l corresponds to the sequence [u + 1],[e₀],[e₁],...,[e_u−1],[0] of USHORT values e_l[i] for i = 0,..., u + 1; one should note that we have added a leading zero.

Let f := , and for i = 0,..., f let s_i := and d_i := ki mod 16. With these settings the following statements hold:

1. There are f + 1 digits in ; that is, n − 1 = f.

2. contains the least-significant bit of the digit .

3. d_i specifies the position of the least-significant bit of in (counting of positions begins with 0). If i < f and d_i > 16 − k, then not all the binary digits of are in ; the remaining (higher-valued) bits of are in +1. The desired digit thus corresponds to the k least-significant binary digits of

   

As a result we have for i Î {0,..., f} the following expression for determining :

(6.5)

Since for the sake of simplicity we set e_l [s_f + 2] ← 0, this expression holds as well for i = f.

We have thus found an efficient method for accessing the digits of the exponent in its CLINT representation, which arise from its representation in a power-of-two base 2^k with k ≤ 16, whereby we are saved an explicit transformation of the exponent. The number of necessary multiplications and squarings for the exponentiation is now

(6.6)

where in comparison to μ₁(k) (see page 86) the expenditure for the precomputations has been reduced by half. The table for determining the favorable values of k (Table 6.4) now has a somewhat different appearance.

Starting with 768 binary digits of the exponent, the favorable values of k are larger by 1 than those given in the table for the previous version of the exponentiation algorithm, while the number of required modular multiplications has easily been reduced. It is to be expected that this procedure is on the whole more favorable than the variant considered previously. Nothing now stands in the way of an implementation of the algorithm.

To demonstrate the implementation of these principles we select an adaptive procedure that uses the appropriate optimal value for k. To accomplish this we rely again on [Cohe] and look for, as is specified there, the smallest integer value k that satisfies the inequality

(6.7)

which comes from the formula μ₂(k) given previously for the number of necessary multiplications based on the condition μ₂(k + 1) − μ₂(k) ≥ 0. The constant number of modular squarings ⌊log₂ e⌋ for all algorithms for exponentiation introduced thus far is eliminated; here only the "real" modular multiplications, that is, those that are not squarings, are considered.

The implementation of exponentiation with variable k requires a large amount of main memory for storing the precomputed powers of a; for k = 8 we require about 64 Kbyte for 127 CLINT variables (this is arrived at via (2⁷ − 1) * sizeof(USHORT) * CLINTMAXSHORT), where two additional automatic CLINT fields were not counted. For applications with processors or memory models with segmented 16-bit architecture this already has reached the limit of what is possible (see in this regard, for example, [Dunc], Chapter 12, or [Petz], Chapter 7).

Depending on the system platform there are thus various strategies appropriate for making memory available. While the necessary memory for the function mexp5_l() is taken from the stack (as automatic CLINT variables), with each call of the following function mexpk_l() memory is allocated from the heap. To save the expenditure associated with this, one may imagine a variant in which the maximum needed memory is reserved during a one-time initialization and is released only at the end of the entire program. In each case it is possible to fit memory management to the concrete requirements and to orient oneself to this in the commentaries on the following code.

One further note for applications: It is recommended always to check whether it suffices to employ the algorithm with the base M = 2⁵. The savings in time that comes with larger values of k is relatively not so large in comparison to the total calculation time so as to justify in all cases the greater demand on memory and the thereby requisite memory management. Typical calculation times for various exponentiation algorithms, on the basis of which one can decide whether to use them, are given in Appendix D.

The algorithm, implemented with M = 2⁵, is contained in the FLINT/C package as the function mexp5_l(). With the macro EXP_L() contained in flint.h one can set the exponentiation function to be used: mexp5_l() or the following function mexpk_l() with variable k.

Function:	modular exponentiation
Syntax:	int mexpk_l (CLINT bas_l, CLINT exp_l, CLINT p_l, CLINT m_l);
Input:	bas_l (base) exp_l (exponent) m_l (modulus)
Output:	p_l (power residue)
Return:	E_CLINT_OK if all is ok E_CLINT_DBZ if division by 0 E_CLINT_MAL if malloc() error

We begin with a segment of the table for representing e_i = 2^tu, u odd, 0 ≤ e_i < 2⁸. The table is represented in the form of two vectors. The first, twotab[], contains the exponents t of the two-factor 2^t, while the second, oddtab[], holds the odd part u of a digit 0 ≤ e_i < 2⁵. The complete table is contained, of course, in the FLINT/C source code.

static int twotab[] =
{0,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5, ...};
static USHORT oddtab[]=
{0,1,1,3,1,5,3,7,1,9,5,11,3,13,7,15,1,17,9,19,5,21,11,23,3,25,13, ...};

int
mexpk_l (CLINT bas_l, CLINT exp_l, CLINT p_l, CLINT m_l)
{

The definitions reserve memory for the exponents plus the leading zero, as well as a pointer clint **aptr_l to the memory still to be allocated, which will take pointers to the powers of bas_l to be precomputed. In acc_l the intermediate results of the exponentiation will be stored.

  CLINT a_l, a2_l;
  clint e_l[CLINTMAXSHORT + 1];
  CLINTD acc_l;
  clint **aptr_l, *ptr_l;
  int noofdigits, s, t, i;
  unsigned int k, lge, bit, digit, fk, word, pow2k, k_mask;

Then comes the usual checking for division by 0 and reduction by 1.

  if (EQZ_L (m_l))
    {
      return E_CLINT_DBZ;
    }

  if (EQONE_L (m_l))
    {
      SETZERO_L (p_l);    /* modulus = 1 ==> residue = 0 */
      return E_CLINT_OK;
    }

Base and exponent are copied to the working variables a_l and e_l, and any leading zeros are purged.

  cpy_l (a_l, bas_l);
  cpy_l (e_l, exp_l);

Now we process the simple cases a⁰ = 1 and 0^e = 0 (e > 0).

  if (EQZ_L (e_l))
    {
      SETONE_L (p_l);
      return E_CLINT_OK;
    }

  if (EQZ_L (a_l))
    {
      SETZERO_L (p_l);
      return E_CLINT_OK;
    }

Next, the optimal value for k is determined; the value 2^k is stored in pow2k, and in k_mask the value 2^k−1. For this the function ld_l() is used, which returns the number of binary digits of its argument.

  lge = ld_l (e_l);
  k = 8;
  while (k > 1 && ((k - 1) * (k << ((k - 1) << 1))/((1 << k ) - k - 1)) >= lge - 1)
    {
      --k;
    }
  pow2k = 1U << k;
  k_mask = pow2k - 1U;

Memory is allocated for the pointers to the powers of a_l to be computed. The base a_l is reduced modulo m_l.

  if ((aptr_l = (clint **) malloc (sizeof(clint *) * pow2k)) == NULL)
    {
      return E_CLINT_MAL;
    }
  mod_l (a_l, m_l, a_l);
  aptr_l[1] = a_l;

If k > 1, then memory is allocated for the powers to be computed. This is not necessary for k = 1, since then no powers have to be precomputed. In the following setting of the pointer aptr_l[i] one should note that in the addition of an offset to a pointer p a scaling of the offset by the compiler takes place, so that it counts objects of the pointer type of p.

We have already mentioned that the allocation of working memory can be carried out alternatively in a one-time initialization. The pointers to the CLINT objects would in this case be contained in global variables outside of the function or in static variables within mexpk_l().

    if (k > 1)
      {
         if ((ptr_l = (clint *) malloc (sizeof(CLINT) * ((pow2k >> 1) - 1))) == NULL)
           {
              return E_CLINT_MAL;
           }

         aptr_l[2] = a2_l;
         for (aptr_l[3] = ptr_l, i = 5; i < (int)pow2k; i+=2)
           {
              aptr_l[i] = aptr_l[i - 2] + CLINTMAXSHORT;
           }

Now comes the precomputation of the powers of the value a stored in a_l. The values a³, a⁵, a⁷,...,a^k − 1 are computed (a² is needed only in an auxiliary role).

    msqr_l (a_l, aptr_l[2], m_l);
    for (i = 3; i < (int)pow2k; i += 2)
       {
         mmul_l (aptr_l[2], aptr_l[i - 2], aptr_l[i], m_l);
       }
   }

This ends the case distinction for k > 1. The exponent is lengthened by the leading zero.

  *(MSDPTR_L (e_l) + 1) = 0;

The determination of the value f (represented by the variable noofdigits).

  noofdigits = (lge - 1)/k;
  fk = noofdigits * k;

Word position s_i and bit position d_i of the digit e_i in the variables word and bit.

  word = fk >> LDBITPERDGT;    /* fk div 16 */
  bit = fk & (BITPERDGT-1U);    /* fk mod 16 */

Calculation of the digit e_n−1 with the above-derived formula; e_n−1 is represented by the variable digit.

  switch (k)
    {
      case 1:
      case 2:
      case 4:
      case 8:
         digit = ((ULONG)(e_l[word + 1] ) >> bit) & k_mask;
         break;
      default:
         digit = ((ULONG)(e_l[word + 1] | ((ULONG)e_l[word + 2]
                                               << BITPERDGT)) >> bit) & k_mask;
    }

First run through step 3 of the algorithm, the case digit = e_n−1 ≠ 0.

    if (digit != 0)    /* k-digit > 0 */
       {
         cpy_l (acc_l, aptr_l[oddtab[digit]]);

Calculation of ; t is set to the two-part of e_n−1 via twotab[e_n−1]; p is represented by acc_l.

         t = twotab[digit];
         for (; t > 0; t--)
            {
              msqr_l (acc_l, acc_l, m_l);
            }
            }
         else    /* k-digit == 0 */
            {
              SETONE_L (acc_l);
            }

Loop over noofdigits beginning with f - 1.

  for (--noofdigits, fk -= k; noofdigits >= 0; noofdigits--, fk -= k)
    {

Word position s_i and bit position d_i of the digit e_i in the variables word and bit.

     word = fk >> LDBITPERDGT;     /* fk div 16 */
     bit = fk & (BITPERDGT - 1U);     /* fk mod 16 */

Computation of the digit e_i with the formula derived above; e_i is represented by the variable digit.

     switch (k)
        {
          case 1:
          case 2:
          case 4:
          case 8:
             digit = ((ULONG)(e_l[word + 1] ) >> bit) & k_mask;
             break;
          default:
             digit = ((ULONG)(e_l[word + 1] | ((ULONG)e_l[word + 2]
                                                   << BITPERDGT)) >> bit) & k_mask;
        }

Step 3 of the algorithm, the case digit = e_i ≠ 0; t is set via the table twotab[e_i] to the two-part of e_i.

     if (digit != 0)    /* k-digit > 0 */
        {
          t = twotab[digit];

Calculation of a^u in acc_l. Access to a^u with the odd part u of e_i is via aptr_l [oddtab [e_i]].

         for (s = k - t; s > 0; s--)
            {
              msqr_l (acc_l, acc_l, m_l);
            }

         mmul_l (acc_l, aptr_l[oddtab[digit]], acc_l, m_l);

Calculation of ; p is still represented by acc_l.

         for (; t > 0; t--)
            {
              msqr_l (acc_l, acc_l, m_l);
            }
       }
    else    /* k-digit == 0 */
       {

Step 3 of the algorithm, case e_i = 0: Calculate .

         for (s = k; s > 0; s--)
            {
              msqr_l (acc_l, acc_l, m_l);
            }
        }
    }

End of the loop; output of acc_l as power residue modulo m_l.

  cpy_l (p_l, acc_l);

At the end, allocated memory is released.

  free (aptr_l);
  if (ptr_l != NULL)
    free (ptr_l);
  return E_CLINT_OK;
}

The various processes of M-ary exponentiation can be clarified with the help of a numerical example. To this end let us examine the calculation of the power 1234⁶⁶⁷ mod 18577, which will be carried out by the function mexpk_l() in the following steps:

1. Precomputations

   The representation of the exponent e = 667 can be expressed to the base 2^k with k = 2 (cf. the algorithm for M-ary exponentiation on page 86), whereby the exponent e has the representation e = .

   The power a³ mod 18577 has the value 17354. Further powers of a do not arise in the precomputation because of the small size of the exponent.

2. Exponentiation loop

   exponent digit ei = 2tu	21 · 1	21 · 1	20 · 1	21 · 1	20 · 3
   p ← p2 mod n	-	14132	13261	17616	13599
   p ← mod n	-	-	4239	-	17343
   p ← pau mod n	1234	13662	10789	3054	4445
   p ← p2 mod n	18019	7125	-	1262	-

3. Result

   p = 1234⁶⁶⁷ mod 18577 = 4445.

As an extension to the general case we shall introduce a special version of exponentiation with a power of two 2^k as exponent. From the above considerations we know that this function can be implemented very easily by means of k-fold exponentiation. The exponent 2^k will be specified by k.

Function:	modular exponentiation with exponent a power of 2
Syntax:	int mexp2_l (CLINT a_l, USHORT k, CLINT p_l, CLINT m_l);
Input:	a_l (base) k (exponent of 2) m_l (modulus)
Output:	p_l (residue of mod m_l)
Return:	E_CLINT_OK if all is ok E_CLINT_DBZ if division by 0

int
mexp2_l (CLINT a_l, USHORT k, CLINT p_l, CLINT m_l)
{
  CLINT tmp_l;

  if (EQZ_L (m_l))
    {
      return E_CLINT_DBZ;
    }

If k > 0, then a_l is squared k times modulo m_l.

  if (k > 0)
    {
      cpy_l (tmp_l, a_l);
      while (k-- > 0)
         {
           msqr_l (tmp_l, tmp_l, m_l);
         }
      cpy_l (p_l, tmp_l);
    }
  else

Otherwise, if k = 0, we need only to reduce modulo m_l.

    {
      mod_l (a_l, m_l, p_l);
    }
  return E_CLINT_OK;
}