Pepin’s Primality Test

A Fermat number is of the form F_n=2^{2^{n}}+1 where n is a nonnegative integer. A Fermat prime is a Fermat number that is a prime number. The first five Fermat numbers F_0 to F_4 are known to be prime. Pierre de Fermat (1601-1665) conjectured that all Fermat numbers are prime. In 1732, Euler discovered a counterexample by showing that F_5 is composite with F_5=641 \times 6700417. The challenge for working with Fermat numbers is that the numbers grow very rapidly. However, Fermat numbers have many interesting properties. See here for a more detailed look. In this post, we discuss Pepin’s test, which is a primality test solely for Fermat numbers, as well as a sequence of integers that are associated with Fermat numbers. Out of this sequence of integers comes the notion of elite prime numbers.

The following are two versions of the Pepin’s test. Theorem 1 is proved in this previous post.

Theorem 1 (Pepin’s Test)
For any n \ge 1, the Fermat number F_n is a prime number if and only if 3^{(F_n-1) / 2} \equiv -1 \ (\text{mod} \ F_n).

Theorem 2 (Pepin’s Test)
For any n \ge 2, the Fermat number F_n is a prime number if and only if 5^{(F_n-1) / 2} \equiv -1 \ (\text{mod} \ F_n).

To determine whether a Fermat number F_n is prime, simply raise 3 (or 5) to the exponent of (F_n-1)/2 modulo F_n. If the result is congruent to -1, then F_n is prime. Otherwise it is composite. The exponentiation can be performed using the fast powering algorithm (also called square-and-multiply). The idea is to use the binary expansion of the exponent to convert the exponentiation into a series of squarings and multiplications. For Pepin’s test, the exponent is a power of 2. Thus the binary expansion is a 1 followed by a number of zeros. As a result, in programming the “square-and-multiply”, we only need to program the squarings.

The above two versions of Pepin’s test use 3 and 5 as bases. An interesting and natural question arises. What are all the numbers that can serve as bases for the Pepin’s test? In fact, there are many other bases. What would be the property that unites all the possible bases? The sequence A129802 in OEIS describes all possible bases for use in Pepin’s test. In this post, we give a formulation of Pepin’s test that takes into account all the possible bases in this sequence. We also discuss the prime numbers in the sequence, which are also captured in the sequence A102742 (these are called elite prime numbers).

___________________________________________________________________

The Periodic Nature of Fermat Numbers

Let a be any positive integer with a \ge 3. Consider all Fermat numbers F_n reduced modulo a, i.e. F_n \ (\text{mod} \ a). We call these values the Fermat remainders modulo a. Since there are only a many such values, i.e. 0, 1, 2, \cdots, a-1 and since Fermat numbers can be derived by the recursive formula

    F_{n+1}=(F_n-1)^2+1 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (0)

the Fermat remainders are periodic. Let’s look at an example.

Let a=41. The following shows the first several Fermat remainders modulo 41. Each value is computed from the previous value via the recursive relation (0), which is then reduced modulo 41.

    F_0 \equiv 3 \ (\text{mod} \ 41)

    F_1 \equiv 5 \ (\text{mod} \ 41)

    F_2 \equiv 17 \ (\text{mod} \ 41)

    F_3 \equiv 257 \equiv 11 \ (\text{mod} \ 41)

    F_4 \equiv 101 \equiv 19 \ (\text{mod} \ 41)

    F_5 \equiv 325 \equiv 38 \ (\text{mod} \ 41)

    F_6 \equiv 1370 \equiv 17 \ (\text{mod} \ 41)

It is clear the Fermat remainders 17, 11, 19, 38 repeat themselves as n increases. For Fermat remainders modulo 41, the periodicity starts at F_2 and the length of the period is 4 (we call this the Fermat period).

Let’s describe idea behind the example more generally. Because any Fermat remainder must be an element of \left\{0,1,\cdots,a-1 \right\}, there exist at least two Fermat numbers sharing the same remainder, say F_h and F_k with h<k. Assume that h and k are the smallest possible choices. In the modulo 41 example, h would be 2 and k would be 6. As a result, F_h is the starting point of the period and L=k-h is the length of the Fermat period.

Based on the discussion in the above paragraph, for any positive integer a \ge 3, the Fermat remainders modulo a are periodic with the repeating remainders starting at some Fermat number F_h. The length of the Fermat period is the smallest positive integer L satisfying the congruence F_{h+j} \equiv F_{h+j+m \cdot L} \ (\text{mod} \ a) for all j=0,1,\cdots,L-1 and for all positive integer m.

The symbol (\frac{p}{q}) is a Legendre symbol if the bottom number is a prime number, in which case the value of the symbol is 1 or -1. If the bottom number is a composite number, then it is a Jacobi symbol with possible values 0, 1 or -1.

We now wish to evaluate the symbol (\frac{a}{F_n}) where the bottom part is a Fermat number and the top part is a possible base for Pepin’s test. Consider the following derivation.

    \displaystyle \biggl(\frac{a}{F_n}\biggr)=\biggl(\frac{F_n}{a}\biggr)=\biggl(\frac{r}{a}\biggr) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1)

