07. Modular Arithmetic, Linear Congruential Generators, and Pseudo-Random Numbers

Inference Theory 1

©2018 Raazesh Sainudiin. Attribution 4.0 International (CC BY 4.0)

  • What's this all about?
  • Modular Arithmetic
  • Linear Congruential Generators
  • More Sophisticated Pseudo-Random Number Generators
  • Accumulating sequences with numpy.cumsum
  • Simulating a Drunkard's Walk

What's this all about?

Question:

How can we produce realisations from Uniform$(0,1)$, the fundamental random variable?

i.e., how can we produce samples $(x_1, x_2, \ldots, x_n)$ from $X_1, X_2, \ldots, X_n$ $\overset{IID}{\thicksim}$ Uniform$(0,1)$?

What is Sage doing when we ask for random()?

In [9]:
random()
Out[9]:
0.4744053394643051

Answer:

Modular arithmetic and number theory gives us pseudo-random number generators.

Question:

What can we do with samples from a Uniform$(0,1)$ RV? Why bother?

Answer:

We can use them to sample or simulate from other, more complex, random variables. These simulations can be used to understand real-world phenomenon such as:

  • modelling traffic queues on land, air or sea for supply chain management
  • estimate missing data in Statistical survey to better manage or administer
  • helping a Hospital to manage critical care for pre-term babies
  • helping NZ's Department of Conservation to minimise the extinction probability of various marine organisms
  • help the Government find if certain fishing boats are illegally under-reporting their catches from NZ's waters
  • find cheaper air tickets for a vacation
  • various physical systems can be modeled. See http://en.wikipedia.org/wiki/Monte_Carlo_method for a bigger picture.

See how the modern history of Monte carlo starts in Los Alamos!

In [1]:
def showURL(url, ht=500):
    """Return an IFrame of the url to show in notebook with height ht"""
    from IPython.display import IFrame
    return IFrame(url, width='95%', height=ht) 
showURL('http://en.wikipedia.org/wiki/Monte_Carlo_method',600)
Out[1]:

The starting point for Monte Carlo methods is modular arithmetic ...

Modular arithmetic

Modular arithmetic (arithmetic modulo $m$) is a central theme in number theory and is crucial for generating random numbers from a computer (i.e., machine-implementation of probabilistic objects). Being able to do this is essential for computational statistical experiments and methods that help do this are called Monte Carlo methods. Such computer-generated random numbers are technically called pseudo-random numbers.

In this notebook we are going to learn to add and multiply modulo $m$ (this part of our notebook is adapted from William Stein's SageMath worksheet on Modular Arithmetic for the purposes of linear congruential generators). If you want a more thorough treatment see http://en.wikipedia.org/wiki/Modular_arithmetic.

In [2]:
showURL("http://en.wikipedia.org/wiki/Modular_arithmetic",500)
Out[2]:

Remember when we talked about the modulus operator %? The modulus operator gives the remainder after division:

In [3]:
14%12 # "14 modulo 12" or just "14 mod 12"
Out[3]:
2
In [4]:
zeroto23Hours = range(0,24,1)  # 0,1,2...,23 hours in the 24-hours clock
print(zeroto23Hours)
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]
In [5]:
[x%12 for x in range(0,24,1)]  # x%12 for the 24 hours in the analog clock
Out[5]:
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
In [6]:
set([x%12 for x in range(0,24,1)]) # unique hours 0,1,2,...,11 in analog clock by making a set out of the list
Out[6]:
{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11}

Arithmetic modulo $m$ is like usual arithmetic, except every time you add or multiply, you also divide by $m$ and return the remainder. For example, working modulo $m=12$, we have:

$$8 + 6 = 14 = 2$$

since $2$ is the remainder of the division of $14$ by $12$.

Think of this as like the hours on a regular analog clock. We already do modular addition on a regular basis when we think about time, for example when answering the question:

Question: If it is 8pm now then what time will it be after I have spent the 6 hours that I am supposed to dedicate this week toward this course?

Answer: <p style="color:#F8F9F9";>2am.</p>

Modular addition and multiplication with numbers modulo $m$ is well defined, and has the following properties:

  • $a + b = b+a$ (addition is commutative)
  • $a \cdot b = b \cdot a$ (multiplication is commutative)
  • $a\cdot (b + c) = a\cdot b + a\cdot c$ (multiplication is distributive over addition)
  • If $a$ is coprime to $m$ (i.e., not divisible by any of the same primes), then there is a unique $b$ (mod $m$) such that $a\cdot b=1$.

Let us make a matrix of results from addition and multiplication modulo 4.

Before that realise that arithmetic over a set $\mathbb{M}$ is merely a function from the Cartesian product $\mathbb{M} \times \mathbb{M}$ of all pairs in $\mathbb{M}$ to $\mathbb{M}$:

$$ \star : \mathbb{M} \times \mathbb{M} \to \mathbb{M} $$

where $\star$ is usually $+$ for "addition" or $*$ for "multiplication".

Thus, we can think of arithmetic as a set of three-tuples: $$(x_1, x_2, x_3) \in \mathbb{M}^3 := \mathbb{M} \times \mathbb{M} \times \mathbb{M}, \quad \text{where} \quad x_3 = x_1 \star x_2. $$

In [7]:
mySageMatrix = matrix(2,3,[1, 2, 3, 4, 5, 6]) # this is how you make a 2X3 matrix in Sage
mySageMatrix
Out[7]:
[1 2 3]
[4 5 6]
In [8]:
type(mySageMatrix) # the type says a lot
Out[8]:
<type 'sage.matrix.matrix_integer_dense.Matrix_integer_dense'>
In [ ]:
mySageMatrix.

Let us list all the three-tuples for addition modulo 4 next.

In [9]:
m=4;
# list (i,j, (i+j) mod m) as (i,j) range in [0,1,2,3]
[(i,j,(i+j)%m) for i in range(m) for j in range(m)]
Out[9]:
[(0, 0, 0),
 (0, 1, 1),
 (0, 2, 2),
 (0, 3, 3),
 (1, 0, 1),
 (1, 1, 2),
 (1, 2, 3),
 (1, 3, 0),
 (2, 0, 2),
 (2, 1, 3),
 (2, 2, 0),
 (2, 3, 1),
 (3, 0, 3),
 (3, 1, 0),
 (3, 2, 1),
 (3, 3, 2)]
In [10]:
[(i+j)%m for i in range(m) for j in range(m)] # just the third element of the three-tuple from above
Out[10]:
[0, 1, 2, 3, 1, 2, 3, 0, 2, 3, 0, 1, 3, 0, 1, 2]

Since the arguments to addition modulo m are fixed in a square array $[0,1,2\ldots,m-1]\times[0,1,2\ldots,m-1]$ we can simply focus on the image of the arithmetic modulo $m$ operation and place it over the square array using the following matrix call:

In [11]:
#addition mod m
matrix(m,m,[(i+j)%m for i in range(m) for j in range(m)])
Out[11]:
[0 1 2 3]
[1 2 3 0]
[2 3 0 1]
[3 0 1 2]
In [12]:
# multiplication mod m
matrix(m,m,[(i*j)%m for i in range(m) for j in range(m)])
Out[12]:
[0 0 0 0]
[0 1 2 3]
[0 2 0 2]
[0 3 2 1]

Visualize Modular Arithmetic

In the following interactive image (created by William Stein) we make an addition and multiplication table modulo $m$, where you can control $m$ with the slider. Try changing $m$ and seeing what effect it has:

In [13]:
import matplotlib.cm
cmaps = ['gray'] + list(sorted(matplotlib.cm.datad.keys()))

@interact
def mult_table(m=(4,(1..200)), cmap=("Color Map",cmaps)):
    Add = matrix(m,m,[(i+j)%m for i in range(m) for j in range(m)])
    Mult = matrix(m,m,[(i*j)%m for i in range(m) for j in range(m)])
    print "Addition and multiplication table modulo %s"%m
    show(graphics_array( (plot(Add, cmap=cmap),plot(Mult, cmap=cmap))),figsize=[10,5])

Look at this for a while. Make $m$ close to 200 slowly and see how the outputs of the two arithmetic operations change.

Answer the following questions (which refer to the default colour map) to be sure you understand the table:

  • Question: Why is the top row and leftmost column of the multiplication table always black?
  • Question: Why is there an anti-diagonal block line in the addition table (the left-most table)?
  • Question: Why are the two halves of each table (the two pictures) symmetric about the diagonal?

You can change the colour map if you want to.

YouTry

You should be able to understand a bit of what is happening here. See that there are two list comprehensions in there. Take the line matrix(m,m,[(i+j)%m for i in range(m) for j in range(m)]). The list comprehension part is [(i+j)%m for i in range(m) for j in range(m)]. Let's pull it out and have a look at it. We have set the default modulo m to 4 but you could change it if you want to:

In [14]:
m = 4
listFormAdd = [(i+j)%m for i in range(m) for j in range(m)]
listFormAdd
Out[14]:
[0, 1, 2, 3, 1, 2, 3, 0, 2, 3, 0, 1, 3, 0, 1, 2]

This list comprehension doing double duty: remember that a list comprehension is like a short form for a for loop to create a list? This one is like a short form for one for-loop nested within another for-loop.

We could re-create what is going on here by making the same list with two for-loops. This one uses modulo m = 4 again.

In [23]:
m = 4
listFormAddTheHardWay = []
for i in range(m):
    for j in range(m):
        listFormAddTheHardWay.append((i+j)%m)
listFormAddTheHardWay
Out[23]:
[0, 1, 2, 3, 1, 2, 3, 0, 2, 3, 0, 1, 3, 0, 1, 2]

Notice that the last statement in the list comprehension, for j in range(m), is the inner loop in the nested for-loop.

The next step that Stein's interactive image is to make a matrix out of the list. We won't be doing matrices in detail in this course (we decided to concentrate on numpy.arrays and multi-dimensional numpy.ndarrays with floating-point numbers for numerical computing instead), but if you are interested, this statement uses the matrix function to make a matrix out of the list named listFormAdd. The dimensions of the matrix are given by the first two arguments to the matrix function: (m, m,...).

In [24]:
matrixForm = matrix(m,m, listFormAdd)
matrixForm
Out[24]:
[0 1 2 3]
[1 2 3 0]
[2 3 0 1]
[3 0 1 2]

Optionally, you can find out more about matrices from the documentation.

