- Programming Notes
- References
- Glossary

This document is a gentle introduction to computational number theory. The plan of the paper is to first give a quick overview of arithmetic in the modular integers. Throughout, we will emphasize computation and practical results rather than delving into the why. Simple programs, generally in JavaScript, are available for all of the algorithms mentioned. At the end of the paper we will introduce the Gaussian Integers and Galois Fields and compare them to the modular integers. Companion papers will examine number theory from a more advanced perspective.

Modular arithmetic is arithmetic using integers modulo some fixed
integer *N*. Thus, some examples of operations modulo 12 are:

- 7 + 7 = 14 = 2 (mod 12)
- 5 * 7 = 35 = 11 (mod 12)

Further examples can be generated and checked out with the following short programs. Note that, as JavaScript cannot compute with integers larger than 53 bits, the largest modulus allowed is 26.5 bits, or about 94,900,000.

Among the basic operations we have missed the division operator. If we were working in the integers we would almost never be able to define a quotient (unless the answer is itself an integer). In the modular integers we can often, but not always, define a quotient:

- 11 / 5 = 7 (mod 12) as 5 * 7 = 11 (mod 12)
- 5
^{(-1)}= 1 / 5 = 17 (mod 21) as 5 * 17 = 85 = 1 (mod 21) - 48 / 31 = 72 (mod 91) as 31 * 72 = 48 (mod 91)
- 9 / 3 = 7 (mod 12) as 3 * 7 = 21 = 9 (mod 12)

As the last example points out, modular division does not always produce a unique result, for other correct answers are 3 and 11 (as 3 * 3 = 9 (mod 12) and 3 * 11 = 33 = 9 (mod 12)). In particular, if the modulus and the divisor share a common factor (in this case 3 divides both 3 and 12), the answer will not be unique if it exists at all. We draw two conclusions from this - the first is that we would like to be able to spot such common factors, and the second is that we would like to avoid such situations. Usually we avoid the situation ever showing up by using only prime moduli. Spotting the situation is one of the results of the next section.

The Euclidean Algorithm solves two problems we have posed:

- Common Factors - Given two numbers
*n*and*m*, find any common factors they may have. - Modular Inverses - Given a number
*n*and a modulus*m*, find the inverse of*n*, ie the number*k*such that*n***k*= 1 (mod*m*).

The largest number which divides two numbers
`n` and `m` is called the *greatest common divisor*
of `n` and `m`, and denoted
`n`, `m`)

The algorithm to find gcd(*n*, *m*) runs something
like this:

- Write down both
*n*and*m* - Reduce the larger modulo the smaller
- Repeat step 2 until the result is zero
- Publish the preceding result as gcd(
*n*,*m*)

An example of this algorithm is the following computation of gcd(120,222):

120 222 120 222-120=102 120-102=18 102 18 102-5*18=12 18-12=6 12 6 12-2*6=0

Thus gcd(120,222)=6.

The following short program will allow you to compute examples of the Euclidean Algorithm:

If you think carefully about the Euclidean algorithm you will see that
at each step both numbers are formed by adding some multiples of the original
numbers. Thus there are some numbers, `a` and `b`, such
that gcd(`n`,`m`)=`a`*`n`+`b`*`m`.
The Extended Euclidean algorithm can recover not only gcd(`n`,`m`),
but also these numbers `a` and `b`.

The Extended Euclidean algorithm not only computes gcd(`n`,`m`),
but also returns the numbers `a` and `b` such that
gcd(`n`,`m`)=`a`*`n`+`b`*`m`.
If gcd(`n`,`m`)=1 this solves the problem of computing modular
inverses.

Assume that we want to compute `n`^{(-1)}(mod `m`)
and further assume that gcd(`n`,`m`)=1. Run the Extended
Euclidean algorithm to get `a` and `b` such that
`a`*`n`+`b`*`m`=1. Rearranging this result,
we see that `a`*`n`=1-`b`*`m`, or
`a`*`n`=1(mod `m`). This solves the problem of
finding the modular inverse of `n`, as this shows that
`n`^{(-1)}=`a`(mod `m`).

The Extended Euclidean algorithm is nothing more than the usual Euclidean
algorithm, with side computations to keep careful track of what combination
of the original numbers `n` and `m` have been added at
each step. The algorithm runs generally like this:

