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.

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. 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]

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.

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)
	>>> GaussInt(-2544,1788).xgcd(GaussInt(2241,5238))
	((21 + -12i), (148 + 168i), (73 + -98i))
	# So (21-12i) is the gcd of (-2544+1788i) and (2241+5238i)
	# and (21-12i) = (-2544+1788i)(148+168i) + (2241+5238i)(73-98i)
	# We check this:
	>>> GaussInt(-2544,1788)*GaussInt(148,168) + GaussInt(2241,5238)*GaussInt(73,-98)
	(21 + -12i)

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, 2005