Primitive and Almost Primitive Trinomials^{∗}

Richard P. Brent Computing Laboratory University of Oxford

rpb@comlab.ox.ac.uk

3 May 2002 Abstract

∗UKCRC Algorithms and Complexity Day, Warwick Copyright \begin{matrix}c \\ ◯ \end{matrix}2002, R. P. Brent. warwacdt

We consider the problem of testing trinomials over GF(2) for irreducibility or primitivity. In particular, we consider trinomials whose degree r is a Mersenne exponent. We describe a new algorithm for testing primitivity of such trinomials. The algorithm has been used to find primitive trinomials of very high degree, e.g. r = 3021377. For certain r, primitive trinomials of degree r do not exist. We show how to overcome this difficulty by finding almost primitive trinomials – polynomials with a primitive factor of degree r and overall degree slightly greater than r. In most applications, almost primitive trinomials are almost as useful as primitive trinomials.

Outline

• Introduction

• Definitions

• Random number generators • Sieving

• The standard algorithm

• The new algorithm

• Performance

• Some new primitive trinomials • Almost primitive trinomials • Acknowledgement

• Bibliography

Introduction

Irreducible and primitive polynomials over finite fields have many applications in cryptography, coding theory, random number generation etc. See, for example, the books by Golomb, Knuth (Vol. 2), and Menezes et al.

In this talk I will describe a new algorithm which has been used to find primitive polynomials of very high (in fact, “world record”) degree over the field GF(2) of two elements.

The results can be generalised to almost primitive polynomials, which are (roughly speaking) polynomials with a large primitive factor.

Polynomials over GF(2)

GF(2) is just the set {0, 1} with operations of addition and multiplication modulo 2.

Equivalently, GF(2) is the set of Boolean values {F, T } with operations ⊕ (exclusive or) and & (and).

We consider polynomials over GF(2), that is, polynomials whose coefficients are in GF(2). For the sake of brevity, we won’t repeat this statement every time!

Note that, for polynomials u, v over GF(2),

2u = 2v = 0 .

This implies that u − v = u + v and

(u + v)^{2} = u^{2} + v^{2} .

Some Definitions

We say that a polynomial P (x) is reducible if it has nontrivial factors; otherwise it is irreducible.

If P (x) is irreducible of degree r > 1, then

GF(2^{r}) ≈ Z_{2}[x]/(P (x)) ,

so we have a representation of the finite field GF(2^{r}) with 2^{r} elements. If x is generator for the multiplicative group of Z_{2}[x]/(P (x)), then we say that P (x) is primitive.

Since the multiplicative group has order 2^{r} − 1, we need to know the complete factorisation of 2^{r} − 1 in order to test if an irreducible polynomial is primitive. However, if r is a Mersenne exponent, i.e. 2^{r} − 1 is prime, then irreducibility implies primitivity.

Random Number Generators

Pseudo-random number generators (RNGs) are widely used in simulation. A program running on a fast computer or cluster of PCs might use 10^{9} random numbers per second for many hours. Small correlations or other deficiencies could easily lead to spurious results. In order to have confidence in the results of simulations, we need to have confidence in the statistical properties of the random numbers used. There is no time to talk about random number generators today. Suffice it to say that primitive trinomials are useful for the construction of high quality random number generators. If you are interested, see the references below, or a longer talk which I gave recently in Bristol,

available from http://www.comlab.ox.ac.uk/
oucl/work/richard.brent/talks.html

Some Well-Known Results

The following results can be found in texts such as Lidl, Menezes et al. Here μ is the Möbius function, and ϕ is Euler’s phi function.

1. x^{2^{n}} + x is the product of all irreducible polynomials of degree d|n. For example,

x^{8} + x = x(1 + x)(1 + x + x^{3})(1 + x^{2} + x^{3}) .

2. Let J_{n} be the number of irreducible polynomials of degree n. Then

∑_{d|n} dJ_{d} = 2^{n} and J_{n} = \frac{1}{n} ∑_{d|n} 2^{d}μ(n/d) .

In particular, if n is prime then J_{n} = (2^{n} − 2)/n.

3. The number of primitive polynomials of degree n is P_{n} = ϕ(2^{n} − 1)/n ≤ J_{n}.

In particular, if n is a Mersenne exponent, then P_{n} = J_{n} = (2^{n} − 2)/n.

