Speeding up modular exponentiation using CRT

This is the fifth post in a series of posts on the Chinese remainder theorem (CRT). When solving a congruence equation with a composite modulus, it is often easier to convert the problem to one of solving several congruence equations with smaller moduli that are primes or powers of primes. Then combine the individual solutions using the Chinese remainder theorem. In this post, we demonstrate this process for modular exponentiation x \equiv c^d \ (\text{mod} \ m) where the exponent d and the modulus m are large. Given x \equiv c^d \ (\text{mod} \ m), the algorithm discussed here is to produce a system of equations with smaller moduli and exponents that give the same answer as for the original problem.

The previous posts in the series on CRT: first post; second post; third post; fourth post

____________________________________________________________________________

Preliminary discussion

The exponentiation c^d is usually programmed using the fast powering algorithm. The CRT method will convert the exponentiation c^d to several exponentiations that involve much smaller exponents and moduli, thus greatly reducing the calculation time, in particular speeding up the fast powering algorithm. One application of the CRT method is to improve the run time of the decryption process in the RSA algorithm (up to four times faster).

As already mentioned, the goal of the algorithm discussed here is to produce an equivalent system of linear congruence equations. Once this system is produced, the Chinese remainder theorem only guarantees a solution and does not actually produce a solution. So we need to know how to Chinese remainder, i.e. using an algorithm for solving a simultaneous linear congruences. We can use the one discussed here (the first method) or here (the second method). The following examples are worked using the second method.

The CRT algorithm discussed here makes use of Euler’s theorem, which states that a^{\phi(m)} \equiv 1 \ (\text{mod} \ m) whenever a and the modulus m are relatively prime where \phi(m) is the phi function evaluated at m. For this reason, the algorithm requires the evaluation of the phi function.

When calculating c^d \ (\text{mod} \ m), the use of the phi function \phi(m) is to reduce the exponent d by the largest multiple of \phi(m). For example, since \phi(17)=16 and 2^{16} \equiv 1 \ (\text{mod} \ 17), the problem of 2^{250} \ (\text{mod} \ 17) is converted to finding 2^{10} \ (\text{mod} \ 17). Here, the original exponent of 250 is reduced to 10 after taking out the largest multiple of 16. Note that 10 \equiv 250 \ (\text{mod} \ \phi(17)=16). In general, we want to replace the original exponent d by a smaller exponent d_1 where d_1 \equiv d \ (\text{mod} \ \phi(m)).

____________________________________________________________________________

Examples

In the following three examples, we use the algorithm discussed here to solve systems of linear congruence equations (this is the iterative approach). These examples are by no means realistic since the numbers used are small. So they are for demonstration of how CRT works.

Example 1
Calculate x \equiv 2^{3163} \ (\text{mod} \ 3969).

First, factor the modulus 3969=3^4 \times 7^2=81 \times 49. Now the problem is converted to solving the following system of two equations:

    x \equiv 2^{3163} \ (\text{mod} \ 81)

    x \equiv 2^{3163} \ (\text{mod} \ 49)

By CRT, any solution to the two equations is also a solution to the original equation. However, the exponent of 3163 should first be reduced. To this end, calculate the phi function, where \phi(3^4)=3^2 \cdot (3-1)=54 and \phi(7^2)=7 \cdot (7-1)=42. We should reduce from 3163 the largest multiple of 54 in the first equation and reduce the largest multiple of 42 in the second equation. In other words, reduce the exponent 3163 modulo the two phi function values:

    3163 \equiv 31 \ (\text{mod} \ 54)

    3163 \equiv 13 \ (\text{mod} \ 42)

As a result, we solve the following two equations:

    x \equiv 2^{3163} \equiv 2^{31} \equiv 65 \ (\text{mod} \ 81)

    x \equiv 2^{3163} \equiv 2^{13} \equiv 9 \ (\text{mod} \ 49)

Note that the original exponentiation 2^{3163} is turned into the easier ones of 2^{31} and 2^{13}. The following gives the solution to the above two equations.

    \displaystyle \begin{aligned} x_0&=65+81 \cdot 23 \cdot (9-65) \ (\text{mod} \ 3969)  \\&\equiv -104263 \ (\text{mod} \ 3969) \\&\equiv 2900 \ (\text{mod} \ 3969) \end{aligned}

    where 23 is obtained by solving for y in 81y \equiv 1 \ (\text{mod} \ 49)

By CRT, the answer to the original problem is 2^{3163} \equiv 2900 \ (\text{mod} \ 3969). \square

Example 2
Calculate x \equiv 3^{3163} \ (\text{mod} \ 3969).

In this example, the number 3 and the modulus 3969 are not relatively prime. The CRT method still applies. As in Example 1, the problem can be reduced in the following way:

    x \equiv 3^{3163} \equiv 3^{31} \equiv 0 \ (\text{mod} \ 81)

    x \equiv 2^{3163} \equiv 3^{13} \equiv 10 \ (\text{mod} \ 49)

The first equation is congruent to 0 since 3^{31} contains 81 as a factor. The following gives the solution to the above two equations:

    \displaystyle \begin{aligned} x_0&=0+81 \cdot 23 \cdot (10-0) \ (\text{mod} \ 3969)  \\&\equiv 18630 \ (\text{mod} \ 3969) \\&\equiv 2754 \ (\text{mod} \ 3969) \end{aligned}

By CRT, the answer to the original problem is 3^{3163} \equiv 2754 \ (\text{mod} \ 3969). \square

Example 3
The above 2 examples use small numbers to illustrate the CRT technique. In this example, we use slightly larger numbers. Calculate x \equiv 355^{d} \ (\text{mod} \ m) where d=\text{1,759,695,794} and m=\text{3,055,933,789}=1277 \cdot 1439 \cdot 1663.

As in the other examples, we break up the problem in three congruences. The three factors of the modulus are prime numbers. Thus we reduce the exponent d by multiples of a prime factor less one.

    \text{1,759,695,794} \equiv 1198 \ (\text{mod} \ 1276)

    \text{1,759,695,794} \equiv 814 \ (\text{mod} \ 1438)

    \text{1,759,695,794} \equiv 110 \ (\text{mod} \ 1662)

Then the original problem is transformed to solving the following three equations.

    x \equiv 355^{d} \equiv 355^{1198} \equiv 189 \ (\text{mod} \ 1277)

    x \equiv 355^{d} \equiv 355^{814} \equiv 1010 \ (\text{mod} \ 1439)

    x \equiv 355^{d} \equiv 355^{110} \equiv 315 \ (\text{mod} \ 1663)

Notice that the original exponentiation is transformed to the smaller ones of 355^{1198}, 355^{814} and 355^{315}. The remaining task is to solve the system of three equations. One way to find to solution to the above three equations is to use the iterative approach, starting with the solution x_1=189 to the first equation. Then find the solution x_2 to the first two equations and then the solution x_3 to all three equations.

    \displaystyle \begin{aligned} x_2&=189+1277 \cdot 151 \cdot (1010-189) \ (\text{mod} \ 1277 \cdot 1439)  \\&\equiv 158311156 \ (\text{mod} \ 1837603) \\&\equiv 277298 \ (\text{mod} \ 1837603) \end{aligned}

    where 151 is obtained by solving for y in 1277y \equiv 1 \ (\text{mod} \ 1439)

    \displaystyle \begin{aligned} x_3&=277298+(1277 \cdot 1439) \cdot 970 \cdot (315-277298) \ (\text{mod} \ 1277 \cdot 1439 \cdot 1663)  \\&\equiv 277298+1782474910 \cdot (-276983) \ (\text{mod} \ m) \\&\equiv 277298+1414954310 \ (\text{mod} \ m) \\&\equiv 1415231608 \ (\text{mod} \ m) \end{aligned}

    where 970 is obtained by solving for y in (1277 \cdot 1439) y \equiv 1 \ (\text{mod} \ 1663)

By CRT, the answer to the original problem is 355^{d} \equiv 1415231608 \ (\text{mod} \ m). \square

____________________________________________________________________________

The CRT algorithm to speed exponentiation

