ExternalVendorCode/tahoe-lafs
1
0
Fork 0
mirror of https://github.com/tahoe-lafs/tahoe-lafs.git synced 2025-01-18 10:46:24 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

import py_ecc, a pure python fec by Emin Martinian, which is under a permissive licence

Browse Source 
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...
This commit is contained in:
Zooko O'Whielacronx 2006-12-29 13:46:45 -07:00
parent dcfddcac41
commit 6eaaf1fe9d
4 changed files with 2097 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

764
src/py_ecc/ffield.py Normal file
View File

@ -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.'

218
src/py_ecc/file_ecc.py Normal file
View File

@ -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'

869
src/py_ecc/genericmatrix.py Normal file
View File

@ -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.'

246
src/py_ecc/rs_code.py Normal file
View File

@ -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'