The first Jacobi symbol can be flipped without changing sign because F_n \equiv 1 \ (\text{mod} \ 4) for all n \ge 1. This is due to the law of quadratic reciprocity of Jacobi symbol (see bullet # 6 in Theorem 6 here). The number r in the top part of the third Jacobi symbol represents all the possible Fermat remainders modulo a. For example, if a=41, r would be 17, 11, 19 or 38. The above relation in Jacobi symbol is how we tie Fermat remainders with the Pepin’s test.

___________________________________________________________________

The Sequences

The following theorem gives the linkage between Fermat remainders and the Pepin’s test.

Theorem 3
Let a be a positive integer where a \ge 3. The following conditions are equivalent.

  1. The Jacobi symbol (\frac{r}{a})=-1 for all Fermat remainders r modulo a.
  2. There is a positive integer h such that the Jacobi symbol (\frac{a}{F_n})=-1 for all n \ge h.

Proof of Theorem 3
Note that if condition 1 is not true, then there is one Fermat remainder r such that the Jacobi symbol (\frac{r}{a})=1. Since the Fermat remainders are periodic, via relation (1) above there are infinitely many F_n with (\frac{a}{F_n})=1. Thus condition 2 is not true. In other words, 2 \rightarrow 1 holds. On the other hand, if condition 2 is not true, then condition 1 is not true. \square

Any positive integer a belongs to the sequence A129802 if and only of one of the conditions in Theorem 1 holds. Note that 2 does not belong to the sequence. To see this, note that F_{n} \equiv 1  \ (\text{mod} \ 2) for all n. So there is only one Fermat remainder modulo 2, namely 1. Since the Legendre symbol (\frac{1}{2})=1, condition 1 in Theorem 1 is not true. Thus all possible members of the sequence A129802 must be 3 or higher. The numbers 3 and 5 both belong to the sequence. Both numbers have only one Fermat remainder, namely 2. Note that F_{n} \equiv 2  \ (\text{mod} \ 3) for all n \ge 1 and F_{n} \equiv 2  \ (\text{mod} \ 5) for all n \ge 2. Also note that the Legendre symbols (\frac{2}{3})=-1 and (\frac{2}{5})=-1.

What if we limit the scope of Theorem 3 to just the prime numbers? The following is the translation of Theorem 3 where p is a prime number (using p to replace a).

Theorem 4a
Let p be a prime number. The following conditions are equivalent.

  1. The Legendre symbol (\frac{r}{p})=-1 for all Fermat remainders r modulo p.
  2. There is a positive integer h such that the Legendre symbol (\frac{F_n}{p})=-1 for all n \ge h.

Condition 2 of Theorem 4a should start with the Jacobi symbol (\frac{a}{F_n}). Since it can be flipped without changing sign, it become the Legendre symbol (\frac{F_n}{a}). Next, the following theorem is obtained by rewriting Theorem 4a using the meaning of the Legendre symbol.

Theorem 4b
Let p be a prime number. The following conditions are equivalent.

  1. For each Fermat remainders r modulo p, r is a quadratic nonresidue modulo p
  2. For each Fermat remainders r modulo p, the congruence r^{\frac{p-1}{2}} \equiv -1 \ (\text{mod} \ p) holds.
  3. There is a positive integer h such that the Fermat number F_n is a quadratic nonresidue modulo p for all n \ge h.

Note that condition 2 in Theorem 4b is due to Euler’s criterion.

Any prime number p that satisfies any one of the conditions in Theorem 4b is said to be an elite prime number. Thus an elite prime can be used as a base for the Pepin’s test and is a member of the sequence A129802. In fact, another sequence A102742 is to capture all the elite prime numbers.

Note that 3 and 5 are elite primes as is the prime number 41. The Fermat remainders modulo 41 are 17, 11, 19 and 38 (as shown above). All of the 4 Fermat remainders are quadratic nonresidue modulo 41. The Legendre symbols have value -1 as shown below:

    \displaystyle \biggl(\frac{17}{41}\biggr)=\biggl(\frac{11}{41}\biggr)=\biggl(\frac{19}{41}\biggr)=\biggl(\frac{38}{41}\biggr)=-1

Before further discussing elite primes, let’s have a more general statement of Pepin’s test.

___________________________________________________________________

Pepin’s Test – General Version

Theorem 5 (Pepin’s Test)
Let a \ge 3 be an integer such that there exists a positive integer h such that the Jacobi symbol (\frac{a}{F_n})=-1 for all n \ge h. In other words, a is a member of the sequence A129802. Then for any n \ge h, the Fermat number F_n is a prime number if and only if a^{(F_n-1) / 2} \equiv -1 \ (\text{mod} \ F_n).

Proof of Theorem 5
Let a \ge 3 be a number that satisfies the requirement in Theorem 5, i.e. it is a member of the sequence A129802. Let n \ge h. Suppose that F_n=2^{2^n}+1 is prime. Then the symbol (\frac{a}{F_n}) is a Legendre symbol. By the Euler’s criterion, we have a^{(F_n-1) / 2} \equiv -1 \ (\text{mod} \ F_n).

Suppose that the congruence a^{(F_n-1) / 2} \equiv -1 \ (\text{mod} \ F_n) holds. Then by Lucas primality test, F_n is prime. For the details, see the proof of Theorem 7 in this previous post. \square

Theorem 5 is a version of Pepin’s test that takes into account of all the possibilities for a base. Any number in the sequence A129802 will work. This version is more interesting from a mathematical standpoint. In practice, to use Pepin’s test, we must first settle on a base. Pick a number from list of A129802 or find some number a that satisfies one condition in Theorem 3. It is also important to determine the least number h for the condition in Theorem 3. The h is the lower bound for the F_n that can be tested by Pepin’s test. Of course, it is handy to use 3 or 5 as a base. Thus the versions such as Theorem 1 and Theorem 2 above are useful in practice. Mathematically Theorem 5 is more interesting. It brings out the concept of the integer sequence A129802 and the notion of elite prime.

___________________________________________________________________

More on the Sequences

There is one natural question about the A129802 (possible bases for Pepin’s test) and the sequence A102742 (possible prime bases for Pepin’s test, the elite primes). Does each sequence contain infinitely many numbers? The answer is not known for the second sequence (the elite primes). It is not known whether there are infinitely many elite primes. However, we can show that the first sequence infinitely many members.

Theorem 6

  • If a is a member of the sequence A129802, then the number 2 \times a is also a member of the sequence.
  • If a is a member of the sequence A129802, then for any integer k \ge 0, the number 2^k \times a is also a member of the sequence.

Proof of Theorem 6
We prove the first bullet point. The second follows from the first. Since a is in the sequence A129802 (i.e. it satisfies Theorem 3), there is a least h such that (\frac{a}{F_n})=-1 for all n \ge h. The following derivation shows that (\frac{2a}{F_n})=-1 for all n \ge h.

    \displaystyle \biggl(\frac{2a}{F_n}\biggr)=\biggl(\frac{2}{F_n}\biggr) \times \biggl(\frac{a}{F_n}\biggr)=\biggl(\frac{F_n}{2}\biggr) \times \biggl(\frac{F_n}{a}\biggr)=\biggl(\frac{1}{2}\biggr) \times -1=-1

In the above derivation, a property of Jacobi symbol is used along with the fact that 1 is a quadratic residue modulo 2 and the fact that (\frac{a}{F_n})=(\frac{F_n}{a})=-1. \square

Since 3 is a base for Pepin’s test, 3, 6, 12, 24, 48, 96, … are also bases for Pepin’s test.

___________________________________________________________________

Elite Primes

From the perspective of Fermat numbers, elite prime numbers are simply prime numbers that can serve as bases for Pepin’s primality test for Fermat numbers. Of course this is not a good working definition. Theorem 5 tells us that any prime number that can work as a base for Pepin’s test satisfies the conditions in Theorem 4b (or 4a). Thus we can define elite primes by Theorem 4b. Any odd positive integer p is an elite prime if it satisfies anyone of the conditions in Theorem 4b. We introduce the notion of elite prime by way of Pepin’s test of primality of Fermat numbers. Of course, it is possible to define the notion by using only the conditions in Theorem 4b.

As of the writing of this post, there are 29 known elite primes. The smallest are 3, 5, 7, 41, 15361, 23041, 26881, 61441, 87041. The largest is 3,580,135,407,617, about 3.58 trillion. According to the prime number theorem, there are approximately 124 billion prime numbers in this range. So elite primes are very sparse in the number line. No one know whether there are infinitely many elite primes.

How does one test whether a prime number is an elite prime? One way is to use Theorem 4b. Identify the Fermat remainders modulo the prime p. Compute the congruence r^{\frac{p-1}{2}} \ (\text{mod} \ p) for each remainder r. If r^{\frac{p-1}{2}} \equiv -1 \ (\text{mod} \ p) for all r, then p is an elite prime. it turns out that if p is indeed elite, the calculation to verify is quite light. If p is not elite, there usually is a Fermat remainder early on that disproves the eliteness. We close with two examples.

Example 1
Verify that the largest known elite prime 3,580,135,407,617 is indeed elite.

Let p = 3,580,135,407,617. Then p-1=2^{16} \times 54,628,531. This means that the period will start at or before F_{16}. The following shows the Fermat remainders.

    F_0 \equiv 3
    F_1 \equiv 5
    F_2 \equiv 17
    F_3 \equiv 257
    F_4 \equiv 65537
    F_5 \equiv 4294967297
    F_6 \equiv 3302442361075
    F_7 \equiv 3290048612191
    F_8 \equiv 2705228890692
    F_9 \equiv 3257180532505
    F_{10} \equiv 517702690191
    F_{11} \equiv 690766084782
    F_{12} \equiv 1306969066177
    F_{13} \equiv 2982034914426
    F_{14} \equiv 3080321286431
    F_{15} \equiv 123155921421
    F_{16} \equiv 3185964351132
    F_{17} \equiv 2462332130355
    F_{18} \equiv 3164234988981
    F_{19} \equiv 395355620647
    F_{20} \equiv 2009257348278
    F_{21} \equiv 2441138263337
    F_{22} \equiv 499814121188
    F_{23} \equiv 123155921421

The Fermat period starts at F_{15}. The length of the period is 8 (the remainders are in bold). Let r_j denote the remainder of F_j. Now raise each remainder to the power of \frac{p-1}{2} modulo p. The following shows the results.

    (r_{15})^{\frac{p-1}{2}} \equiv 3580135407616 \equiv -1 \ (\text{mod} \ p)

    (r_{16})^{\frac{p-1}{2}} \equiv 3580135407616 \equiv -1 \ (\text{mod} \ p)

    (r_{17})^{\frac{p-1}{2}} \equiv 3580135407616 \equiv -1 \ (\text{mod} \ p)

    (r_{18})^{\frac{p-1}{2}} \equiv 3580135407616 \equiv -1 \ (\text{mod} \ p)

    (r_{19})^{\frac{p-1}{2}} \equiv 3580135407616 \equiv -1 \ (\text{mod} \ p)

    (r_{20})^{\frac{p-1}{2}} \equiv 3580135407616 \equiv -1 \ (\text{mod} \ p)

    (r_{21})^{\frac{p-1}{2}} \equiv 3580135407616 \equiv -1 \ (\text{mod} \ p)

    (r_{22})^{\frac{p-1}{2}} \equiv 3580135407616 \equiv -1 \ (\text{mod} \ p)

The calculation shows that each Fermat remainder is a quadratic nonresidue modulo 3,580,135,407,617. Thus it is an elite prime. \square.

Example 2
The next prime number beyond the elite prime in Example 1 is p = 3,580,135,407,647. Show that it is not elite.

Note that p-1 = 2 x 1,790,067,703,823. Because the exponent in the power of 2 is 1, the period would begin at F_1 or F_0. For some reason, the length of the Fermat period for a non-elite prime can be very long. The following shows the first 13 or 14 remainders.

    F_0 \equiv 3
    F_1 \equiv 5
    F_2 \equiv 17
    F_3 \equiv 257
    F_4 \equiv 65537
    F_5 \equiv 4294967297
    F_6 \equiv 3302287785295
    F_7 \equiv 2157778759259
    F_8 \equiv 2823039261272
    F_9 \equiv 1745925786363
    F_{10} \equiv 1237739799756
    F_{11} \equiv 2960576514326
    F_{12} \equiv 258269474223
    F_{13} \equiv 2361028595996

Fortunately one of the early remainders disproves the eliteness. The following shows that calculation.

    (r_{0})^{\frac{p-1}{2}} \equiv 1 \ (\text{mod} \ p)
    (r_{1})^{\frac{p-1}{2}} \equiv 3580135407646 \equiv -1 \ (\text{mod} \ p)
    (r_{2})^{\frac{p-1}{2}} \equiv 3580135407646 \equiv -1 \ (\text{mod} \ p)
    (r_{3})^{\frac{p-1}{2}} \equiv 1 \ (\text{mod} \ p)
    (r_{4})^{\frac{p-1}{2}} \equiv 3580135407646 \equiv -1 \ (\text{mod} \ p)
    (r_{5})^{\frac{p-1}{2}} \equiv 3580135407646 \equiv -1 \ (\text{mod} \ p)
    (r_{6})^{\frac{p-1}{2}} \equiv 1 \ (\text{mod} \ p)
    (r_{7})^{\frac{p-1}{2}} \equiv 3580135407646 \equiv -1 \ (\text{mod} \ p)
    (r_{8})^{\frac{p-1}{2}} \equiv 3580135407646 \equiv -1 \ (\text{mod} \ p)
    (r_{9})^{\frac{p-1}{2}} \equiv 3580135407646 \equiv -1 \ (\text{mod} \ p)
    (r_{10})^{\frac{p-1}{2}} \equiv 3580135407646 \equiv -1 \ (\text{mod} \ p)
    (r_{11})^{\frac{p-1}{2}} \equiv 3580135407646 \equiv -1 \ (\text{mod} \ p)
    (r_{12})^{\frac{p-1}{2}} \equiv 1 \ (\text{mod} \ p)
    (r_{13})^{\frac{p-1}{2}} \equiv 1 \ (\text{mod} \ p)

Even though the Fermat period may be very long in this case, there are several Fermat remainders early that show r^{\frac{p-1}{2}} \equiv 1 \ (\text{mod} \ p). This shows that the prime in this example cannot be elite. \square

The notion of elite primes is an interesting one. Even though the initial motivation is Pein’s test – a primality test for Fermat numbers, the notion of elite prime is interesting in its own right. There had been a lot of research done, both computationally and mathematically. See the references listed in the two OEIS sequences A129802 and A102742.

___________________________________________________________________
\copyright \ \ 2016 \ \text{Dan Ma} Let

The first prime number after the 8th Fermat number

In this post, we discuss a primality testing exercise involving the eighth Fermat number. A Fermat number is of the form F_n=2^{2^n}+1 where n is any nonnegative integer. We search for the first prime number that is greater than F_8. The basic idea is to search for the first probable prime base 2 among the odd numbers after F_8. Once the first probable prime base 2 is identified, we apply the Miller-Rabin primality test to confirm that it is a prime number. At the outset of this exercise, we did not know how many numbers we had to check before reaching the first prime number.

The first five Fermat numbers F_0, F_1, F_2, F_3 and F_4 are the only Fermat numbers that are known to be prime (it was conjectured by Fermat that all Fermat numbers are prime). It is unknown whether there exists prime Fermat number beyond F_4. What is clear, however, is that all the higher Fermat numbers that were studied turn out to be composite. The 8th Fermat number F_8 has 78 decimal digits with two factors with 16 and 62 digits (it was factored in 1961). The largest Fermat number that has been completely factored (as of the writing of this post) is F_{11} which has 617 decimal digits. Many Fermat numbers larger than F_{11} have been partially factored.

___________________________________________________________________

The basic approach

The following is the number 2^{256}, which has 78 decimal digits.

    2^{256}=
    11579208923731619542357098500868790785326998466564
    0564039457584007913129639936

Define P_j=2^{256}+j where j is an odd positive integer, i.e., j=1,3,5,7,\cdots. The exercise is to find the smallest j such that P_j is a prime number. According to Euclid’s proof that there are infinitely many prime numbers, such a P_j is sure to exist. Just that we do not know at the outset how far we have to go to find it. Of course, P_1 is the 8th Fermat number, which is a composite number with two prime factors with 16 and 62 decimal digits. So the search starts with j=3.

The key is to do the following two quick checks to eliminate composite numbers so that we can reach a probable prime as quickly as possible.

  • For any given P_j, the first step is to look for small prime factors, i.e., to factor P_j using prime numbers less than a bound B. If a small prime factor is found, then we increase j by 2 and start over. Note that we skip any P_j where the sum of digits is divisible by 3. We also skip any P_j that ends with the digit 5.
  • If no small factors are found, then compute the congruence 2^{P_j-1} \ (\text{mod} \ P_j). If the answer is not congruent to 1, then we know P_j is composite and work on the next number. If 2^{P_j-1} \equiv 1 \ (\text{mod} \ P_j), then P_j is said to be a probable prime base 2. Once we know that a particular P_j is a probable prime base 2, it is likely a prime number. To further confirm, we apply the Miller-Rabin primality test on that P_j.

In the first check, we check for prime factors among the first 100 odd prime numbers (i.e. all odd primes up to and including 547).

___________________________________________________________________

Searching the first probable prime

At the outset, we did not know how many numbers we will have to check. Since there can be a long gap between two successive prime numbers, the worse fear is that the number range we are dealing with is situated in such a long gap, in which case we may have to check thousands of numbers (or even tens of thousands). Luckily the search settles on a probable prime rather quickly. The magic number is 297. In other words, for the number

    P_{297}=2^{256}+297

, we find that 2^{P_{297}-1} \equiv 1 \ (\text{mod} \ P_{297}). Thus P_{297} is a probable prime in base 2. The following shows the decimal digits of P_{297}.

    P_{297}=
    11579208923731619542357098500868790785326998466564
    0564039457584007913129640233

To further give a sense of how the magic number P_{297} is reached, the following table lists the 25 calculations leading to the magic number.

    \left[\begin{array}{rrrrrrr}      j & \text{ } & \text{last 5 digits of } P_j & \text{ }  & \text{least factor of } P_j & \text{ } & 2^{P_j-1} \ \text{mod} \ P_j \\           \text{ } & \text{ } & \text{ }  \\                      259 & \text{ } & 40195 & \text{ } & 5 & \text{ } & \text{ }  \\      261 & \text{ } & 40197 & \text{ } & * & \text{ } & \not \equiv 1  \\      263 & \text{ } & 40199 & \text{ } & 3 & \text{ } & \text{ }  \\      265 & \text{ } & 40201 & \text{ } & * & \text{ } & \not \equiv 1  \\      267 & \text{ } & 40203 & \text{ } & * & \text{ } & \not \equiv 1  \\        269 & \text{ } & 40205 & \text{ } & 3 & \text{ } & \text{ }  \\      271 & \text{ } & 40207 & \text{ } & 7 & \text{ } & \text{ }  \\      273 & \text{ } & 40209 & \text{ } & * & \text{ } & \not \equiv 1  \\      275 & \text{ } & 40211 & \text{ } & 3 & \text{ } & \text{ }  \\      277 & \text{ } & 40213 & \text{ } & 11 & \text{ } & \text{ }  \\        279 & \text{ } & 40215 & \text{ } & 5 & \text{ } & \text{ }  \\      281 & \text{ } & 40217 & \text{ } & 3 & \text{ } & \text{ }  \\      283 & \text{ } & 40219 & \text{ } & 13 & \text{ } & \text{ }  \\      285 & \text{ } & 40221 & \text{ } & 7 & \text{ } & \text{ }  \\      287 & \text{ } & 40223 & \text{ } & 3 & \text{ } & \text{ }  \\        289 & \text{ } & 40225 & \text{ } & 5 & \text{ } & \text{ }  \\      291 & \text{ } & 40227 & \text{ } & 23 & \text{ } & \text{ }  \\      293 & \text{ } & 40229 & \text{ } & 3 & \text{ } & \text{ }  \\      295 & \text{ } & 40231 & \text{ } & 71 & \text{ } & \text{ }  \\      297 & \text{ } & 40233 & \text{ } & * & \text{ } & \equiv 1      \end{array}\right]

The first number is in the table P_{259} ends in a 5 and is thus composite. The third number P_{263} is composite since the sum of the digits of its is divisible by 3. The third column of the above table shows the least prime factor below 547 (if one is found). An asterisk in the third column means that none of the prime numbers below 547 is a factor. For such numbers, we compute the modular exponentiation 2^{P_j-1} \ (\text{mod} \ P_j).

In the above table, 4 of the asterisks lead to the result 2^{P_j-1} \not \equiv 1 \ (\text{mod} \ P_j). These numbers P_j are thus composite. For example, for P_{273}, the following is the result:

    2^{P_{273}-1} \ (\text{mod} \ P_{273}) \equiv
    55365573520609500639906523255562025480037454102798
    631593548187358338340281435

The last number P_{297} in the table is a probable prime base 2 since our calculation shows that 2^{P_{297}-1} \equiv 1 \ (\text{mod} \ P_{297}). Being a probable prime to base 2 is actually very strong evidence that the number is a prime number. We want even stronger evidence that P_{297} is a prime. For example, we can carry out the Miller-Rabin test in such a way that the probability of mistaking a composite number as prime is at most one in a septillion! A septillion is the square of a trillion. A trillion is 10^{12}. Thus a septillion is 10^{24}. One in a septillion is for all practical purposes zero. But if one wants more reassurance, one can always run the Miller-Rabin test with more bases.

___________________________________________________________________

The Miller-Rabin primality test

The Miller-Rabin test is a variant of the Fermat test because Miller-Rabin still relies on Fermat’s little theorem. But Miller-Rabin uses Fermat’s little theorem in such a way that it eliminates the issue of the Fermat test mistakenly identifying Carmichael numbers as prime.

Given an odd positive integer whose “prime or composite” status is not known, the Miller-Rabin test will output “composite” or “probable prime”. Like the Fermat test, the Miller-Rabin test calculates a^{n-1} \ (\text{mod} \ n) for several values of a. But the test organizes the congruence a^{n-1} \ (\text{mod} \ n) a little differently to capture additional information about prime numbers.

Here’s how to set up the calculation for Miller-Rabin. Because n is odd, n-1 is even. We can factor n-1 as a product of a power of 2 and an odd number. So we have n-1=2^k \cdot q where k \ge 1 and q is odd (q may not be prime). Then we calculate the following sequence:

    a^q, \ a^{2 \cdot q}, \ a^{2^2 \cdot q}, \cdots, a^{2^{k-1} \cdot q}, \ a^{2^{k} \cdot q} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1)