Suppose we wish to evaluate x \equiv c^{d} \ (\text{mod} \ m) where the prime factorization of m is m=p_1^{n_1} \cdot p_2^{n_2} \cdots p_t^{n_t}. The numbers p_i are distinct prime numbers and each exponent n_i \ge 1. To prepare for the calculation, do the following:

    Let m_i=p_i^{n_i} for each i.

    Calculate \phi(m_i) for each i.

Case 1. The base c and the modulus m are relatively prime.
Then the answer to the original exponentiation problem is identical to the solution to the following system of t congruence equations:

    x \equiv c^{d_1} \ (\text{mod} \ m_1)

    x \equiv c^{d_2} \ (\text{mod} \ m_2)

    \cdots

    x \equiv c^{d_t} \ (\text{mod} \ m_t)

where d_i \equiv d \ (\text{mod} \ \phi(m_i)) for each i. If possible, each c^{d_i} should be reduced modulo m_i.

Case 2. The base c and the modulus m are not relatively prime.
In this case, c and m have prime factors in common (at least one p_i). The idea here is that for any p_i that is a prime factor of c, the equation x \equiv c^{d_i} \ (\text{mod} \ m_i) in Case 1 is replaced by x \equiv 0 \ (\text{mod} \ m_i). Then solve the resulting system of equations (see Example 2). Essentially, Case 2 can fall under Case 1 with c^{d_i} being congruent to zero. We call out Case 2 for the sake of clarity.

The original exponentiation c^d boils down to solving an appropriate system of CRT congruences as described above. Once the equivalent system of congruences is set up, use the algorithm discussed here or here to do Chinese remaindering.

Comment
Both Case 1 and Case 2 produce a system of linear congruence equations that have identical solution to the original equation. This is a result of using CRT (see Theorem D and Theorem G here). The savings in the calculation come in the form of the smaller exponentiations in the resulting congruence equations.

In Case 2, some of the congruence equations are x \equiv 0. This is because the base c and some moduli m_i are not relatively prime. For these moduli, c^{d_i} would contain p_i^{n_i} as a factor (assuming that d_i \ge n_i). Hence, x \equiv c^{d_i} \equiv 0 \ (\text{mod} \ m_i).

The two-equation case
When the modulus is the product of two factors that are relatively prime, the CRT algorithm involves only two equations. We write out the solution explicitly for this case. To evaluate x \equiv c^{d} \ (\text{mod} \ m) where m=m_1 \cdot m_2 and m_1 and m_2 are relatively prime, solve the following system of two equations,

    x \equiv c^{d_1} \ (\text{mod} \ m_1)

    x \equiv c^{d_2} \ (\text{mod} \ m_2)

The solution is x \equiv  c^{d_1}+m_1 \cdot v_1 \cdot (c^{d_2}-c^{d_1}) \ (\text{mod} \ m) where v_1 is the multiplicative inverse of m_1 modulo m_2. If possible, each c^{d_i} should be reduced modulo m_i.

____________________________________________________________________________

RSA application

The algorithm to speed the exponentiation is possible because the factorization of the modulus m is known (as a result, the values of \phi(m_i) are known). Knowing the values of \phi(m_i) makes it possible to reduce the large exponent d. For this reason, the decryption process in the RSA algorithm is a perfect place to apply the CRT technique described here.

With the RSA cryptosystem, a public key consists of N and e, where N is a large modulus that is a product of two large primes p and q (the two primes are not published) and e is the encryption exponent. Say Bob is the originator of the RSA public key. Bob also generates a private key, which is a number d that is used for decrypting any messages that he receives. The number d must be kept private. The prime factors p and q must also be kept secret since knowing p and q can derive d.

Suppose that Alice has a message to send to Bob. She can do so using the published key of N and e through the exponentiation c \equiv m^e \ (\text{mod} \ N). Here, m is the plaintext (the message to be sent) and c is the ciphertext (the encrypted message). Upon receiving the ciphertext c, Bob can then decrypt through the exponentiation m \equiv c^d \ (\text{mod} \ N) where d is the decryption exponent. In realistic RSA calculation, the public modulus N and the private decryption exponent d are large integers (N is at minimum a 2048-bit number). With CRT, the decryption can be reduced to two much smaller exponentiations. The effect can be at least four times faster.

We illustrate this with an example. This is a toy example since the numbers used are small. It is only intended as an illustration.

Example 4
Suppose the public key consists of N=\text{17,086,049} and e=65537. Bob has the additional information of N=p \cdot q where p=3863 and q=4423, which are kept secret. Knowing p and q allows Bob to compute d=\text{5,731,241}. Suppose that Bob receives a message c= 4831984 from Alice. Use the CRT approach to find the plaintext m.

The exponentiation is m \equiv 4831984^{5731241} \ (\text{mod} \ N), which is equivalent to the following two equations by CRT:

    m \equiv 4831984^{5731241} \ (\text{mod} \ 3863)

    m \equiv 4831984^{5731241} \ (\text{mod} \ 4423)

which is further simplified to:

    m \equiv 4831984^{5731241} \equiv 4831984^{33} \equiv 3084 \ (\text{mod} \ 3863)

    m \equiv 4831984^{5731241} \equiv 4831984^{329} \equiv 1436 \ (\text{mod} \ 4423)

where 33 \equiv 5731241 \ (\text{mod} \ 3862) and 329 \equiv 5731241 \ (\text{mod} \ 4422). Then the following is the plaintext (the original message).

    \displaystyle \begin{aligned} m&=3084+3863 \cdot 1319 \cdot (1436-3084) \ (\text{mod} \ 3863 \cdot 4423)  \\&\equiv 3084+5095297 \cdot (-1648) \ (\text{mod} \ N) \\&\equiv 9289736 \ (\text{mod} \ N) \end{aligned}

    where 1319 is obtained by solving for y in 3683y \equiv 1 \ (\text{mod} \ 4423)

The answer is 4831984^{5731241} \equiv 9289736 \ (\text{mod} \ N). \square

____________________________________________________________________________

Closing comment

In conclusion, we state the explcit formula for providing the CRT answer to the RSA decryption.

    m \equiv  c^{d_p}+p \cdot p_{inv} \cdot (c^{d_q}-c^{d_p}) \ (\text{mod} \ p \cdot q) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1)

where d_p \equiv d \ (\text{mod} \ p-1), d_q \equiv d \ (\text{mod} \ q-1) and p_{inv} is the multiplicative inverse of p modulo q.

The decryption formula of (1) represends tremendous saving in calculation (up to four times faster). It is possible only for the holder of the RSA private key. It requires knowledge of the decryption exponent d, which is calculated from the factors p and q of the modulus N.
____________________________________________________________________________
\copyright \ 2015 \text{ by Dan Ma}

Advertisements

More about checking for primitive roots

Finding primitive roots modulo a number is of great interest in number theory, both in a theoretical standpoint and in a computational standpoint. In this post we compare and contrast three different ways of checking for primitive roots, continuing a discussion in an earlier post An elementary algorithm for finding primitive roots.

___________________________________________________________________________________________________________________

Background

Let m be a positive integer. Let a be a positive integer that is relative prime to m. Let \phi be Euler’s phi function, which counts the number of least residues that are relatively prime to the modulus. For example, \phi(6)=2 as 1 and 5 are the only numbers relatively prime to 6 (out of the numbers 0,1,2,3,4,5). Furthermore, \phi(p)=p-1 for any prime number p. Previous posts on the phi function: Euler’s phi function, part 1 and Euler’s phi function, part 2.

Euler’s theorem tells us that a^{\phi(m)} \equiv 1 \ (\text{mod} \ m). By the order of a modulo m we mean the least positive exponent k such that a^{k} \equiv 1 \ (\text{mod} \ m) (Euler’s theorem indicates that this notion of order is well defined). The number a is said to be a primitive root modulo m if the order of a modulo m is \phi(m).

___________________________________________________________________________________________________________________

Three Checks

Let m be a positive integer. Let a be a positive integer that is relative prime to m. How can we determine whether the number a is a primitive root modulo m? We discuss three ways of answering this question.

    ____________________________________________________________________________________________
    Check # 1

      Check a^j \ (\text{mod} \ m) for all positive integers j<\phi(m).

      If each such congruence \not \equiv 1, then the number a is a primitive root modulo m.

    ____________________________________________________________________________________________