The Reciprocal Polynomial

If P (x) = ∑_{j=0}^{r} a_{j}x^{j} is a polynomial of degree r, with a_{0} ≠ 0, then

P_{R}(x) = x^{r}P (1/x) = ∑_{j=0}^{r} a_{j}x^{r−j}

is the reciprocal polynomial. Clearly P (x) is irreducible (or primitive) iff P_{R}(x) is irreducible (or primitive).

In particular, if

P (x) = 1 + x^{s} + x^{r} , 0 < s < r

is a trinomial, then the reciprocal trinomial is

P_{R}(x) = 1 + x^{r−s} + x^{r} .

If it is convenient, we can assume that s ≤ r/2 (else consider the reciprocal trinomial).

Similarly, if r is odd, we can assume that s is odd (this will be useful later).

Searching for Irreducible Polynomials

The irreducible polynomials (over GF(2), as usual) of degree r are analogous to primes with r digits. When searching for large primes we can quickly eliminate most candidates by sieving out multiples of small primes. Similarly, when searching for irreducible polynomials, we can eliminate candidates by checking if they are divisible by irreducible polynomials of low degree.

Sieving

Given a candidate P (x), we check if GCD(P (x), x^{2^{n}} + x) ≠ 1, for n = 1, 2, 3, …, N where N < deg(P) is some bound. Since all irreducible polynomials of degree n are divisors of x^{2^{n}} + x, P (x) passes these checks iff it has no irreducible factors of degree ≤ N.

If 2^{n} > r = deg(P), we can save time and space by computing x^{2^{n}} mod P (x) by repeated squaring and reduction, before computing the GCD.

Irreducible Trinomials

For applications such as random number generation, we want irreducible (or better, primitive) polynomials of high degree r and with a small number of nonzero terms. Hence, we restrict attention to trinomials of the form

P (x) = P_{r,s}(x) = 1 + x^{s} + x^{r} , 0 < s < r .

For simplicity, assume r is an odd prime.

Swan’s Theorem

Swan (1962) determines the parity of the number of irreducible factors by an argument involving the discriminant (actually, Swan’s main result is due to Stickelberger (1897)).

If r is an odd prime, then Swan’s theorem implies that P_{r,s}(x) has an even number of irreducible factors (and hence is reducible) if r = ±3 mod 8 and s ≠ 2 or r − 2.

The condition on s can not be omitted, e.g. x^{29} + x^{2} + 1 is irreducible.

Expectation of Success

The probability that a randomly chosen polynomial of degree r is irreducible is of order 1/r. Empirically, it seems that the same holds for trinomials of prime degree r = ±1 mod 8 (this condition implies that Swan’s theorem is not applicable).

Thus, if we consider all s in the range 0 < s < r/2, we expect a small constant number c of irreducible trinomials of degree r. Empirical evidence suggests that c ≈ 3.2 .

For example, considering the 523 prime r ∈ [1000, 10000] such that r = ±1 mod 8, we find exactly 1683 irreducible trinomials, giving an estimate c = 3.22 ± 0.08 .

Open Question. What is this constant c ?

Searching for Irreducible Trinomials

Suppose r is an odd prime, r = ±1 mod 8, and sieving has failed to show that P (x) = P_{r,s}(x) is reducible. The standard algorithm computes

x^{2^{r}} mod P (x)

by r steps of squaring and reduction, then uses the result that P (x) is irreducible iff

x^{2^{r}} = x mod P (x) .

All authors of papers which give tables of irreducible trinomials (Watson, Rodemich, Zierler, Kurita, Heringa, Kumada, …) seem to have used essentially this algorithm, which is why we call it the standard algorithm.

Complexity of the Standard Algorithm

Since we are working over GF(2),

(∑_{j} a_{j}x^{j})^{2} = ∑_{j} a_{j}x^{2j} .

Thus, each squaring step takes O(r) operations.

Each reduction step also takes O(r) operations, since P (x) is a trinomial and we can apply

x^{j+r} = x^{j+s} + x^{j} mod P (x)

for j = r − 2, r − 3, …, 0 to reduce the result of squaring to a polynomial of degree less than r.

Thus, the complete test for reducibility of P_{r,s} takes O(r^{2}) operations, and to test all s takes O(r^{3}) operations (assuming that sieving leaves a constant fraction of trinomials to test).