In [26]:
?matrix
In [27]:
%%sh
cat ~/all/software/sage/SageMath/src/sage/matrix/constructor.pyx
"""
General matrix Constructor
"""

#*****************************************************************************
#       Copyright (C) 2005 William Stein <wstein@gmail.com>
#       Copyright (C) 2016 Jeroen Demeyer <jdemeyer@cage.ugent.be>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#                  http://www.gnu.org/licenses/
#*****************************************************************************
from __future__ import absolute_import

import types
from .matrix_space import MatrixSpace
from sage.rings.ring import is_Ring
from sage.rings.all import ZZ, RDF, CDF
from sage.arith.srange import srange
from sage.structure.coerce cimport is_numpy_type, py_scalar_parent
from sage.structure.element cimport Vector
from sage.structure.sequence import Sequence
from sage.misc.long cimport pyobject_to_long


class MatrixFactory(object):
    """
    Create a matrix.

    This implements the ``matrix`` constructor::

        sage: matrix([[1,2],[3,4]])
        [1 2]
        [3 4]

    It also contains methods to create special types of matrices, see
    ``matrix.[tab]`` for more options. For example::

        sage: matrix.identity(2)
        [1 0]
        [0 1]

    INPUT:

    The matrix command takes the entries of a matrix, optionally
    preceded by a ring and the dimensions of the matrix, and returns a
    matrix.

    The entries of a matrix can be specified as a flat list of
    elements, a list of lists (i.e., a list of rows), a list of Sage
    vectors, a callable object, or a dictionary having positions as
    keys and matrix entries as values (see the examples). If you pass
    in a callable object, then you must specify the number of rows and
    columns. You can create a matrix of zeros by passing an empty list
    or the integer zero for the entries.  To construct a multiple of
    the identity (`cI`), you can specify square dimensions and pass in
    `c`. Calling matrix() with a Sage object may return something that
    makes sense. Calling matrix() with a NumPy array will convert the
    array to a matrix.

    The ring, number of rows, and number of columns of the matrix can
    be specified by setting the ``ring``, ``nrows``, or ``ncols``
    keyword parameters or by passing them as the first arguments to the
    function in specified order. The ring defaults to ``ZZ`` if it is
    not specified and cannot be determined from the entries. If the
    number of rows and columns are not specified and cannot be
    determined, then an empty 0x0 matrix is returned.

    INPUT:

    -  ``ring`` -- the base ring for the entries of the
       matrix.

    -  ``nrows`` -- the number of rows in the matrix.

    -  ``ncols`` -- the number of columns in the matrix.

    -  ``sparse`` -- create a sparse matrix. This defaults to ``True``
       when the entries are given as a dictionary, otherwise defaults to
       ``False``.

    - ``entries`` -- see examples below.

    OUTPUT:

    a matrix

    EXAMPLES::

        sage: m = matrix(2); m; m.parent()
        [0 0]
        [0 0]
        Full MatrixSpace of 2 by 2 dense matrices over Integer Ring

    ::

        sage: m = matrix(2,3); m; m.parent()
        [0 0 0]
        [0 0 0]
        Full MatrixSpace of 2 by 3 dense matrices over Integer Ring

    ::

        sage: m = matrix(QQ,[[1,2,3],[4,5,6]]); m; m.parent()
        [1 2 3]
        [4 5 6]
        Full MatrixSpace of 2 by 3 dense matrices over Rational Field

    ::

        sage: m = matrix(QQ, 3, 3, lambda i, j: i+j); m
        [0 1 2]
        [1 2 3]
        [2 3 4]
        sage: m = matrix(3, lambda i,j: i-j); m
        [ 0 -1 -2]
        [ 1  0 -1]
        [ 2  1  0]

    ::

        sage: matrix(QQ, 2, 3, lambda x, y: x+y)
        [0 1 2]
        [1 2 3]
        sage: matrix(QQ, 5, 5, lambda x, y: (x+1) / (y+1))
        [  1 1/2 1/3 1/4 1/5]
        [  2   1 2/3 1/2 2/5]
        [  3 3/2   1 3/4 3/5]
        [  4   2 4/3   1 4/5]
        [  5 5/2 5/3 5/4   1]

    ::

        sage: v1=vector((1,2,3))
        sage: v2=vector((4,5,6))
        sage: m = matrix([v1,v2]); m; m.parent()
        [1 2 3]
        [4 5 6]
        Full MatrixSpace of 2 by 3 dense matrices over Integer Ring

    ::

        sage: m = matrix(QQ,2,[1,2,3,4,5,6]); m; m.parent()
        [1 2 3]
        [4 5 6]
        Full MatrixSpace of 2 by 3 dense matrices over Rational Field

    ::

        sage: m = matrix(QQ,2,3,[1,2,3,4,5,6]); m; m.parent()
        [1 2 3]
        [4 5 6]
        Full MatrixSpace of 2 by 3 dense matrices over Rational Field

    ::

        sage: m = matrix({(0,1): 2, (1,1):2/5}); m; m.parent()
        [  0   2]
        [  0 2/5]
        Full MatrixSpace of 2 by 2 sparse matrices over Rational Field

    ::

        sage: m = matrix(QQ,2,3,{(1,1): 2}); m; m.parent()
        [0 0 0]
        [0 2 0]
        Full MatrixSpace of 2 by 3 sparse matrices over Rational Field

    ::

        sage: import numpy
        sage: n = numpy.array([[1,2],[3,4]],float)
        sage: m = matrix(n); m; m.parent()
        [1.0 2.0]
        [3.0 4.0]
        Full MatrixSpace of 2 by 2 dense matrices over Real Double Field

    ::

        sage: v = vector(ZZ, [1, 10, 100])
        sage: m = matrix(v); m; m.parent()
        [  1  10 100]
        Full MatrixSpace of 1 by 3 dense matrices over Integer Ring
        sage: m = matrix(GF(7), v); m; m.parent()
        [1 3 2]
        Full MatrixSpace of 1 by 3 dense matrices over Finite Field of size 7

    ::

        sage: g = graphs.PetersenGraph()
        sage: m = matrix(g); m; m.parent()
        [0 1 0 0 1 1 0 0 0 0]
        [1 0 1 0 0 0 1 0 0 0]
        [0 1 0 1 0 0 0 1 0 0]
        [0 0 1 0 1 0 0 0 1 0]
        [1 0 0 1 0 0 0 0 0 1]
        [1 0 0 0 0 0 0 1 1 0]
        [0 1 0 0 0 0 0 0 1 1]
        [0 0 1 0 0 1 0 0 0 1]
        [0 0 0 1 0 1 1 0 0 0]
        [0 0 0 0 1 0 1 1 0 0]
        Full MatrixSpace of 10 by 10 dense matrices over Integer Ring

    ::

        sage: matrix(ZZ, 10, 10, range(100), sparse=True).parent()
        Full MatrixSpace of 10 by 10 sparse matrices over Integer Ring

    ::

        sage: R = PolynomialRing(QQ, 9, 'x')
        sage: A = matrix(R, 3, 3, R.gens()); A
        [x0 x1 x2]
        [x3 x4 x5]
        [x6 x7 x8]
        sage: det(A)
        -x2*x4*x6 + x1*x5*x6 + x2*x3*x7 - x0*x5*x7 - x1*x3*x8 + x0*x4*x8

    TESTS:

    There are many ways to create an empty matrix::

        sage: m = matrix(); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 dense matrices over Integer Ring
        sage: m = matrix(sparse=True); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 sparse matrices over Integer Ring
        sage: m = matrix(QQ); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 dense matrices over Rational Field
        sage: m = matrix(ring=QQ); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 dense matrices over Rational Field
        sage: m = matrix(0); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 dense matrices over Integer Ring
        sage: m = matrix(0, 0, ring=QQ); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 dense matrices over Rational Field
        sage: m = matrix([]); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 dense matrices over Integer Ring
        sage: m = matrix(QQ, []); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 dense matrices over Rational Field
        sage: m = matrix(QQ, {}); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 sparse matrices over Rational Field

    Only a ring and dimensions::

        sage: m = matrix(2); m; m.parent()
        [0 0]
        [0 0]
        Full MatrixSpace of 2 by 2 dense matrices over Integer Ring
        sage: m = matrix(QQ,2); m; m.parent()
        [0 0]
        [0 0]
        Full MatrixSpace of 2 by 2 dense matrices over Rational Field
        sage: m = matrix(QQ,2,3); m; m.parent()
        [0 0 0]
        [0 0 0]
        Full MatrixSpace of 2 by 3 dense matrices over Rational Field

    A ring, dimensions and a scalar::

        sage: m = matrix(2,2,1); m; m.parent()
        [1 0]
        [0 1]
        Full MatrixSpace of 2 by 2 dense matrices over Integer Ring
        sage: m = matrix(QQ,2,2,5); m; m.parent()
        [5 0]
        [0 5]
        Full MatrixSpace of 2 by 2 dense matrices over Rational Field

    For non-square matrices, only zero works::

        sage: m = matrix(2,3,0); m; m.parent()
        [0 0 0]
        [0 0 0]
        Full MatrixSpace of 2 by 3 dense matrices over Integer Ring
        sage: m = matrix(QQ,2,3,0); m; m.parent()
        [0 0 0]
        [0 0 0]
        Full MatrixSpace of 2 by 3 dense matrices over Rational Field
        sage: matrix(QQ,2,3,1)
        Traceback (most recent call last):
        ...
        TypeError: identity matrix must be square
        sage: matrix(QQ,2,3,5)
        Traceback (most recent call last):
        ...
        TypeError: nonzero scalar matrix must be square

    Matrices specified by a list of entries::

        sage: m = matrix([[1,2,3],[4,5,6]]); m; m.parent()
        [1 2 3]
        [4 5 6]
        Full MatrixSpace of 2 by 3 dense matrices over Integer Ring
        sage: m = matrix(QQ,2,[[1,2,3],[4,5,6]]); m; m.parent()
        [1 2 3]
        [4 5 6]
        Full MatrixSpace of 2 by 3 dense matrices over Rational Field
        sage: m = matrix(QQ,3,[[1,2,3],[4,5,6]]); m; m.parent()
        Traceback (most recent call last):
        ...
        ValueError: number of rows does not match up with specified number
        sage: m = matrix(QQ,2,3,[[1,2,3],[4,5,6]]); m; m.parent()
        [1 2 3]
        [4 5 6]
        Full MatrixSpace of 2 by 3 dense matrices over Rational Field
        sage: m = matrix(QQ,2,4,[[1,2,3],[4,5,6]]); m; m.parent()
        Traceback (most recent call last):
        ...
        ValueError: number of columns does not match up with specified number
        sage: m = matrix([(1,2,3),(4,5,6)]); m; m.parent()
        [1 2 3]
        [4 5 6]
        Full MatrixSpace of 2 by 3 dense matrices over Integer Ring
        sage: m = matrix([1,2,3,4,5,6]); m; m.parent()
        [1 2 3 4 5 6]
        Full MatrixSpace of 1 by 6 dense matrices over Integer Ring
        sage: m = matrix((1,2,3,4,5,6)); m; m.parent()
        [1 2 3 4 5 6]
        Full MatrixSpace of 1 by 6 dense matrices over Integer Ring
        sage: m = matrix(QQ,[1,2,3,4,5,6]); m; m.parent()
        [1 2 3 4 5 6]
        Full MatrixSpace of 1 by 6 dense matrices over Rational Field
        sage: m = matrix(QQ,3,2,[1,2,3,4,5,6]); m; m.parent()
        [1 2]
        [3 4]
        [5 6]
        Full MatrixSpace of 3 by 2 dense matrices over Rational Field
        sage: m = matrix(QQ,2,4,[1,2,3,4,5,6]); m; m.parent()
        Traceback (most recent call last):
        ...
        ValueError: entries has the wrong length
        sage: m = matrix(QQ,5,[1,2,3,4,5,6]); m; m.parent()
        Traceback (most recent call last):
        ...
        TypeError: cannot construct an element of
        Full MatrixSpace of 5 by 1 dense matrices over Rational Field
        from [1, 2, 3, 4, 5, 6]!

    Matrices specified by a dict of entries::

        sage: m = matrix({(1,1): 2}); m; m.parent()
        [0 0]
        [0 2]
        Full MatrixSpace of 2 by 2 sparse matrices over Integer Ring
        sage: m = matrix({(1,1): 2}, sparse=False); m; m.parent()
        [0 0]
        [0 2]
        Full MatrixSpace of 2 by 2 dense matrices over Integer Ring
        sage: m = matrix(QQ,{(1,1): 2}); m; m.parent()
        [0 0]
        [0 2]
        Full MatrixSpace of 2 by 2 sparse matrices over Rational Field
        sage: m = matrix(QQ,3,{(1,1): 2}); m; m.parent()
        [0 0 0]
        [0 2 0]
        [0 0 0]
        Full MatrixSpace of 3 by 3 sparse matrices over Rational Field
        sage: m = matrix(QQ,3,4,{(1,1): 2}); m; m.parent()
        [0 0 0 0]
        [0 2 0 0]
        [0 0 0 0]
        Full MatrixSpace of 3 by 4 sparse matrices over Rational Field
        sage: m = matrix(QQ,2,{(1,1): 2}); m; m.parent()
        [0 0]
        [0 2]
        Full MatrixSpace of 2 by 2 sparse matrices over Rational Field
        sage: m = matrix(QQ,1,{(1,1): 2}); m; m.parent()
        Traceback (most recent call last):
        ...
        IndexError: invalid entries list
        sage: m = matrix({}); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 sparse matrices over Integer Ring
        sage: m = matrix(QQ,{}); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 sparse matrices over Rational Field
        sage: m = matrix(QQ,2,{}); m; m.parent()
        [0 0]
        [0 0]
        Full MatrixSpace of 2 by 2 sparse matrices over Rational Field
        sage: m = matrix(QQ,2,3,{}); m; m.parent()
        [0 0 0]
        [0 0 0]
        Full MatrixSpace of 2 by 3 sparse matrices over Rational Field
        sage: m = matrix(2,{}); m; m.parent()
        [0 0]
        [0 0]
        Full MatrixSpace of 2 by 2 sparse matrices over Integer Ring
        sage: m = matrix(2,3,{}); m; m.parent()
        [0 0 0]
        [0 0 0]
        Full MatrixSpace of 2 by 3 sparse matrices over Integer Ring

    Matrices with zero rows or columns::

        sage: m = matrix(0,2); m; m.parent()
        []
        Full MatrixSpace of 0 by 2 dense matrices over Integer Ring
        sage: m = matrix(2,0); m; m.parent()
        []
        Full MatrixSpace of 2 by 0 dense matrices over Integer Ring
        sage: m = matrix(0,[1]); m; m.parent()
        Traceback (most recent call last):
        ...
        ValueError: entries has the wrong length
        sage: m = matrix(1,0,[]); m; m.parent()
        []
        Full MatrixSpace of 1 by 0 dense matrices over Integer Ring
        sage: m = matrix(0,1,[]); m; m.parent()
        []
        Full MatrixSpace of 0 by 1 dense matrices over Integer Ring
        sage: m = matrix(0,[]); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 dense matrices over Integer Ring
        sage: m = matrix(0,{}); m; m.parent()
        []
        Full MatrixSpace of 0 by 0 sparse matrices over Integer Ring
        sage: m = matrix(0,{(1,1):2}); m; m.parent()
        Traceback (most recent call last):
        ...
        IndexError: invalid entries list
        sage: m = matrix(2,0,{(1,1):2}); m; m.parent()
        Traceback (most recent call last):
        ...
        IndexError: invalid entries list

    Check conversion from numpy::

        sage: import numpy
        sage: n = numpy.array([[numpy.complex(0,1),numpy.complex(0,2)],[3,4]],complex)
        sage: m = matrix(n); m; m.parent()
        [1.0*I 2.0*I]
        [  3.0   4.0]
        Full MatrixSpace of 2 by 2 dense matrices over Complex Double Field
        sage: n = numpy.array([[1,2],[3,4]],'int32')
        sage: m = matrix(n); m; m.parent()
        [1 2]
        [3 4]
        Full MatrixSpace of 2 by 2 dense matrices over Integer Ring
        sage: n = numpy.array([[1,2,3],[4,5,6],[7,8,9]],'float32')
        sage: m = matrix(n); m; m.parent()
        [1.0 2.0 3.0]
        [4.0 5.0 6.0]
        [7.0 8.0 9.0]
        Full MatrixSpace of 3 by 3 dense matrices over Real Double Field
        sage: n = numpy.matrix([[1,2,3],[4,5,6],[7,8,9]],'float64')
        sage: m = matrix(n); m; m.parent()
        [1.0 2.0 3.0]
        [4.0 5.0 6.0]
        [7.0 8.0 9.0]
        Full MatrixSpace of 3 by 3 dense matrices over Real Double Field
        sage: n = numpy.array([[1,2,3],[4,5,6],[7,8,9]],'complex64')
        sage: m = matrix(n); m; m.parent()
        [1.0 2.0 3.0]
        [4.0 5.0 6.0]
        [7.0 8.0 9.0]
        Full MatrixSpace of 3 by 3 dense matrices over Complex Double Field
        sage: n = numpy.matrix([[1,2,3],[4,5,6],[7,8,9]],'complex128')
        sage: m = matrix(n); m; m.parent()
        [1.0 2.0 3.0]
        [4.0 5.0 6.0]
        [7.0 8.0 9.0]
        Full MatrixSpace of 3 by 3 dense matrices over Complex Double Field
        sage: a = matrix([[1,2],[3,4]])
        sage: b = matrix(a.numpy()); b
        [1 2]
        [3 4]
        sage: a == b
        True
        sage: c = matrix(a.numpy('float32')); c
        [1.0 2.0]
        [3.0 4.0]
        sage: matrix(numpy.array([[5]]))
        [5]
        sage: matrix(numpy.matrix([[5]]))
        [5]

    A ring and a numpy array::

        sage: n = numpy.array([[1,2,3],[4,5,6],[7,8,9]],'float32')
        sage: m = matrix(ZZ, n); m; m.parent()
        [1 2 3]
        [4 5 6]
        [7 8 9]
        Full MatrixSpace of 3 by 3 dense matrices over Integer Ring
        sage: n = matrix(QQ, 2, 2, [1, 1/2, 1/3, 1/4]).numpy(); n
        array([[ 1.        ,  0.5       ],
               [ 0.33333333,  0.25      ]])
        sage: matrix(QQ, n)
        [  1 1/2]
        [1/3 1/4]

    The dimensions of a matrix may be given as numpy types::

        sage: matrix(numpy.int32(2), ncols=numpy.int32(3))
        [0 0 0]
        [0 0 0]

    The dimensions of a matrix must have an integral type::

        sage: matrix(RR, 2.0, 2.0)
        Traceback (most recent call last):
        ...
        TypeError: invalid matrix constructor: type matrix? for help

    More tests::

        sage: v = vector(ZZ, [1, 10, 100])
        sage: m = matrix(ZZ['x'], v); m; m.parent()
        [  1  10 100]
        Full MatrixSpace of 1 by 3 dense matrices over Univariate Polynomial Ring in x over Integer Ring
        sage: matrix(ZZ, 10, 10, range(100)).parent()
        Full MatrixSpace of 10 by 10 dense matrices over Integer Ring
        sage: m = matrix(GF(7), [[1/3,2/3,1/2], [3/4,4/5,7]]); m; m.parent()
        [5 3 4]
        [6 5 0]
        Full MatrixSpace of 2 by 3 dense matrices over Finite Field of size 7
        sage: m = matrix([[1,2,3], [RDF(2), CDF(1,2), 3]]); m; m.parent()
        [        1.0         2.0         3.0]
        [        2.0 1.0 + 2.0*I         3.0]
        Full MatrixSpace of 2 by 3 dense matrices over Complex Double Field
        sage: m = matrix(3,3,1/2); m; m.parent()
        [1/2   0   0]
        [  0 1/2   0]
        [  0   0 1/2]
        Full MatrixSpace of 3 by 3 dense matrices over Rational Field
        sage: matrix([[1],[2,3]])
        Traceback (most recent call last):
        ...
        ValueError: list of rows is not valid (rows are wrong types or lengths)
        sage: matrix([[1],2])
        Traceback (most recent call last):
        ...
        ValueError: list of rows is not valid (rows are wrong types or lengths)
        sage: matrix(vector(RR,[1,2,3])).parent()
        Full MatrixSpace of 1 by 3 dense matrices over Real Field with 53 bits of precision

    Check :trac:`10158`::

        sage: matrix(ZZ, [[0] for i in range(10^5)]).is_zero()
        True

    Test conversion using a ``_matrix_`` method::

        sage: A = gap(MatrixSpace(QQ, 2, 2)(range(4)))
        sage: matrix(QQ, A)
        [0 1]
        [2 3]
        sage: matrix(A, ring=QQ)
        [0 1]
        [2 3]

    A redundant ``ring`` argument::

        sage: matrix(ZZ, 3, 3, ring=ZZ)
        Traceback (most recent call last):
        ...
        TypeError: invalid matrix constructor: type matrix? for help

    TESTS:

    Some calls using an iterator (note that xrange is no longer available
    in Python 3)::

        sage: from six.moves import range
        sage: matrix(QQ, 3, 6, range(18), sparse=true)
        [ 0  1  2  3  4  5]
        [ 6  7  8  9 10 11]
        [12 13 14 15 16 17]
        sage: matrix(4, 4, range(16))
        [ 0  1  2  3]
        [ 4  5  6  7]
        [ 8  9 10 11]
        [12 13 14 15]

    AUTHORS:

    - William Stein: Initial implementation

    - Jason Grout (2008-03): almost a complete rewrite, with bits and
      pieces from the original implementation

    - Jeroen Demeyer (2016-02-05): major clean up, see :trac:`20015`
      and :trac:`20016`
    """
    def __call__(self, *Args, ring=None, nrows=None, ncols=None, sparse=None):
        cdef list args = list(Args)

        # ring argument
        if ring is None and args and is_Ring(args[0]):
            ring = args.pop(0)

        # object with _matrix_ method
        if args:
            try:
                makematrix = args[0]._matrix_
            except AttributeError:
                pass
            else:
                if ring is None:
                    args.pop(0)
                else:
                    args[0] = ring
                if sparse is None:
                    return makematrix(*args)
                else:
                    return makematrix(*args, sparse=sparse)

        # nrows argument
        if nrows is None and args:
            arg = args[0]
            try:
                if is_numpy_type(type(arg)):
                    import numpy
                    if isinstance(arg, numpy.ndarray):
                        raise TypeError
                nrows = pyobject_to_long(arg)
            except TypeError:
                pass
            else:
                args.pop(0)

        # ncols argument
        if ncols is None and args:
            arg = args[0]
            try:
                if is_numpy_type(type(arg)):
                    import numpy
                    if isinstance(arg, numpy.ndarray):
                        raise TypeError
                ncols = pyobject_to_long(arg)
            except TypeError:
                pass
            else:
                args.pop(0)

        # Now we've taken care of initial ring, nrows, and ncols arguments.
        # We've also taken care of the Sage object case.

        # Now the rest of the arguments are a list of rows, a flat list of
        # entries, a callable, a dict, a numpy array, or a single value.
        entry_ring = ZZ
        if not args:
            # If no entries are specified, pass back a zero matrix
            entries = 0
        elif len(args) == 1:
            arg = args[0]
            if isinstance(arg, (types.FunctionType, types.LambdaType, types.MethodType)):
                if ncols is None and nrows is None:
                    raise TypeError("when passing in a callable, the dimensions of the matrix must be specified")
                if ncols is None:
                    ncols = nrows
                elif nrows is None:
                    nrows = ncols

                irange = srange(nrows)
                jrange = srange(ncols)
                arg = [[arg(i, j) for j in jrange] for i in irange]

            if isinstance(arg, xrange):
                arg = list(arg)
            if isinstance(arg, (list, tuple)):
                if not arg:
                    # no entries are specified, pass back the zero matrix
                    entries = 0
                elif isinstance(arg[0], (list, tuple)) or isinstance(arg[0], Vector):
                    # Ensure we have a list of lists, each inner list having the same number of elements
                    first_len = len(arg[0])
                    if not all( (isinstance(v, (list, tuple)) or isinstance(v, Vector)) and len(v) == first_len for v in arg):
                        raise ValueError("list of rows is not valid (rows are wrong types or lengths)")
                    # We have a list of rows or vectors
                    if nrows is None:
                        nrows = len(arg)
                    elif nrows != len(arg):
                        raise ValueError("number of rows does not match up with specified number")
                    if ncols is None:
                        ncols = len(arg[0])
                    elif ncols != len(arg[0]):
                        raise ValueError("number of columns does not match up with specified number")

                    entries = []
                    for v in arg:
                        entries.extend(v)

                else:
                    # We have a flat list; figure out nrows and ncols
                    if nrows is None:
                        nrows = 1

                    if nrows > 0:
                        if ncols is None:
                            ncols = len(arg) // nrows
                        elif ncols != len(arg) // nrows:
                            raise ValueError("entries has the wrong length")
                    elif len(arg) > 0:
                        raise ValueError("entries has the wrong length")

                    entries = arg

                if nrows > 0 and ncols > 0 and ring is None:
                    entries, ring = prepare(entries)

            elif isinstance(arg, dict):
                # We have a dictionary: default to sparse
                if sparse is None:
                    sparse = True
                if not arg:
                    # no entries are specified, pass back the zero matrix
                    entries = 0
                else:
                    entries, entry_ring = prepare_dict(arg)
                    if nrows is None:
                        nrows = nrows_from_dict(entries)
                        ncols = ncols_from_dict(entries)
                    # note that ncols can still be None if nrows is set --
                    # it will be assigned nrows down below.

                # See the construction after the numpy case below.
            else:
                if is_numpy_type(type(arg)):
                    import numpy
                    if isinstance(arg, numpy.ndarray):
                        # Convert to a numpy array if it was a matrix.
                        if type(arg) is not numpy.ndarray:
                            arg = numpy.array(arg)

                        str_dtype = str(arg.dtype)

                        if not (arg.flags.c_contiguous is True or arg.flags.f_contiguous is True):
                            raise TypeError('numpy matrix must be either c_contiguous or f_contiguous')

                        if str_dtype.count('float32') == 1:
                            m = matrix(RDF, arg.shape[0], arg.shape[1], 0)
                            m._replace_self_with_numpy32(arg)
                        elif str_dtype.count('float64') == 1:
                            m = matrix(RDF, arg.shape[0], arg.shape[1], 0)
                            m._replace_self_with_numpy(arg)
                        elif str_dtype.count('complex64') == 1:
                            m = matrix(CDF, arg.shape[0], arg.shape[1], 0)
                            m._replace_self_with_numpy32(arg)
                        elif str_dtype.count('complex128') == 1:
                            m = matrix(CDF, arg.shape[0], arg.shape[1], 0)
                            m._replace_self_with_numpy(arg)
                        elif str_dtype.count('int') == 1:
                            m = matrix(ZZ, [list(row) for row in list(arg)])
                        elif str_dtype.count('object') == 1:
                            # Get the raw nested list from the numpy array
                            # and feed it back into matrix
                            m = matrix([list(row) for row in list(arg)])
                        else:
                            raise TypeError("cannot convert NumPy matrix to Sage matrix")

                        if ring is not None and m.base_ring() is not ring:
                            m = m.change_ring(ring)

                        return m
                elif nrows is not None and ncols is not None:
                    # assume that we should just pass the thing into the
                    # MatrixSpace constructor and hope for the best
                    # This is not documented, but it is doctested
                    if ring is None:
                        entry_ring = arg.parent()
                    entries = arg
                else:
                    raise TypeError("invalid matrix constructor: type matrix? for help")
        else:  # len(args) >= 2
            raise TypeError("invalid matrix constructor: type matrix? for help")

        if ring is None:
            ring = entry_ring
        if nrows is None:
            nrows = 0
        if ncols is None:
            ncols = nrows

        return MatrixSpace(ring, nrows, ncols, sparse=sparse)(entries)

