From: Zooko O'Whielacronx Date: Fri, 29 Dec 2006 20:46:45 +0000 (-0700) Subject: import py_ecc, a pure python fec by Emin Martinian, which is under a permissive licence X-Git-Tag: tahoe_v0.1.0-0-UNSTABLE~418 X-Git-Url: https://git.rkrishnan.org/pf/content/en.html?a=commitdiff_plain;h=6eaaf1fe9d2cdaa9a9dd39ca5867675bdfb1deea;p=tahoe-lafs%2Ftahoe-lafs.git import py_ecc, a pure python fec by Emin Martinian, which is under a permissive licence 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... --- diff --git a/src/py_ecc/ffield.py b/src/py_ecc/ffield.py new file mode 100644 index 00000000..a17a3594 --- /dev/null +++ b/src/py_ecc/ffield.py @@ -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<= 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 index 00000000..66d9b65b --- /dev/null +++ b/src/py_ecc/file_ecc.py @@ -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 index 00000000..1b38b4f0 --- /dev/null +++ b/src/py_ecc/genericmatrix.py @@ -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 + +>>> vi = v.Inverse() +>>> vi + +>>> (vi * v) - v.MakeSimilarMatrix(v.Size(),'i') + + +# See what happens when we try to invert a non-invertible matrix + +>>> v[0,1] = 0.0 +>>> v + +>>> 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 + +>>> u + +>>> p + +>>> p * v - l * u + + +# 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 + +>>> 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 = '' + 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 + +>>> 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 index 00000000..e0703143 --- /dev/null +++ b/src/py_ecc/rs_code.py @@ -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 = ('') + 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<