Check # 1 is merely a restatement of the definition of primitive root. It is a dumb test as it requires too much calculation. For large moduli, it would be an inefficient method of checking for primitive roots. The following is a much better test.

    ____________________________________________________________________________________________
    Check # 2

      Check a^j \ (\text{mod} \ m) for all positive divisors j of \phi(m) with j<\phi(m).

      If each such congruence \not \equiv 1, then the number a is a primitive root modulo m.

    ____________________________________________________________________________________________

Check # 2 narrows down the checking by quite a bit – simply checking a^j among the divisors of \phi(m). This works because the only possible numbers for the order modulo m of the number a are the divisors of \phi(m). So we can skip all j that are not divisors of \phi(m). The following lemma shows why this is so. Actually, the lemma proves something more general. It shows that if a^n \equiv 1 \ (\text{mod} \ m), then the order of a must be a divisor of n. Euler’s theorem says that a^{\phi(m)} \equiv 1 \ (\text{mod} \ m). So the order of a must be a divisor of \phi(m).

    Lemma 1

      Let k be the order of a modulo m. If a^n \equiv 1 \ (\text{mod} \ m), then k \ \lvert \ n.

Proof of Lemma 1
We have k \le n since k is least with the property a^k \equiv 1 \ (\text{mod} \ m). By the division algorithm, we have n=q \cdot k+r where q is some integer and 0 \le r <k. We have the following:

    1 \equiv a^{n} \equiv a^{q \cdot k+r} \equiv (a^k)^q \cdot a^r \equiv a^r \ (\text{mod} \ m)

With a^r \equiv 1 \ (\text{mod} \ m) and r < k, it follows that r=0 and n=q \cdot k. Thus k is a divisor of n. \blacksquare

Though Check # 2 is definitely an improvement over Check # 1, the following further narrows the list of exponents to check.

    ____________________________________________________________________________________________
    Check # 3

      Find all prime divisors q of \phi(m). Then compute \displaystyle j=\frac{\phi(m)}{q} over all q.

      Check a^j \ (\text{mod} \ m) for all j calculated above.

      If each such congruence \not \equiv 1, then the number a is a primitive root modulo m.

    ____________________________________________________________________________________________

Check # 3 further eliminates the exponents to try when we check a^j \ (\text{mod} \ m). Instead of checking over all the divisors of \phi(m), we only need to try the divisors of the form \displaystyle \frac{\phi(m)}{q} where q is a prime divisor of \phi(m). The following lemma shows why this works.

    Lemma 2

      The number a is a primitive root modulo m if and only if \displaystyle a^{\frac{\phi(m)}{q}} \not \equiv 1 \ (\text{mod} \ m) for all prime divisors q of \phi(m).

Proof of Lemma 2
The direction \Longrightarrow is clear.

To show \Longleftarrow, suppose a is not a primitive root modulo m. Then
\displaystyle a^{t} \equiv 1 \ (\text{mod} \ m) for some t that is a divisor of \phi(m). We have \phi(m)=t \cdot y for some integer y. Let q be a prime factor of y. Then \phi(m)=t \cdot q \cdot b for some integer b. Consider the following derivation.

    \displaystyle 1 \equiv (a^t)^b =(a^{\frac{\phi(m)}{qb}})^b \equiv a^{\frac{\phi(m)}{q}} \ (\text{mod} \ m)

Thus if \displaystyle a^{\frac{\phi(m)}{q}} \not \equiv 1 \ (\text{mod} \ m) for all prime divisors q of \phi(m), then a must be a primitive root modulo m. \blacksquare

___________________________________________________________________________________________________________________

Examples

We now work some examples using Check # 3. The modular arithmetic is done using an online calculator. It can also be done using the fast powering algorithm (discussed in the post Congruence Arithmetic and Fast Powering Algorithm).

Example 1
Consider m=37. Find all primitive roots modulo m=37.

First \phi(37)=36. The divisors of 36 are:

    1,2,3,4,6,9,12,18,36

To use Check # 2, in order to find out if a is a primitive root, we would need to calculate a^j nine times, one for each of the above divisors of \phi(37)=36.

To use Check # 3, only two of these nine divisors are needed. There are two prime divisors of 36, namely 2 and 3. We use \displaystyle \frac{36}{2}=18 and \displaystyle \frac{36}{3}=12. So we check a^{12} and a^{18} modulo 37. The calculation is presented in the following tables.

    \displaystyle \begin{bmatrix} a&\text{ }&a^{12}&\text{ }&a^{18}  \\\text{ }&\text{ }&\text{ } \\ 1&\text{ }&1&\text{ }&1  \\ 2&\text{ }&26&\text{ }&36  \\ 3&\text{ }&10&\text{ }&1  \\ 4&\text{ }&10&\text{ }&1  \\ 5&\text{ }&10&\text{ }&36  \\ 6&\text{ }&1&\text{ }&36 \\ 7&\text{ }&10&\text{ }&1 \\ 8&\text{ }&1&\text{ }&36 \\ 9&\text{ }&26&\text{ }&1 \\ 10&\text{ }&1&\text{ }&1   \end{bmatrix} \ \ \ \ \ \displaystyle \begin{bmatrix} a&\text{ }&a^{12}&\text{ }&a^{18}  \\\text{ }&\text{ }&\text{ } \\ 11&\text{ }&1&\text{ }&1  \\ 12&\text{ }&26&\text{ }&1  \\ 13&\text{ }&10&\text{ }&36  \\ 14&\text{ }&1&\text{ }&36  \\ 15&\text{ }&26&\text{ }&36  \\ 16&\text{ }&26&\text{ }&1 \\ 17&\text{ }&26&\text{ }&36 \\ 18&\text{ }&10&\text{ }&36 \\ 19&\text{ }&10&\text{ }&36 \\ 20&\text{ }&26&\text{ }&36   \end{bmatrix}

    \displaystyle \begin{bmatrix} a&\text{ }&a^{12}&\text{ }&a^{18}  \\\text{ }&\text{ }&\text{ } \\ 21&\text{ }&26&\text{ }&1  \\ 22&\text{ }&26&\text{ }&36  \\ 23&\text{ }&1&\text{ }&36  \\ 24&\text{ }&10&\text{ }&36  \\ 25&\text{ }&26&\text{ }&1  \\ 26&\text{ }&1&\text{ }&1 \\ 27&\text{ }&1&\text{ }&1 \\ 28&\text{ }&26&\text{ }&1 \\ 29&\text{ }&1&\text{ }&36 \\ 30&\text{ }&10&\text{ }&1   \end{bmatrix} \ \ \ \ \ \displaystyle \begin{bmatrix} a&\text{ }&a^{12}&\text{ }&a^{18}  \\\text{ }&\text{ }&\text{ } \\ 31&\text{ }&1&\text{ }&36  \\ 32&\text{ }&10&\text{ }&36  \\ 33&\text{ }&10&\text{ }&1  \\ 34&\text{ }&10&\text{ }&1  \\ 35&\text{ }&26&\text{ }&36  \\ 36&\text{ }&1&\text{ }&1 \\ \text{ }&\text{ }&\text{ }&\text{ }&\text{ } \\ \text{ }&\text{ }&\text{ }&\text{ }&\text{ } \\ \text{ }&\text{ }&\text{ }&\text{ }&\text{ } \\ \text{ }&\text{ }&\text{ }&\text{ }&\text{ }  \end{bmatrix}

The primitive roots are the rows with both congruences \not \equiv 1. They are:

    2, 5, 13, 15, 17, 18, 19, 20, 22, 24, 32, 35

One comment about the above tables. The non-one values in the above table seem to follow a pattern. In the columns for the calculation for a^{18}, the values are either 1 or 36. The non-one value is 36. It turns out that it has order 2 modulo 37. The non-one values in the columns for a^{12} are 10 and 26. It turns out that they have order 3 modulo 37. See the exercise stated below.

Example 2
Consider m=17. Find all primitive roots modulo m=17.