If the sieve limit is N and 2^{N} ≈ r^{c} for some c < 1, then N ≈ c log_{2} r and we expect O(r/N) trinomials to survive sieving, so the overall complexity might be reduced from O(r^{3}) to O(r^{3}/ log r).

Improving the Standard Algorithm

The standard algorithm uses 2r bits of memory for squaring, and 2r + O(1) ⊕ operations for each reduction (we count bit-operations but in practice we perform 32 or 64 bit-operations in parallel using word-operations; this also applies to our new algorithm). Many of these operations are on bits which are necessarily zero. There is a better algorithm which avoids these redundant operations.

Both algorithms represent a polynomial A(x) = ∑_{j=0}^{r−1} a_{j}x^{j} as a bit-vector a_{0} … a_{r−1}.

Since r is odd and we can consider either P_{r,s} or its reciprocal P_{r,r−s}, we can assume that s is odd.

The New Algorithm – Squaring

The first point is that there is no need to actually perform the squaring step! The standard algorithm would replace the bit vector

a_{0}a_{1}a_{2} … a_{r−2}a_{r−1}

by its “square”

a_{0}0a_{1}0a_{2} … 0a_{r−2}0a_{r−1} .

However, we can simply keep the bit vector

a_{0}a_{1}a_{2} … a_{r−2}a_{r−1}

and regard it as implicitly representing a square (in other words, we do not store the coefficients of odd-degree terms, since they are known to be zero).

The New Algorithm – Reduction

To see how the reduction can be performed after our “implicit squaring”, consider the example r = 7, s = 3. We initialise A(x) ← x, i.e. a_{0} … a_{6} ← 0100000. The “squaring” operation is implicit: we keep the bit-vector 0100000 and regard this as representing

a_{0}a_{2}a_{4}a_{6}a_{8}a_{10}a_{12}

We now reduce mod P (x) = 1 + x^{3} + x^{7}. Observe that x^{12} = x^{5} + x^{8} mod P (x), so we should replace a_{8} by a_{8} ⊕ a_{12} and a_{5} by a_{5} ⊕ a_{12}, but a_{5} is currently zero, so we can simply regard the rightmost bit as representing a_{5} rather than a_{12}. Thus, after the first step of the reduction we have a bit-vector representing

a_{0}a_{2}a_{4}a_{6}\underline{a_{8}}a_{10}a_{5} .

The only bit(s) which could have changed, because they depend on the result of an ⊕ operation, are underlined.

Example of Reduction continued

Proceeding in a similar fashion, we observe that x^{10} = x^{3} + x^{6} mod P (x), but a_{3} = 0 , so we replace a_{6} by a_{6} ⊕ a_{10} and implicitly regard the second bit from the right as representing a_{3} rather than a_{10}. Thus, after the reduction we have a bit-vector representing

a_{0}a_{2}a_{4}\underline{a_{6}}a_{8}a_{3}a_{5} .

One more step of reduction gives a bit-vector representing

a_{0}a_{2}\underline{a_{4}}a_{6}a_{1}a_{3}a_{5} .

Observe that this bit-vector contains the coefficients of A(x)^{2} mod P (x), but they are in a shuffled order. We need to apply an interleave permutation to get back to the natural order

a_{0}a_{1}a_{2}a_{3}a_{4}a_{5}a_{6} .

Interleaving

Interleaving is closely related to squaring. In fact, if we square a_{0}a_{2}a_{4}a_{6}:

a_{0}a_{2}a_{4}a_{6} → a_{0}0a_{2}0a_{4}0a_{6}0 ,

square and rightshift a_{1}a_{3}a_{5}0:

a_{1}a_{3}a_{5}0 → 0a_{1}0a_{3}0a_{5}00 ,

and apply a bitwise ∨ operation, we obtain

a_{0}a_{1}a_{2}a_{3}a_{4}a_{5}a_{6}0 .

Thus, interleaving can be implemented by squaring and a few additional operations (shifting and ∨-ing). Although two squarings are necessary, the bit-vectors are only half as long as before, so the work involved is almost the same.

A Complete Example

Consider the example r = 7, s = 3. The k-th operation of (implicitly) squaring and reducing mod P (x) is denoted by S_{k}, and the k-th operation of interleaving by I_{k}.