The first term in (1) can be calculated using the fast powering algorithm (using the binary expansion of q to convert the calculation of a^q into a series of squarings and multiplications). Each subsequent term is then the square of the preceding term. The last term is of course a^{n-1}. Each squaring or multiplication is reduced modulo n. The Miller-Rabin test is based on the following property of prime numbers:

    Theorem 1
    Let n be an odd prime number such that n-1=2^k \cdot q where k \ge 1 and q is odd. Let a be a positive integer not divisible by n. Then the following two conditions are true about the sequence (1).

    • At least one term in the sequence (1) is congruent to 1 modulo n.
    • Either the first term in (1) is congruent to 1 modulo n or the term preceding the first 1 is congruent to -1 modulo n.

How the Miller-Rabin test works
Suppose that the “prime or composite” status of an odd integer n is not known. If both conditions in the above theorem are satisfied with respect to the number a, then n is said to be a strong probable prime in base a. If a strong probable prime in base a happens to be composite, then it is said to be a strong pseudoprime in base a. In other words, a strong pseudoprime is a composite number that possesses a prime-like property, namely it satisfies the two conditions in Theorem 1 with respect to one base a.

The test procedure of Miller-Rabin is to check whether n is a strong probable prime to several bases that are randomly chosen. The following determines the outcome of the test:

  • If n is not a strong probable prime in one of the chosen bases, then n is proved to be composite.
  • If n is shown to be a strong probable prime in all the chosen bases (say there are k of them), then n is “probably prime” with an error probability of at most 0.25^k.

