Making a random number generator

A psuedo-random number sequence is a sequence of numbers that appears random to the casual observer, but is generated by a deterministic, repeatable algorithm. One of the simplest algorithms used to make a pseudo-random number sequence is the linear congruential random number generator. This is a simple mathematical formula that generates a sequence of integers. To form the sequence we pick an arbitrary integer as the seed for the sequence. We then substitute that starting value into a simple linear function to generate the next number in the sequence. Most often that simple linear function uses two large prime numbers as its coefficients.

x1 = 12553 x0 + 15401

We repeat that process to make a whole sequence of numbers. If we did just that, the sequence would typically march off toward positive or negative infinity. To prevent that, we mod each number by a large prime integer b on each round.

x1 = (12553 x0 + 15401) mod 6133

This guarantees that the numbers stay trapped in the range from 0 to b-1. Often, the sequence of integers generated by this process is not immediately useful to us, because, say, we may be interested in generating a sequence of numbers in the range from 0 to c-1 for some other integer c. That is easy to fix. If c < b, we can compute x % c for each number we generate and save the list of x % c values as we generate them.

Here is a simple example program in Python that uses this process to generate and print a list of integers in the range from 0-99:

x = 52 # The seed value for the pseudo-random sequence

for n in range(100):
    print(x % 100)
    x = (12553*x+15401)%6133

Next, we would like some way to encapsulate the random number generation process so we can start and stop the process on demand. The key to doing this is to place the current value of x in a safe location so we can access it on demand to generate the next x in the sequence. This is a very good application for the class concept in Python. Python classes store state information in the form of member variables. In this case, we are going to make a random sequence class that stores the current value of the x variable as a member variable. Next, we will also need to equip our class with at least one method that clients can use to fetch a new number from the sequence.

Here is the code for a Python class that can do these things.

class RndSeq():
    def __init__(self,seed = 0):
      self.x = seed
      
    def nextInt(self,N):
        """Return a random integer in the range from 0 to N-1."""
        if N > 5000:
            print("RndSeq.nextInt requires its argument to be less than or equal to 5000.")
            return 0
        self.x = (12553*self.x+15401)%6133
        return self.x % N

This code shows the bare minimum set of things needed to make a useful class. The class contains an __init__ function that is used to initialize objects. In this case, we need to initialize the member variable x with a seed value. Since classes store useful information in member variables, we also need member functions to access and modify those values. In this case, we only need one member function, nextInt(), that will compute and store the next x in the sequence and return a result based on the new value of the x hidden in the object.

This code is stored in a file named 'rndseq.py'. We can import and use this class in a simple test program.

from rndseq import RndSeq

seed = input("Enter a seed value for the sequence: ")
s = RndSeq(int(seed))
for n in range(100):
    print(s.nextInt(100))

A complex number class

For our next example of a Python class I am going to construct a class that represents complex numbers. This is a very simple class with only two data members, the real and imaginary parts of the number. The class contains methods to read those data members and also compute and return the modulus of the complex number.

Because we will want to be able to do arithmetic with our complex number objects, I have also supplied this class with a number of operator overloads to implement the various operations.

When you try to do an operation with a couple of data items, such as adding two Complex objects, Python will try to make sense of the operation by converting the code into something that looks like a method call. Here is an example. Suppose we write the code

a = Complex(2,3)
b = Complex(3,4)
c = a*b

When the Python interpret reaches the last statement it try to make sense of the code there by first converting the code into an alternative form:

c = a.__mul__(b)

This has the form of a method call being applied to the object a with a parameter of b. To make this work, Python will inspect the object stored in a to see if it has a method with that name. If it does, it will call that method.

Here is the code for the __mul__ method in the Complex class.

def __mul__(self,other):
  r = self.real*other.getReal() - self.imag*other.getImag()
  i = self.real*other.getImag() + self.imag*other.getReal()
  return Complex(r,i)

This method implements complex multiplication using the standard mathematical rule for that operation. After computing the real and imaginary parts of the product, the function constructs and returns a Complex object with those parts.

Another handy method to provide in a class is the __str__ method. This method will get invoked when you try to use str() to convert an object to text. Here is an example.

x = Complex(2,1)
print(str(x))

When the Python interpreter encounters the second statement it will attempt to rewrite it as

print(x.__str__())

This will succeed if x is an object and the class that x belongs to implements the __str__ method.

Here is the code for the __str__ method for the Complex class.

def __str__(self):
  if self.imag < 0:
    return str(self.real)+"-"+str(-self.imag)+"i"
  elif self.imag > 0:
    return str(self.real)+"+"+str(self.imag)+"i"
  else:
    return str(self.real)

Code for the class

Here now is the complete source code for the Complex class. This code lives in a separate source file named 'complex.py'

import math

class Complex():
    def __init__(self,real,imag = 0):
        self.real = real
        self.imag = imag
        
    def getReal(self):
        return self.real
        
    def getImag(self):
        return self.imag
        
    def getAbs(self):
        return math.sqrt(self.real*self.real+self.imag*self.imag)
        
    def __add__(self,other):
        return Complex(self.real+other.getReal(),self.imag+other.getImag())
        
    def __sub__(self,other):
        return Complex(self.real-other.getReal(),self.imag-other.getImag())
        
    def __mul__(self,other):
        r = self.real*other.getReal() - self.imag*other.getImag()
        i = self.real*other.getImag() + self.imag*other.getReal()
        return Complex(r,i)
        
    def __truediv__(self,other):
        d = other.getReal()*other.getReal() + other.getImag()*other.getImag()
        r = (self.real*other.getReal() + self.imag*other.getImag()) / d
        i = (self.imag*other.getReal() - self.real*other.getImag()) / d
        return Complex(r,i)
        
  def __str__(self):
    if self.imag < 0:
      return str(self.real)+"-"+str(-self.imag)+"i"
    elif self.imag > 0:
      return str(self.real)+"+"+str(self.imag)+"i"
    else:
      return str(self.real)

Application: Newton's method with complex numbers

Here now is an application program that makes use of the complex number class. This program is an implementation of Newton's method that takes advantage of the fact that Newton's method works just fine in the presence of complex numbers.

This program attempts to find the roots of the polynomial

x4 + 2 x3 - 6 x2 + 8 x + 80

This polynomial has no real roots: all four of its roots lie in the complex plane. If you run this program with a complex number as the starting guess it will converge to one of those four complex roots. If you start the program with a real number as the starting guess, Newton's method falls into an infinite loop and does not converge. When that happens you will have to stop the program by choosing the stop command in your IDE or by hitting the control-C key combination in the terminal.

# Program to find the roots of f(x) = x^4 + 2 x^3 - 6 x^2 + 8 x + 80
# by Newton's method. This polynomial has four complex roots.
from complex import Complex

# Define the function and its derivative
def f(x):
    return (((x+Complex(2))*x+Complex(-6))*x+Complex(8))*x + Complex(80)
    
def fp(x):
    return ((Complex(4)*x+Complex(6))*x+Complex(-12))*x+Complex(8)
    
r = input("Enter the real part for your starting guess: ")
i = input("Enter the imaginary part for your starting guess: ")
x = Complex(float(r),float(i))

y = f(x)
tolerance = 10e-6

while y.getAbs() > tolerance:
    # Compute f'(x)
    yp = fp(x) 
    # Do one round of Newton's method
    x = x - y/yp
    y = f(x)
    
print("The root estimate is " + str(x))

An important thing to note here is that all of the numbers that appear in the definitions for f(x) and its derivative have to be complex numbers. We need to do this because the operator functions that allow us to do arithmetic with complex number objects require that both operands in the operation by complex number objects.