If we start with A(x) = x and perform operations S_{1}, I_{1}, S_{2}, I_{2}, …, S_{7}, I_{7} we obtain the following:

S_{1} → 0100000, I_{1} → 0010000 ≡ x^{2}
S_{2} → 0010000, I_{2} → 0000100 ≡ x^{4}
S_{3} → 0010100, I_{3} → 0100100 ≡ x + x^{4}
S_{4} → 0110100, I_{4} → 0110100 ≡ x + x^{2} + x^{4}
S_{5} → 0100100, I_{5} → 0110000 ≡ x + x^{2}
S_{6} → 0110000, I_{6} → 0010100 ≡ x^{2} + x^{4}
S_{7} → 0000100, I_{7} → 0100000 ≡ x

Since the final result is x, we deduce that P (x) = 1 + x^{3} + x^{7} is irreducible.

The New Algorithm

We now describe the new algorithm formally, in terms of bit-operations. As before, assume that r and s are odd.

To avoid confusion, we denote the working bit-array by b_{0}b_{1} ⋯ b_{r−1}. This bit-array is used to represent the coefficients a_{0}a_{1} ⋯ a_{r−1} of the polynomial A(x), but not necessarily in the natural order.

Let α = (r − 1)/2 and δ = (r − s)/2. Since r and s are odd, α and δ are integers. Initially we set b_{1} ← 1 and the other b_{j} ← 0 to represent A(x) = x.

Squaring and Reduction

Each step S_{k} is implemented by

for j ← r − 1 downto α + 1 do b_{j−δ} ← b_{j−δ} ⊕ b_{j}.

Squaring and reduction step S_{k}

Note that there are only r/2 + O(1) “⊕” bit-operations in the loop, which is a 75% reduction over the 2r + O(1) for the reduction step of the standard algorithm.

Interleaving

The obvious implementation of the interleaving step I_{k} requires a temporary bit-array (say c_{0}c_{1} ⋯ c_{r−1}). For example:

c_{0} ← b_{0};

for j ← 1 to α do {forward interleave} begin c_{2j−1} ← b_{j+α}; c_{2j} ← b_{j}; end;

for j ← 0 to r − 1 do b_{j} ← c_{j}.

Forward interleave I_{k} with copy

We call this a “forward interleave” because the first loop index j increases.

We can avoid the final loop (copying the c array to b) by alternately using the array b and the array c (or by interchanging pointers appropriately). However, the space required is still 2r + O(1) bits, the same as for the standard algorithm.

A Refinement: Overlapping Arrays

We can interleave in the backward direction (replace “for j ← 1 to α” by “for j ← α downto 1” above). Suppose we also interchange the roles of b and c to avoid the final copy.

The point of interleaving alternately in the forward and backward directions is that we can save space by using a single working array of size 3r/2 + O(1) bits. The b and c arrays can partially overlap – in fact b_{j} can occupy the same memory as c_{j+α} (j = 0, 1, …), as shown:

b_{0} b_{1} ⋯ b_{α} ⋯ b_{r−1}

c_{0} c_{1} ⋯ c_{α} ⋯ c_{r−1}

Note that the “forward interleave” transmits data from b to c (i.e. to the left) and the “backward interleave” transmits data from c to b (i.e. to the right)!

Effect of the Refinement

Partially overlapping the arrays b and c can improve performance dramatically on machines with memory hierarchies and cache sizes of less than 2r bits, because the working set is reduced in size by 25%. It has little effect on machines with much larger caches.

Generalisation of the New Algorithm

If we replace 1 + x^{s} + x^{r} by

1 + x^{s_{1}} + ⋯ + x^{s_{k}} + x^{r} ,

then an obvious generalisation of our new algorithm is applicable provided that r is odd and the s_{i} all have the same parity (all odd or all even).

Comparison of the Algorithms

The new algorithm has 75% fewer ⊕ operations than the standard algorithm.

Perhaps more significant than the number of operations is the number of memory references, which is reduced by 56%, from 8r/w + O(1) loads/stores to \frac{7r}{2w} + O(1) loads/stores, on a machine with wordlength w bits.

Also significant on some machines is that the working set size is reduced by 25%, so memory references are more likely to be in the cache.