To prove the integer n is composite, we look for a base a for which n is not a strong probable prime. Such a value of a is also called a Miller-Rabin witness for the compositeness of n. For primality, the Miller-Rabin test does not give a mathematical proof that a number is prime. The Miller-Rabin test is a probable prime test. It gives strong evidence that n is a prime number, with an error probability that can be made arbitrarily small by using a large random sample of values of a.

Take the prime candidate P_{297} that is discussed above. We plan to run the Miller-Rabin test on P_{297} using 40 random values of a where 1<a<P_{297}. If P_{297} is shown to be a strong probable prime in all 40 bases, then the prime candidate P_{297} is likely a prime number with an error probability of at most 0.25^{40}. This probability works out to be less than 1 in 10 raised to 24 (hence the one in a septillion that is mentioned earlier). If one wants stronger evidence, we can compute for more values of a. Thus if P_{297} is in actuality a composite number, there is at most a one in septillion chance that the Miller-Rabin test will declare P_{297} is a prime number.

How can the Miller-Rabin test make the claim of having such a small error probability? The fact the the error probability of Miller-Rabin can be made arbitrarily small stems from the following fact.

    Theorem 2
    Suppose that n is a composite odd number. At most 25% of the numbers in the interval 1<a<n are bases in which n is a strong pseudoprime. Putting it in another way, at least 75% of the numbers in 1<a<n are bases in which n is not a strong pseudoprime.