- Write down
`n`,`m`, and the two-vectors (1,0) and (0,1) - Divide the larger of the two numbers by the smaller - call this
quotient
`q` - Subtract
`q`times the smaller from the larger (ie reduce the larger modulo the smaller) - Subtract
`q`times the vector corresponding to the smaller from the vector corresponding to the larger - Repeat steps 2 through 4 until the result is zero
- Publish the preceding result as gcd(
`n`,`m`)

An example of this algorithm is the following computation of
30^{(-1)}(mod 53):

53 30 (1,0) (0,1) 53-1*30=23 30 (1,0)-1*(0,1)=(1,-1) (0,1) 23 30-1*23=7 (1,-1) (0,1)-1*(1,-1)=(-1,2) 23-3*7=2 7 (1,-1)-3*(-1,2)=(4,-7) (-1,2) 2 7-3*2=1 (4,-7) (-1,2)-3*(4,7)=(-13,23) 2-2*1=0 1 (4,-7)-2*(-13,23)=(30,-53) (-13,23)

From this we see that gcd(30,53)=1 and, rearranging terms, we see
that ^{(-1)}=23(mod 53).
(This can be confirmed by checking that 23*30=1(mod 53).)

The following short program will allow you to compute examples of the extended Euclidean Algorithm:

When computing a power of a number with a finite modulus there are efficient ways to do it and inefficient ways to do it. In this section we will outline a commonly used efficient method which has the curious name "the Russian Peasant algorithm".

The most obvious way to compute 12^{10} (mod 23) is to
multiply 12 a total of nine times, reducing the result mod 23 at
each step. A more efficient method which takes only four multiplications
is accomplished by first noting that:

```
12
```^{2}=6 (mod 23)

12^{4}=6^{2}=13 (mod 23)

12^{8}=13^{2}=8 (mod 23)

We have now performed three squarings and, by noting that the exponent breaks into powers of 2 as 10=8+2, we can rewrite our computation:

```
12
```^{10}=12^{(8+2)}

=12^{8}*12^{2}

=8*6=2(mod 23)

So our algorithm consists of writing the exponent
as powers of two. (This can be done by writing it as
a number base 2 and reading off successive digits -
eg 10_{10}=1010_{2}.) Now we multiply
successive squares of the base number for each digit
of the exponent which is a "1".

The following short program will allow you to compute examples of the Russian Peasant method for exponentiation:

The (slightly cryptic) output of this program can be read as:

- The base two digits of the exponent
- The successive squarings of the base (ie
`b`^{2},`b`^{4},`b`^{8}, ...) - The product of the appropriate squares

A *prime* is a number which has no divisors
other than 1 and the number itself.

Thus, 13 is a prime but 12 is not - it has divisors 2, 3, and 6. A number which is not prime is said to be composite.

The idea of primes raises several questions - how does one find a prime and how do you reduce a number into a product of its prime factors (ie how do you factor it)?

There are two ways to find primes. The oldest is to search a large number of numbers, sieving out the non-prime (ie composite) numbers, leaving only the primes. The second method is to test an individual number for primality.

The *Sieve of Eratosthenes* is a recent method
(dating from the 3rd century BC). The idea is this - write
down all the integers up to some point:

2 3 4 5 6 7 8 9 10 11 12 13 14 15.

Now cross out all the even ones larger than 2:

2 3 ~~4~~ 5 ~~6~~ 7 ~~8~~
9 ~~10~~ 11 ~~12~~ 13
~~14~~ 15.

Now cross out all those divisible by 3 (other than 3):

2 3 ~~4~~ 5 ~~6~~ 7 ~~8~~
~~9~~ ~~10~~ 11 ~~12~~ 13
~~14~~ ~~15~~.

At each step we find the next prime (number not crossed out)
publish it as a prime and cross out all of its multiples.
So, in this example we see that the numbers 2, 3, 5, 7, 11,
and 13 are prime.

While the sieve is an efficient way to find a large number of successive primes it cannot be used to test an arbitrary number for primality. This can be done with a pseudoprime test.

A *pseudoprime test* would be more accurately called
a compositeness test as when it is applied to an integer
`N` it returns one of two possible results:

- Succeeds:
`N`is definitely composite. - Fails:
`N`*might*be prime (or might be composite)