Matrix = matrix = MatrixFactory()


def prepare(w):
    """
    Given a list w of numbers, find a common ring that they all
    canonically map to, and return the list of images of the elements
    of w in that ring along with the ring.

    This is for internal use by the matrix function.

    INPUT:

    - ``w`` - list

    OUTPUT:

    list, ring

    EXAMPLES::

        sage: sage.matrix.constructor.prepare([-2, Mod(1,7)])
        ([5, 1], Ring of integers modulo 7)

    Notice that the elements must all canonically coerce to a common
    ring (since Sequence is called)::

        sage: sage.matrix.constructor.prepare([2/1, Mod(1,7)])
        Traceback (most recent call last):
        ...
        TypeError: unable to find a common ring for all elements

    TESTS:

    Check that :trac:`19920` is fixed::

        sage: import numpy
        sage: matrix([[numpy.int8(1)]])
        [1]
    """
    if not w:
        return Sequence([], ZZ), ZZ
    entries = Sequence(w)
    ring = entries.universe()
    if isinstance(ring,type):
        ring = py_scalar_parent(ring)
    if not is_Ring(ring):
        raise TypeError("unable to find a common ring for all elements")
    return entries, ring

def prepare_dict(w):
    """
    Given a dictionary w of numbers, find a common ring that they all
    canonically map to, and return the dictionary of images of the
    elements of w in that ring along with the ring.

    This is for internal use by the matrix function.

    INPUT:

    - ``w`` -- dict

    OUTPUT:

    dict, ring

    EXAMPLES::

        sage: sage.matrix.constructor.prepare_dict({(0,1):2, (4,10):Mod(1,7)})
        ({(0, 1): 2, (4, 10): 1}, Ring of integers modulo 7)
    """
    Z = list(w.items())
    X = [x for _, x in Z]
    entries, ring = prepare(X)
    return {Z[i][0]: ent for i, ent in enumerate(entries)}, ring