To paraphrase Theorem 2, if n is composite to begin with, at least 75% of the numbers in 1<a<n will prove its compositeness. That means that at most 25% of the numbers a will exhibit the prime-like property described in Theorem 1. The power of Miller-Rabin comes from the fact that for composite numbers there are more values of a that will give a correct result (in fact, at least 3 times more).

Thus if you apply the Miller-Rabin test on a composite number n, you will bound to stumble on a base a that will prove its compositeness, especially if the bases are randomly chosen. Any random choice of a where 1<a<n has at least a 75% chance of being correct on the composite number n. In a series of 100 random choices of a, it will be hard to miss such values of a. The only way that Miller-Rabin can make a mistake by declaring a composite number as prime is to pick all the values of a from the (at most) 25% of the pool of values of a that are strong pseudoprime prime. This probability is bounded by 0.25^k (if k is the number of selections of a).

___________________________________________________________________

Applying Miller-Rabin on the prime candidate

The first task is to factor P_{297}-1. We find that P_{297}-1=2^3 \times q where q is the following odd number:

    q=
    14474011154664524427946373126085988481658748083205
    070504932198000989141205029

For each randomly selected a, we calculate the following sequence:

    a^q, \ a^{2q}, \ a^{4q}, \ a^{8q} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (2)

The first term is calculated using the fast powering algorithm (a series of squarings and multiplications). Each subsequent term is the square of the preceding term. Each term in the sequence is reduced modulo P_{297}. The goal is to see if the two conditions in Theorem 1 are satisfied. One is that one of the 4 values in (2) is a 1. The other is that the term preceding the first 1 in (2) has to be a -1.