The following is the simplest pseudoprime test. In the literature it is normally not given a name but as it is based on Fermat's Little Theorem (stated in the next section) I will call it the Fermat test. This test runs like this:

- Choose a base
`b`for the test. - If
`b`^{(N-1)}=1(mod`N`) then`N`might be prime, if it is any other value then`N`is definitely composite.

Lets try a few examples:

- Test 35=5*7 for primality:
- Choose base 2
- Compute 2
^{34}(mod 35) = 9 - As the result is not 1 we conclude that 35 is composite.

- Test 31 (known to be prime) for primality:
- Choose base 2
- Compute 2
^{30}(mod 31) = 1

So 31 is a pseudoprime for base 2. - Also compute:

3^{30}(mod 31) = 1

5^{30}(mod 31) = 1 - From this see that 31 is a pseudoprime for bases
2, 3, and 5, and we conclude that 31 is
*probably*prime.

- Test 341=11*31 for primality:
- Choose base 2
- Compute 2
^{340}(mod 341) = 1

So 341 is a pseudoprime for base 2. - Choose base 3
- Compute 3
^{340}(mod 341) = 56

So 341 is definitely composite.

The last example should be somewhat disturbing as it shows that a composite number can masquerade as a prime under the Fermat pseudoprime test. For this there is both good news and bad news.

- The first piece of good news is that the vast majority of pseudoprimes are actually prime.
- The next piece of good news is that most composite pseudoprimes will be revealed as composite after very few bases are tried.
- The bad news is that there are some composite numbers which are pseudoprimes (masquerading as primes) for all bases. These rare numbers are called Carmichael numbers.

More examples of this test can be computed by using the previously mentioned method of fast exponentiation:

Why did the primality test work? The Fermat pseudoprime test,
as well as a number of
others, depends on a simple property of primes first noticed
by Pierre Fermat in the 17^{th} century:

