Number Theory

with Python


Contents

  1. Integer Operations
  2. Number Theory Operations
  3. Gaussian Integer Operations
  4. Polynomial Operations
  5. Finite Field Operations
  1. Packages for Number Theory in Python
  2. References & Links

[Prev] [Python Intro] [Next]

The Python programming language has basic commands which implement integer arithmetic. More involved number theory will require us to write short programs and modules in Python. (Given the option, the best way to do number theory in Python is to use SAGE, a Python-based symbolic algebra system.) A general introduction to Python use and where it can be found/installed at UMBC can be found in a separate document.

A particular strength of Python for number theory is its native support for arbitrary sized integers. These numbers are entered by writing the number followed by the letter "L" (for example 1234512L). Starting with Python 2.x there is an automatic conversion from regular integers to long integers when the size of the number is large enough.

Integer Operations

The basic arithmetic operations, +, -, *, /, are available from the keyboard. Note that when used with integers, the division operator returns the greatest integer less than the exact result. (Note that this is done consistently, even for negative values - this differs from many programming languages.) The modular reduction operator is represented by the % operator (eg 7 % 2 returns 1. Again, this differs from many languages in that, if the modulus is positive the result is always positive.) Exponentiation is represented with the ** symbol. The pow() function can be used for exponentiation if called with two parameters, or efficient modular exponentiation if called with three parameters. (pow(a,b,c) returns the same result as (a**b % c), but is much more efficient.)

Examples:

[There is a fine point with the division operation which may become important with Python 3.0. Currently both the "/" and the the "//" operation on integers return the floor of the quotient. In future it may be possible to use "/" to produce the floating point quotient, and "//" to produce the floor of the quotient. A full treatment is in PEP 238.]

Number Theory Operations

Beyond the basic arithmetic operations Python has few natively implemented number theory functions. It is easy to write Python code which implements these functions though. Python source code for these operations is here.

The gcd function can be computed by hand as a succession of modular reductions. An example is this computation of the gcd of the two integers 12345678 and 87654321, whose gcd is 9. >>> a = 12345678; b = 87654321 >>> b = b%a; b 1234575 >>> a = a%b; a 1234503 >>> b = b%a; b 72 >>> a = a%b; a 63 >>> b = b%a; b 9 >>> a = a%b; a 0

A simple recursive implementation of the gcd function is:

	  def gcd(a,b):
	     if a == 0:
	        return b
	     return gcd(b % a, a)