Since \phi(17)=16=2^4, the only prime divisor of $latex \phi(17) is 2. We use \displaystyle \frac{16}{2}=8. For any a, we only need to calculate a^8.

    \displaystyle \begin{bmatrix} a&\text{ }&a^{8}&\text{ }&a&\text{ }&a^{8}  \\\text{ }&\text{ }&\text{ } \\ 1&\text{ }&1&\text{ }&11&\text{ }&16  \\ 2&\text{ }&1&\text{ }&12&\text{ }&16  \\ 3&\text{ }&16&\text{ }&13&\text{ }&1  \\ 4&\text{ }&1&\text{ }&14&\text{ }&16  \\ 5&\text{ }&16&\text{ }&15&\text{ }&1  \\ 6&\text{ }&16&\text{ }&16&\text{ }&1 \\ 7&\text{ }&16&\text{ }&\text{ } \\ 8&\text{ }&1&\text{ }&\text{ } \\ 9&\text{ }&1&\text{ }&\text{ } \\ 10&\text{ }&16&\text{ }&\text{ }     \end{bmatrix}

The primitive roots modulo 17 are:

    3, 5, 6, 7, 10, 11, 12, 14

Note that the non-one value 16 in the above table has order 2 modulo 17. See the exercise below.

___________________________________________________________________________________________________________________

Special Case

Based on Example 2, the following is a special case for Check # 3.

    ____________________________________________________________________________________________
    Check # 3 (A Special Case)

      Let p be a prime such that p-1=2^n for some positive integer n.

      Note that 2 is the only prime divisor of \phi(p)=p-1.

      Check a^j \ (\text{mod} \ p) where \displaystyle j=\frac{p-1}{2}.

      If a^j \not \equiv 1, then the number a is a primitive root modulo p.

    ____________________________________________________________________________________________

___________________________________________________________________________________________________________________

Exercise

This is the exercise mentioned at the end of Example 1.

Let p be a prime number. Let q be a prime divisor of p-1. Let a be an integer with 1 \le a \le p-1. Show that the number \displaystyle a^{\frac{p-1}{q}} is either \equiv 1 or has order q modulo p.

___________________________________________________________________________________________________________________

\copyright \ 2013 \text{ by Dan Ma}

Solving Quadratic Congruences

This post is an introductory discussion on the congruence equations of the form x^2 \equiv a \ (\text{mod} \ p) where the modulus p is an odd prime and the number a is relatively prime to p. A discussion on the related notion of quadratic residues is found here. Specific algorithms for solving quadratic congruence eqautions with odd prime moduli are discussed in this subsequent post.

____________________________________________________________________________

Simple Example

We start off with a simple example. Calculate x^2 modulo m=11 for x=0,1,2,\cdots,10.

    0^2 \equiv 0 \ (\text{mod} \ 11)

    1^2 \equiv 1 \ (\text{mod} \ 11)

    2^2 \equiv 4 \ (\text{mod} \ 11)

    3^2 \equiv 9 \ (\text{mod} \ 11)

    4^2 \equiv 5 \ (\text{mod} \ 11)

    5^2 \equiv 3 \ (\text{mod} \ 11)

    6^2 \equiv 3 \ (\text{mod} \ 11)

    7^2 \equiv 5 \ (\text{mod} \ 11)

    8^2 \equiv 9 \ (\text{mod} \ 11)

    9^2 \equiv 4 \ (\text{mod} \ 11)

    10^2 \equiv 1 \ (\text{mod} \ 11)

The above calculation shows that the values of x^2 modulo m=11 can only be 1,3,4,5,9. So equations such as x^2 \equiv a \ (\text{mod} \ 11) for a=1,3,4,5,9 have solutions. For example, the solutions for the equation x^2 \equiv 5 \ (\text{mod} \ 11) are x=4 and x=7.

On the other hand, the equations x^2 \equiv b \ (\text{mod} \ 11) for b=2,6,7,8,10 have no solutions.

Also note that whenever a \ne 0 and the equation x^2 \equiv a \ (\text{mod} \ 11) has a solution, the solutions come in pairs.

____________________________________________________________________________

Quadratic Congruences

Let m be an odd prime number. Let a be an integer that is not divisible by p (equivalently relatively prime to p). The object of interest here is the quadratic congruence equation x^2 \equiv a \ (\text{mod} \ p). It turns out that each such equation has exactly two solutions whenever the number a and the modulus p are relatively prime (as demonstrated in the above simple example). The following lemma and corollary confirm what we see in the above example.

    Lemma 1

      Let p be an odd prime number. Let a be an integer that is not divisible by p. Then the equation

        x^2 \equiv a \ (\text{mod} \ p) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1)

      has no solutions or exactly two solutions.

Proof of Lemma 1
If equation (1) has no solutions, then we are done. Suppose it has at least one solution, say x=r. We have r^2 \equiv a \ (\text{mod} \ p). It follows that x=-r \equiv p-r \ (\text{mod} \ p) is a solution of equation (1) too.

We claim x=r and x=p-r are distinct modulo p. To see this, suppose p-r \equiv r \ (\text{mod} \ p). Then p \ \lvert \ (p-2r). Because p is an odd prime, p \not \lvert \ 2. So p \ \lvert \ r. This implies that p \ \lvert \ r^2. Because p \ \lvert \ (r^2-a), p \ \lvert \ a, contradicting the assumption that p \not \lvert \ a. Thus p-r \not \equiv r \ (\text{mod} \ p), demonstrating that equation (1) has at least two solutions.

It remains to be shown that any solution of equation (1) must be congruent to one of x=r and x=p-r. Suppose t^2 \equiv a \ (\text{mod} \ p). Then t^2 \equiv r^2 \ (\text{mod} \ p). It follows that p \ \lvert \ (t-r)(t+r). Thus p must divide one of the two factors (Euclid’s lemma). The case p \ \lvert \ (t-r) implies t \equiv r \ (\text{mod} \ p). The case p \ \lvert \ (t+r) implies t \equiv -r \ (\text{mod} \ p). \blacksquare

    Corollary 2

      Let p be an odd prime number. The equation x^2 \equiv a \ (\text{mod} \ p) has exactly two solutions for \displaystyle \frac{p-1}{2} many numbers a \in \left\{1,2,\cdots,p-1 \right\} and has no solutions for the other \displaystyle \frac{p-1}{2} numbers a \in \left\{1,2,\cdots,p-1 \right\}.

Remark
For the even prime p=2, the equation x^2 \equiv a \ (\text{mod} \ 2) is not an interesting one. For x^2 \equiv 0 \ (\text{mod} \ 2), x=0 is the only solution. For x^2 \equiv 1 \ (\text{mod} \ 2), x=1 is the only solution.

For composite moduli, the quadratic equation x^2 \equiv a \ (\text{mod} \ m) can have more than two solutions. For example, x^2 \equiv 1 \ (\text{mod} \ 8) has four solutions x=1,3,5,7.

For these reasons, we only work with odd prime moduli for quadratic congruences.

____________________________________________________________________________

General Case

What about the general case of the quadratic congruence equation of the following form?

    \alpha y^2+\beta y+\gamma \equiv 0 \ (\text{mod} \ p)  \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (2)

Of course, we only consider the equations where \alpha \not \equiv 0 \ (\text{mod} \ p) and p is an odd prime. It turns out that equation (2) can be replaced by an equivalent congruence equation of the same form as equation (1) above. So the general case of equation (2) presents no new problem. We just convert equation (2) to its equivalence and solve it accordingly. We now discuss how this is done.

The coefficient \alpha, the coefficient of the y^2 term, has a multiplicative inverse modulo p. So multiplying equation (2) by \alpha^{-1} gives the following equation.

    y^2+\alpha^{-1} \beta y+\alpha^{-1} \gamma \equiv 0 \ (\text{mod} \ p)  \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (3)

So we can now focus on solving equation (3), which has the same solutions as equation (2). Consider the coefficient of the y term. If the coefficient \alpha^{-1} \beta is even, we can complete the square and obtain an equivalent equation of the same form as equation (1). If the coefficient \alpha^{-1} \beta is odd, then we can add p to it and obtain an even coefficient. We can then proceed to complete the square as in the even case. We demonstrate with two examples.