**Theorem:** (Fermat's Little Theorem) If `p` is prime
and 0<`b`<`p` then
`b`^{(p-1)}
=1(mod `p`)

In order to winnow primes from non-primes we make use of the
fact that for a prime `p` this is true for any base
`b`, while for a non-prime `q` and a chosen
base `b` the value of
`b`^{(q-1)}(mod `q`)

The smallest non-zero exponent `e` such that
`a`^{e} = 1 (mod `N`) is the
*order* of `a` mod `N`. This is generally
denoted o(`a`) (mod `N`).

A consequence of Fermat's Little Theorem is that, for prime `p`
and 0<`a`<`p`, the order of `a` mod
`p` divides (`p`-1).

The simplest way to factor an integer is to trial divide it by all of the primes less than its square root. For large numbers this is a very inefficient method, though. A number of better methods have been developed during the mid and late 20th Century, though. We will describe Pollard's p-1 algorithm.

To use Pollard's p-1 factoring method to factor a composite number
`N`:

- Choose a base
`b`. - Start with an exponent of
`e`=2.- Replace
`b`with`b`^{e}(mod`N`) - If gcd(
`N`,`b`-1) is not 1, publish it as a factor and**STOP**. - Add 1 to
`e`.

- Replace

The principle behind this method is again Fermat's Little Theorem.
Assume for simplicity that the composite number `N` has two
prime factors, `p` and `q`, so `N`=`pq`.
Assume that `b` does not divide `p` (an unlikely
circumstance). Now we know that
`b`^{(p-1)}=1(mod `p`)`E` such that
(`p`-1) divides `E`, say
`E`=`k`(`p`-1)`b`^{E}
=`b`^{k}
=1(mod `p`)`p`
in advance. We get around that by *assuming* that the largest
prime factor of (`p`-1) is some bound `B`. If this
is true, then by the `B`^{th} step of the algorithm
we have computed
gcd(`N`,`b`^{2*3*...*B}-1).
We expect that
`b`^{2*3*...*B}=1(mod `p`)`b`^{2*3*...*B}-1=0(mod `p`)`b`^{2*3*...*B}-1(mod `N`)`q`. Thus we expect
gcd(`N`,`b`^{2*3*...*B}-1)
will be `q`.

This method usually works if you are willing to run enough
to large enough values of the exponent. Sometimes you are lucky,
(`p`-1) has only small factors, and factorization is
easy. Sometimes you are unlucky, (`p`-1) has a large
prime factor and it takes a lot of work to factor `N`.

A number `g` is *primitive* mod `p` if
the order of `g` mod `p` is (`p`-1).

If `p` is prime, Fermat's Little Theorem that, for any `g`
not divisible by `p`, `g`^{(p-1)} = 1 (mod `p`).
If `r` is not prime, say `r` = `pq`, then there are no
primitive elements mod `r`.
Conversely, it is (fairly) simple to prove that there are primitive elements mod any
prime `p`.

When `p` is prime there are many elements which are primitive mod `p`,
but there is no way to tell which elements are primitive short of finding their order.
Thus 2 is not primitive mod 7 (2^{1}=2, 2^{2}=4 and 2^{3}=8=1 (mod 7),
so the order of 2 mod 7 is 3), but 3 is primitive mod 7 (3^{1}=3, 3^{2}=9=2,
3^{3}=27=6, 3^{4}=81=4, 3^{5}=243=5
so the order of 3 mod 7 is 6).

(Strictly speaking, we have not given the correct definition of primitivity. A more
correct definition is that `g` is primitive mod `N` if the successive
powers of `g` mod `N` produce all integers less than `N` and
coprime to `N`. This sharper definition allows us to state that there are
primitive elements mod `N` if and only if `N` has the form `p`,
`p`^{n}, or 2`p`^{n}, where
`p` is prime.)

Now think about a number which is not prime, say `N`=`pq`, where both
`p` and `q` are prime. Consider some `a` which is coprime to
both `p` and `q` (in other words 0<`a`(mod `p`)<`p`,
and similarly for `q`). Thus from Fermat's Little Theorem we know that
`a`^{(p-1)(q-1)} = (`a`^{(p-1)})^{(q-1)} = 1^{(q-1)} = 1 (mod `p`)`a`^{(p-1)(q-1)} = 1 (mod `q`)`a`^{(p-1)(q-1)} = 1`a` mod `pq` must divide (`p`-1)(`q`-1).
In particular, as (`p`-1)(`q`-1) < (`pq`-1), we know that
if `N`=`pq` is not prime, then no element has order `N`-1.
We will use this to prove that a number is prime.

Among other things, we have shown the following, which is a sub-case of what is called Euler's Theorem. (One of the many theorems called "Euler's Theorem".)

**Theorem:** If `N`=`pq`, where both `p`
and `q` are prime, and if `b` is chosen so that gcd(`b`,`N`)=1,
then `b`^{(p-1)(q-1)} = 1 (mod `N`)

The Miller-Rabin primality test is, like the Fermat test, a pseudoprime test. It
uses a consequence of the structure of the group of units - that the only square
roots of 1 mod a prime `p` are 1 and -1, but mod any odd composite `N`
with `n` distinct prime factors, there are 2^{n} square
roots of 1. Thus, the Miller-Rabin test looks for square roots of 1, declaring
`N` composite if it finds one other than 1 or -1.

Thus far, every test has been a pseudoprime test - it positively identifies composite numbers but can identify primes at best with some likelihood of error. We can use what we now know to produce a test which proves that some number is prime. We know that:

- If
`N`is prime, there is some 0<`a`<`N`which is primitive, in other words has order exactly (`N`-1) mod`N`. - If
`N`is composite, then for any 0<`a`<`N`, either gcd(`b`,`N`)=1 (easily checked) or has order strictly less than (`N`-1) mod`N`.

Thus, if we can find an element `a` such that it has order exactly
(`N`-1) mod `N`, then we have proven that `N` is prime.
Our approach is brute force - we try successive values of `a` and see if
each is primitive. The work is not as bad as it sounds, though. For any `a`
we first confirm that `a`^{(N-1)}=1 (mod `N`).
Now we know that the true order of `a` is a divisor of `N`-1.
Checking `a`^{n} for every divisor `n` of
`N`-1 is not bad. There is a second trick we can apply though. If the
order of `a` is strictly less than `N`-1 then there must be some
prime divisor `p` of `N`-1 such that
`a`^{(N-1)/p}=1 (mod `N`)

- Factor (
`N`-1) =`p`_{1}^{e1}...`p`_{k}^{ek}, where every`p`_{i}is prime - For some sequence of
`a`'s (`a`= 2,3,4,5,... is probably good enough)- For each
`p`_{i}, if`a`^{(N-1)/p}=1 (mod`N`), reject`a`

`a`^{(N-1)/p}≠1 (mod`N`) for all`p`_{i}, then`a`is primitive and`N`is prime. - For each

Sometimes easier said than done, particularly the part where you factor (`N`-1).

Number theory is the study of more than just the integers. This is true both because various other algebraic structures are natural generalizations of the integers and also because they often give us information about the integers.

The first generalization of the integers we will touch on
is the Gaussian Integers, an
extension of the integers by *i*, the square root
of negative one.

The elements of the Gaussian Integers are all numbers of
the form `n`+`m`*i*, where both
`n` and `m` are integers. Addition and
multiplication are done in the usual way.

Things start getting interesting when we try to decide
what we mean by primes in the Gaussian Integers. The simple
observation that *i*)(2-*i*) = 4-(-1) = 5