The following shows the 40 numbers that are randomly chosen in the interval 1<a<P_{297}.

    a_1=
    03006957708701194503849170682264647623506815369915
    7798209693214442533348380872
    a_2=
    02223067440101780765895379553626469438082041828085
    0568523022714143509352911267
    a_3=
    04531895131849504635258523281146698909008537921009
    6337435091877410129499153591
    a_4=
    05434508269932993379745836263818598804800824522102
    0278113825689716192402178622
    a_5=
    08799241442673378780142202326330306495270149563840
    3866810486309815815031353521
    a_6=
    02638607393577034288802880492058261281940769238659
    8928068666401909247319838064
    a_7=
    04283430251977183138176255955338099404217762991191
    9192783003754562986178981473
    a_8=
    09773398144692973692102006868849010147546139698798
    3443958657834362269077067224
    a_9=
    05504666974469005713839308880951115507992521746498
    7157086751623602877205126361
    a_{10}=
    11369425784373951812019794994427515082375862595853
    6524984616385315102874812557
    a_{11}=
    11280428157869817083329641054154150272024966029283
    2165114734540900026838117128
    a_{12}=
    11208322317253928483879618989535357346499197200982
    7728283667193655956607063861
    a_{13}=
    05585951853297694372636067012444311272073854408338
    4421611399136081624631900538
    a_{14}=
    06831924581003106427566658433259804779354874917795
    9811865334330929987281859876
    a_{15}=
    07339174229323952008915772840377019251465052264221
    1294344032116313026124007734
    a_{16}=
    05117387267263929174559713098463717229625661656017
    7194611080485470890280573816
    a_{17}=
    06599941646668915168578091934085890873056463577356
    8090503454939353325803291530
    a_{18}=
    07545265152740184887140788322673806569482388835389
    5577110370797470603035554930
    a_{19}=
    02591621894664804222839429868664505564743756550515
    2520842332602724614579447809
    a_{20}=
    04791002227899384351266879075743764807247161403811
    8767378458621521760044966007
    a_{21}=
    03251071871924939761772100645669847224066002842238
    6690935371046248267119967874
    a_{22}=
    07211128555514235391448579740428274673170438137060
    9390617781010839144521896079
    a_{23}=
    02839820419745979344283855308465698534375525126267
    1701870835230228506944995955
    a_{24}=
    06304631891686637702274634195264042846471748931602
    4893381338158934204519928855
    a_{25}=
    06492095235781034422561843267711627481401158404402
    2978856782776323231230432687
    a_{26}=
    11078868891712009912929762366314190797941038596568
    5459274315695355251764942151
    a_{27}=
    05795069944009506186885816367149671702413127414386
    2708093175566185349033983346
    a_{28}=
    01712922833914010148104423892201355622294341143990
    7524285008693345292476544524
    a_{29}=
    09743541325262594740093734822046739122734773994479
    9814337973200740861495044676
    a_{30}=
    02503872375817370838455279068302037475992008315394
    2976462871038003917493744995
    a_{31}=
    06980677383898331402575574511880992071872803011356
    6498794763450065008785347168
    a_{32}=
    01507075889390134242331585173319278262699562685820
    7121480322563439665642035394
    a_{33}=
    02471785068822350832987019936892052187736451275830
    5372059292781558599916131031
    a_{34}=
    10950891460180297156465120507537244257810396062906
    9207306297501015755045004254
    a_{35}=
    11052976297188507170707306917942099264941855478856
    2965936913589165233381674539
    a_{36}=
    03911878231948499128291863266472008604449261315172
    1053813631612297577166335941
    a_{37}=
    06903294587603383022211116535092146484651980588002
    9291840261276683214113088012
    a_{38}=
    03942020579038616658412018517396703874933208670283
    3087287933190554281896471934
    a_{39}=
    04338728160253711124705740270085271024911573570055
    1690460857511205663297661796
    a_{40}=
    06707597137792150532106913489524457238449067437061
    7211249957355483821516113140