def nrows_from_dict(d):
    """
    Given a dictionary that defines a sparse matrix, return the number
    of rows that matrix should have.

    This is for internal use by the matrix function.

    INPUT:

    - ``d`` - dict

    OUTPUT:

    integer

    EXAMPLES::

        sage: sage.matrix.constructor.nrows_from_dict({})
        0

    Here the answer is 301 not 300, since there is a 0-th row.

    ::

        sage: sage.matrix.constructor.nrows_from_dict({(300,4):10})
        301
    """
    if 0 == len(d):
        return 0
    return max([0] + [ij[0] for ij in d.keys()]) + 1

def ncols_from_dict(d):
    """
    Given a dictionary that defines a sparse matrix, return the number
    of columns that matrix should have.

    This is for internal use by the matrix function.

    INPUT:

    - ``d`` - dict

    OUTPUT:

    integer

    EXAMPLES::

        sage: sage.matrix.constructor.ncols_from_dict({})
        0

    Here the answer is 301 not 300, since there is a 0-th row.

    ::

        sage: sage.matrix.constructor.ncols_from_dict({(4,300):10})
        301
    """
    if 0 == len(d):
        return 0
    return max([0] + [ij[1] for ij in d.keys()]) + 1


from .special import *

YouTry

Try recreating the matrix for multiplication, just as we have just recreated the one for addition.