Consider the equation 3 y^2+4y+1 \equiv 0 \ (\text{mod} \ 11). Since 4 \cdot 3 \equiv 1 \ (\text{mod} \ 11). The multiplicative inverse of 4 is 3. So we multiply 4 across and obtain y^2+16y+4 \equiv 0 \ (\text{mod} \ 11). The coefficient of the y term is even. We complete the square as follows.

    y^2+16y+4 \equiv 0 \ (\text{mod} \ 11)

    y^2+16y+64 \equiv 64-4 \ (\text{mod} \ 11)

    (y+8)^2 \equiv 60 \ (\text{mod} \ 11)

    (y+8)^2 \equiv 5 \ (\text{mod} \ 11)

The last equation in the above derivation is of the form x^2 \equiv 5 \ (\text{mod} \ 11) where the solutions are x=4 and x=7. Thus we have y+8 \equiv 4 \ (\text{mod} \ 11) and y+8 \equiv 7 \ (\text{mod} \ 11). These two congruences give y \equiv 7 \ (\text{mod} \ 11) and y \equiv 10 \ (\text{mod} \ 11).

For the odd case, consider the equation 5 y^2+y+8 \equiv 0 \ (\text{mod} \ 11). The multiplicative inverse of 5 is 9 as 5 \cdot 9 \equiv 1 \ (\text{mod} \ 11). After multiplying by the inverse, we obtain y^2+9y+72 \equiv 0 \ (\text{mod} \ 11). We further reduce 72 modulo 11 to get y^2+9y+6 \equiv 0 \ (\text{mod} \ 11). Note that the coefficient of the y term is odd. So we add modulus to that coefficient to obtain the equation y^2+20y+6 \equiv 0 \ (\text{mod} \ 11). We then proceed to complete the square as follows.

    y^2+20y+100 \equiv 100-6 \ (\text{mod} \ 11)

    (y+10)^2 \equiv 94 \ (\text{mod} \ 11)

    (y+10)^2 \equiv 6 \ (\text{mod} \ 11)

The last equation in the above derivation is of the form x^2 \equiv 6 \ (\text{mod} \ 11), which has no solutions (based on the simple example above). Thus the original equation 5 y^2+y+8 \equiv 0 \ (\text{mod} \ 11) has no solutions.

____________________________________________________________________________

Examples

To solve the quadratic congruence x^2 \equiv a \ (\text{mod} \ p), one way is to compute the entire table of values for x^2 modulo p. For very small prime such as the simple example above, this approach is workable. For large primes, this requires a lot of computational resources.

To further illustrate the quadratic congruences, we work three examples with help from Euler’s Criterion and from using Excel to do some of the calculations.

According to Euler’s Criterion, the equation x^2 \equiv a \ (\text{mod} \ p) has solutions if and only if \displaystyle a^{\frac{p-1}{2}} \equiv 1 \ (\text{mod} \ p). Equivalently, the equation x^2 \equiv a \ (\text{mod} \ p) has no solutions if and only if \displaystyle a^{\frac{p-1}{2}} \equiv -1 \ (\text{mod} \ p). So the solvability of the quadratic congruence equation can be translated as a modular exponentiation calculation.

The computation for \displaystyle a^{\frac{p-1}{2}} \ (\text{mod} \ p) can be done directly using an online modular arithmetic calculator or using the fast-powering algorithm (discussed in the post Congruence Arithmetic and Fast Powering Algorithm). For a discussion and a proof of Euler’s Criterion, see the post Euler’s Criterion.

When Euler’s Criterion indicates there are solutions, how do we find the solutions? We demonstrate using the following examples.

Example 1
Solve x^2 \equiv 5 \ (\text{mod} \ 61).

According to Euler’s Criterion, the equation x^2 \equiv 5 \ (\text{mod} \ 61) has solutions since 5^{30} \equiv 1 \ (\text{mod} \ 61). To find the solutions, we keep adding the modulus to a=5 until we get a perfect square.

    \displaystyle x^2 \equiv 5 \equiv 5+61 \equiv 5+2(61) \equiv \cdots \equiv 5+20(61)=1225=35^2 \ (\text{mod} \ 61)

So we have x^2 \equiv 35^2 \ (\text{mod} \ 61), which gives x=35 and x=-35. The solutions are x \equiv -35 \equiv 26 \ (\text{mod} \ 61) and x \equiv 35 \ (\text{mod} \ 61).

Example 2
Solve x^2 \equiv 899 \ (\text{mod} \ 50261).

Since 899^{25130} \equiv 1 \ (\text{mod} \ 50261), the equation has solutions. We then add the modulus repeatedly to 899 until we get a perfect square (with the aid of an Excel spreadsheet).

    \displaystyle x^2 \equiv 899 \equiv 899+50261 \equiv 899+2(50261) \equiv \cdots \equiv 899+4297(50261)=215972416=14696^2 \ (\text{mod} \ 50261)

So we have x^2 \equiv 14696^2 \ (\text{mod} \ 50261), which gives x=14696 and x=-14696. The solutions are x \equiv 14696 \ (\text{mod} \ 50261) and x \equiv -14696 \equiv 35565 \ (\text{mod} \ 50261).

Example 3
Solve x^2 \equiv 13961 \ (\text{mod} \ 50261).

Since 13961^{25130} \equiv -1 \ (\text{mod} \ 50261), the equation has no solutions according to Euler’s Criterion.

____________________________________________________________________________

\copyright \ 2013 - 2015 \text{ by Dan Ma}
Revised December 9, 2015

Using Fermat’s Little Theorem to Test Primality

Fermat’s little theorem describes a property that is common to all prime numbers. This property can be used as a way to detect the “prime or composite” status of an integer. Primality testing using Fermat’s little theorem is called the Fermat primality test. In this post, we explain how to use this test and to discuss some issues surrounding the Fermat test.

___________________________________________________________________

Describing the test

The Fermat primality test, as mentioned above, is based on Fermat’s little theorem. The following is the statement of the theorem.

    Fermat’s little theorem
    If n is a prime number and if a is an integer that is relatively prime to n, then the following congruence relationship holds:

      a^{n-1} \equiv 1 (\text{mod} \ n) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1)

The above theorem indicates that all prime numbers possess a certain property. Therefore if a given positive integer does not possess this property, we know for certain that this integer is not prime. Suppose that the primality of an integer n is not known. If we can find an integer a that is relatively prime to n such that a^{n-1} \not \equiv 1 \ (\text{mod} \ n), then we have conclusive proof that n is composite. Such a number a is said to be a Fermat witness for (the compositeness of) n.

The Fermat test is closedly linked to the notations of probable primes and pseudoprimes. If the congruence relation (1) is true for n and a, then n is said to be a probable prime to base a. Furthermore, if n happens to be a composite number, then n is said to be a pseudoprime to base a. Pseudoprime prime is a composite number that possesses the prime-like property as indicated by (1) for one base a.

The Fermat primality test from a compositeness perspective is about looking for Fermat witnesses. If a Fermat witness is found, the number being tested is proved to be composite. On the other hand, the Fermat primality test, from a primality perspective, consists of checking the congruence relation (1) for several bases that are randomly selected. If the number n is found to be a probable prime to all the randomly chosen bases, then n is likely a prime number.

If the number n is in reality a prime number, then the Fermat test will always give the correct result (as a result of Fermat’s little theorem). If the number n is in reality a composite number, the Fermat test can make the mistake of identifying the composite number n as prime (i.e. identifying a pseudoprime as a prime). For most composite numbers this error probability can be made arbitrarily small (by testing a large number of bases a). But there are rare composite numbers that evade the Fermat test. Such composite numbers are called Carmichael numbers. No matter how many bases you test on a Carmichael number, the Fermat test will always output Probably Prime. Carmichael numbers may be rare but there are infinitely many of them over the entire number line. More about Carmichael numbers below.

The following describes the steps of the Fermat primality test.

    Fermat primality test
    The test is to determine whether a large positive integer n is prime or composite. The test will output one of two results: n is Composite or n is Probably Prime.

  • Step 1. Choose a random integer a \in \left\{2,3,\cdots,n-1 \right\}.
  • Step 2. Compute \text{GCD}(a,n). If it is greater than 1, then stop and output n is Composite. Otherwise go to the next step.
  • Step 3. Compute a^{n-1} \ (\text{mod} \ n).
    • If a^{n-1} \not \equiv 1 \ (\text{mod} \ n), then stop and output n is Composite.
    • If a^{n-1} \equiv 1 \ (\text{mod} \ n), then n may be a prime number. Do one of the following:
      • Return to Step 1 and repeat the process with a new a.
      • Output n is Probably Prime and stop.

    \text{ }

The exponentiation in Step 3 can be done by the fast powering algorithm. This involves a series of squarings and multiplications. Even for numbers that have hundreds of digits, the fast powering algorithm is efficient.

One comment about Step 2 in the algorithm. Step 2 could be called the GCD test for primality. If you can find an integer a such that 1<a<n and such that \text{GCD}(a,n) \ne 1, then the integer n is certainly composite. Such a number a is called a GCD witness for the compositeness of n. So the Fermat test as described above combines the GCD test and the Fermat test. We can use the Euclidean algorithm to find the GCD. If we happen to stumble upon a GCD witness, then we can try another n for a candidate of a prime number. For most composite numbers, it is not likely to stumble upon a GCD witness. Thus when using the Fermat test, it is likely that Step 3 in the algorithm is used.

An example of Fermat primality testing is the post called A primality testing exercise from RSA-100.

____________________________________________________________________________

More about the test

When using the Fermat test, what is the probability of the test giving the correct result? Or what is the probability of making an error? Because the Fermat test is not a true probabilistic primality test, the answers to these questions are conditional. In one scenario which covers most of the cases, the test works like an efficient probabilistic test. In another scenario which occurs very rarely, the Fermat test fails miserably.

As with most diagnostic tests, the Fermat test can make two types of mistakes – false positives or false negatives. For primality testing discussed in this post, we define a positive result as the outcome that says the number being tested is a prime number and a negative result as the outcome that says the number being tested is a composite number. Thus a false positive is identifying a composite number as a prime number and a false negative is identifying a prime number as a composite number.

For the Fermat test, there is no false negative. If n is a prime number in reality, the statement of Fermat’s little theorem does not allow the possibility that n be declared a composite number. Thus if the Fermat test gives a negative result, it would be a true negative. In other words, finding a Fermat witness for n is an irrefutable proof that n is composite.

However, there can be false positives for the Fermat test. This is where things can get a little tricky. A composite number n is said to be a Carmichael number if the above congruence relationship (1) holds for all bases a relatively prime to n. In other words, n is a Carmichael number if a^{n-1} \equiv 1 (\text{mod} \ n) for all a that are relatively prime to n. Saying it in another way, n is a Carmichael number if there exists no Fermat witness for n.

The smallest Carmichael number is 561. Carmichael numbers are rare but there are infinitely many of them. The existence of such numbers poses a challenge for the Fermat test. If you apply the Fermat test on a Carmichael number, the outcome will always be Probably Prime. So the Fermat test will always give a false positive when it is applied on a Carmichael number. To put it in another way, with respect to Carmichael numbers, the error probability of the Fermat test is virtually 100%!

So should a primality tester do? To keep things in perspective, Carmichael numbers are rare (see this post). If the primality testing is done on randomly chosen numbers, choosing a Carmichael number is not likely. So the Fermat test will often give the correct results. For those who are bothered by the nagging fear of working with Carmichael numbers, they can always switch to a Carmichael neutral test such as the Miller-Rabin test.

___________________________________________________________________

One bright spot about the Fermat test

There is one bright spot about the Fermat test. When applying the Fermat test on numbers that are not Carmichael numbers, the error probability can be made arbitrarily small. In this sense the Fermat test works like a true probabilistic primality test. Consider the following theorem.

    Theorem 1
    Let n be a composite integer such that it is not a pseudoprime to at least one base (i.e. n has a Fermat witness). In other words, n is not a Carmichael number. Then n is not a pseudoprime to at least half of the bases a (1<a<n) that are relatively prime to n. In other words, n is a pseudoprime to at most half of the bases a (1<a<n) that are relatively prime to n.

Theorem 1 means that the Fermat test can be very accurate on composite numbers that are not Carmichael numbers. As long as there is one base to which the composite number is not a pseudoprime (i.e. as long as there is a Fermat witness for the composite number in question), there will be enough of such bases (at least 50% of the possible bases). As a result, it is likely that the Fermat test will find a witness, especially if the tester is willing to use enough bases to test and if the bases are randomly chosen. When a base is randomly chosen, there is at least a 50% chance that the number n is not a pseudoprime to that base (i.e. the Fermat test will detect the compositeness) or putting it in another way, there is at most a 50% chance that the Fermat test will not detect the compositeness of the composite number n. So if k values of a are randomly selected, there is at most 0.5^k probability that the Fermat test will not detect the compositeness of the composite number n (i.e. making a mistake). So the probability of a false positive is at most 0.5^k. For a large enough k, this probability is practically zero.

Proof of Theorem 1
A base to which n is a pseudoprime or not a pseudoprime should be a number in the interval 1<a<n that is relatively prime to n. If n is a pseudoprime to base a, then a raised to some power is congruent to 1 modulo n. For this to happen, a must be relatively prime to the modulus n. For this reason, when we consider a base, it must be a number that is relatively prime to the composite integer n (see the post on Euler’s phi function).

Let a be a base to which n is not a pseudoprime. We make the following claim.

    Claim
    If b is a number such that 1<b<n and such that n is a pseudoprime to base b, then n is not a pseudoprime to base a \cdot b.

Since both integers a and b are assumed to be relatively prime to n, the product a \cdot b is also relatively prime to n (see Lemma 4 in this post). Now consider the congruence (ab)^{n-1} \ (\text{mod} \ n), which is derived as follows:

    (ab)^{n-1} \equiv a^{n-1} \cdot b^{n-1} \equiv a^{n-1} \not \equiv 1 \ (\text{mod} \ n)

In the above derivation, we use the fact that n is not a pseudoprime to base a and n is a pseudoprime to base b. The above derivation shows that n is not a pseudoprime to base ab.

If n is not a pseudoprime to all bases in 1<w<n, then we are done. So assume that n is a pseudoprime to at least one base. Let b_1,b_2,\cdots,b_k enumerate all bases to which n is a pseudoprime. We assume that the b_j are all distinct. So b_i \not \equiv b_j \ (\text{mod} \ n) for all i \ne j. By the above claim, the composite number n is not a pseudoprime to all the following k numbers:

    a \cdot b_1, \ a \cdot b_2, \cdots, \ a \cdot b_k

It is also clear that a \cdot b_i \not \equiv a \cdot b_j \ (\text{mod} \ n) for i \ne j. What we have just shown is that there are at least as many bases to which n is not a pseudoprime as there are bases to which n is a pseudoprime. This means that n is not a pseudoprime to at least 50% of the bases that are relatively prime to n. In other words, as long as there exists one Fermat witness for n, at least 50% of the bases are Fermat witnesses for n. It then follows that n is a pseudoprime to no more than 50% of the bases relatively prime to n. \blacksquare

There is another way to state Theorem 1. Recall that Euler’s phi function \phi(n) is defined to be the number of integers a in the interval 1<a<n that are relatively prime to n. With this in mind, Theorem 1 can be restated as the following:

    Corollary 2
    Let n be a composite integer such that it is not a pseudoprime to at least one base. Then n is not a pseudoprime to at least \displaystyle \frac{\phi(n)}{2} many bases in the interval 1<a<n.

___________________________________________________________________

Concluding remarks

Of course, Theorem 1 works only for the composite numbers that are not pseudoprime to at least one base (i.e. they are not Carmichael numbers). When you test the compositeness of a number, you do not know in advance if it is Carmichael or not. On the other hand, if the testing is done on randomly chosen numbers, it is not likely to randomly stumble upon Carmichael numbers. The Fermat test works well for the most part and often give the correct results. If one is concerned about the rare chance of a false positive in the form of a Carmichael number, then the Miller-Rabin test will be a good alternative.