For each random number a_j, we calculated the 4 numbers indicated in sequence (2). The following 3 tables show the results of the calculation.

    \left[\begin{array}{rrrrrrrrr}      j & \text{ } & a_j^q & \text{ }  & a_j^{2q} & \text{ } & a_j^{4q} & \text{ } & a_j^{8q} \\           \text{ } & \text{ } & \text{ }  \\                  1 & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1 & \text{ } & 1  \\      2 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      3 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1  \\      4 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      5 & \text{ } & -1 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1  \\        6 & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1 & \text{ } & 1  \\      7 & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1 & \text{ } & 1  \\      8 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      9 & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1 & \text{ } & 1  \\      10 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\        11 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      12 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      13 & \text{ } & -1 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1  \\      14 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      15 & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1 & \text{ } & 1        \end{array}\right]

    \left[\begin{array}{rrrrrrrrr}      j & \text{ } & a_j^q & \text{ }  & a_j^{2q} & \text{ } & a_j^{4q} & \text{ } & a_j^{8q} \\           \text{ } & \text{ } & \text{ }  \\                  16 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1  \\      17 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      18 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      19 & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1 & \text{ } & 1  \\      20 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1  \\        21 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      22 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      23 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      24 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1  \\      25 & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1 & \text{ } & 1  \\        26 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      27 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      28 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1  \\      29 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1  \\      30 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1       \end{array}\right]

    \left[\begin{array}{rrrrrrrrr}      j & \text{ } & a_j^q & \text{ }  & a_j^{2q} & \text{ } & a_j^{4q} & \text{ } & a_j^{8q} \\           \text{ } & \text{ } & \text{ }  \\                  31 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      32 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      33 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      34 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      35 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1  \\        36 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      37 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      38 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1  \\      39 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1 & \text{ } & 1  \\      40 & \text{ } & \text{ } & \text{ } & \text{ } & \text{ } & -1 & \text{ } & 1             \end{array}\right]

There are 4 columns of calculation results, one for each term in sequence (2). If a calculation result is a blank in the above tables, it means that the result is a number that is not 1 or -1 modulo P_{297}. For example, a_1^q \ (\text{mod} \ P_{297}) and a_2^q \ (\text{mod} \ P_{297}) are congruent to the following two numbers:

    a_1^q \equiv
    86168678768024029811437552745042076645410792873480
    629834883948094184848812907

    a_2^q \equiv
    10235477176842589582260882228891913141693105976929
    7597880545619812030150151760

In the above 3 tables, all results match the conditions of Theorem 1. For each number a_j, the calculated results are eventually 1. On some of the rows, the first result is a 1. In all the other rows, the term right before the first 1 is a -1. For example, in the first row where j=1, the first 1 a_1^{4q} and the term preceding that is a -1.

The results in the above 3 tables show that the number P_{297} is a strong probable prime in all 40 of the randomly chosen bases. We have very strong evidence that the number P_{297} is a prime number. The probability that it is a composite number but we mistakenly identify it as prime is at most one in a septillion!

___________________________________________________________________

Exercise

In our search for probable primes larger than the 8th Fermat number, we also find that the number P_{301}=2^{301}+301 is also a probable prime base 2. The following shows the decimal digits:

    P_{301}=
    11579208923731619542357098500868790785326998466564
    0564039457584007913129640237

Is it a prime number? Perform the Miller-Rabin test on this number.

___________________________________________________________________
\copyright \ \ 2014 \ \text{Dan Ma}

A primality testing exercise from the 9th Fermat number

In this post, we discuss a primality testing exercise on a number derived from a prime factor of the 9th Fermat number. The idea of the exercise is that we test the primality of a number generated from a known large prime number by tagging or inserting a digit into the known prime number. This is an excellent way to generate exercises for primality testing. This approach of altering a known prime number is not meant to be a prime number generation process.

A Fermat number is a positive integer of the form F_n=2^{2^n}+1 where n=0,1,2,\cdots. The 9th Fermat number is F_9=2^{2^9}+1=2^{512}+1. This is a 155-digit number that was completely factored in 1993 [1]. Its three prime factors have 7 digits, 49 digits and 99 digits. The following is the 49-digit factor of F_9:

    7,455,602,825,647,884,208,337,395,736,200,454,918,783,366,342,657