In practice the improvement provided by the new algorithm depends on many factors: the values of r and (to a lesser extent) s, the cache size, the compiler and compiler options used, whether inner loops are written in assembler, etc, but it is generally at least a factor of two.

Performance of the New Algorithm

Table 1 gives normalised times for the standard and new algorithms on various processors, for r = 3021377. The third column is the “normalised time” c = time(nsec)/r^{2}.

processor  algorithm  c
300 Mhz P-II  standard  6.31
”  new (no overlap) 2.60
”  new (overlap) 1.64
500 Mhz P-III  ”  0.77
833 Mhz P-III  ”  1.66
300 Mhz SGI R12000  ”  1.16
667 Mhz DEC Alpha  ”  0.60

Table 1: Normalised time to test reducibility

Note that 3r/2 bits is 553KB. The L2 cache size was 512KB on the P-II and P-III machines except only 256KB for the 833 Mhz P-III. The program was written in C, except that on PCs the inner loops were written in assembler to use the 64-bit MMX registers.

Times for Various Degrees

In Table 2 we show the time for a full reducibility test with our new algorithm and various degrees r on a machine (300 Mhz Pentium P-II) with 512KB L2 cache.

\begin{matrix}r & time T (sec) & c = 10^{9}T/r^{2} \\ 19937 & 0.42 & 1.06 \\ 44497 & 2.10 & 1.06 \\ 110503 & 14.4 & 1.18 \\ 132049 & 21.7 & 1.24 \\ 756839 & 812 & 1.42 \\ 859433 & 1027 & 1.39 \\ 3021377 & 15010 & 1.64 \\ 6972593 & 198000 & 4.10 \end{matrix}

Table 2: Time to test reducibility on a P-II

New Primitive Trinomials

In Table 3 we give a table of primitive trinomials x^{r} + x^{s} + 1 where r is a Mersenne exponent (i.e. 2^{r} − 1 is prime). We assume that 0 < 2s ≤ r (so x^{r} + x^{r−s} + 1 is not listed).

Results for r < 756839 are given by Heringa et al. [9]. We have confirmed these results.

The entries for r < 3021377 have been checked by running at least two different programs on different machines.

During this checking process, the entry with

r = 859433, s = 170340

was found. This was surprising, because Kumada et al. [11] claimed to have searched the whole range for r = 859433. It turns out that Kumada et al. missed this entry because of a bug in the sieving routine!

New Primitive Trinomials cont.

The three entries for r = 756839 are new (Kumada et al. did not search for this r), as are the two entries for r = 3021377.

r  s  Notes
756839 215747 BLZ, 14 June 2000
267428 BLZ, 11 June 2000
279695 BLZ, 9 June 2000
859433 170340 BLZ, 26 June 2000
288477 Kumada et al. [11]
3021377 361604 BLZ, 8 August 2000
1010202 BLZ, 17 Dec 2000
6972593 none yet search 36% complete

Table 3: Primitive trinomials

The new entries are from Brent, Larvala and Zimmermann (BLZ) [5].

“Almost primitive” trinomials

For about half the Mersenne exponents r (those with r = ±3 mod 8, r > 5), primitive trinomials of degree r probably do not exist. Examples are r = 13, 19, 61, . . ..

In applications it would be almost as good to find trinomials of slightly higher degree, say r + δ, having a primitive polynomial of degree r as a factor. Thus the period of the associated linear recurrence would be a small multiple of 2^{r} − 1 (except for certain exceptional initial conditions which are easy to avoid). Blake, Gao and Lambert recently found some such trinomials of degree up to 500.

We can use a slight modification of our searching algorithm to find such trinomials. The sieving phase needs some modifications, and we have to allow the possibility of trinomials of even degree. Some examples are given in Table 4.

Some Almost Primitive Trinomials

For the values of r, δ and s given in Table 4,

x^{r+δ} + x^{s} + 1

has a primitive factor of degree exactly r and period

(2^{r} − 1)f > 2^{r+δ−1} .

The values of r are all the Mersenne exponents for which primitive trinomials of degree r do not exist, 107 < r < 2976221.