In [31]:
m = 4
listFormMultiply = [] # make the list non-empty!
In [32]:
listFormMultiply
Out[32]:
[]

Modular arithmetic in SageMath

The simplest way to create a number modulo $m$ in Sage is to use the Mod(a,m) command. We illustrate this below.

In [15]:
Mod(8, 12)
Out[15]:
8

Let's assign it to a variable so that we can explore it further:

In [16]:
myModInt = Mod(8, 12)
myModInt
Out[16]:
8
In [17]:
type(myModInt)
Out[17]:
<type 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>
In [18]:
parent(myModInt)
Out[18]:
Ring of integers modulo 12

We will compare myModInt to a "normal" SageMath integer:

In [19]:
myInt = 8
myInt
Out[19]:
8
In [20]:
type(myInt)
Out[20]:
<type 'sage.rings.integer.Integer'>
In [21]:
parent(myInt) # ZZ
Out[21]:
Integer Ring

We can see that myModInt and myInt are different types, but what does this mean? How do they behave?

Try addition:

In [22]:
myModInt + 6
Out[22]:
2
In [23]:
myInt + 6
Out[23]:
14

Was this what you already expected?

What about multiplication?

In [24]:
myModInt * 6
Out[24]:
0
In [25]:
myInt * 6
Out[25]:
48

What's going on here? As we said above, arithmetic modulo mm is like usual arithmetic, except every time you add or multiply, you also divide by mm and return the remainder. 8 x 6 is 48, and the remainder of 48 divided by 12 (12 is our modulo) is 0.

What about this one? What's happening here?

In [26]:
Mod(-45,2018) # Raaz was born in the year
Out[26]:
1973
In [27]:
Mod(-1,12)  # what's an hour before mid-night
Out[27]:
11

YouTry

Can you create the number giving your year of birth (±1 year) in a similar way. For example, if you are 19 years old now then find the number -19 modulo 2018.

In [ ]:
 

Let us assign 10 modulo 12 to a variable named now:

In [49]:
now = Mod(10, 12)
now
Out[49]:
10

<img src="images/Week6NowEquals10.png" width=200>

YouTrys

Add -1 to now: <img src="images/Week6NowEquals9.png" width=200>

Put in the expression to do this in SageMath into the cell below:

In [50]:
now -1
Out[50]:
9

And subtract 13 from the previous expression. <img src="images/Week6NowEquals8.png" width=250>

Put in the expression to do this in SageMath into the cell below:

In [51]:
now - 1 -13
Out[51]:
8

Also try adding 14 to the previous expression -- the new postition is not shown on this clock but you should be able to see what it should be.

<img src="images/Week6NowWas8.png" width=200>

Put in the expression to do this in SageMath into the cell below:

In [ ]:
 

And multiplying 2 by now (or now by 2)

In [ ]:
 

Try making and using some modular integers of your own. Make sure that you can predict the results of simple addition and multiplication operations on them: this will confirm that you understand what's happening.

In [ ]:
 
In [ ]:
 

YouTry Later!

What happens if you try arithmetic with two modular integers? Does it matter if the moduli of both operands a and b in say a+b are the same if a and b are both modular integers? Does this make sense to you?

In [ ]:
 
In [ ]:
 
In [ ]:
 

(end of YouTrys)



Linear Congruential Generators

A linear congruential generator (LCG) is a simple pseudo-random number generator (PRNG) - a simple way of imitating samples or realizations from the Uniform$(0,1)$.

"Pseudo-random" means that the numbers are not really random [From Ancient Greek ψευδής (pseudḗs, “false, lying”)].

We will look at what we mean by that as we find out about linear congruential generators.

The theory behind LCGs is easy to understand, and they are easily implemented and fast.

To make a LCG we need:

  • a modulus $m$ ($m > 0$)
  • an integer multiplier $a$ ($0 \le a < m$)
  • an integer increment $c$ ($0 \le c < m$)
  • an integer seed $x_0$ ($0 \le x_0 < m$)
  • an integer sequence length $n$

Using these inputs, the LCG generates numbers $x_1, x_2, \ldots x_{n-1} $ where $x_{i+1}$ is calculated from $x_i$ as defined by the following recurrence relation:

$$x_{i+1} \gets mod \left(a x_i + c , m \right)$$

In the above expression $\gets$ denotes "gets" as in variable assignment.

$x_0,x_1,\ldots,x_{n-1}$ is the sequence of pseudo-random numbers called the linear congruential sequence.

We can define a function parameterised by $(m,a,c,x_0,n)$ to give us a linear congruential sequence in the form of a list.