Galois fields (also referred to as finite fields)
are extensions
of the integers modulo some prime `p`.

We start with the integers modulo `p` and some
polynomial which cannot be solved (ie factored) in this
field. Lets continue with a two tangible examples:

Our example will start with the integers
modulo 2 (containing only the two numbers 0 and 1). In keeping
with the usual notation we will call this **F**_{2}.
Given this field we now pull out of our hat the polynomial
`x`^{3}+`x`+1`a` and add this symbol to our original field.
The thing we now have is called **F**_{2}[`a`].
The elements of this thing have the form
`n`+`m``a`+`k``a`^{2}`a` as
`a`^{3} can be replaced by `a`-1`a` was presumed to be a root of the original
polynomial. In fact, as we are working modulo 2, so +1 and -1 are the
same, we can replace `a`^{3} by `a`+1.
Enumerating the 2^{3}=8 elements we have
`a`, `a`+1, `a`^{2},
`a`^{2}+1, `a`^{2}+`a`,
`a`^{2}+`a`+1}

Starting with the naive assumption that the defining polynomial
`x`^{4}+2`x`^{3}+3`x`^{2}+`x`+1

A bit more interesting is computing powers of elements. Just as in the case of computing with the integers, we can efficiently compute using the Russian Peasant algorithm of squaring and multiplying.

An example of this efficient multiplication is computing `a`^{26}:

- As 26 = 16+8+2, we can write in base 2, 26 = 11010
_{2} - Thus
`a`^{26}= (`a`^{16})(`a`^{8})(`a`^{2}) - Computing:
`a`^{2}:`a`^{2}=`a`^{2}`a`^{4}: (`a`^{2})^{2}= 3`a`^{3}+ 2`a`^{2}+ 4`a`+ 4`a`^{8}: (`a`^{4})^{2}= (3`a`^{3}+ 2`a`^{2}+ 4`a`+ 4)^{2}= 3`a`^{3}+ 3`a`^{16}: (`a`^{8})^{2}= (3`a`^{3}+ 3)^{2}= 2`a`^{2}+ 4`a` - Thus
`a`^{26}= (`a`^{16})(`a`^{8})(`a`^{2}) = (2`a`^{2}+ 4`a`)(3`a`^{3}+ 3)(`a`^{2}) = 3`a`^{3}+ 2`a`+ 4

There are a number of different ways in which we can use the basic Russian Peasant algorithm. The following slightly different approach is a little more difficult to understand but cleaner to program as it doesn't require temporary variables:

- As 26 = 11010
_{2}we can write 26 = 0+2(1+2(0+2(1+2(1)))) - and
`a`^{26}=`a`^{(0+2(1+2(0+2(1+2(1)))))}= (`a`((`a`(`a`)^{2})^{2})^{2})^{2} - Computing:

(`a`)^{2}=`a`^{2}`a`(`a`)^{2}=`a`^{3}(which is`a`^{3}=`a`(`a`)^{2})

(`a`^{3})^{2}= 3`a`^{3}+ 3`a`^{2}+`a`+ 4 (which is`a`^{6}= (`a`(`a`)^{2})^{2})

(3`a`^{3}+ 3`a`^{2}+`a`+ 4)^{2}= 2`a`^{2}+ 3 (which is`a`^{12}= ((`a`(`a`)^{2})^{2})^{2})`a`(2`a`^{2}+ 3) = (2`a`^{3}+ 3`a`) (which is`a`^{13}=`a`((`a`(`a`)^{2})^{2})^{2})