\begin{matrix}r & δ & s & f \\ 2203 & 3 & 355 & 7 \\ 4253 & 8 & 1806 & 255 \\ 9941 & 3 & 1077 & 7 \\ 11213 & 6 & 227 & 63 \\ 21701 & 3 & 6999 & 7 \\ 86243 & 2 & 2288 & 3 \\ 216091 & 12 & 42930 & 3937 \\ 1257787 & 3 & 74343 & 7 \\ 1398269 & 5 & 417719 & 21 \end{matrix}

Table 4: Some Almost Primitive Trinomials

Acknowledgement

This is joint work with Samuli Larvala (HUT) and Paul Zimmermann (INRIA Lorraine).

Thanks are also due to Nate Begeman, Nicolas Daminelli, Shuhong Gao, Robert Hedges, Brendan McKay, Barry Mead, Mark Rodenkirch, Juan Varona, and Mike Yoder for their assistance in various ways.

References

[1] S. L. Anderson, Random number generators on vector supercomputers and other advanced architectures, SIAM Rev. 32 (1990), 221–251.

[2] I. F. Blake, S. Gao and R. Lambert, Construction and distribution problems for irreducible trinomials over finite fields, preprint, July 2001.

[3] R. P. Brent, On the periods of generalized Fibonacci recurrences, Math. Comp. 63 (1994),

389–401. http://www.comlab.ox.ac.uk/oucl/
work/richard.brent/pub/pub133.html

[4] R. P. Brent, Random number generation and simulation on vector and parallel computers, Proc. Fourth Euro-Par Conference, LNCS 1470, Springer-Verlag, Berlin, 1998, 1–20.

· · ·/pub185.html

[5] R. P. Brent, S. Larvala and P. Zimmermann, A fast algorithm for testing reducibility of trinomials mod 2 and some new primitive trinomials of degree 3021377, Math. Comp., to appear. Also “A fast algorithm for testing irreducibility of trinomials mod 2 (preliminary report)”, Report PRG TR-13-00, 30 December

2000. · · ·/pub199.html

[6] D. Coppersmith, H. Krawczyk and Y. Mansour, The shrinking generator, Proc. CRYPTO’93, LNCS 773 (1994), 22–39.

[7] GIMPS, The Great Internet Prime Search,

http://www.mersenne.org/

[8] S. W. Golomb, Shift register sequences, Aegean Park Press, revised edition, 1982.

[9] J. R. Heringa, H. W. J. Blöte and A. Compagner. New primitive trinomials of Mersenne-exponent degrees for random-number generation, International J. of Modern Physics C 3 (1992), 561–564.

[10] D. E. Knuth, The art of computer programming, Volume 2: Seminumerical algorithms (third ed.), Addison-Wesley, Menlo Park, CA, 1998.

[11] T. Kumada, H. Leeb, Y. Kurita and M.

Matsumoto, New primitive t-nomials (t = 3, 5) over GF(2) whose degree is a Mersenne exponent, Math. Comp. 69 (2000), 811–814.

[12] Y. Kurita and M. Matsumoto, Primitive t-nomials (t = 3, 5) over GF(2) whose degree is a Mersenne exponent ≤ 44497, Math. Comp. 56 (1991), 817–821.

[13] R. Lidl and H. Niederreiter, Introduction to Finite Fields and their Applications, Cambridge Univ. Press, Cambridge, second edition, 1994.

[14] G. Marsaglia, A current view of random number generators, in Computer Science and Statistics: The Interface, Elsevier Science Publishers B. V.,1985, 3–10.

[15] A. J. Menezes, P. C. van Oorschot and S. A. Vanstone, Handbook of Applied Cryptography, CRC Press, New York, 1997.

http://cacr.math.uwaterloo.ca/hac/

[16] R. G. Swan, Factorization of polynomials over finite fields, Pacific J. Math. 12 (1962), 1099–1106.

[17] S. Tezuka, Efficient and portable combined Tausworthe random number generators, ACM Trans. on Modeling and Computer Simulation 1 (1991), 99–112.

[18] I. Vattulainen, T. Ala-Nissila and K. Kankaala, Physical tests for random numbers in simulations, Phys. Rev. Lett. 73 (1994), 2513–2516.

[19] N. Zierler and J. Brillhart, On primitive trinomials (mod 2), Inform. and Control 13 (1968), 541–554 and 14 (1969), 566–569.

[20] N. Zierler, Primitive trinomials whose degree is a Mersenne exponent, Inform. and Control 15 (1969), 67–69.