(Remember about function parameters? The function parameters here are m, a, c, x0, and n. Note that the return value of the function is a list.

Implementing the LCG

Let's see how we can define a function for our LCG next.'

In [28]:
def linConGen(m, a, c, x0, n):
    '''A linear congruential sequence generator.
    
    Param m is the integer modulus to use in the generator.
    Param a is the integer multiplier.
    Param c is the integer increment.
    Param x0 is the integer seed.
    Param n is the integer number of desired pseudo-random numbers.
    
    Returns a list of n pseudo-random integer modulo m numbers.'''
    
    x = x0 # the seed
    retValue = [Mod(x,m)]  # start the list with x=x0
    for i in range(2, n+1, 1):
        x = mod(a * x + c, m) # the generator, using modular arithmetic
        retValue.append(x) # append the new x to the list
    return retValue

You know enough SageMath/Python to understand what every line in the above function linConGen is doing!

The function is merely implementing the pseudocode or algorithm of the linear congruential generator using a for-loop and modular arithmetic: note that the generator produces integers modulo $m$.

Linear Congruential Generators: the Good, the Bad, and the Ugly

Are all linear congruential generators as good as each other? What makes a good LCG?

One desirable property of a LCG is to have the longest possible period.

The period is the length of sequence we can get before we get a repeat. The longest possible period a LCG can have is $m$. Lets look at an example.

In [29]:
?mod
In [30]:
# first assign values to some variables to pass as arguments to the function
m  = 17  # set modulus to 17
a  = 2   # set multiplier to 2
c  = 7   # set increment to 7
x0 = 1   # set seed to 1
n  = 18  # set length of sequence to 18 = 1 + maximal period 
L1 = linConGen(m, a, c, x0, n) # use our linConGren function to make the sequence
L1                    # this sequence repeats itself with period 8
Out[30]:
[1, 9, 8, 6, 2, 11, 12, 14, 1, 9, 8, 6, 2, 11, 12, 14, 1, 9]

You should be able to see the repeating pattern 9, 8, 6, 2, 11, 12, 14. If you can't you can see that the sequence actually contains 8 unique numbers by making a set out of it:

In [32]:
Set(L1)                        # make a Sage Set out of the sequence in the list L1
Out[32]:
{1, 2, 6, 8, 9, 11, 12, 14}

Changing the seed $x_0$ will, at best, make the sequence repeat itself over other numbers but with the same period:

In [33]:
x0 = 3             # leave m, a, c as it is but just change the seed from 1 to 3 here
L2 = linConGen(m, a, c, x0, n)
L2                 # changing the seed makes this sequence repeat itself over other numbers but also with period 8
Out[33]:
[3, 13, 16, 5, 0, 7, 4, 15, 3, 13, 16, 5, 0, 7, 4, 15, 3, 13]
In [35]:
Set(L2)                        # make a Sage Set out of the sequence in the list L2
Out[35]:
{0, 3, 4, 5, 7, 13, 15, 16}

At worst, a different seed specified by x0 might make the sequence get stuck immediately:

In [36]:
x0 = 10            # leave m, a, c as it is but just change the seed to 10
L3 = linConGen(m, a, c, x0, n)
L3
Out[36]:
[10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10]
In [37]:
Set(L3)                        # make a Sage Set out of the sequence in the list L3
Out[37]:
{10}

An LCG is just a discrete determinstic dynamical system and each seed acts as the initial condition with different behaviours. The behaviour above can be seen as a fixed point at 10 for this LCG initialized at 10.

In [38]:
Set(L3)+Set(L2)+Set(L1) # + on Sage Sets gives the union of the three Sets and it is the set of integers modulo 17
Out[38]:
{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}

Thus, the three sequences with different initial conditions or seeds cover all the points in our Ring of integers modulo 17.

What about changing the multiplier $a$?

In [39]:
m  = 17 # same as before
a  = 3  # change the multiplier from 2 to 3 here
c  = 7  # same as before
x0 = 1  # set seed at 1 
n  = 18 # set length of sequence to 18 = 1 + maximal period
L4 = linConGen(m, a, c, x0, n)
L4                             # now we get a longer period of 16 but with the number 5 missing
Out[39]:
[1, 10, 3, 16, 4, 2, 13, 12, 9, 0, 7, 11, 6, 8, 14, 15, 1, 10]
In [41]:
Set(L4)
Out[41]:
{0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}
In [42]:
x0 = 5  # just change the seed to 5
linConGen(m, a, c, x0, n)
Out[42]:
[5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5]

We want an LCG with a full period of $m$ so that we can use it with any seed and not get stuck at fixed points or short periodic sequences. This is a minimal requirement for simulation purposes that we will employ such sequences for. Is there anyway to know what makes a LCG with a full period?

It turns out that an LCG will have a full period if and only if:

  • $c$ and $m$ are relatively prime or coprime. i.e. the greatest common divisor (gcd) of $c$ and $m$ is $1$; and
  • $a-1$ is divisible by all prime factors of $m$; and
  • $a-1$ is a multiple of 4 if $m$ is a multiple of 4.

(Different conditions apply when $c=0$. The Proposition and Proof for this are in Knuth, The Art of Computer Programming, vol. 2, section 3.3).

SageMath has a function gcd which can calculate the greatest common divisor of two numbers:

In [43]:
?gcd
In [46]:
gcd(7,17)
Out[46]:
1

SageMath can also help us to calculate prime factors with the prime_factors function:

In [47]:
prime_factors(m)
Out[47]:
[17]
In [48]:
prime_factors(7634887623876428376482746)
Out[48]:
[2, 25504980841, 149674443424853]

$(m, a, c, x0, n) = (256, 137, 123, 13, 256)$ gives a linear congruential sequence with the longest possible period of 256. Let us see how these parameters satisfy the above three requirements while those earlier with $m=17$, $a=2$, and $c=7$ do not.

In [49]:
(m,a,c,x0,n)=(256, 137, 123, 13, 256) # faster way to assign a bunch of parameters
In [50]:
gcd(c,m)         # checking if the greatest common divisor of c and m is indeed 1
Out[50]:
1
In [51]:
prime_factors(m)   # it is easy to get a list of all the prime factors of m
Out[51]:
[2]
In [52]:
(a-1) % 2             # checking if a-1=136 is divisible by 2, the only prime factors of m
Out[52]:
0
In [53]:
[ (a-1)%x for x in prime_factors(m) ]   # a list comprehension check for an m with more prime factors
Out[53]:
[0]
In [54]:
m % 4         # m is a multiple of 4 check
Out[54]:
0
In [55]:
(a-1) % 4       # if m is a multiple of 4 then a-1 is also a multiple of 4
Out[55]:
0
In [56]:
(m, a, c, x0, n)  # therefore these parameter values satisfy all conditions to have maximum period length m
Out[56]:
(256, 137, 123, 13, 256)

Thus, the parameters $(m, a, c, x_0, n) = (256, 137, 123, 13, 256)$ do indeed satisfy the three conditions to guarantee the longest possible period of 256.

In contrast, for the LCG example earlier with $m=17$, $a=2$, and $c=7$, although $m$ and $c$ are relatively prime, i.e., $gcd(m,c)=1$, we have the violation that $a-1=2-1=1$ is not divisible by the only prime factor $17$ of $m=17$. Thus, we cannot get a period of maximal length $17$ in that example.

In [57]:
[ (2-1)%x for x in prime_factors(17) ]
Out[57]:
[1]

Let us see if the parameters $(m, a, c, x_0, n) = (256, 137, 123, 13, 256)$ that satisfy the three conditions to guarantee the longest possible period of 256 do indeed produce such a sequence:

In [58]:
(m,a,c,x0,n)=(256, 137, 123, 13, 256) # faster way to assign a bunch of parameters
ourLcg = linConGen(m,a,c,x0,n)
print(ourLcg)
[13, 112, 107, 190, 41, 108, 71, 122, 197, 232, 163, 182, 225, 228, 127, 114, 125, 96, 219, 174, 153, 92, 183, 106, 53, 216, 19, 166, 81, 212, 239, 98, 237, 80, 75, 158, 9, 76, 39, 90, 165, 200, 131, 150, 193, 196, 95, 82, 93, 64, 187, 142, 121, 60, 151, 74, 21, 184, 243, 134, 49, 180, 207, 66, 205, 48, 43, 126, 233, 44, 7, 58, 133, 168, 99, 118, 161, 164, 63, 50, 61, 32, 155, 110, 89, 28, 119, 42, 245, 152, 211, 102, 17, 148, 175, 34, 173, 16, 11, 94, 201, 12, 231, 26, 101, 136, 67, 86, 129, 132, 31, 18, 29, 0, 123, 78, 57, 252, 87, 10, 213, 120, 179, 70, 241, 116, 143, 2, 141, 240, 235, 62, 169, 236, 199, 250, 69, 104, 35, 54, 97, 100, 255, 242, 253, 224, 91, 46, 25, 220, 55, 234, 181, 88, 147, 38, 209, 84, 111, 226, 109, 208, 203, 30, 137, 204, 167, 218, 37, 72, 3, 22, 65, 68, 223, 210, 221, 192, 59, 14, 249, 188, 23, 202, 149, 56, 115, 6, 177, 52, 79, 194, 77, 176, 171, 254, 105, 172, 135, 186, 5, 40, 227, 246, 33, 36, 191, 178, 189, 160, 27, 238, 217, 156, 247, 170, 117, 24, 83, 230, 145, 20, 47, 162, 45, 144, 139, 222, 73, 140, 103, 154, 229, 8, 195, 214, 1, 4, 159, 146, 157, 128, 251, 206, 185, 124, 215, 138, 85, 248, 51, 198, 113, 244, 15, 130]
In [59]:
S = Set(ourLcg)  # sort it in a set to see if it indeed has maximal period of 256
S
Out[59]:
{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255}

ourLcg is an example of a good LCG.

The next example will demonstrate what can go wrong.

Consider $m= 256$, $a = 136$, $c = 3$, $x_0 = 0$, $n = 15$.

In [60]:
m,a,c,x0,n = 256, 136, 3, 0, 15
gcd(m,c)
Out[60]:
1
In [61]:
prime_factors(m)
Out[61]:
[2]

But, since $a-1=$135 is not divisible by 2, the only prime factor of $m=$256, we get into the fixed point 91 no matter where we start from.

YouTry

See if changing the seed $x_0$ makes a difference?

In [62]:
linConGen(m, a, c, x0, n)
Out[62]:
[0, 3, 155, 91, 91, 91, 91, 91, 91, 91, 91, 91, 91, 91, 91]
In [ ]:
 

We can look at the linear congruential sequence generated by $(m,a,c,x_0,n)=(256,137,0,123,256)$ from Knuth's classic [The Art of Computer Programming, vol. 2, Sec. 3.3.4, Table 1, Line 5] and compare it with ourLcg.

In [63]:
m, a, c, x0, n = 256, 137, 0, 123, 256
lcgKnuth334T1L5 = linConGen(m, a, c, x0, n)
print(lcgKnuth334T1L5)
[123, 211, 235, 195, 91, 179, 203, 163, 59, 147, 171, 131, 27, 115, 139, 99, 251, 83, 107, 67, 219, 51, 75, 35, 187, 19, 43, 3, 155, 243, 11, 227, 123, 211, 235, 195, 91, 179, 203, 163, 59, 147, 171, 131, 27, 115, 139, 99, 251, 83, 107, 67, 219, 51, 75, 35, 187, 19, 43, 3, 155, 243, 11, 227, 123, 211, 235, 195, 91, 179, 203, 163, 59, 147, 171, 131, 27, 115, 139, 99, 251, 83, 107, 67, 219, 51, 75, 35, 187, 19, 43, 3, 155, 243, 11, 227, 123, 211, 235, 195, 91, 179, 203, 163, 59, 147, 171, 131, 27, 115, 139, 99, 251, 83, 107, 67, 219, 51, 75, 35, 187, 19, 43, 3, 155, 243, 11, 227, 123, 211, 235, 195, 91, 179, 203, 163, 59, 147, 171, 131, 27, 115, 139, 99, 251, 83, 107, 67, 219, 51, 75, 35, 187, 19, 43, 3, 155, 243, 11, 227, 123, 211, 235, 195, 91, 179, 203, 163, 59, 147, 171, 131, 27, 115, 139, 99, 251, 83, 107, 67, 219, 51, 75, 35, 187, 19, 43, 3, 155, 243, 11, 227, 123, 211, 235, 195, 91, 179, 203, 163, 59, 147, 171, 131, 27, 115, 139, 99, 251, 83, 107, 67, 219, 51, 75, 35, 187, 19, 43, 3, 155, 243, 11, 227, 123, 211, 235, 195, 91, 179, 203, 163, 59, 147, 171, 131, 27, 115, 139, 99, 251, 83, 107, 67, 219, 51, 75, 35, 187, 19, 43, 3, 155, 243, 11, 227]
In [64]:
Set(lcgKnuth334T1L5)
Out[64]:
{51, 131, 139, 115, 147, 171, 27, 35, 43, 99, 91, 179, 155, 59, 19, 123, 83, 195, 211, 203, 11, 75, 163, 219, 3, 107, 67, 227, 235, 243, 187, 251}

Note that although ourLcg has maximal period of 256, lcgKnuth334T1L5 has period of 32. Let us look at them as points.

We can plot each number in the sequence, say ourLcg, against its index in the sequence -- i.e. plot the first number 13 against 0 as the tuple (0, 13), the second number against 1 as the tuple (1, 112), etc.

To do this, we need to make a list of the index values, which is simple using range(256) which as you know will give us a list of numbers from 0 to 255 going up in steps of 1. Then we can zip this list with the sequence itself to make the list of the desired tuples.

In [65]:
ourPointsToPlot = zip(range(256), ourLcg)
knuPointsToPlot = zip(range(256), lcgKnuth334T1L5)

Then we plot using points for ourLcg and for lcgKnuth334T1L5.

In [66]:
?text
In [67]:
p1 = points(ourPointsToPlot, pointsize='1')
t1 = text('ourLcg', (150,300), rgbcolor='blue',fontsize=10) 
p2 = points(knuPointsToPlot, pointsize='1',color='black')
t2 = text(' lcgKnuth334T1L5', (150,300), rgbcolor='black',fontsize=10) 
show(graphics_array((p1+t1,p2+t2)),figsize=[6,3])

We can see that in section 3.3.4, Table 1, line 5 in The Art of Computer Programming, Knuth is giving an example of a particularly bad LCG.

When we introducted LCGs, we said that using an LCG was a simple way to imiate independent samples from the Uniform$(0, 1)$ RV, but clearly so far we have been generating sequences of integers. How does that help?

To get a simple pseudo-random Uniform$(0,1)$ generator, we scale the linear congruential sequence over [0, 1]. We can do this by dividing each element by the largest number in the sequence (256 in the case of ourLcg).

Important note: The numbers in the list returned by our linConGen function are integers modulo $m$.

In [68]:
type(ourLcg[0])
Out[68]:
<type 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>

You will need to convert these to Sage integers or Sage multi-precision floating point numbers using int() or RR() command before dividing by $m$. Otherwise you will be doing modular arithmetic when you really want the usual division.

In [69]:
type(RR(ourLcg[0]))
Out[69]:
<type 'sage.rings.real_mpfr.RealNumber'>

Having sorted that out, we can use a list comprehension as a nice neat way to do our scaling:

In [70]:
#ourLcgScaled = [RR(x)/256 for x in ourLcg]    # convert to mpfr real numbers before division
ourLcgScaled = [(RR(x)/256).n(digits=8) for x in ourLcg]    # convert to mpfr real numbers with non-zero digits
#ourLcgScaled = [int(x)/256 for x in ourLcg]  # or convert x to usual integer before division
#ourLcgScaled = [x/256 for x in ourLcg]       # a very bad idea
print(ourLcgScaled)
[0.050781250, 0.43750000, 0.41796875, 0.74218750, 0.16015625, 0.42187500, 0.27734375, 0.47656250, 0.76953125, 0.90625000, 0.63671875, 0.71093750, 0.87890625, 0.89062500, 0.49609375, 0.44531250, 0.48828125, 0.37500000, 0.85546875, 0.67968750, 0.59765625, 0.35937500, 0.71484375, 0.41406250, 0.20703125, 0.84375000, 0.074218750, 0.64843750, 0.31640625, 0.82812500, 0.93359375, 0.38281250, 0.92578125, 0.31250000, 0.29296875, 0.61718750, 0.035156250, 0.29687500, 0.15234375, 0.35156250, 0.64453125, 0.78125000, 0.51171875, 0.58593750, 0.75390625, 0.76562500, 0.37109375, 0.32031250, 0.36328125, 0.25000000, 0.73046875, 0.55468750, 0.47265625, 0.23437500, 0.58984375, 0.28906250, 0.082031250, 0.71875000, 0.94921875, 0.52343750, 0.19140625, 0.70312500, 0.80859375, 0.25781250, 0.80078125, 0.18750000, 0.16796875, 0.49218750, 0.91015625, 0.17187500, 0.027343750, 0.22656250, 0.51953125, 0.65625000, 0.38671875, 0.46093750, 0.62890625, 0.64062500, 0.24609375, 0.19531250, 0.23828125, 0.12500000, 0.60546875, 0.42968750, 0.34765625, 0.10937500, 0.46484375, 0.16406250, 0.95703125, 0.59375000, 0.82421875, 0.39843750, 0.066406250, 0.57812500, 0.68359375, 0.13281250, 0.67578125, 0.062500000, 0.042968750, 0.36718750, 0.78515625, 0.046875000, 0.90234375, 0.10156250, 0.39453125, 0.53125000, 0.26171875, 0.33593750, 0.50390625, 0.51562500, 0.12109375, 0.070312500, 0.11328125, 0.00000000, 0.48046875, 0.30468750, 0.22265625, 0.98437500, 0.33984375, 0.039062500, 0.83203125, 0.46875000, 0.69921875, 0.27343750, 0.94140625, 0.45312500, 0.55859375, 0.0078125000, 0.55078125, 0.93750000, 0.91796875, 0.24218750, 0.66015625, 0.92187500, 0.77734375, 0.97656250, 0.26953125, 0.40625000, 0.13671875, 0.21093750, 0.37890625, 0.39062500, 0.99609375, 0.94531250, 0.98828125, 0.87500000, 0.35546875, 0.17968750, 0.097656250, 0.85937500, 0.21484375, 0.91406250, 0.70703125, 0.34375000, 0.57421875, 0.14843750, 0.81640625, 0.32812500, 0.43359375, 0.88281250, 0.42578125, 0.81250000, 0.79296875, 0.11718750, 0.53515625, 0.79687500, 0.65234375, 0.85156250, 0.14453125, 0.28125000, 0.011718750, 0.085937500, 0.25390625, 0.26562500, 0.87109375, 0.82031250, 0.86328125, 0.75000000, 0.23046875, 0.054687500, 0.97265625, 0.73437500, 0.089843750, 0.78906250, 0.58203125, 0.21875000, 0.44921875, 0.023437500, 0.69140625, 0.20312500, 0.30859375, 0.75781250, 0.30078125, 0.68750000, 0.66796875, 0.99218750, 0.41015625, 0.67187500, 0.52734375, 0.72656250, 0.019531250, 0.15625000, 0.88671875, 0.96093750, 0.12890625, 0.14062500, 0.74609375, 0.69531250, 0.73828125, 0.62500000, 0.10546875, 0.92968750, 0.84765625, 0.60937500, 0.96484375, 0.66406250, 0.45703125, 0.093750000, 0.32421875, 0.89843750, 0.56640625, 0.078125000, 0.18359375, 0.63281250, 0.17578125, 0.56250000, 0.54296875, 0.86718750, 0.28515625, 0.54687500, 0.40234375, 0.60156250, 0.89453125, 0.031250000, 0.76171875, 0.83593750, 0.0039062500, 0.015625000, 0.62109375, 0.57031250, 0.61328125, 0.50000000, 0.98046875, 0.80468750, 0.72265625, 0.48437500, 0.83984375, 0.53906250, 0.33203125, 0.96875000, 0.19921875, 0.77343750, 0.44140625, 0.95312500, 0.058593750, 0.50781250]

This is more like it! We could have a look at this on a plot. Again again want tuples (index, element in scaled sequence at index position), which we can get using range(256) (to get the indexes 0, .., 255) and zip:

In [71]:
ourPointsToPlotScaled = zip(range(256), ourLcgScaled)
p = points(ourPointsToPlotScaled, pointsize='1')
show(p, figsize = (3,3))

Now we have points on the real line. We could use a histogram to look at their distribution. If we are hoping that our LCG, once we have scaled the results, is imitating the Uniform$(0,1)$, what kind of shape would we want our histogram to be?

In [72]:
import pylab # This is a nice Python numerics module on top of mumpy
pylab.clf() # clear current figure
n, bins, patches = pylab.hist(ourLcgScaled, 40) # make the histogram (don't have to have n, bins, patches = ...)
pylab.xlabel('pseudo-random number') # use pyplot methods to set labels, titles etc similar to as in matlab
pylab.ylabel('count')
pylab.title('Count histogram for linear congruential prns')
pylab.savefig('myHist', dpi=(50)) # save figure (dpi) to display the figure
pylab.show() # and finally show it

Is this roughly what you expected?

(Please note: we are using the histogram plots to help you to see the data, but you are not expected to be able to make one yourself.)

We could repeat this for the Knuth bad LCG example:

In [73]:
%%sh
ls
00.html
00.ipynb
00.md
01.html
01.ipynb
01.md
02.html
02.ipynb
02.md
03.html
03.ipynb
03.md
04.html
04.ipynb
04.md
05.html
05.ipynb
05.md
06.html
06.ipynb
06.md
07.html
07.ipynb
07.md
08.html
08.ipynb
08.md
09.html
09.ipynb
09.md
10.html
10.ipynb
10.md
11.html
11.ipynb
11.md
12.html
12.ipynb
12.md
data
images
mydir
myHist.png
pdf
sageMathIpynbArchive
In [74]:
knuLcgScaled = [RR(x)/256 for x in lcgKnuth334T1L5] 
knuPointsToPlotScaled = zip(range(256), knuLcgScaled)
p = points(knuPointsToPlotScaled, pointsize='1', color='black')
show(p, figsize = (3,3))

And show it as a histogram. Given the pattern above, what would you expect the histogram to look like?

In [75]:
import pylab
pylab.clf() # clear current figure
n, bins, patches = pylab.hist(knuLcgScaled, 40) # make the histogram (don't have to have n, bins, patches = ...)
pylab.xlabel('pseudo-random number [Knuth 3.3.4 Table 1, Line 5]') # use pyplot methods to set labels, titles etc similar to as in matlab
pylab.ylabel('count')
pylab.title('Count histogram for linear congruential prns')
pylab.savefig('myHist', dpi=(50)) # save figure (dpi) to display the figure
pylab.show() # and finally show it

Larger LCGs

The above generators are cute but not useful for simulating Lotto draws with 40 outcomes. Minimally, we need to increase the period length with a larger modulus $m$.

But remember that the quality of the pseudo-random numbers obtained from a LCG is extremely sensitive to the choice of $m$, $a$, and $c$.

To illustrate that having a large $m$ alone is not enough we next look at RANDU, an infamous LCG which generates sequences with strong correlations between 3 consecutive points, which can been seen if we manipulate the sequence to make 3-dimensional tuples out of groups of 3 consecutive points.

In the cell below we make our scaled sequence in one step, using a list comprehension which contains the expression to generate the LCG and the scaling.

In [76]:
m, a, c, x0, n = 2147483648, 65539, 0, 1, 5010
RANDU = [RR(x)/m for x in linConGen(m, a, c, x0, n)]

Have a look at the results as a histogram:

In [77]:
import pylab
pylab.clf() # clear current figure
n, bins, patches = pylab.hist(RANDU, 40) # make the histogram (don't have to have n, bins, patches = ...)
pylab.xlabel('pseudo-random number of RANDU') # use pyplot methods to set labels, titles etc similar to as in matlab
pylab.ylabel('count')
pylab.title('Count histogram for linear congruential prns')
pylab.savefig('myHist', dpi=(50)) # save figure (dpi) to display the figure
pylab.show() # and finally show it

Now we are going to use some of the array techniques we have learned about to resize the sequence from the RANDU LCG to an array with two columns. We can then zip the two columns together to make tuples (just pairs or two-tuples).

In [78]:
import pylab
m, a, c, x0, n = 2147483648, 65539, 0, 1, 5010
randu = pylab.array(linConGen(m, a, c, x0, n))
In [79]:
pylab.shape(randu)
Out[79]:
(5010,)
In [80]:
randu.resize(5010/2, 2) # resize the randu array to 2 columns
In [81]:
pylab.shape(randu)
Out[81]:
(2505, 2)
In [83]:
seqs = zip(randu[:, 0], randu[:, 1]) # zip to make tuples from columns
p = points(seqs, pointsize='1', color='black')
show(p, figsize = (3,3))

Let us resize the LCG to an array with three columns. We can then zip the three columns together to make tuples. The effect will be that if our original sequence was $x_0, x_1, x_2, x_3, \ldots, x_{n-3}, x_{n-2}, x_{n-1}$, we will get a list of triplets, tuples of length three, $(x_0, x_3, x_6), (x_1, x_4, x_7), \ldots, (x_{n-1-6}, x_{n-1-3}, x_{n-1})$. Unlike the pairs in 2D which seem well-scattered and random, triplets from the RANDU LCG are not very random at all! They all lie on parallel planes in 3D.

In [84]:
import pylab
m, a, c, x0, n = 2147483648, 65539, 0, 1, 5010
randu = pylab.array(linConGen(m, a, c, x0, n))
randu.resize(5010/3, 3) # resize the randu array to 3 columns
seqs = zip(randu[:, 0], randu[:, 1], randu[:, 2]) # zip to make tuples from columns
point3d(seqs, size=3)
Out[84]:

You can alter your perspective on this image using the mouse. From a particular perspective you can see that something has gone horribly wrong ... RANDU is a really ugly LCG.

The above generators are of low quality for producing pseudo-random numbers to drive statistical simulations. We end with a positive note with a LCG that is in use in the Gnu Compiler Collection. It does not have obvious problems as in small periods or as high a correlation as RANDU.

In [85]:
#jmol error JmolInitCheck is not defined
import pylab
glibCGCCLcg = pylab.array([RR(x)/(2^32) for x in linConGen(2^32, 1103515245,12345,13,5010)])
glibCGCCLcg.resize(1670, 3) # resize the randu array to 3 columns
seqs = zip(glibCGCCLcg[:, 0], glibCGCCLcg[:, 1], glibCGCCLcg[:, 2]) # zip to make tuples from columns
point3d(seqs, size=3)
Out[85]:
In [107]:
import pylab
glibCGCCLcg = pylab.array([RR(x)/(2^32) for x in linConGen(2^32, 1103515245,12345,13,5010)])
pylab.clf() # clear current figure
n, bins, patches = pylab.hist(glibCGCCLcg, 40) # make the histogram (don't have to have n, bins, patches = ...)
pylab.xlabel('pseudo-random number of glibc used by Gnu Compiler Collection') # use pyplot methods to set labels, titles etc similar to as in matlab
pylab.ylabel('count')
pylab.title('Count histogram for linear congruential prns')
pylab.savefig('myHist',dpi=(50)) # seem to need to have this to be able to actually display the figure
pylab.show() # and finally show it

Even good LCG are not suited for realistic statistical simulation problems.

This is because of the strong correlation between successive numbers in the sequence. For instance, if an LCG is used to choose points in an n-dimensional space, the points will lie on, at most, $m^{1/n}$ hyper-planes. There are various statistical tests that one can use to test the quality of a pseudo-random number generator. For example, the spectral test checks if the points are not on a few hyper-planes. Of course, the Sample Mean, Sample Variance, etc. should be as expected. Let us check those quickly:

Recall that the population mean for a Uniform$(0,1)$ RV is $\frac{1}{2}$ and the population variance is $\frac{1}{12}$.

In [86]:
glibCGCCLcg.mean()  # check that the mean is close to the population mean of 0.5 for Uniform(0,1) RV
Out[86]:
0.49233286396638415
In [87]:
glibCGCCLcg.var()   # how about the variance
Out[87]:
0.082447854409686105
In [88]:
1/12.0
Out[88]:
0.0833333333333333

To go into this topic in detail is clearly beyond the scope of this course. You should just remember that using computers to generate pseudo-random numbers is not a trivial problem and use care when employing them especially in higher dimensional or less smooth problems. The mathematics behind this has one of the most beautiful roots of the subject called number theory.

More Sophisticated Pseudo-Random Number Generators

We will use a pseudo-random number generator (PRNG) called the Mersenne Twister for simulation purposes in this course. It is based on more sophisticated theory than that of LCG but the basic principles of recurrence relations are the same.

(The Mersenne Twister is a variant of the recursive relation known as a twisted generalised feedback register. See [Makato Matsumoto and Takuji Nishimura, "Mersenne Twister: A 623-dimensionally equidistributed uniform pseudo-random number generator, ACM Transactions on Modelling and Computer Simulation, vol. 8, no. 1, Jan. 1998, pp. 3-20.], or at http://en.wikipedia.org/wiki/Mersenne_twister.)

The Mersenne Twister has a period of $2^{19937}-1 \approx 10^{6000}$ (which is essentially a Very Big Number) and is currently widely used by researchers interested in statistical simulation.

In [89]:
showURL("http://en.wikipedia.org/wiki/Mersenne_twister",500)
Out[89]:
In [90]:
?random # find out the source file for this function
In [92]:
%%sh
# you can just cat or concatenate the source file to see what is really under the hoood! Power of SageMath/Python!!
#cat ~/all/software/sage/SageMath/local/lib/python2.7/site-packages/sage/misc/prandom.py

Some Warmup before Simulating a Drunkard's Walk

Accumulating sequences with numpy.cumsum

Our example using prngs is going to use a very useful function from the numpy module called cumsum. This calculates a cumulative sum from a sequence of numeric values. What do we mean by a cumulative sum? If we have a sequence $x_0, x_1, x_2, x_3, \ldots, x_n$, then we can derive a sequence that is the cumulative sum,

$$x_0, \,x_0 + x_1,\, x_0+ x_1+x_2, \,\dots, \sum_{i=0}^nx_i$$

Try evaluating the next cell to get the cumulative sum of the sequence 0, 1, 2, ..., 9

In [93]:
from numpy import cumsum, array # make sure we have the numpy classes and methods we need
cumsum(range(10))
Out[93]:
array([ 0,  1,  3,  6, 10, 15, 21, 28, 36, 45])

Since pylab is a superset of numpy, we can also do the following imports instead:

In [94]:
from pylab import cumsum, array # make sure we have the pylab classes and methods we need
cumsum(range(10))
Out[94]:
array([ 0,  1,  3,  6, 10, 15, 21, 28, 36, 45])

You will see that the result is in a pylab.array or numpy.array. This can be useful, but all we need is a list, we can convert this back to a list with the array's tolist method:

In [95]:
cumsum(range(10)).tolist()
Out[95]:
[0, 1, 3, 6, 10, 15, 21, 28, 36, 45]

Using the list function to make a list will do the same thing:

In [96]:
list(cumsum(range(10)))
Out[96]:
[0, 1, 3, 6, 10, 15, 21, 28, 36, 45]

Now, say we have some probabilities in a list and we want to get the cumulative sum of them.

We can use cumsum as we did above. We are going to start with just two probabilties which (for reasons that will shortly become clear), we will assign as values to variables pLeft and pRight.

In [97]:
pLeft, pRight = 0.25, 0.25 # assign values to variables
cumulativeProbs = cumsum((pLeft, pRight)).tolist() # cumsum works on a numeric tuple
cumulativeProbs
Out[97]:
[0.25, 0.5]

Note that cumsum works on a tuple as well as a list, providing that the elements of the tuple are some type of number.

Simulating a Drunkard's Walk

You are now going to use some simulated Uniform$(0,1)$ samples to simulate a Drunkard's Walk.

The idea of the Drunkard's Walk is that the Drunkard has no idea where s/he is going: at each decision point s/he makes a random decision about which way to go. We simulate this random decison using our Uniform$(0,1)$ pseudo-random samples.

We are going to have quite a limited version of this. At each decision point the Drunkard can either go one unit left, right, up or down (pretend we are in the middle of Manhattan with its streets and avenues).

Effectively, s/he is moving around on a $(x, y)$ coordinate grid. The points on the grid will be tuples. Each tuple will have two elements and will represent a point in 2-d space. For example, $(0,0)$ is a tuple which we could use as a starting point.

First, we recall useful feature of tuples: we can unpack a tuple and assign its elements as values to specified variables:

In [98]:
some_point = (3,2) # make a tuple
some_point_x, some_point_y = some_point # unpack the tuple and assign values to variables
some_point_x # disclose one of the variables
Out[98]:
3
In [99]:
some_point_y
Out[99]:
2

We use this useful feature in the functions we have defined below.

You now know what is happening in each step of these functions and should be able to understand what is happening.

In [101]:
def makeLeftTurnPoint(listOfPoints):
    '''Function to make a point representing the destination of a left turn from the current path.
    
    Param listOfPoints is a list of tuples representing the path so far.
    Returns a new point to the immediate left of the last point in listOfPoints.

    Tuples in the list represent points in 2d space.
    The destination of a left turn is a tuple representing a point which on a 2d axis would
       be to the immediate left of the last point in the list representing the current path.'''
    
    newPoint = (0,0)  # a default return value
    if len(listOfPoints) > 0:  # check there is at least one point in the list
        lastPoint_x, lastPoint_y = listOfPoints[len(listOfPoints)-1]  # unpack the last point in the list
        new_x = lastPoint_x - 1  # a new x coordinate, one unit to the left of the last one in the list
        newPoint = (new_x, lastPoint_y)  # a new point one unit to the left of the last one in the list
    return newPoint
In [102]:
def makeRightTurnPoint(listOfPoints):
    '''Function to make a point representing the destination of a right turn from the current path.
    
    Param listOfPoints is a list of tuples representing the path so far.
    Returns a new point to the immediate right of the last point in listOfPoints.

    Tuples in the list represent points in 2d space.
    The destination of a right turn is a tuple representing a point which on a 2d axis would
       be to the immediate right of the last point in the list representing the current path.'''
    
    newPoint = (0,0)  # a default return value
    if len(listOfPoints) > 0:  # check there is a least one point in the list
        lastPoint_x, lastPoint_y = listOfPoints[len(listOfPoints)-1]  # the last point in the list
        new_x = lastPoint_x + 1  # a new x coordinate one unit to the right of the last one in the list
        newPoint = (new_x, lastPoint_y) # a new point one unit to the right of the last one in the list
    return newPoint

You should be thinking "Hey that is a lot of duplicated code: I thought that functions should help us not to duplicate code!". You are completely right, but this example will be easier to understand if we just live with some bad programming for the sake of clarity.

Now, experience the flexibility of Python: We can have lists of functions!

In [104]:
movementFunctions = [makeLeftTurnPoint, makeRightTurnPoint]
type(movementFunctions[0])
Out[104]:
<type 'function'>

To demonstrate how we can use our list of functions, take some path which is represented by a list of points (tuples):

In [105]:
somePath = [(0,0), (1,0), (2,0)]
# find the point that is the rightTurn from the last point in the path

If you want to, see if you can find the tuple which is the last point in the path so far using len(somePath) and the indexing operator [ ].

We put the functions into the list in order, left to right, so we know that the first function in the list (index 0) is makeLeftTurnPoint, and the second function in the list (index 1) is makeRightTurnPoint.

We can now call makeLeftTurnPoint by calling movementFunctions[0], and passing it the arguement somePath which is our example path so far. We should get back a new point (tuple) which is the point to the immediate left of the last point in the list somePath: