import py_ecc, a pure python fec by Emin Martinian, which is under a permissive licence
authorZooko O'Whielacronx <zooko@zooko.com>
Fri, 29 Dec 2006 20:46:45 +0000 (13:46 -0700)
committerZooko O'Whielacronx <zooko@zooko.com>
Fri, 29 Dec 2006 20:46:45 +0000 (13:46 -0700)
It is too slow for a real product, but is a quick way to get a working prototype, and also is freely redistributable by us...

src/py_ecc/ffield.py [new file with mode: 0644]
src/py_ecc/file_ecc.py [new file with mode: 0644]
src/py_ecc/genericmatrix.py [new file with mode: 0644]
src/py_ecc/rs_code.py [new file with mode: 0644]

diff --git a/src/py_ecc/ffield.py b/src/py_ecc/ffield.py
new file mode 100644 (file)
index 0000000..a17a359
--- /dev/null
@@ -0,0 +1,764 @@
+
+# Copyright Emin Martinian 2002.  See below for license terms.
+# Version Control Info: $Id: ffield.py,v 1.10 2003/10/28 21:19:43 emin Exp $
+
+"""
+This package contains the FField class designed to perform calculations
+in finite fields of characteristic two.  The following docstrings provide
+detailed information on various topics:
+
+  FField.__doc__    Describes the methods of the FField class and how
+                    to use them.
+
+  FElement.__doc__  Describes the FElement class and how to use it.
+
+  fields_doc        Briefly describes what a finite field is and
+                    establishes notation for further documentation.
+
+  design_doc        Discusses the design of the FField class and attempts
+                    to justify why certain decisions were made.
+
+  license_doc       Describes the license and lack of warranty for
+                    this code.
+
+  testing_doc       Describes some tests to make sure the code is working
+                    as well as some of the built in testing routines.
+  
+"""
+
+import string, random, os, os.path, cPickle
+
+
+# The following list of primitive polynomials are the Conway Polynomials
+# from the list at
+# http://www.math.rwth-aachen.de/~Frank.Luebeck/ConwayPol/cp2.html
+
+gPrimitivePolys = {}
+gPrimitivePolysCondensed = {
+    1  : (1,0),
+    2  : (2,1,0),
+    3  : (3,1,0),
+    4  : (4,1,0),
+    5  : (5,2,0),
+    6  : (6,4,3,1,0),
+    7  : (7,1,0),
+    8  : (8,4,3,2,0),
+    9  : (9,4,0),
+    10 : (10,6,5,3,2,1,0),
+    11 : (11,2,0),
+    12 : (12,7,6,5,3,1,0),
+    13 : (13,4,3,1,0),
+    14 : (14,7,5,3,0),
+    15 : (15,5,4,2,0),
+    16 : (16,5,3,2,0),
+    17 : (17,3,0),
+    18 : (18,12,10,1,0),
+    19 : (19,5,2,1,0),
+    20 : (20,10,9,7,6,5,4,1,0),
+    21 : (21,6,5,2,0),
+    22 : (22,12,11,10,9,8,6,5,0),
+    23 : (23,5,0),
+    24 : (24,16,15,14,13,10,9,7,5,3,0),
+    25 : (25,8,6,2,0),
+    26 : (26,14,10,8,7,6,4,1,0),
+    27 : (27,12,10,9,7,5,3,2,0),
+    28 : (28,13,7,6,5,2,0),
+    29 : (29,2,0),
+    30 : (30,17,16,13,11,7,5,3,2,1,0),
+    31 : (31,3,0),
+    32 : (32,15,9,7,4,3,0),
+    33 : (33,13,12,11,10,8,6,3,0),
+    34 : (34,16,15,12,11,8,7,6,5,4,2,1,0),
+    35 : (35, 11, 10, 7, 5, 2, 0),
+    36 : (36, 23, 22, 20, 19, 17, 14, 13, 8, 6, 5, 1, 0),
+    37 : (37, 5, 4, 3, 2, 1, 0),
+    38 : (38, 14, 10, 9, 8, 5, 2, 1, 0),
+    39 : (39, 15, 12, 11, 10, 9, 7, 6, 5, 2 , 0),
+    40 : (40, 23, 21, 18, 16, 15, 13, 12, 8, 5, 3, 1, 0),
+    97 : (97,6,0),
+    100 : (100,15,0)
+    }
+
+for n in gPrimitivePolysCondensed.keys():
+    gPrimitivePolys[n] = [0]*(n+1)
+    if (n < 16):
+        unity = 1
+    else:
+        unity = long(1)
+    for index in gPrimitivePolysCondensed[n]:
+        gPrimitivePolys[n][index] = unity
+                
+
+class FField:
+    """
+    The FField class implements a finite field calculator.  The
+    following functions are provided:
+
+    __init__
+    Add
+    Subtract
+    Multiply
+    Inverse
+    Divide
+    FindDegree
+    MultiplyWithoutReducing
+    ExtendedEuclid
+    FullDivision
+    ShowCoefficients
+    ShowPolynomial
+    GetRandomElement
+    ConvertListToElement
+    TestFullDivision
+    TestInverse
+
+    Most of these methods take integers or longs representing field 
+    elements as arguments and return integers representing the desired 
+    field elements as output.  See ffield.fields_doc for an explanation
+    of the integer representation of field elements.
+
+    Example of how to use the FField class:
+    
+>>> import ffield
+>>> F = ffield.FField(5) # create the field GF(2^5) 
+>>> a = 7 # field elements are denoted as integers from 0 to 2^5-1
+>>> b = 15
+>>> F.ShowPolynomial(a) # show the polynomial representation of a
+'x^2 + x^1 + 1'
+>>> F.ShowPolynomial(b)
+'x^3 + x^2 + x^1 + 1'
+>>> c = F.Multiply(a,b) # multiply a and b modulo the field generator
+>>> c
+4
+>>> F.ShowPolynomial(c)
+'x^2'
+>>> F.Multiply(c,F.Inverse(a)) == b # verify multiplication works
+1
+>>> F.Multiply(c,F.Inverse(b)) == a # verify multiplication works
+1
+>>> d = F.Divide(c,b) # since c = F.Multiply(a,b), d should give a
+>>> d
+7
+
+    See documentation on the appropriate method for further details.
+    """
+    
+    def __init__(self,n,gen=0,useLUT=-1):
+        """
+        This method constructs the field GF(2^p).  It takes one
+        required argument, n = p, and two optional arguments, gen,
+        representing the coefficients of the generator polynomial
+        (of degree n) to use and useLUT describing whether to use
+        a lookup table.  If no gen argument is provided, the
+        Conway Polynomial of degree n is obtained from the table
+        gPrimitivePolys.
+
+        If useLUT = 1 then a lookup table is used for
+        computing finite field multiplies and divides.
+        If useLUT = 0 then no lookup table is used.
+        If useLUT = -1 (the default), then the code
+        decides when a lookup table should be used.
+
+        Note that you can look at the generator for the field object
+        F by looking at F.generator.
+        """
+        
+        self.n = n 
+        if (gen):
+            self.generator = gen
+        else:
+            self.generator = self.ConvertListToElement(gPrimitivePolys[n])
+
+
+        if (useLUT == 1 or (useLUT == -1 and self.n < 10)): # use lookup table
+            self.unity = 1
+            self.Inverse = self.DoInverseForSmallField            
+            self.PrepareLUT()
+            self.Multiply = self.LUTMultiply
+            self.Divide = self.LUTDivide
+            self.Inverse = lambda x: self.LUTDivide(1,x)            
+        elif (self.n < 15):
+            self.unity = 1
+            self.Inverse = self.DoInverseForSmallField            
+            self.Multiply = self.DoMultiply
+            self.Divide = self.DoDivide
+        else: # Need to use longs for larger fields
+            self.unity = long(1)
+            self.Inverse = self.DoInverseForBigField            
+            self.Multiply = lambda a,b: self.DoMultiply(long(a),long(b))
+            self.Divide = lambda a,b: self.DoDivide(long(a),long(b))
+
+
+
+    def PrepareLUT(self):
+        fieldSize = 1 << self.n
+        lutName = 'ffield.lut.' + `self.n`
+        if (os.path.exists(lutName)):
+            fd = open(lutName,'r')
+            self.lut = cPickle.load(fd)
+            fd.close()
+        else:
+            self.lut = LUT()
+            self.lut.mulLUT = range(fieldSize)
+            self.lut.divLUT = range(fieldSize)
+            self.lut.mulLUT[0] = [0]*fieldSize
+            self.lut.divLUT[0] = ['NaN']*fieldSize
+            for i in range(1,fieldSize):
+                self.lut.mulLUT[i] = map(lambda x: self.DoMultiply(i,x),
+                                         range(fieldSize))
+                self.lut.divLUT[i] = map(lambda x: self.DoDivide(i,x),
+                                         range(fieldSize))
+            fd = open(lutName,'w')
+            cPickle.dump(self.lut,fd)
+            fd.close()
+
+            
+    def LUTMultiply(self,i,j):
+        return self.lut.mulLUT[i][j]
+
+    def LUTDivide(self,i,j):
+        return self.lut.divLUT[i][j]
+        
+    def Add(self,x,y):
+        """
+        Adds two field elements and returns the result.
+        """
+        
+        return x ^ y
+
+    def Subtract(self,x,y):
+        """
+        Subtracts the second argument from the first and returns
+        the result.  In fields of characteristic two this is the same
+        as the Add method.
+        """
+        return self.Add(x,y)
+
+    def DoMultiply(self,f,v):
+        """
+        Multiplies two field elements (modulo the generator
+        self.generator) and returns the result.
+
+        See MultiplyWithoutReducing if you don't want multiplication
+        modulo self.generator.
+        """
+        m = self.MultiplyWithoutReducing(f,v)
+        return self.FullDivision(m,self.generator,
+                                 self.FindDegree(m),self.n)[1]
+            
+    def DoInverseForSmallField(self,f):
+        """
+        Computes the multiplicative inverse of its argument and
+        returns the result.
+        """
+        return self.ExtendedEuclid(1,f,self.generator,
+                                   self.FindDegree(f),self.n)[1]
+
+    def DoInverseForBigField(self,f):
+        """
+        Computes the multiplicative inverse of its argument and
+        returns the result.
+        """
+        return self.ExtendedEuclid(self.unity,long(f),self.generator,
+                                   self.FindDegree(long(f)),self.n)[1]
+
+    def DoDivide(self,f,v):
+        """
+        Divide(f,v) returns f * v^-1.
+        """
+        return self.DoMultiply(f,self.Inverse(v))
+
+    def FindDegree(self,v):
+        """
+        Find the degree of the polynomial representing the input field
+        element v.  This takes O(degree(v)) operations.
+
+        A faster version requiring only O(log(degree(v)))
+        could be written using binary search...
+        """
+
+        if (v):
+            result = -1
+            while(v):
+                v = v >> 1
+                result = result + 1
+            return result
+        else:
+            return 0
+
+    def MultiplyWithoutReducing(self,f,v):
+        """
+        Multiplies two field elements and does not take the result
+        modulo self.generator.  You probably should not use this
+        unless you know what you are doing; look at Multiply instead.
+
+        NOTE: If you are using fields larger than GF(2^15), you should
+        make sure that f and v are longs not integers.
+        """
+        
+        result = 0
+        mask = self.unity
+        i = 0
+        while (i <= self.n):
+            if (mask & v): 
+                result = result ^ f
+            f = f << 1
+            mask = mask << 1
+            i = i + 1
+        return result
+
+
+    def ExtendedEuclid(self,d,a,b,aDegree,bDegree):
+        """
+        Takes arguments (d,a,b,aDegree,bDegree) where d = gcd(a,b)
+        and returns the result of the extended Euclid algorithm
+        on (d,a,b).
+        """
+        if (b == 0):
+            return (a,self.unity,0)
+        else:
+            (floorADivB, aModB) = self.FullDivision(a,b,aDegree,bDegree)
+            (d,x,y) = self.ExtendedEuclid(d, b, aModB, bDegree,
+                                          self.FindDegree(aModB))
+            return (d,y,self.Subtract(x,self.DoMultiply(floorADivB,y)))
+
+    def FullDivision(self,f,v,fDegree,vDegree):
+        """
+        Takes four arguments, f, v, fDegree, and vDegree where
+        fDegree and vDegree are the degrees of the field elements
+        f and v represented as a polynomials.
+        This method returns the field elements a and b such that
+
+            f(x) = a(x) * v(x) + b(x).  
+
+        That is, a is the divisor and b is the remainder, or in
+        other words a is like floor(f/v) and b is like f modulo v.
+        """
+
+        result = 0
+        i = fDegree
+        mask = self.unity << i
+        while (i >= vDegree):
+            if (mask & f):
+                result = result ^ (self.unity << (i - vDegree))
+                f = self.Subtract(f, v << (i - vDegree))
+            i = i - 1
+            mask = mask >> self.unity
+        return (result,f)
+
+
+    def ShowCoefficients(self,f):
+        """
+        Show coefficients of input field element represented as a
+        polynomial in decreasing order.
+        """
+
+        fDegree = self.n
+
+        result = []
+        for i in range(fDegree,-1,-1):
+            if ((self.unity << i) & f):
+                result.append(1)
+            else:
+                result.append(0)
+            
+        return result
+
+    def ShowPolynomial(self,f):
+        """
+        Show input field element represented as a polynomial.
+        """
+
+        fDegree = self.FindDegree(f)
+        result = ''
+
+        if (f == 0):
+            return '0'
+        
+        for i in range(fDegree,0,-1):
+            if ((1 << i) & f):
+                result = result + (' x^' + `i`)
+        if (1 & f):
+            result = result + ' ' + `1`
+        return string.replace(string.strip(result),' ',' + ')
+
+    def GetRandomElement(self,nonZero=0,maxDegree=None):
+        """
+        Return an element from the field chosen uniformly at random
+        or, if the optional argument nonZero is true, chosen uniformly
+        at random from the non-zero elements, or, if the optional argument
+        maxDegree is provided, ensure that the result has degree less
+        than maxDegree.
+        """
+
+        if (None == maxDegree):
+            maxDegree = self.n
+        if (maxDegree <= 1 and nonZero):
+            return 1
+        if (maxDegree < 31):
+            return random.randint(nonZero != 0,(1<<maxDegree)-1)
+        else:
+            result = 0L
+            for i in range(0,maxDegree):
+                result = result ^ (random.randint(0,1) << long(i))
+            if (nonZero and result == 0):
+                return self.GetRandomElement(1)
+            else:
+                return result
+                    
+
+
+    def ConvertListToElement(self,l):
+        """
+        This method takes as input a binary list (e.g. [1, 0, 1, 1])
+        and converts it to a decimal representation of a field element.
+        For example, [1, 0, 1, 1] is mapped to 8 | 2 | 1 = 11.
+
+        Note if the input list is of degree >= to the degree of the
+        generator for the field, then you will have to call take the
+        result modulo the generator to get a proper element in the
+        field.
+        """
+        
+        temp = map(lambda a, b: a << b, l, range(len(l)-1,-1,-1))
+        return reduce(lambda a, b: a | b, temp)
+
+    def TestFullDivision(self):
+        """
+        Test the FullDivision function by generating random polynomials
+        a(x) and b(x) and checking whether (c,d) == FullDivision(a,b)
+        satsifies b*c + d == a
+        """
+        f = 0
+
+        a = self.GetRandomElement(nonZero=1)
+        b = self.GetRandomElement(nonZero=1)
+        aDegree = self.FindDegree(a)
+        bDegree = self.FindDegree(b)
+
+        (c,d) = self.FullDivision(a,b,aDegree,bDegree)
+        recon = self.Add(d, self.Multiply(c,b))
+        assert (recon == a), ('TestFullDivision failed: a='
+                              + `a` + ', b=' + `b` + ', c='
+                              + `c` + ', d=' + `d` + ', recon=', recon)
+            
+    def TestInverse(self):
+        """
+        This function tests the Inverse function by generating
+        a random non-zero polynomials a(x) and checking if
+        a * Inverse(a) == 1.
+        """
+
+        a = self.GetRandomElement(nonZero=1)
+        aInv = self.Inverse(a)
+        prod = self.Multiply(a,aInv)
+        assert 1 == prod, ('TestInverse failed:' + 'a=' + `a` + ', aInv='
+                           + `aInv` + ', prod=' + `prod`)
+
+class LUT:
+    """
+    Lookup table used to speed up some finite field operations.
+    """
+    pass
+
+
+class FElement:
+    """
+    This class provides field elements which overload the
+    +,-,*,%,//,/ operators to be the appropriate field operation.
+    Note that before creating FElement objects you must first
+    create an FField object.  For example,
+    
+>>> import ffield
+>>> F = FField(5)
+>>> e1 = FElement(F,7)
+>>> e1
+x^2 + x^1 + 1
+>>> e2 = FElement(F,19)
+>>> e2
+x^4 + x^1 + 1
+>>> e3 = e1 + e2
+>>> e3
+x^4 + x^2
+>>> e4 = e3 / e2
+>>> e4
+x^4 + x^3 + x^2 + x^1 + 1
+>>> e4 * e2 == (e3)
+1
+    
+    """
+    
+    def __init__(self,field,e):
+        """
+        The constructor takes two arguments, field, and e where
+        field is an FField object and e is an integer representing
+        an element in FField.
+
+        The result is a new FElement instance.
+        """
+        self.f = e
+        self.field = field
+        
+    def __add__(self,other):
+        assert self.field == other.field
+        return FElement(self.field,self.field.Add(self.f,other.f))
+
+    def __mul__(self,other):
+        assert self.field == other.field
+        return FElement(self.field,self.field.Multiply(self.f,other.f))
+
+    def __mod__(self,o):
+        assert self.field == o.field
+        return FElement(self.field,
+                        self.field.FullDivision(self.f,o.f,
+                                                self.field.FindDegree(self.f),
+                                                self.field.FindDegree(o.f))[1])
+
+    def __floordiv__(self,o):
+        assert self.field == o.field
+        return FElement(self.field,
+                        self.field.FullDivision(self.f,o.f,
+                                                self.field.FindDegree(self.f),
+                                                self.field.FindDegree(o.f))[0])
+
+    def __div__(self,other):
+        assert self.field == other.field
+        return FElement(self.field,self.field.Divide(self.f,other.f))
+
+    def __str__(self):
+        return self.field.ShowPolynomial(self.f)
+
+    def __repr__(self):
+        return self.__str__()
+
+    def __eq__(self,other):
+        assert self.field == other.field
+        return self.f == other.f
+        
+def FullTest(testsPerField=10,sizeList=None):
+    """
+    This function runs TestInverse and TestFullDivision for testsPerField
+    random field elements for each field size in sizeList.  For example,
+    if sizeList = (1,5,7), then thests are run on GF(2), GF(2^5), and
+    GF(2^7).  If sizeList == None (which is the default), then every
+    field is tested.
+    """
+    
+    if (None == sizeList):
+        sizeList = gPrimitivePolys.keys()
+    for i in sizeList:
+        F = FField(i)
+        for j in range(testsPerField):
+            F.TestInverse()
+            F.TestFullDivision()
+
+
+fields_doc = """
+Roughly speaking a finite field is a finite collection of elements
+where most of the familiar rules of math work.  Specifically, you
+can add, subtract, multiply, and divide elements of a field and
+continue to get elements in the field.  This is useful because
+computers usually store and send information in fixed size chunks.
+Thus many useful algorithms can be described as elementary operations
+(e.g. addition, subtract, multiplication, and division) of these chunks.
+
+Currently this package only deals with fields of characteristic 2.  That
+is all fields we consider have exactly 2^p elements for some integer p.
+We denote such fields as GF(2^p) and work with the elements represented
+as p-1 degree polynomials in the indeterminate x.  That is an element of
+the field GF(2^p) looks something like
+
+     f(x) = c_{p-1} x^{p-1} + c_{p-2} x^{p-2} + ... + c_0
+
+where the coefficients c_i are in binary.
+
+Addition is performed by simply adding coefficients of degree i
+modulo 2.  For example, if we have two field elements f and v
+represented as f(x) = x^2 + 1 and v(x) = x + 1 then s = f + v
+is given by (x^2 + 1) + (x + 1) = x^2 + x.  Multiplication is
+performed modulo a p degree generator polynomial g(x).
+For example, if f and v are as in the above example, then s = s * v
+is given by (x^2 + 1) + (x + 1) mod g(x).  Subtraction turns out
+to be the same as addition for fields of characteristic 2.  Division
+is defined as f / v = f * v^-1 where v^-1 is the multiplicative
+inverse of v.  Multiplicative inverses in groups and fields
+can be calculated using the extended Euclid algorithm.
+
+Roughly speaking the intuition for why multiplication is
+performed modulo g(x), is because we want to make sure s * v
+returns an element in the field.  Elements of the field are
+polynomials of degree p-1, but regular multiplication could
+yield terms of degree greater than p-1.  Therefore we need a
+rule for 'reducing' terms of degree p or greater back down
+to terms of degree at most p-1.  The 'reduction rule' is
+taking things modulo g(x).
+
+For another way to think of
+taking things modulo g(x) as a 'reduction rule', imagine
+g(x) = x^7 + x + 1 and we want to take some polynomial,
+f(x) = x^8 + x^3 + x, modulo g(x).  We can think of g(x)
+as telling us that we can replace every occurence of
+x^7 with x + 1.  Thus f(x) becomes x * x^7 + x^3 + x which
+becomes x * (x + 1) + x^3 + x = x^3 + x^2 .  Essentially, taking
+polynomials mod x^7 by replacing all x^7 terms with x + 1 will 
+force down the degree of f(x) until it is below 7 (the leading power
+of g(x).  See a book on abstract algebra for more details.
+"""
+
+design_doc = """
+The FField class implements a finite field calculator for fields of
+characteristic two.  This uses a representation of field elements
+as integers and has various methods to calculate the result of
+adding, subtracting, multiplying, dividing, etc. field elements
+represented AS INTEGERS OR LONGS.
+
+The FElement class provides objects which act like a new kind of
+numeric type (i.e. they overload the +,-,*,%,//,/ operators, and
+print themselves as polynomials instead of integers).
+
+Use the FField class for efficient storage and calculation.
+Use the FElement class if you want to play around with finite
+field math the way you would in something like Matlab or
+Mathematica.
+
+--------------------------------------------------------------------
+                           WHY PYTHON?
+
+You may wonder why a finite field calculator written in Python would
+be useful considering all the C/C++/Java code already written to do
+the same thing (and probably faster too).  The goals of this project
+are as follows, please keep them in mind if you make changes:
+
+o  Provide an easy to understand implementation of field operations.
+   Python lends itself well to comments and documentation.  Hence,
+   we hope that in addition to being useful by itself, this project
+   will make it easier for people to implement finite field
+   computations in other languages.  If you've ever looked at some
+   of the highly optimized finite field code written in C, you will
+   understand the need for a clear reference implementation of such
+   operations.
+
+o  Provide easy access to a finite field calculator.
+   Since you can just start up the Python interpreter and do
+   computations, a finite field calculator in Python lets you try
+   things out, check your work for other algorithms, etc.
+   Furthermore since a wealth of numerical packages exist for python,
+   you can easily write simulations or algorithms which draw upon
+   such routines with finite fields.
+
+o  Provide a platform independent framework for coding in Python.
+   Many useful error control codes can be implemented based on
+   finite fields.  Some examples include error/erasure correction,
+   cyclic redundancy checks (CRCs), and secret sharing.  Since
+   Python has a number of other useful Internet features being able
+   to implement these kinds of codes makes Python a better framework
+   for network programming.
+
+o  Leverages Python arbitrary precision code for large fields.
+   If you want to do computations over very large fields, for example
+   GF(2^p) with p > 31 you have to write lots of ugly bit field
+   code in most languages.  Since Python has built in support for
+   arbitrary precision integers, you can make this code work for
+   arbitrary field sizes provided you operate on longs instead of
+   ints.  That is if you give as input numbers like
+   0L, 1L, 1L << 55, etc., most of the code should work.
+
+--------------------------------------------------------------------
+                            BASIC DESIGN
+
+
+The basic idea is to index entries in the finite field of interest
+using integers and design the class methods to work properly on this
+representation.  Using integers is efficient since integers are easy
+to store and manipulate and allows us to handle arbitrary field sizes
+without changing the code if we instead switch to using longs.
+
+Specifically, an integer represents a bit string
+
+  c = c_{p-1} c_{p-2} ... c_0.
+
+which we interpret as the coefficients of a polynomial representing a
+field element
+
+  f(x) = c_{p-1} x^{p-1} + c_{p-2} x^{p-2} + ... + c_0.
+
+--------------------------------------------------------------------
+                             FUTURE
+In the future, support for fields of other
+characteristic may be added (if people want them).  Since computers
+have built in parallelized operations for fields of characteristic
+two (i.e. bitwise and, or, xor, etc.), this implementation uses
+such operations to make most of the computations efficient.
+
+"""
+  
+
+license_doc = """
+  This code was originally written by Emin Martinian (emin@allegro.mit.edu).
+  You may copy, modify, redistribute in source or binary form as long
+  as credit is given to the original author.  Specifically, please
+  include some kind of comment or docstring saying that Emin Martinian
+  was one of the original authors.  Also, if you publish anything based
+  on this work, it would be nice to cite the original author and any
+  other contributers.
+
+  There is NO WARRANTY for this software just as there is no warranty
+  for GNU software (although this is not GNU software).  Specifically
+  we adopt the same policy towards warranties as the GNU project:
+
+  BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM 'AS IS' WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+"""
+
+testing_doc = """
+The FField class has a number of built in testing functions such as
+TestFullDivision, TestInverse.  The simplest thing to
+do is to call the FullTest method.  
+
+>>> import ffield
+>>> ffield.FullTest(sizeList=None,testsPerField=100)
+
+# To decrease the testing time you can either decrease the testsPerField
+# or you can only test the field sizes you care about by doing something
+# like sizeList = [2,7,20] in the ffield.FullTest command above.
+
+If any problems occur, assertion errors are raised.  Otherwise
+nothing is returned.  Note that you can also use the doctest
+package to test all the python examples in the documentation
+by typing 'python ffield.py' or 'python -v ffield.py' at the
+command line.
+"""
+
+
+# The following code is used to make the doctest package
+# check examples in docstrings.
+
+__test__ = {
+    'testing_doc' : testing_doc
+}
+
+def _test():
+    import doctest, ffield
+    return doctest.testmod(ffield)
+
+if __name__ == "__main__":
+    print 'Starting automated tests (this may take a while)'
+    _test()
+    print 'Tests passed.'
+
diff --git a/src/py_ecc/file_ecc.py b/src/py_ecc/file_ecc.py
new file mode 100644 (file)
index 0000000..66d9b65
--- /dev/null
@@ -0,0 +1,218 @@
+
+# Copyright Emin Martinian 2002.  See below for license terms.
+# Version Control Info: $Id: file_ecc.py,v 1.4 2003/10/28 21:28:29 emin Exp $
+
+__doc__ = """
+This package implements an erasure correction code for files.
+Specifically it lets you take a file F and break it into N
+pieces (which are named F.p_0, F.p_1, ..., F.p_N-1) such that
+F can be recovered from any K pieces.  Since the size of each
+piece is F/K (plus some small header information).
+
+How is this better than simply repeating copies of a file?
+
+Firstly, this package lets you get finer grained
+redunancy control since producing a duplicate copy of a file
+requires at least 100% redundancy while this package lets you
+expand the redunancy by n/k (e.g. if n=11, k=10 only 10%
+redundancy is added).
+
+Secondly, using a Reed-Solomon code as is done in this package,
+allows better loss resistance.  For example, assume you just
+divided a file F into 4 pieces, F.1, F.2, ..., F.4, duplicated
+each piece and placed them each on a different disk.  If the
+two disks containing a copy of F.1 go down then you can no longer
+recover F.
+
+With the Reed-Solomon code used in this package, if you use n=8, k=4
+you divide F into 8 pieces such that as long as at least 4 pieces are
+available recovery can occur.  Thus if you placed each piece on a
+seprate disk, you could recover data as if any combination of 4 or
+less disks fail.
+
+The docstrings for the functions EncodeFile and DecodeFiles
+provide detailed information on usage and the docstring
+license_doc describes the license and lack of warranty.
+
+The following is an example of how to use this file:
+
+>>> import file_ecc
+>>> testFile = '/bin/ls'      # A reasonable size file for testing.
+>>> prefix = '/tmp/ls_backup' # Prefix for shares of file.
+>>> names = file_ecc.EncodeFile(testFile,prefix,15,11) # break into N=15 pieces
+
+# Imagine that only pieces [0,1,5,4,13,8,9,10,11,12,14] are available.
+>>> decList = map(lambda x: prefix + '.p_' + `x`,[0,1,5,4,13,8,9,10,11,12,14])
+
+>>> decodedFile = '/tmp/ls.r' # Choose where we want reconstruction to go.
+>>> file_ecc.DecodeFiles(decList,decodedFile)
+>>> fd1 = open(testFile,'rb')
+>>> fd2 = open(decodedFile,'rb')
+>>> fd1.read() == fd2.read()
+1
+"""
+
+
+
+from rs_code import RSCode
+from array import array
+
+import os, struct, string
+
+headerSep = '|'
+
+def GetFileSize(fname):
+    return os.stat(fname)[6]
+
+def MakeHeader(fname,n,k,size):
+    return string.join(['RS_PARITY_PIECE_HEADER','FILE',fname,
+                        'n',`n`,'k',`k`,'size',`size`,'piece'],
+                       headerSep) + headerSep
+
+def ParseHeader(header):
+    return string.split(header,headerSep)
+
+def ReadEncodeAndWriteBlock(readSize,inFD,outFD,code):
+    buffer = array('B')
+    buffer.fromfile(inFD,readSize)
+    for i in range(readSize,code.k):
+        buffer.append(0)
+    codeVec = code.Encode(buffer)
+    for j in range(code.n):
+        outFD[j].write(struct.pack('B',codeVec[j]))
+
+def EncodeFile(fname,prefix,n,k):
+    """
+    Function:     EncodeFile(fname,prefix,n,k)
+    Description:  Encodes the file named by fname into n pieces named
+                  prefix.p_0, prefix.p_1, ..., prefix.p_n-1.  At least
+                  k of these pieces are needed for recovering fname.
+                  Each piece is roughly the size of fname / k (there
+                  is a very small overhead due to some header information).
+
+                  Returns a list containing names of files for the pieces.
+
+                  Note n and k must satisfy 0 < k < n < 257.
+                  Use the DecodeFiles function for decoding.
+    """
+    fileList = []
+    if (n > 256 or k >= n or k <= 0):
+        raise Exception, 'Invalid (n,k), need 0 < k < n < 257.'
+    inFD = open(fname,'rb')
+    inSize = GetFileSize(fname)
+    header = MakeHeader(fname,n,k,inSize)
+    code = RSCode(n,k,8,shouldUseLUT=-(k!=1))
+    outFD = range(n)
+    for i in range(n):
+        outFileName = prefix + '.p_' + `i`
+        fileList.append(outFileName)
+        outFD[i] = open(outFileName,'wb')
+        outFD[i].write(header + `i` + '\n')
+
+    if (k == 1): # just doing repetition coding
+        str = inFD.read(1024)
+        while (str):
+            map( lambda x: x.write(str), outFD)
+            str = inFD.read(256)
+    else: # do the full blown RS encodding
+        for i in range(0, (inSize/k)*k,k):
+            ReadEncodeAndWriteBlock(k,inFD,outFD,code)
+        
+        if ((inSize % k) > 0):
+            ReadEncodeAndWriteBlock(inSize % k,inFD,outFD,code)
+
+    return fileList
+
+def ExtractPieceNums(fnames,headers):
+    l = range(len(fnames))
+    pieceNums = range(len(fnames))
+    for i in range(len(fnames)):
+        l[i] = ParseHeader(headers[i])
+    for i in range(len(fnames)):
+        if (l[i][0] != 'RS_PARITY_PIECE_HEADER' or
+            l[i][2] != l[0][2] or l[i][4] != l[0][4] or
+            l[i][6] != l[0][6] or l[i][8] != l[0][8]):
+            raise Exception, 'File ' + `fnames[i]` + ' has incorrect header.'
+        pieceNums[i] = int(l[i][10])
+    (n,k,size) = (int(l[0][4]),int(l[0][6]),long(l[0][8]))
+    if (len(pieceNums) < k):
+        raise Exception, ('Not enough parity for decoding; needed '
+                          + `l[0][6]` + ' got ' + `len(fnames)` + '.')
+    return (n,k,size,pieceNums)
+
+def ReadDecodeAndWriteBlock(writeSize,inFDs,outFD,code):
+    buffer = array('B')
+    for j in range(code.k):
+        buffer.fromfile(inFDs[j],1)
+    result = code.Decode(buffer.tolist())
+    for j in range(writeSize):
+        outFD.write(struct.pack('B',result[j]))
+
+
+def DecodeFiles(fnames,outName):
+    """
+    Function:     DecodeFiles(fnames,outName)
+    Description:  Takes pieces of a file created using EncodeFiles and
+                  recovers the original file placing it in outName.
+                  The argument fnames must be a list of at least k
+                  file names generated using EncodeFiles.
+    """
+    inFDs = range(len(fnames))
+    headers = range(len(fnames))
+    for i in range(len(fnames)):
+        inFDs[i] = open(fnames[i],'rb')
+        headers[i] = inFDs[i].readline()
+    (n,k,inSize,pieceNums) = ExtractPieceNums(fnames,headers)
+    outFD = open(outName,'wb')
+    code = RSCode(n,k,8)
+    decList = pieceNums[0:k]
+    code.PrepareDecoder(decList)
+    for i in range(0, (inSize/k)*k,k):
+        ReadDecodeAndWriteBlock(k,inFDs,outFD,code)
+    if ((inSize%k)>0):
+        ReadDecodeAndWriteBlock(inSize%k,inFDs,outFD,code)
+
+license_doc = """
+  This code was originally written by Emin Martinian (emin@allegro.mit.edu).
+  You may copy, modify, redistribute in source or binary form as long
+  as credit is given to the original author.  Specifically, please
+  include some kind of comment or docstring saying that Emin Martinian
+  was one of the original authors.  Also, if you publish anything based
+  on this work, it would be nice to cite the original author and any
+  other contributers.
+
+  There is NO WARRANTY for this software just as there is no warranty
+  for GNU software (although this is not GNU software).  Specifically
+  we adopt the same policy towards warranties as the GNU project:
+
+  BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM 'AS IS' WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+"""
+
+# The following code is used to make the doctest package
+# check examples in docstrings.
+
+def _test():
+    import doctest, file_ecc
+    return doctest.testmod(file_ecc)
+
+if __name__ == "__main__":
+    _test()
+    print 'Tests passed'
diff --git a/src/py_ecc/genericmatrix.py b/src/py_ecc/genericmatrix.py
new file mode 100644 (file)
index 0000000..1b38b4f
--- /dev/null
@@ -0,0 +1,869 @@
+
+# Copyright Emin Martinian 2002.  See below for license terms.
+# Version Control Info: $Id: genericmatrix.py,v 1.7 2003/10/28 21:18:41 emin Exp $
+
+"""
+This package implements the GenericMatrix class to provide matrix
+operations for any type that supports the multiply, add, subtract,
+and divide operators.  For example, this package can be used to
+do matrix calculations over finite fields using the ffield package
+available at http://martinian.com.
+
+The following docstrings provide detailed information on various topics:
+
+  GenericMatrix.__doc__   Describes the methods of the GenericMatrix
+                          class and how to use them.
+
+  license_doc             Describes the license and lack of warranty
+                          for this code.
+
+  testing_doc             Describes some tests to make sure the code works.
+
+"""
+
+import operator
+
+class GenericMatrix:
+
+    """
+    The GenericMatrix class implements a matrix with works with
+    any generic type supporting addition, subtraction, multiplication,
+    and division.  Matrix multiplication, addition, and subtraction
+    are implemented as are methods for finding inverses,
+    LU (actually LUP) decompositions, and determinants.  A complete
+    list of user callable methods is:
+
+    __init__
+    __repr__
+    __mul__
+    __add__
+    __sub__
+    __setitem__
+    __getitem__
+    Size
+    SetRow
+    GetRow
+    GetColumn
+    Copy
+    MakeSimilarMatrix
+    SwapRows
+    MulRow
+    AddRow
+    AddCol
+    MulAddRow
+    LeftMulColumnVec
+    LowerGaussianElim
+    Inverse
+    Determinant
+    LUP
+
+    A quick and dirty example of how to use the GenericMatrix class
+    for matricies of floats is provided below.
+    
+>>> import genericmatrix
+>>> v = genericmatrix.GenericMatrix((3,3))
+>>> v.SetRow(0,[0.0, -1.0, 1.0])
+>>> v.SetRow(1,[1.0, 1.0, 1.0])
+>>> v.SetRow(2,[1.0, 1.0, -1.0])
+>>> v
+<matrix
+  0.0 -1.0  1.0
+  1.0  1.0  1.0
+  1.0  1.0 -1.0>
+>>> vi = v.Inverse()
+>>> vi
+<matrix
+  1.0  0.0  1.0
+ -1.0  0.5 -0.5
+ -0.0  0.5 -0.5>
+>>> (vi * v) - v.MakeSimilarMatrix(v.Size(),'i')
+<matrix
+ 0.0 0.0 0.0
+ 0.0 0.0 0.0
+ 0.0 0.0 0.0>
+
+# See what happens when we try to invert a non-invertible matrix
+
+>>> v[0,1] = 0.0
+>>> v
+<matrix
+  0.0  0.0  1.0
+  1.0  1.0  1.0
+  1.0  1.0 -1.0>
+>>> abs(v.Determinant())
+0.0
+>>> v.Inverse()
+Traceback (most recent call last):
+     ...
+ValueError: matrix not invertible
+
+# LUP decomposition will still work even if Inverse() won't.
+
+>>> (l,u,p) = v.LUP()
+>>> l
+<matrix
+ 1.0 0.0 0.0
+ 0.0 1.0 0.0
+ 1.0 0.0 1.0>
+>>> u
+<matrix
+  1.0  1.0  1.0
+  0.0  0.0  1.0
+  0.0  0.0 -2.0>
+>>> p
+<matrix
+ 0.0 1.0 0.0
+ 1.0 0.0 0.0
+ 0.0 0.0 1.0>
+>>> p * v - l * u
+<matrix
+ 0.0 0.0 0.0
+ 0.0 0.0 0.0
+ 0.0 0.0 0.0>
+
+# Operate on some column vectors using v.
+# The LeftMulColumnVec methods lets us do this without having
+# to construct a new GenericMatrix to represent each column vector.
+>>> v.LeftMulColumnVec([1.0,2.0,3.0])
+[3.0, 6.0, 0.0]
+>>> v.LeftMulColumnVec([1.0,-2.0,1.0])
+[1.0, 0.0, -2.0]
+
+# Most of the stuff above could be done with something like matlab.
+# But, with this package you can do matrix ops for finite fields.
+>>> XOR = lambda x,y: x^y
+>>> AND = lambda x,y: x&y
+>>> DIV = lambda x,y: x  
+>>> m = GenericMatrix(size=(3,4),zeroElement=0,identityElement=1,add=XOR,mul=AND,sub=XOR,div=DIV)
+>>> m.SetRow(0,[0,1,0,0])
+>>> m.SetRow(1,[0,1,0,1])
+>>> m.SetRow(2,[0,0,1,0])
+>>> # You can't invert m since it isn't square, but you can still 
+>>> # get the LUP decomposition or solve a system of equations.
+>>> (l,u,p) = v.LUP()
+>>> p*v-l*u
+<matrix
+ 0.0 0.0 0.0
+ 0.0 0.0 0.0
+ 0.0 0.0 0.0>
+>>> b = [1,0,1]
+>>> x = m.Solve(b)
+>>> b == m.LeftMulColumnVec(x)
+1
+
+    """
+
+   
+    def __init__(self, size=(2,2), zeroElement=0.0, identityElement=1.0,
+                 add=operator.__add__, sub=operator.__sub__,
+                 mul=operator.__mul__, div = operator.__div__,
+                 eq = operator.__eq__, str=lambda x:`x`,
+                 equalsZero = None,fillMode='z'):
+        """
+        Function:     __init__(size,zeroElement,identityElement,
+                               add,sub,mul,div,eq,str,equalsZero fillMode)
+
+        Description:  This is the constructor for the GenericMatrix
+                      class.  All arguments are optional and default
+                      to producing a 2-by-2 zero matrix for floats.
+                      A detailed description of arguments follows:
+
+             size: A tuple of the form (numRows, numColumns)
+                   zeroElement: An object representing the additive
+                   identity (i.e. 'zero') for the data
+                   type of interest.
+                   
+             identityElement: An object representing the multiplicative
+                              identity (i.e. 'one') for the data
+                              type of interest.
+
+             add,sub,mul,div: Functions implementing basic arithmetic
+                              operations for the type of interest.
+
+             eq: A function such that eq(x,y) == 1 if and only if x == y.
+
+             str: A function used to produce a string representation of
+                  the type of interest.
+
+             equalsZero: A function used to decide if an element is
+                         essentially zero.  For floats, you could use
+                         lambda x: abs(x) < 1e-6.
+
+             fillMode: This can either be 'e' in which case the contents
+                       of the matrix are left empty, 'z', in which case
+                       the matrix is filled with zeros, 'i' in which
+                       case an identity matrix is created, or a two
+                       argument function which is called with the row
+                       and column of each index and produces the value
+                       for that entry.  Default is 'z'.
+        """
+        if (None == equalsZero):
+            equalsZero = lambda x: self.eq(self.zeroElement,x)
+
+        self.equalsZero = equalsZero
+        self.add = add
+        self.sub = sub
+        self.mul = mul
+        self.div = div
+        self.eq = eq
+        self.str = str
+        self.zeroElement = zeroElement
+        self.identityElement = identityElement
+        self.rows, self.cols = size
+        self.data = []
+
+
+        def q(x,y,z):
+            if (x):
+                return y
+            else:
+                return z
+
+        if (fillMode == 'e'):
+            return
+        elif (fillMode == 'z'):
+            fillMode = lambda x,y: self.zeroElement
+        elif (fillMode == 'i'):
+            fillMode = lambda x,y: q(self.eq(x,y),self.identityElement,
+                                     self.zeroElement)
+
+        for i in range(self.rows):
+            self.data.append(map(fillMode,[i]*self.cols,range(self.cols)))
+
+    def MakeSimilarMatrix(self,size,fillMode):
+        """
+        MakeSimilarMatrix(self,size,fillMode)
+
+        Return a matrix of the given size filled according to fillMode
+        with the same zeroElement, identityElement, add, sub, etc.
+        as self.
+
+        For example, self.MakeSimilarMatrix(self.Size(),'i') returns
+        an identity matrix of the same shape as self.
+        """
+        return GenericMatrix(size=size,zeroElement=self.zeroElement,
+                               identityElement=self.identityElement,
+                               add=self.add,sub=self.sub,
+                               mul=self.mul,div=self.div,eq=self.eq,
+                               str=self.str,equalsZero=self.equalsZero,
+                               fillMode=fillMode)
+        
+
+    def __repr__(self):
+        m = 0
+        # find the fattest element
+        for r in self.data:
+            for c in r:
+                l = len(self.str(c))
+                if l > m:
+                    m = l
+        f = '%%%ds' % (m+1)
+        s = '<matrix'
+        for r in self.data:
+            s = s + '\n'
+            for c in r:
+                s = s + (f % self.str(c))
+        s = s + '>'
+        return s
+
+    def __mul__(self,other):
+        if (self.cols != other.rows):
+            raise ValueError, "dimension mismatch"
+        result = self.MakeSimilarMatrix((self.rows,other.cols),'z')
+                               
+        for i in range(self.rows):
+            for j in range(other.cols):
+                result.data[i][j] = reduce(self.add,
+                                           map(self.mul,self.data[i],
+                                               other.GetColumn(j)))
+        return result
+
+    def __add__(self,other):
+        if (self.cols != other.rows):
+            raise ValueError, "dimension mismatch"
+        result = self.MakeSimilarMatrix(size=self.Size(),fillMode='z')
+        for i in range(self.rows):
+            for j in range(other.cols):
+                result.data[i][j] = self.add(self.data[i][j],other.data[i][j])
+        return result
+    
+    def __sub__(self,other):
+        if (self.cols != other.cols or self.rows != other.rows):
+            raise ValueError, "dimension mismatch"
+        result = self.MakeSimilarMatrix(size=self.Size(),fillMode='z')
+        for i in range(self.rows):
+            for j in range(other.cols):
+                result.data[i][j] = self.sub(self.data[i][j],
+                                             other.data[i][j])
+        return result
+
+    def __setitem__ (self, (x,y), data):
+        "__setitem__((x,y),data) sets item row x and column y to data."
+        self.data[x][y] = data
+
+    def __getitem__ (self, (x,y)):
+        "__getitem__((x,y)) gets item at row x and column y."
+        return self.data[x][y]
+
+    def Size (self):
+        "returns (rows, columns)"
+        return (len(self.data), len(self.data[0]))
+
+    def SetRow(self,r,result):
+        "SetRow(r,result) sets row r to result."
+        
+        assert len(result) == self.cols, ('Wrong # columns in row: ' +
+                                          'expected ' + `self.cols` + ', got '
+                                          + `len(result)`)
+        self.data[r] = list(result)
+
+    def GetRow(self,r):
+        "GetRow(r) returns a copy of row r."
+        return list(self.data[r])
+
+    def GetColumn(self,c):
+        "GetColumn(c) returns a copy of column c."
+        if (c >= self.cols):
+            raise ValueError, 'matrix does not have that many columns'
+        result = []
+        for r in self.data:
+            result.append(r[c])
+        return result
+
+    def Transpose(self):
+        oldData = self.data
+        self.data = []
+        for r in range(self.cols):
+            self.data.append([])
+            for c in range(self.rows):
+                self.data[r].append(oldData[c][r])
+        rows = self.rows
+        self.rows = self.cols
+        self.cols = rows
+
+    def Copy(self):
+        result = self.MakeSimilarMatrix(size=self.Size(),fillMode='e')
+
+        for r in self.data:
+            result.data.append(list(r))
+        return result
+
+    def SubMatrix(self,rowStart,rowEnd,colStart=0,colEnd=None):
+        """
+        SubMatrix(self,rowStart,rowEnd,colStart,colEnd)
+        Create and return a sub matrix containg rows
+        rowStart through rowEnd (inclusive) and columns
+        colStart through colEnd (inclusive).
+        """
+        if (not colEnd):
+            colEnd = self.cols-1
+        if (rowEnd >= self.rows):
+            raise ValueError, 'rowEnd too big: rowEnd >= self.rows'
+        result = self.MakeSimilarMatrix((rowEnd-rowStart+1,colEnd-colStart+1),
+                                        'e')
+
+        for i in range(rowStart,rowEnd+1):
+            result.data.append(list(self.data[i][colStart:(colEnd+1)]))
+
+        return result
+    def UnSubMatrix(self,rowStart,rowEnd,colStart,colEnd):
+        """
+        UnSubMatrix(self,rowStart,rowEnd,colStart,colEnd)
+        Create and return a sub matrix containg everything except
+        rows rowStart through rowEnd (inclusive) 
+        and columns colStart through colEnd (inclusive).
+        """
+        result = self.MakeSimilarMatrix((self.rows-(rowEnd-rowStart),
+                                         self.cols-(colEnd-colStart)),'e')
+
+        for i in range(0,rowStart) + range(rowEnd,self.rows):
+            result.data.append(list(self.data[i][0:colStart] +
+                                    self.data[i][colEnd:]))
+
+        return result
+
+
+    def SwapRows(self,i,j):
+        temp = list(self.data[i])
+        self.data[i] = list(self.data[j])
+        self.data[j] = temp
+
+    def MulRow(self,r,m,start=0):
+        """
+        Function: MulRow(r,m,start=0)
+        Multiply row r by m starting at optional column start (default 0).
+        """
+        row = self.data[r]
+        for i in range(start,self.cols):
+            row[i] = self.mul(row[i],m)
+
+    def AddRow(self,i,j):
+        """
+        Add row i to row j.
+        """
+        self.data[j] = map(self.add,self.data[i],self.data[j])
+
+    def AddCol(self,i,j):
+        """
+        Add column i to column j.
+        """
+        for r in range(self.rows):
+            self.data[r][j] = self.add(self.data[r][i],self.data[r][j])
+
+    def MulAddRow(self,m,i,j):
+        """
+        Multiply row i by m and add to row j.
+        """
+        self.data[j] = map(self.add,
+                           map(self.mul,[m]*self.cols,self.data[i]),
+                           self.data[j])
+
+    def LeftMulColumnVec(self,colVec):
+        """
+        Function:       LeftMulColumnVec(c)
+        Purpose:        Compute the result of self * c.
+        Description:    This function taks as input a list c,
+                        computes the desired result and returns it
+                        as a list.  This is sometimes more convenient
+                        than constructed a new GenericMatrix to represent
+                        c, computing the result and extracting c to a list.
+        """
+        if (self.cols != len(colVec)):
+            raise ValueError, 'dimension mismatch'
+        result = range(self.rows)
+        for r in range(self.rows):
+            result[r] = reduce(self.add,map(self.mul,self.data[r],colVec))
+        return result
+
+    def FindRowLeader(self,startRow,c):
+        for r in range(startRow,self.rows):
+            if (not self.eq(self.zeroElement,self.data[r][c])):
+                return r
+        return -1
+
+    def FindColLeader(self,r,startCol):
+        for c in range(startCol,self.cols):
+            if (not self.equalsZero(self.data[r][c])):
+                return c
+        return -1    
+
+    def PartialLowerGaussElim(self,rowIndex,colIndex,resultInv):
+        """
+        Function: PartialLowerGaussElim(rowIndex,colIndex,resultInv)
+        
+        This function does partial Gaussian elimination on the part of
+        the matrix on and below the main diagonal starting from
+        rowIndex.  In addition to modifying self, this function
+        applies the required elmentary row operations to the input
+        matrix resultInv.
+
+        By partial, what we mean is that if this function encounters
+        an element on the diagonal which is 0, it stops and returns
+        the corresponding rowIndex.  The caller can then permute
+        self or apply some other operation to eliminate the zero
+        and recall PartialLowerGaussElim.
+        
+        This function is meant to be combined with UpperInverse
+        to compute inverses and LU decompositions.
+        """
+
+        lastRow = self.rows-1
+        while (rowIndex < lastRow):
+            if (colIndex >= self.cols):
+                return (rowIndex, colIndex)
+            if (self.eq(self.zeroElement,self.data[rowIndex][colIndex])):
+                # self[rowIndex,colIndex] = 0 so quit.
+                return (rowIndex, colIndex)
+            divisor = self.div(self.identityElement,
+                               self.data[rowIndex][colIndex])
+            for k in range(rowIndex+1,self.rows):
+                nextTerm = self.data[k][colIndex]
+                if (self.zeroElement != nextTerm):
+                    multiple = self.mul(divisor,self.sub(self.zeroElement,
+                                                         nextTerm))
+                    self.MulAddRow(multiple,rowIndex,k)
+                    resultInv.MulAddRow(multiple,rowIndex,k)
+            rowIndex = rowIndex + 1
+            colIndex = colIndex + 1
+        return (rowIndex, colIndex)
+
+    def LowerGaussianElim(self,resultInv=''):
+        """
+        Function:       LowerGaussianElim(r)
+        Purpose:        Perform Gaussian elimination on self to eliminate
+                        all terms below the diagonal.
+        Description:    This method modifies self via Gaussian elimination
+                        and applies the elementary row operations used in
+                        this transformation to the input matrix, r
+                        (if one is provided, otherwise a matrix with
+                         identity elements on the main diagonal is
+                         created to serve the role of r).
+
+                        Thus if the input, r, is an identity matrix, after
+                        the call it will represent the transformation
+                        made to perform Gaussian elimination.
+
+                        The matrix r is returned.
+        """
+        if (resultInv == ''):
+            resultInv = self.MakeSimilarMatrix(self.Size(),'i')
+            
+        (rowIndex,colIndex) = (0,0)
+        lastRow = min(self.rows - 1,self.cols)
+        lastCol = self.cols - 1
+        while( rowIndex < lastRow and colIndex < lastCol):
+            leader = self.FindRowLeader(rowIndex,colIndex)
+            if (leader < 0):
+                colIndex = colIndex + 1
+                continue
+            if (leader != rowIndex):
+                resultInv.AddRow(leader,rowIndex)
+                self.AddRow(leader,rowIndex)
+            (rowIndex,colIndex) = (
+                self.PartialLowerGaussElim(rowIndex,colIndex,resultInv))
+        return resultInv
+
+    def UpperInverse(self,resultInv=''):
+        """
+        Function: UpperInverse(resultInv)
+        
+        Assumes that self is an upper triangular matrix like
+
+          [a b c ... ]
+          [0 d e ... ]
+          [0 0 f ... ]
+          [.     .   ]
+          [.      .  ]
+          [.       . ]
+
+        and performs Gaussian elimination to transform self into
+        the identity matrix.  The required elementary row operations
+        are applied to the matrix resultInv passed as input.  For
+        example, if the identity matrix is passed as input, then the
+        value returned is the inverse of self before the function
+        was called.
+
+        If no matrix, resultInv, is provided as input then one is
+        created with identity elements along the main diagonal.
+        In either case, resultInv is returned as output.
+        """
+        if (resultInv == ''):
+            resultInv = self.MakeSimilarMatrix(self.Size(),'i')
+        lastCol = min(self.rows,self.cols)
+        for colIndex in range(0,lastCol):
+            if (self.zeroElement == self.data[colIndex][colIndex]):
+                raise ValueError, 'matrix not invertible'
+            divisor = self.div(self.identityElement,
+                               self.data[colIndex][colIndex])
+            if (self.identityElement != divisor):
+                self.MulRow(colIndex,divisor,colIndex)
+                resultInv.MulRow(colIndex,divisor)
+            for rowToElim in range(0,colIndex):
+                multiple = self.sub(self.zeroElement,
+                                    self.data[rowToElim][colIndex])
+                self.MulAddRow(multiple,colIndex,rowToElim)
+                resultInv.MulAddRow(multiple,colIndex,rowToElim)
+        return resultInv
+    
+    def Inverse(self):
+        """
+        Function:       Inverse
+        Description:    Returns the inverse of self without modifying
+                        self.  An exception is raised if the matrix
+                        is not invertable.
+        """
+
+        workingCopy = self.Copy()
+        result = self.MakeSimilarMatrix(self.Size(),'i')
+        workingCopy.LowerGaussianElim(result)
+        workingCopy.UpperInverse(result)
+        return result
+
+    def Determinant(self):
+        """
+        Function:       Determinant
+        Description:    Returns the determinant of the matrix or raises
+                        a ValueError if the matrix is not square.
+        """
+        if (self.rows != self.cols):
+            raise ValueError, 'matrix not square'
+        workingCopy = self.Copy()
+        result = self.MakeSimilarMatrix(self.Size(),'i')
+        workingCopy.LowerGaussianElim(result)
+        det = self.identityElement
+        for i in range(self.rows):
+            det = det * workingCopy.data[i][i]
+        return det
+
+    def LUP(self):
+        """
+        Function:       (l,u,p) = self.LUP()
+        Purpose:        Compute the LUP decomposition of self.
+        Description:    This function returns three matrices
+                        l, u, and p such that p * self = l * u
+                        where l, u, and p have the following properties:
+
+                        l is lower triangular with ones on the diagonal
+                        u is upper triangular 
+                        p is a permutation matrix.
+
+                        The idea behind the algorithm is to first
+                        do Gaussian elimination to obtain an upper
+                        triangular matrix u and lower triangular matrix
+                        r such that r * self = u, then by inverting r to
+                        get l = r ^-1 we obtain self = r^-1 * u = l * u.
+                        Note tha since r is lower triangular its
+                        inverse must also be lower triangular.
+
+                        Where does the p come in?  Well, with some
+                        matrices our technique doesn't work due to
+                        zeros appearing on the diagonal of r.  So we
+                        apply some permutations to the orginal to
+                        prevent this.
+                        
+        """
+        upper = self.Copy()
+        resultInv = self.MakeSimilarMatrix(self.Size(),'i')
+        perm = self.MakeSimilarMatrix((self.rows,self.rows),'i')
+
+        (rowIndex,colIndex) = (0,0)
+        lastRow = self.rows - 1
+        lastCol = self.cols - 1
+        while( rowIndex < lastRow and colIndex < lastCol ):
+            leader = upper.FindRowLeader(rowIndex,colIndex)
+            if (leader < 0):
+                colIndex = colIndex+1
+                continue            
+            if (leader != rowIndex):
+                upper.SwapRows(leader,rowIndex)
+                resultInv.SwapRows(leader,rowIndex)
+                perm.SwapRows(leader,rowIndex)
+            (rowIndex,colIndex) = (
+                upper.PartialLowerGaussElim(rowIndex,colIndex,resultInv))
+
+        lower = self.MakeSimilarMatrix((self.rows,self.rows),'i')
+        resultInv.LowerGaussianElim(lower)
+        resultInv.UpperInverse(lower)
+        # possible optimization: due perm*lower explicitly without
+        # relying on the * operator.
+        return (perm*lower, upper, perm)
+    
+    def Solve(self,b):
+        """
+        Solve(self,b):
+
+        b:   A list.
+        
+        Returns the values of x such that Ax = b.
+
+        This is done using the LUP decomposition by 
+        noting that Ax = b implies PAx = Pb implies LUx = Pb.
+        First we solve for Ly = Pb and then we solve Ux = y.
+        The following is an example of how to use Solve:
+
+>>> # Floating point example
+>>> import genericmatrix
+>>> A = genericmatrix.GenericMatrix(size=(2,5),str=lambda x: '%.4f' % x)
+>>> A.SetRow(0,[0.0, 0.0, 0.160, 0.550, 0.280])
+>>> A.SetRow(1,[0.0, 0.0, 0.745, 0.610, 0.190])
+>>> A
+<matrix
+ 0.0000 0.0000 0.1600 0.5500 0.2800
+ 0.0000 0.0000 0.7450 0.6100 0.1900>
+>>> b = [0.975, 0.350]
+>>> x = A.Solve(b)
+>>> z = A.LeftMulColumnVec(x)
+>>> diff = reduce(lambda xx,yy: xx+yy,map(lambda aa,bb: abs(aa-bb),b,z))
+>>> diff > 1e-6
+0
+>>> # Boolean example
+>>> XOR = lambda x,y: x^y
+>>> AND = lambda x,y: x&y
+>>> DIV = lambda x,y: x  
+>>> m=GenericMatrix(size=(3,6),zeroElement=0,identityElement=1,add=XOR,mul=AND,sub=XOR,div=DIV)
+>>> m.SetRow(0,[1,0,0,1,0,1])
+>>> m.SetRow(1,[0,1,1,0,1,0])
+>>> m.SetRow(2,[0,1,0,1,1,0])
+>>> b = [0, 1, 1]
+>>> x = m.Solve(b)
+>>> z = m.LeftMulColumnVec(x)
+>>> z
+[0, 1, 1]
+
+        """
+        assert self.cols >= self.rows
+        
+        (L,U,P) = self.LUP()
+        Pb = P.LeftMulColumnVec(b)
+        y = [0]*len(Pb)
+        for row in range(L.rows):
+            y[row] = Pb[row]
+            for i in range(row+1,L.rows):
+                Pb[i] = L.sub(Pb[i],L.mul(L[i,row],Pb[row]))
+        x = [0]*self.cols
+        curRow = self.rows-1
+
+        for curRow in range(len(y)-1,-1,-1):
+            col = U.FindColLeader(curRow,0)
+            assert col > -1
+            x[col] = U.div(y[curRow],U[curRow,col])
+            y[curRow] = x[col]
+            for i in range(0,curRow):
+                y[i] = U.sub(y[i],U.mul(U[i,col],y[curRow]))
+        return x
+
+
+def DotProduct(mul,add,x,y):
+    """
+    Function:    DotProduct(mul,add,x,y)
+    Description: Return the dot product of lists x and y using mul and
+                 add as the multiplication and addition operations.
+    """
+    assert len(x) == len(y), 'sizes do not match'
+    return reduce(add,map(mul,x,y))
+
+class GenericMatrixTester:
+    def DoTests(self,numTests,sizeList):
+        """
+        Function:       DoTests(numTests,sizeList)
+
+        Description:    For each test, run numTests tests for square
+                        matrices with the sizes in sizeList.
+        """
+
+        for size in sizeList:
+            self.RandomInverseTest(size,numTests)
+            self.RandomLUPTest(size,numTests)
+            self.RandomSolveTest(size,numTests)
+            self.RandomDetTest(size,numTests)
+
+
+    def MakeRandom(self,s):
+        import random 
+        r = GenericMatrix(size=s,fillMode=lambda x,y: random.random(),
+                          equalsZero = lambda x: abs(x) < 1e-6)
+        return r
+
+    def MatAbs(self,m):
+        r = -1
+        (N,M) = m.Size()
+        for i in range(0,N):
+            for j in range(0,M):
+                if (abs(m[i,j]) > r):
+                    r = abs(m[i,j])
+        return r
+
+    def RandomInverseTest(self,s,n):
+        ident = GenericMatrix(size=(s,s),fillMode='i')
+        for i in range(n):
+            m = self.MakeRandom((s,s))
+            assert self.MatAbs(ident - m * m.Inverse()) < 1e-6, (
+                'offender = ' + `m`)
+
+    def RandomLUPTest(self,s,n):
+        ident = GenericMatrix(size=(s,s),fillMode='i')
+        for i in range(n):
+            m = self.MakeRandom((s,s))
+            (l,u,p) = m.LUP()
+            assert self.MatAbs(p*m - l*u) < 1e-6, 'offender = ' + `m`
+
+    def RandomSolveTest(self,s,n):
+        import random
+        if (s <= 1):
+            return
+        extraEquations=3
+        
+        for i in range(n):
+            m = self.MakeRandom((s,s+extraEquations))
+            for j in range(extraEquations):
+                colToKill = random.randrange(s+extraEquations)
+                for r in range(m.rows):
+                    m[r,colToKill] = 0.0
+            b = map(lambda x: random.random(), range(s))
+            x = m.Solve(b)
+            z = m.LeftMulColumnVec(x)
+            diff = reduce(lambda xx,yy:xx+yy, map(lambda aa,bb:abs(aa-bb),b,z))
+            assert diff < 1e-6, ('offenders: m = ' + `m` + '\nx = ' + `x`
+                                 + '\nb = ' + `b` + '\ndiff = ' + `diff`)
+
+    def RandomDetTest(self,s,n):
+        for i in range(n):
+            m1 = self.MakeRandom((s,s))
+            m2 = self.MakeRandom((s,s))
+            prod = m1 * m2
+            assert (abs(m1.Determinant() * m2.Determinant()
+                        - prod.Determinant() )
+                    < 1e-6), 'offenders = ' + `m1` + `m2`
+
+
+license_doc = """
+  This code was originally written by Emin Martinian (emin@allegro.mit.edu).
+  You may copy, modify, redistribute in source or binary form as long
+  as credit is given to the original author.  Specifically, please
+  include some kind of comment or docstring saying that Emin Martinian
+  was one of the original authors.  Also, if you publish anything based
+  on this work, it would be nice to cite the original author and any
+  other contributers.
+
+  There is NO WARRANTY for this software just as there is no warranty
+  for GNU software (although this is not GNU software).  Specifically
+  we adopt the same policy towards warranties as the GNU project:
+
+  BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM 'AS IS' WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+"""
+
+testing_doc = """
+The GenericMatrixTester class contains some simple
+testing functions such as RandomInverseTest, RandomLUPTest,
+RandomSolveTest, and RandomDetTest which generate random floating
+point values and test the appropriate routines.  The simplest way to
+run these tests is via
+
+>>> import genericmatrix
+>>> t = genericmatrix.GenericMatrixTester()
+>>> t.DoTests(100,[1,2,3,4,5,10])
+
+# runs 100 tests each for sizes 1-5, 10
+# note this may take a few minutes
+
+If any problems occur, assertion errors are raised.  Otherwise
+nothing is returned.  Note that you can also use the doctest
+package to test all the python examples in the documentation
+by typing 'python genericmatrix.py' or 'python -v genericmatrix.py' at the
+command line.
+"""
+
+
+# The following code is used to make the doctest package
+# check examples in docstrings when you enter
+
+__test__ = {
+    'testing_doc' : testing_doc
+}
+
+def _test():
+    import doctest, genericmatrix
+    return doctest.testmod(genericmatrix)
+
+if __name__ == "__main__":
+    _test()
+    print 'Tests Passed.'
diff --git a/src/py_ecc/rs_code.py b/src/py_ecc/rs_code.py
new file mode 100644 (file)
index 0000000..e070314
--- /dev/null
@@ -0,0 +1,246 @@
+
+# Copyright Emin Martinian 2002.  See below for license terms.
+# Version Control Info: $Id: rs_code.py,v 1.5 2003/07/04 01:30:05 emin Exp $
+
+"""
+This package implements the RSCode class designed to do
+Reed-Solomon encoding and (erasure) decoding.  The following
+docstrings provide detailed information on various topics.
+
+  RSCode.__doc__   Describes the RSCode class and how to use it.
+
+  license_doc      Describes the license and lack of warranty.
+
+"""
+
+import ffield
+import genericmatrix
+import math
+
+
+class RSCode:
+    """
+    The RSCode class implements a Reed-Solomon code
+    (currently, only erasure decoding not error decoding is
+    implemented).  The relevant methods are:
+
+    __init__
+    Encode
+    DecodeImmediate
+    Decode
+    PrepareDecoder
+    RandomTest
+
+    A breif example of how to use the code follows:
+    
+>>> import rs_code
+
+# Create a coder for an (n,k) = (16,8) code and test 
+# decoding for a simple erasure pattern.
+
+>>> C = rs_code.RSCode(16,8) 
+>>> inVec = range(8)         
+>>> codedVec = C.Encode(inVec)
+>>> receivedVec = list(codedVec)
+
+# now erase some entries in the encoded vector by setting them to None
+>>> receivedVec[3] = None; receivedVec[9] = None; receivedVec[12] = None
+>>> receivedVec
+[0, 1, 2, None, 4, 5, 6, 7, 8, None, 10, 11, None, 13, 14, 15]
+>>> decVec = C.DecodeImmediate(receivedVec)
+>>> decVec
+[0, 1, 2, 3, 4, 5, 6, 7]
+
+# Now try the random testing method for more complete coverage.
+# Note this will take a while.
+>>> for k in range(1,8):
+...     for p in range(1,12):
+...         C = rs_code.RSCode(k+p,k)
+...         C.RandomTest(25)
+>>> for k in range(1,8):
+...     for p in range(1,12):
+...         C = rs_code.RSCode(k+p,k,systematic=0)
+...         C.RandomTest(25)
+"""
+
+    def __init__(self,n,k,log2FieldSize=-1,systematic=1,shouldUseLUT=-1):
+        """
+        Function:   __init__(n,k,log2FieldSize,systematic,shouldUseLUT)
+        Purpose:    Create a Reed-Solomon coder for an (n,k) code.
+        Notes:      The last parameters, log2FieldSize, systematic
+                    and shouldUseLUT are optional.
+
+                    The log2FieldSize parameter 
+                    represents the base 2 logarithm of the field size.
+                    If it is omitted, the field GF(2^p) is used where
+                    p is the smalles integer where 2^p >= n.
+
+                    If systematic is true then a systematic encoder
+                    is created (i.e. one where the first k symbols
+                    of the encoded result always match the data).
+
+                    If shouldUseLUT = 1 then a lookup table is used for
+                    computing finite field multiplies and divides.
+                    If shouldUseLUT = 0 then no lookup table is used.
+                    If shouldUseLUT = -1 (the default), then the code
+                    decides when a lookup table should be used.
+        """
+        if (log2FieldSize < 0):
+            log2FieldSize = int(math.ceil(math.log(n)/math.log(2)))
+        self.field = ffield.FField(log2FieldSize,useLUT=shouldUseLUT)
+        self.n = n
+        self.k = k
+        self.fieldSize = 1 << log2FieldSize
+        self.CreateEncoderMatrix()
+        if (systematic):
+            self.encoderMatrix.Transpose()
+            self.encoderMatrix.LowerGaussianElim()
+            self.encoderMatrix.UpperInverse()
+            self.encoderMatrix.Transpose()
+
+    def __repr__(self):
+        rep = ('<RSCode (n,k) = (' + `self.n` +', ' + `self.k` + ')'
+               + '  over GF(2^' + `self.field.n` + ')\n' +
+               `self.encoderMatrix` + '\n' + '>')
+        return rep
+               
+    def CreateEncoderMatrix(self):                   
+        self.encoderMatrix = genericmatrix.GenericMatrix(
+            (self.n,self.k),0,1,self.field.Add,self.field.Subtract,
+            self.field.Multiply,self.field.Divide)
+        self.encoderMatrix[0,0] = 1
+        for i in range(0,self.n):
+            term = 1
+            for j in range(0, self.k):
+                self.encoderMatrix[i,j] = term
+                term = self.field.Multiply(term,i)
+
+    
+    def Encode(self,data):
+        """
+        Function:       Encode(data)
+        Purpose:        Encode a list of length k into length n.
+        """
+        assert len(data)==self.k, 'Encode: input data must be size k list.'
+        
+        return self.encoderMatrix.LeftMulColumnVec(data)
+
+    def PrepareDecoder(self,unErasedLocations):
+        """
+        Function:       PrepareDecoder(erasedTerms)
+        Description:    The input unErasedLocations is a list of the first
+                        self.k elements of the codeword which were 
+                        NOT erased.  For example, if the 0th, 5th,
+                        and 7th symbols of a (16,5) code were erased,
+                        then PrepareDecoder([1,2,3,4,6]) would
+                        properly prepare for decoding.
+        """
+        if (len(unErasedLocations) != self.k):
+            raise ValueError, 'input must be exactly length k'
+        
+        limitedEncoder = genericmatrix.GenericMatrix(
+            (self.k,self.k),0,1,self.field.Add,self.field.Subtract,
+            self.field.Multiply,self.field.Divide)
+        for i in range(0,self.k):
+            limitedEncoder.SetRow(
+                i,self.encoderMatrix.GetRow(unErasedLocations[i]))
+        self.decoderMatrix = limitedEncoder.Inverse()
+
+    def Decode(self,unErasedTerms):
+        """
+        Function:       Decode(unErasedTerms)
+        Purpose:        Use the
+        Description:
+        """
+        return self.decoderMatrix.LeftMulColumnVec(unErasedTerms)
+
+    def DecodeImmediate(self,data):
+        """
+        Function:       DecodeImmediate(data)
+        Description:    Takes as input a data vector of length self.n
+                        where erased symbols are set to None and
+                        returns the decoded result provided that
+                        at least self.k symbols are not None.
+
+                        For example, for an (n,k) = (6,4) code, a
+                        decodable input vector would be
+                        [2, 0, None, 1, 2, None].
+        """
+
+        if (len(data) != self.n):
+            raise ValueError, 'input must be a length n list'
+
+        unErasedLocations = []
+        unErasedTerms = []
+        for i in range(self.n):
+            if (None != data[i]):
+                unErasedLocations.append(i)
+                unErasedTerms.append(data[i])
+        self.PrepareDecoder(unErasedLocations[0:self.k])
+        return self.Decode(unErasedTerms[0:self.k])
+        
+    def RandomTest(self,numTests):
+        import random
+        
+        maxErasures = self.n-self.k
+        for i in range(numTests):
+            inVec = range(self.k)
+            for j in range(self.k):
+                inVec[j] = random.randint(0, (1<<self.field.n)-1)
+            codedVec = self.Encode(inVec)
+            numErasures = random.randint(0,maxErasures)
+            for j in range(numErasures):
+                j = random.randint(0,self.n-1)
+                while(codedVec[j] == None):
+                    j = random.randint(0,self.n-1)
+                codedVec[j] = None
+            decVec = self.DecodeImmediate(codedVec)
+            assert decVec == inVec, ('inVec = ' + `inVec`
+                                     + '\ncodedVec = ' + `codedVec`
+                                     + '\ndecVec = ' + `decVec`)
+
+license_doc = """
+  This code was originally written by Emin Martinian (emin@allegro.mit.edu).
+  You may copy, modify, redistribute in source or binary form as long
+  as credit is given to the original author.  Specifically, please
+  include some kind of comment or docstring saying that Emin Martinian
+  was one of the original authors.  Also, if you publish anything based
+  on this work, it would be nice to cite the original author and any
+  other contributers.
+
+  There is NO WARRANTY for this software just as there is no warranty
+  for GNU software (although this is not GNU software).  Specifically
+  we adopt the same policy towards warranties as the GNU project:
+
+  BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM 'AS IS' WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+"""
+
+
+# The following code is used to make the doctest package
+# check examples in docstrings.
+
+def _test():
+    import doctest, rs_code
+    return doctest.testmod(rs_code)
+
+if __name__ == "__main__":
+    _test()
+    print 'Tests passed'