Note:
The original post was written in August 10, 2013. On March 29, 2015, this post is replaced with a blog post called The Fermat primality test from the companion math crypto blog.

___________________________________________________________________
\copyright \ \ 2013 - 2015 \ \text{Dan Ma}

Congruence Arithmetic and Fast Powering Algorithm

In some cryptography applications such as RSA algorithm, it is necessary to compute \displaystyle a^w modulo m where the power w and the modulus m are very large numbers. We discuss and demonstrate an efficient algorithm that can handle such calculations. This general algorithm has various names such as fast powering algorithm, square-and-multiply algorithm and exponentiation by squaring.

The problem at hand is to compute a^w \ (\text{mod} \ m). The naïve approach is to compute by repeatedly multiplying by a and reducing modulo m. When the power w is large (e.g. an integer with hundreds of digits), this approach is difficult or even impossible (given the current technology). In this post we discuss an alternative that is known as the fast powering algorithm.

___________________________________________________________________________________________________________________

An Example

Compute 1286^{1171} modulo 1363.

Using the naïve approach described earlier, we would do something like the following:

    \displaystyle \begin{aligned} 1286^{1171} &=(1286 \cdot 1286) \cdot 1286^{1269} \equiv 477 \cdot 1286^{1269} \ \ \ \ \ \ \ \ \ \ \text{mod} \ 1363 \\&=(477 \cdot 1286) \cdot 1286^{1268} \equiv 72 \cdot 1286^{1268} \  \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 1363 \\&=(72 \cdot 1286) \cdot 1286^{1267} \equiv 1271 \cdot 1286^{1267} \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 1363 \\&\text{     }\cdots \\&\text{     }\cdots \\&\text{     }\cdots \end{aligned}

In this naïve approach, we would multiply two numbers at a time and then reduce the result modulo 1363 so that the numbers do not get too large. The above example would involve 1170 multiplications and then 1170 divisions for the reduction modulo 1363. Great difficulty comes when the power is not 1171 and instead is an integer with hundreds or even thousands of digits.

Note that in the above naïve approach, the power is reduced by one at each step. In the fast power alternative, the power comes down by an exponent of two in each step. The idea is to use the binary expansions of the exponent 1171 to transform the computation of 1286^{1171} into a series of squarings and multiplications. To this end, we write 1171 as a sum of powers of two as follows:

    (1) \ \ \ \ \ \ \ \ \ 1171=2^0+2^1+2^4+2^7+2^{10}

Next we compute 1286^{2^0},1286^{2^1},1286^{2^2},1286^{2^3},\cdots,1286^{2^{10}} modulo 1363. Note that each term is the square of the preceding term, hence the word square in the name “square-and-multiply”. To make it easier to see, we put the results in the following table.

    \displaystyle (2) \ \ \ \ \ \ \ \ \ \begin{bmatrix} \text{ i }&\text{ }&1286^{2^i}&\text{ }&\text{squaring}&\text{ }&\text{modulo } 1363&\text{ }&\text{ }  \\\text{ }&\text{ }&\text{ } \\ 0&\text{ }&1286^{2^0}&\text{ }&\text{ }&\text{ }&\equiv 1286&\text{ }&*  \\ 1&\text{ }&1286^{2^1}&\text{ }&\equiv 1286^2&\text{ }&\equiv 477&\text{ }&*  \\ 2&\text{ }&1286^{2^2}&\text{ }&\equiv 477^2&\text{ }&\equiv 1271&\text{ }&\text{ }  \\ 3&\text{ }&1286^{2^3}&\text{ }&\equiv 1271^2&\text{ }&\equiv 286&\text{ }&\text{ }  \\ 4&\text{ }&1286^{2^4}&\text{ }&\equiv 286^2&\text{ }&\equiv 16&\text{ }&*  \\ 5&\text{ }&1286^{2^5}&\text{ }&\equiv 16^2&\text{ }&\equiv 256&\text{ }&\text{ } \\ 6&\text{ }&1286^{2^6}&\text{ }&\equiv 256^2&\text{ }&\equiv 112&\text{ }&\text{ } \\ 7&\text{ }&1286^{2^7}&\text{ }&\equiv 112^2&\text{ }&\equiv 277&\text{ }&* \\ 8&\text{ }&1286^{2^8}&\text{ }&\equiv 277^2&\text{ }&\equiv 401&\text{ }&\text{ } \\ 9&\text{ }&1286^{2^9}&\text{ }&\equiv 401^2&\text{ }&\equiv 1330&\text{ }&\text{ } \\ 10&\text{ }&1286^{2^{10}}&\text{ }&\equiv 1330^2&\text{ }&\equiv 1089&\text{ }&*  \end{bmatrix}

Note that the rows marked by * in the above table are the results that we need. In the above table, there are 10 multiplications for the squarings and 10 divisions for the reduction modulo 1363.

Now 1286^{1171} is calculated as follows:

    \displaystyle \begin{aligned} (3) \ \ \ \ \ \ \ \ \ 1286^{1171} &=1286^{2^0} \cdot 1286^{2^1} \cdot 1286^{2^4} \cdot1286^{2^7} \cdot 1286^{2^{10}} \\&\equiv 1286 \ \ \cdot 477 \ \ \ \ \cdot 16 \ \ \ \ \ \cdot 277 \ \ \ \ \cdot 1089 \ \ \text{mod} \ 1363 \\&\equiv 72 \ \ \ \ \cdot 16 \ \ \ \ \ \cdot 277 \ \ \ \ \cdot 1089 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 1363 \\&\equiv 1152 \ \ \cdot 277 \ \ \cdot 1089 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 1363 \\&\equiv 162 \ \ \cdot 1089 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 1363 \\&\equiv 591 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 1363 \end{aligned}

We have the answer 1286^{1171} \equiv 591 \ (\text{mod} \ 1363). The calculation in (3) is the step that gives the word “multiply” in the name “square-and-multiply”. In this step, we multiply the results obtained from the previous step.

We now tally up the amount of work that is done. The calculation in table (2) requires 10 multiplications for the squaring and 10 divisions for the reduction modulo 1363. The calculation in (3) requires 4 multiplications and 4 divisions for the reduction modulo 1363. Together, there are 14 multiplications and 14 divisions. In contrast, the naïve approach would require 1170 multiplications and 1170 divisions!

___________________________________________________________________________________________________________________

Description of Fast Powering Algorithm

To compute a^w \equiv (\text{mod} \ m), use the following steps. The following steps correspond with the steps in the above example.

Step (1)

Compute the binary expansions of the power w.

    \displaystyle w=C_0 +C_1 \cdot 2^{1}+C_2 \cdot 2^{2}+\cdots+C_{k-1} \cdot 2^{k-1}+C_k \cdot 2^{k}

where each j, C_j=0 or C_j=1. In particular, we assume that C_k=1.

Step (2)

For each j=0,1,2,\cdots,k, compute \displaystyle a^{2^j} \equiv A_j modulo m. Note that each \displaystyle a^{2^j} \equiv A_j is the result of squaring the previous term \displaystyle a^{2^{j-1}} \equiv A_{j-1}. We arrange the calculation in the following table.

    \displaystyle (2) \ \ \ \ \ \ \ \ \ \begin{bmatrix} \text{ i }&\text{ }&a^{2^i}&\text{ }&\text{squaring}&\text{ }&\text{modulo } m&\text{ }&\text{ }  \\\text{ }&\text{ }&\text{ } \\ 0&\text{ }&a^{2^0}&\text{ }&\text{ }&\text{ }&A_0\equiv a&\text{ }&\text{ }  \\ 1&\text{ }&a^{2^1}&\text{ }&\equiv A_0^2&\text{ }&\equiv A_1&\text{ }&\text{ }  \\ 2&\text{ }&a^{2^2}&\text{ }&\equiv A_1^2&\text{ }&\equiv A_2&\text{ }&\text{ }  \\ 3&\text{ }&a^{2^3}&\text{ }&\equiv A_2^2&\text{ }&\equiv A_3&\text{ }&\text{ }  \\ \text{ }&\text{ }&\cdots&\text{ }&\cdots&\text{ }&\cdots&\text{ }&\text{ }  \\ \text{ }&\text{ }&\cdots&\text{ }&\cdots&\text{ }&\cdots&\text{ }&\text{ } \\ \text{ }&\text{ }&\cdots&\text{ }&\cdots&\text{ }&\cdots&\text{ }&\text{ } \\ k-1&\text{ }&a^{2^{k-1}}&\text{ }&\equiv A_{k-2}^2&\text{ }&\equiv A_{k-1}&\text{ }&\text{ } \\ k&\text{ }&a^{2^k}&\text{ }&\equiv A_{k-1}^2&\text{ }&\equiv A_k&\text{ }&\text{ }   \end{bmatrix}