The above 49-digit number is a prime number that was resulted from a published effort of factorizing F_9. The other two factors of F_9 can be found in [1] or here. If we tag a digit in front of this 49-digit number, say the digit 6, is the resulting number a prime number? In this post, we discuss the primality testing of the following 50-digit number:

    A= 67,455,602,825,647,884,208,337,395,736,200,454,918,783,366,342,657

Primality testing is an important aspect of public key cryptography. For example, the modulus in the public key in RSA cryptosystem is a product of two distinct large prime numbers. To generate the key, it is crucial to generate large numbers at random and to be able to efficiently test whether the numbers are prime.

___________________________________________________________________

Pre-checks before primality testing

Primality tests such as Fermat test and Miller-Rabin test are popular and are easy to use. These tests check for the primality of a given integer n by randomly choosing a set of auxiliary numbers and then checking whether any one of these auxiliary values is a witness to the compositeness of the number n. We would like to emphasize the point that before applying a randomized primality testing procedure, it will be advantageous to perform some easy checks to quickly rule out primality. If these quick checks show that the number being tested is composite, we can then continue the search by moving on to another candidate prime number. With the above number A as example, we perform the following pre-checks.

  • Check small prime numbers just in case the number being tested has small prime factors.
  • If the previous step produces no factors, look for small Fermat witnesses.

The idea of the first check is clear. A Fermat witness for an integer n is an integer a that is relatively prime to n such that a^{n-1} \not \equiv 1 \ (\text{mod} \ n). If we can find such a number a, we can say conclusively that n is not a prime number. According to Fermat’s little theorem, if n is a prime number, then the congruence would be a one, i.e., a^{n-1} \equiv 1 \ (\text{mod} \ n). This idea is the basis for the Fermat primality test, which is performed on a random set of values of a. Here in the pre-check, we calculate a^{n-1} \ (\text{mod} \ n) for some small values of a such as 2 or 3. The idea is to quickly rule out the candidate prime numbers that happen to have a small Fermat witness.

If either one of the two checks indicate the number being tested is composite, then we pick another number to test. If both of the above steps fail to indicate compositeness, we can then proceed to use the “probabilistic” part of the Fermat test or another probabilistic primality test.

We like to point out that at the outset, the author of the blog did not know the primality status of the 50-digit number A that is obtained by tagging a 6 in front of the 49-digit factor of the 9th Fermat number. Though this number is far too small to work in the RSA algorithm, it serves as an excellent exercise for primality testing.

___________________________________________________________________

Pre-checking the number A

Very quickly, 7 is found to be a factor of A. We have A=7 \times B where B is:

    B= 9636514689378269172619627962314350702683338048951

We then try to find small factors of B. The number 283 is found to be a factor of B. Then we have A=7 \times 283 \times C where C is the following 47-digit number:

    C= 34051288655046887535758402693690285168492360597

By tagging a 6 in front of the 49-digit factor of the 9th Fermat number, we obtain a composite number. The prime or composite status for the number A is settled. Now the next question: is the remaining factor C prime or composite?

___________________________________________________________________

Pre-checking the number C

We perform the two pre-checks as described above. In searching for small factors, we find that C is not divisible by all prime numbers less than 200,000. So if C has a factor, it will have to be a large one. Or it could be that C is prime. Instead of continuing to search for factors, we can do the second check. We calculate 2^{C-1} \ (\text{mod} \ C) and find that:

    2^{C-1} \equiv t \ (\text{mod} \ C)

    where t= 7426390354305013563302537374271875139618265902

Since t clearly is not 1, this tells us that C is composite.

___________________________________________________________________

Remarks

Instead of doing the randomized primality testing right away, it makes sense to do the pre-checks as described above so that we do not waste time and effort on composite numbers that have small factors or have small Fermat witnesses. As soon as these two checks indicate the compositeness of a number, we then move on to another candidate prime number. If a number passes these two checks, then we can apply a primality test such as the Fermat test or the Miller-Rabin test.

Another remark is that primality testing is easy while factoring is hard. The compositeness of the 47-digit number C can be determined by just one modular exponentiation calculation, which can be computed using the fast powering algorithm. This algorithm uses the binary expansion of the exponent to convert the modular exponentiation into a series of squarings and multiplications. Because the exponent C-1 is a 155-bit number, computing 2^{C-1} \ (\text{mod} \ C) only requires 154 squarings and at most 154 multiplications (the actual count of multiplications is 76). Implemented in a modern computer, this calculation takes less than one second.

On the other hand, factoring the above 47-digit number C by the brute force approach (checking every possible prime factor) would require checking approximately 2 \times 10^{20} many prime numbers as potential factors. The age of the universe is believed to be about 13 billions years, which is about 4 \times 10^{17} seconds. That’s why we stop checking for prime factors of C at some point. The authors of [1] factored a 155-digit number in 1993. So the current computing technology is certainly capable of factoring a 47-digit number. But the effort will take time and will require using more advanced factoring algorithm.

___________________________________________________________________

Exercises

The idea of the exercise is that we test the primality of a number generated from a known large prime number by tagging or inserting a digit into the known prime number. Using the same 49-digit factor of F_9, we can tag another digit in front or insert a digit somewhere in the middle. Another possibility is to try the other two factors of F_9.

There are many other large primes to be found on the Internet that can be used for this exercise. For example, factors of other Fermat numbers and RSA numbers.

___________________________________________________________________

Reference

  1. Lenstra A. K., Lenstra H. W., Manasse M. S., Pollard J. M. The Factorization of the Ninth Fermat Number, Mathematiics of Computation, Volume 61, Number 203, July 1993, Pages 319-349.

___________________________________________________________________
\copyright \ \ 2014 \ \text{Dan Ma}