The extended gcd function returns both gcd(a,b) and those two integers u and v such that gcd(a,b) = ua+vb. An example of computing this by hand with Python is: >>> a = [12345678,1,0]; b = [87654321,0,1]; >>> q = (b[0]//a[0]); b = [b[0]%a[0], b[1]-q*a[1], b[2]-q*a[2]]; b [1234575, -7, 1] >>> q = (a[0]//b[0]); a = [a[0]%b[0], a[1]-q*b[1], a[2]-q*b[2]]; a [1234503, 64, -9] >>> q = (b[0]//a[0]); b = [b[0]%a[0], b[1]-q*a[1], b[2]-q*a[2]]; b [72, -71, 10] >>> q = (a[0]//b[0]); a = [a[0]%b[0], a[1]-q*b[1], a[2]-q*b[2]]; a [63, 1217359, -171459] >>> q = (b[0]//a[0]); b = [b[0]%a[0], b[1]-q*a[1], b[2]-q*a[2]]; b [9, -1217430, 171469] >>> q = (a[0]//b[0]); a = [a[0]%b[0], a[1]-q*b[1], a[2]-q*b[2]]; a [0, 9739369, -1371742] >>> 12345678*(-1217430)+87654321*171469 # Confirm that gcd = u*a + v*b 9 This algorithm is designed so that at any step of the computation you can confirm that for either of the two vectors a and b we have the property v[0] = 12345678*v[1] + 87654321*v[2].

For positive input values this can be implemented by:

	  def xgcd(a,b):
	     a1=1; b1=0; a2=0; b2=1
	     while (1):
	        quot = -(a // b)
	        a = a % b
	        a1 = a1 + quot*a2; b1 = b1 + quot*b2
	        if(a == 0):
	           return [b, a2, b2]
	        quot = -(b // a)
	        b = b % a;
	        a2 = a2 + quot*a1; b2 = b2 + quot*b1
	        if(b == 0):
	           return [a, a1, b1]

A classical way to find small primes is the Sieve of Eratosthenes. We start with a list of the integers starting with 2. Pull off 2 as the first prime and remove all the multiples of 2 from the list. The first element of the list is now 3. Pull off the prime 3 and remove all the multiples of 3 from the list. The first element of the list is now 5 (as 4 was removed as a multiple of 2). Continue this process until the elements of the list are exhausted, recovering all the primes on the list. This can be implemented using Python's list manipulation routines: >>> sieve = range(2,200); primes = [] # Initialize the sieve, starting at 2 >>> nextprime = sieve[0] # The first element of the sieve is prime >>> sieve = [sieve[i] for i in range(len(sieve)) if (sieve[i]%nextprime)!=0] >>> print nextprime, sieve 2 [3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, ... >>> nextprime = sieve[0] >>> sieve = [sieve[i] for i in range(len(sieve)) if (sieve[i]%nextprime)!=0] >>> print nextprime, sieve 3 [5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, 37, 41, 43, ...

Recall that modular exponentiation is implemented by the command pow(b, e, n), which efficiently computes the value of be (mod n). We can use this to implement a simple primality test (Fermat's test). Choosing an arbitrary base b (2 will work quite well), compute pow(b,n-1, n). If the result is equal to 1 it supports a conclusion that n is prime, but any other result is proof that n is composite. We can use this to test the primality of several large integers: >>> pow(2,100,101) # Evidence that 101 is prime 1 >>> pow(3,100,101) # More evidence that 101 is prime 1 >>> pow(2,110,111) # Proof that 111 is composite 4

A more powerful primality test is Euler's test, which tests whether some arbitrary base b has modular powers with a specific property. We first take its n-1 power and confirm that we get a result of 1 (the Fermat test). We then divide this exponent repeatedly by 2 (essentially taking square roots). If we ever observe a square root of 1 other then 1 or -1 we know that n is composite. This can be implemented as:

	    def isprimeE(n,b):
	       if (not pow(b,n-1,n)==1): return False
	       r = n-1
	       while (r % 2 == 0): r //= 2
	       c = pow(b,r,n)
	       if (c == 1): return True
	       while (1):
	          if (c == 1): return False
	          if (c == n-1): return True
	          c = pow(c,2,n)

Finally, Pollard's rho algorithm for factorization can be implemented as:

	    def factorPR(n):
	       for slow in [2,3,4,6]:
	          numsteps=2*math.floor(math.sqrt(math.sqrt(n))); fast=slow; i=1
	          while i<numsteps:
		     slow = (slow*slow + 1) % n
		     i = i + 1
		     fast = (fast*fast + 1) % n
		     fast = (fast*fast + 1) % n
		     g = gcd(fast-slow,n)
		     if (g != 1):
		        if (g == n):
		           break # try a new base
		        else:
		           return g
	       return 1 # fail to factor

Pointers to several number theory modules written in Python are available in the section Packages for Number Theory in Python.

Gaussian Integer Operations

Python implements a complex class. This class is implemented based on floating point values, so we re-implement it as a Python class, requiring at least Python 2.x. The Python source for this class is available here.

An example of the use of the Gaussian Integers class is here: >>> from gaussint import * >>> a = GaussInt(1,2) >>> m = GaussInt(7,2) >>> a + m (8 + 4i) >>> a * m (3 + 16i) >>> 5*a # Multiply an integer times a GaussInt (5 + 10i) >>> a.powmod(12,m) # Compute a**12 mod m (-2 + 2i)

As the Gaussian Integers are a Euclidean ring it they have all sorts of lovely things like GCDs and unique factorization. >>> GaussInt(-2544,1788).xgcd(GaussInt(2241,5238)) ((21 + -12i), (148 + 168i), (73 + -98i)) >>> GaussInt(-2544,1788)*GaussInt(148,168) + GaussInt(2241,5238)*GaussInt(73,-98) (21 + -12i)

Here we computed that (21-12i) is the gcd of (-2544+1788i) and (2241+5238i) and confirmed that (21-12i) = (-2544+1788i)(148+168i) + (2241+5238i)(73-98i).

Finally, we factor a Gaussian integer and check the primality of one of the factors. >>> thefacts = GaussInt(12345,23145).factors(); thefacts [GaussInt(360, -1183), GaussInt(3, 0), GaussInt(1, 2), GaussInt(2, 1), GaussInt(1, 1)] >>> thefacts[0].isprime() # Confirm that the first factor is prime True >>> map(str,thefacts) # Make the format more human-readable ['(360 + -1183i)', '(3 + 0i)', '(1 + 2i)', '(2 + 1i)', '(1 + 1i)']

Find all the prime Gaussian integers in the first quadrant (real and imaginary components positive) with both components less than 10: >>> for i in range(0,20): ... for j in range(0,i): ... n = GaussInt(i,j) ... if n.isprime(): print n

Polynomial Operations

Finite Field Operations

Python does not directly support finite fields. My simple implementation of finite fields is the class FiniteField, and the source for this class is available here.

Here is an extended example of the use of my FiniteField class, working in the finite field GF(35), constructed as Z3(x)/<2+2x+x+x4+x5>. Note that the polynomial must be monic, with a leading coefficient of 1 (so this leading coefficient is assumed, and not part of the input list of coefficients).

	>>> from finitefield import *
	>>> GF243 = FiniteField(3,[2,2,1,0,1])  # GF(3^5) = Z_3(x)/<2+2x+x^2+x^4+x^5>
	>>> str(GF243)                          # Print friendly format
	'GF(3^5)'
	>>> repr(GF243)                         # Internal representation
	'FiniteField(3, [2, 2, 1, 0, 1])'
	

We check to see if 2+2x+x+x4+x5 is a primitive polynomial. (This is the same as asking if "x" has full order, i.e. order (243-1) = 242 = (2)(112).)

	>>> x = FiniteFieldElt(GF243,[0,1])     # Construct x = 0+1*x in GF243
	>>> str(x**2)                           # Various powers of x
	'[0, 0, 1, 0, 0]'
	>>> str(x**3)
	'[0, 0, 0, 1, 0]'
	>>> str(x**4)
	'[0, 0, 0, 0, 1]'
	>>> str(x**5)
	'[1, 1, 2, 0, 2]'
	>>> str(x**242)                         # So x(243-1) = 1 as expected
	'[1, 0, 0, 0, 0]'
	>>> str(x**(242/2))                     # x(242/2) = 1, so not primitive
	'[1, 0, 0, 0, 0]'
	>>> a = 1+2*x+x**3                      # But a = 1+2x+x3 has full order
	>>> a
	FiniteFieldElt(FiniteField(3, [2, 2, 1, 0, 1]),[1, 2, 0, 1, 0])
	>>> str(a**(242/2))
	'[2, 0, 0, 0, 0]'
	>>> str(a**(242/11))
	'[1, 2, 1, 2, 0]'
	>>> str(a**242)
	'[1, 0, 0, 0, 0]'
	

In addition to my finite field class, there are a number of other classes implementing finite fields. There are no guarantees, as I haven't used any of these other classes. (Again, the most complete support for finite fields is found in the Python based SAGE language.)

Appendices:

Packages for Number Theory in Python

A number of authors have implemented packages for number theory operations in Python. A partial list is:

References & Links


[Prev] [Python Intro] [Next]
Robert Campbell
28 Feb, 2013