(2`a`^{3}+ 3`a`)^{2}= 3`a`^{3}+ 2`a`+ 4 (which is`a`^{26}= (`a`((`a`(`a`)^{2})^{2})^{2})^{2}) - Thus
`a`^{26}= 3`a`^{3}+ 2`a`+ 4

Finally, we are in a position to see if our polynomial is irreducible mod 5
(without actually factoring it - a problem for another day). The structure of
finite fields is similar to that of **Z**_{p} in a number
of ways. An analog of Fermat's Little Theorem exists, so for any element
`b`∈GF(5^{4}), we must have
`b`^{(54-1)}=1.
Another similarity is that any Galois field is guaranteed to have primitive elements.
In this case there are 5^{4}=625
elements, of which (5^{4}-1)=624 are non-zero. Thus, a primitive element
will have order exactly 624. As in the integer case, we check this by noting
that 624 has prime factorization 624 = (2^{4})(3)(13).

- Try
`a`∈GF(5^{4})`a`^{(54-1)}= 1`a`^{(54-1)/2}= 1 (so`a`is not primitive - in fact`a`has order 39)

- Try (
`a`+1)∈GF(5^{4})- (
`a`+1)^{(54-1)}= 1 - (
`a`+1)^{(54-1)/2}= 4 ≠ 1 - (
`a`+1)^{(54-1)/3}= 3`a`+2`a`^{3}≠ 1 - (
`a`+1)^{(54-1)/13}=`a`+`a`^{2}+`a`^{3}≠ 1

- (

So (`a`+1) has order exactly 624 and is primitive. Thus we have
shown that the polynomial `x`^{4}+2`x`^{3}+3`x`^{2}+`x`+1
is irreducible.

Some of the more interesting questions in computational
number theory involve large numbers. This can be a problem
as most languages and machines only support integers up to
a certain fixed size, commonly 2^{64} bits (about
1.6×10^{19}) or 2^{32} bits (about
4×10^{9}).

This limit is further reduced by the fact that most of the
algorithms require the computation of intermediate results
twice the size of the number of interest. For example, the
computation `mod(`

is commonly found, where each of `a`*`b`,`m`)`a` and `b`
are roughly the same size as `m`. For this computation
to work, the intermediate value `a`*`b` must
be less than the integer limit of the machine. Thus, on a 32-bit
machine, the largest value of `m` for which this can
be done is 2^{16}=65536. The examples in this paper
are written in JavaScript, which allows 53-bit computations
on any machine, thus allowing `m` to take on values
as large as 2^{26.5}=94,906,265.

The following are simple interpreted languages in which I have written basic number theory programs, suitable for classroom use and modification.

- JavaScript - JavaScript is an interpreted scripting language supported by modern web browsers. The JavaScript/ECMAscript standard supports the use of positive integers up to 53 bits in size.
- Python - The Python language is freely available for any architecture and is already installed on most OSX/Linux/UNIX systems. It allows the use of any size integer and is a very good fit for simple number theory. Programs implementing many of the algorithms are available here.
- HyperCard - The HyperCard program is available on Macintosh machines and allows the use of positive integers up to 63 bits in size. Programs implementing many of the algorithms are available here.
- QBASIC - The QBASIC interpreter is available on MS-DOS and Windows machines shipped during the last few years and allows the use of positive integers up to 31 bits in size. Programs implementing many of the algorithms are available here.
- TI Calculator - The Texas Instruments family of programmable calculators has a simple programming language and allows the use of integers up to 46 bits in size. Programs implementing many of the algorithms are available here.

More details on programming for number theory can be found in a related document here.

*A Concrete Introduction to Higher Algebra*, Lindsay Childs, Springer-Verlag, 1979. A good coverage of the same material as this document with emphasis on algorithms but no example code.*Primes and Programming*, Peter Giblin, Cambridge Univ Press, 1993. A good introduction to number theory with a strong emphasis on algorithms - contains PASCAL code implementing most algorithms.*Factorization and Primality Testing*, David Bressoud, Springer-Verlag, 1989. Covers most current factoring and primality testing algorithms, as well as those elements of number theory needed for them. Contains pseudocode implementations of most algorithms.

Robert Campbell Last modified: Jan 29, 2010