Step (3)

Compute a^w \equiv (\text{mod} \ m) using the following derivation.

    \displaystyle \begin{aligned}(3) \ \ \ \ \ \ \ \ \  a^{w}&=a^{C_0 +C_1 \cdot 2^{1}+C_2 \cdot 2^{2}+\cdots+C_{k-1} \cdot 2^{k-1}+C_k \cdot 2^{k}} \\&=a^{C_0} \cdot a^{C_1 \cdot 2^{1}} \cdot a^{C_2 \cdot 2^{2}} \cdots a^{C_{k-1} \cdot 2^{k-1}} \cdot a^{C_k \cdot 2^{k}} \\&=a^{C_0} \cdot (a^{2^{1}})^{C_1} \cdot (a^{2^{2}})^{C_2} \cdots (a^{2^{k-1}})^{C_{k-1}} \cdot (a^{2^{k}})^{C_k} \\&\equiv A_0^{C_0} \cdot (A_1)^{C_1} \cdot (A_2)^{C_2} \cdots (A_{k-1})^{C_{k-1}} \cdot (A_k)^{C_k} \ \ \ \ \ (\text{mod} \ m) \end{aligned}

The last line in (3) is to be further reduced modulo m. In the actual calculation, only the terms with C_j=1 need to be used.

We now establish an upper bound for the number multiplications. Step (2) requires k multiplications and k divisions to reduce modulo m. Step (3) requires at most k many multiplications since some of the C_j many be zero. Step (3) also requires at most k many divisions to reduce modulo m. So altogether, the algorithm requires at most 2k multiplications and 2k divisions.

From Step (1), we know that \displaystyle 2^k \le w. Take natural log of both sides, we have \displaystyle k \le \frac{\text{ln}(w)}{\text{ln}(2)} and \displaystyle 2 \cdot k \le \frac{2 \cdot \text{ln}(w)}{\text{ln}(2)}. So the fast powering algorithm requires at most

    \displaystyle \frac{2 \cdot \text{ln}(w)}{\text{ln}(2)}

many multiplications and at most that many divisions to compute the congruence calculation a^w \equiv (\text{mod} \ m).

For example, when the power w=2^{127}-1, which is a Mersenne prime, which has 39 digits. Now w \approx 2^{127}. By the above calculation, the fast powering algorithm would take at most 254 multiplications and at most 254 divisions to do the power congruence computation.

The fast powering calculations demonstrated in this post can be done by hand (using a hand-held calculator). In real applications, such calculations should of course be done in a computer.

___________________________________________________________________________________________________________________

Another Example

Use the fast power algorithm to show that

    4030^{2657} \equiv 21144 \ (\text{mod} \ 55049)

    21144^{79081} \equiv 4030 \ (\text{mod} \ 55049)

Note that one congruence is encryption and the other one is decryption. We demonstrate the second calculation.

In doing the second calculation, we use a little bit of help from Fermat’s little theorem. The modulus 55049 is a prime number. So 21144^{55048} \equiv 1 \ (\text{mod} \ 55049). Thus we have:

    21144^{79081}=21144^{24033} \cdot 21144^{55048} \equiv 21144^{24033} \ (\text{mod} \ 55049)

Step (1)

The binary expansions of 24033 are:

    24033=2^0+2^5+2^6+2^7+2^8+2^{10}+2^{11}+2^{12}+2^{14}

Step (2)

Compute 21144^{2^j} modulo 55049 for each j. The computation is displayed in the following table. The rows with * are the results that we need for Step (3).

    \displaystyle \begin{bmatrix} \text{ i }&\text{ }&21144^{2^i}&\text{ }&\text{squaring}&\text{ }&\text{modulo } 55049&\text{ }&\text{ }  \\\text{ }&\text{ }&\text{ } \\ 0&\text{ }&21144^{2^0}&\text{ }&\text{ }&\text{ }&\equiv 21144&\text{ }&*  \\ 1&\text{ }&21144^{2^1}&\text{ }&\equiv 21144^2&\text{ }&\equiv 15807&\text{ }&\text{ }  \\ 2&\text{ }&21144^{2^2}&\text{ }&\equiv 15807^2&\text{ }&\equiv 48887&\text{ }&\text{ }  \\ 3&\text{ }&21144^{2^3}&\text{ }&\equiv 48887^2&\text{ }&\equiv 41483&\text{ }&\text{ }  \\ 4&\text{ }&21144^{2^4}&\text{ }&\equiv 41483^2&\text{ }&\equiv 7549&\text{ }&\text{ }  \\ 5&\text{ }&21144^{2^5}&\text{ }&\equiv 7549^2&\text{ }&\equiv 11686&\text{ }&* \\ 6&\text{ }&21144^{2^6}&\text{ }&\equiv 11686^2&\text{ }&\equiv 41076&\text{ }&* \\ 7&\text{ }&21144^{2^7}&\text{ }&\equiv 41076^2&\text{ }&\equiv 40975&\text{ }&* \\ 8&\text{ }&21144^{2^8}&\text{ }&\equiv 40975^2&\text{ }&\equiv 11174&\text{ }&* \\ 9&\text{ }&21144^{2^9}&\text{ }&\equiv 11174^2&\text{ }&\equiv 7144&\text{ }&\text{ } \\ 10&\text{ }&21144^{2^{10}}&\text{ }&\equiv 7144^2&\text{ }&\equiv 6313&\text{ }&* \\ 11&\text{ }&21144^{2^{11}}&\text{ }&\equiv 6313^2&\text{ }&\equiv 53542&\text{ }&* \\ 12&\text{ }&21144^{2^{12}}&\text{ }&\equiv 53542^2&\text{ }&\equiv 14040&\text{ }&* \\ 13&\text{ }&21144^{2^{13}}&\text{ }&\equiv 14040^2&\text{ }&\equiv 46180&\text{ }&\text{ } \\ 14&\text{ }&21144^{2^{14}}&\text{ }&\equiv 46180^2&\text{ }&\equiv 49189&\text{ }&* \end{bmatrix}

Step (3)

Compute 21144^{79081} \ (\text{mod} \ 55049).

    \displaystyle \begin{aligned} 21144^{24033} &=21144^{2^0} \cdot 21144^{2^5} \cdot 21144^{2^6} \cdot 21144^{2^7} \cdot 21144^{2^{8}} \cdot 21144^{2^{10}} \cdot 21144^{2^{11}} \cdot 21144^{2^{12}} \cdot 21144^{2^{14}} \\&\equiv 21144 \cdot 11686 \cdot 41076 \cdot 40975 \cdot 11174 \cdot 6313 \cdot 53542 \cdot 14040 \cdot 49189 \ \ \text{mod} \ 55049 \\&\equiv 25665 \cdot 40975 \cdot 11174 \cdot 6313 \cdot 53542 \cdot 14040 \cdot 49189 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 55049 \\&\equiv 22328 \cdot 11174 \cdot 6313 \cdot 53542 \cdot 14040 \cdot 49189 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 55049 \\&\equiv 11004 \cdot 6313 \cdot 53542 \cdot 14040 \cdot 49189 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 55049 \\&\equiv 51463 \cdot 53542 \cdot 14040 \cdot 49189 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 55049 \\&\equiv 9300 \cdot 14040 \cdot 49189 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 55049 \\&\equiv 50821 \cdot 49189 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 55049 \\&\equiv 4030 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{mod} \ 55049   \end{aligned}

___________________________________________________________________________________________________________________

\copyright \ 2013 \text{ by Dan Ma}