first commit

.gitignore

@@ -0,0 +1 @@
+*.pyc__pycache__

README.md

@@ -0,0 +1,4 @@
+welch-berlekamp
+===============
+
+Python code and tests for the post ["The Welch-Berlekamp Algorithm for Correcting Errors in Data"](http://jeremykun.com/)

example.py

@@ -0,0 +1,63 @@
+
+from reedsolomon import makeEncoderDecoder
+import random
+
+def simpleTest():
+   integerMessage = [2,3,2]
+   k = 3
+   n = 7
+   p = 11
+
+   enc, dec, solveSystem = makeEncoderDecoder(n, k, p)
+   encoded = enc(integerMessage)
+
+   print("plain message is: %r" % (integerMessage,))
+   print("encoded message is: %r" % (encoded,)) # cleaner output
+
+   e=1
+   corrupted = encoded[:]
+   corrupted[n//2][1] = corrupted[n//2][1] + 1
+   print("corrupted message is: %r" % (corrupted,))
+
+   Q,E = solveSystem(corrupted)
+   P, remainder = (Q.__divmod__(E))
+
+   print("P(x) = %r" % P)
+   print("r(x) = %r" % remainder)
+
+
+def corrupt(message, numErrors, minVal=0, maxVal=131):
+   indices = random.sample(list(range(len(message))), numErrors)
+
+   for i in indices:
+      message[i][1] = random.randint(minVal, maxVal)
+
+   return message
+
+
+def tkamTest():
+   message = '''When he was nearly thirteen, my brother Jem got his arm badly broken at the elbow. When it healed, and Jem's fears of never being able to play football were assuaged, he was seldom self-conscious about his injury. His left arm was somewhat shorter than his right; when he stood or walked, the back of his hand was at right angles to his body, his thumb parallel to his thigh. He couldn't have cared less, so long as he could pass and punt.
+
+When enough years had gone by to enable us to look back on them, we sometimes discussed the events leading to his accident. I maintain that the Ewells started it all, but Jem, who was four years my senior, said it started long before that. He said it began the summer Dill came to us, when Dill first gave us the idea of making Boo Radley come out.'''
+
+   k = len(message)
+   n = len(message) * 2
+   p = 2087
+   integerMessage = [ord(x) for x in message]
+
+   enc, dec, solveSystem = makeEncoderDecoder(n, k, p)
+   print("encoding...")
+   encoded = enc(integerMessage)
+
+   e = int(k/2)
+   print("corrupting...")
+   corrupted = corrupt(encoded[:], e, 0, p)
+
+   print("decoding...")
+   Q,E = solveSystem(corrupted)
+   P, remainder = (Q.__divmod__(E))
+
+   recovered = [chr(x) for x in P.coefficients]
+   print(recovered)
+
+tkamTest()

example.txt

@@ -0,0 +1,55 @@
+plain message is: [2, 3, 2]
+encoded message is: [[0, 2], [1, 0], [2, 2], [3, 1], [4, 4]]
+corrupted message is: [[0, 2], [1, 0], [2, 3], [3, 1], [4, 4]]
+
+e is 2
+
+system is:
+
+
+	[2, 0, 0, 6, 0, 0, 0, 0, 0]
+	[0, 0, 0, 6, 6, 6, 6, 6, 0]
+	[3, 6, 5, 6, 5, 3, 6, 5, 0]
+	[1, 3, 2, 6, 4, 5, 1, 3, 0]
+	[4, 2, 1, 6, 3, 5, 6, 3, 0]
+	[0, 0, 1, 0, 0, 0, 0, 0, 1]
+
+reduced system is:
+
+
+	[1, 0, 0, 0, 0, 0, 1, 4, 0]
+	[0, 1, 0, 0, 0, 0, 3, 3, 1]
+	[0, 0, 1, 0, 0, 0, 0, 0, 1]
+	[0, 0, 0, 1, 0, 0, 2, 1, 0]
+	[0, 0, 0, 0, 1, 0, 2, 2, 5]
+	[0, 0, 0, 0, 0, 1, 4, 5, 2]
+solution is [2, 2, 1, 4, 1, 0, 1, 1]
+Q is 4 + 1 t^1 + 0 t^2 + 1 t^3 + 1 t^4
+E is 2 + 2 t^1 + 1 t^2
+
+e is 1
+
+system is:
+
+
+	[2, 0, 6, 0, 0, 0, 0]
+	[0, 0, 6, 6, 6, 6, 0]
+	[3, 6, 6, 5, 3, 6, 0]
+	[1, 3, 6, 4, 5, 1, 0]
+	[4, 2, 6, 3, 5, 6, 0]
+	[0, 1, 0, 0, 0, 0, 1]
+
+reduced system is:
+
+
+	[1, 0, 0, 0, 0, 0, 5]
+	[0, 1, 0, 0, 0, 0, 1]
+	[0, 0, 1, 0, 0, 0, 3]
+	[0, 0, 0, 1, 0, 0, 3]
+	[0, 0, 0, 0, 1, 0, 6]
+	[0, 0, 0, 0, 0, 1, 2]
+solution is [5, 1, 3, 3, 6, 2]
+Q is 3 + 3 t^1 + 6 t^2 + 2 t^3
+E is 5 + 1 t^1
+P(x) = 2 + 3 t^1 + 2 t^2
+r(x) = 0

finitefield/euclidean.py

@@ -0,0 +1,35 @@
+
+# a general Euclidean algorithm for any number type with
+# a divmod and a valuation abs() whose minimum value is zero
+def gcd(a, b):
+   if abs(a) < abs(b):
+      return gcd(b, a)
+
+   while abs(b) > 0:
+      _,r = divmod(a,b)
+      a,b = b,r
+
+   return a
+
+
+# extendedEuclideanAlgorithm: int, int -> int, int, int
+# input (a,b) and output three numbers x,y,d such that ax + by = d = gcd(a,b).
+# Works for any number type with a divmod and a valuation abs()
+# whose minimum value is zero
+def extendedEuclideanAlgorithm(a, b):
+   if abs(b) > abs(a):
+      (x,y,d) = extendedEuclideanAlgorithm(b, a)
+      return (y,x,d)
+
+   if abs(b) == 0:
+      return (1, 0, a)
+
+   x1, x2, y1, y2 = 0, 1, 1, 0
+   while abs(b) > 0:
+      q, r = divmod(a,b)
+      x = x2 - q*x1
+      y = y2 - q*y1
+      a, b, x2, x1, y2, y1 = b, r, x1, x, y1, y
+
+   return (x2, y2, a)
+

finitefield/finitefield.py

@@ -0,0 +1,118 @@
+import random
+from .polynomial import polynomialsOver
+from .modp import *
+
+
+
+# isIrreducible: Polynomial, int -> bool
+# determine if the given monic polynomial with coefficients in Z/p is
+# irreducible over Z/p where p is the given integer
+# Algorithm 4.69 in the Handbook of Applied Cryptography
+def isIrreducible(polynomial, p):
+   ZmodP = IntegersModP(p)
+   if polynomial.field is not ZmodP:
+      raise TypeError("Given a polynomial that's not over %s, but instead %r" %
+                        (ZmodP.__name__, polynomial.field.__name__))
+
+   poly = polynomialsOver(ZmodP).factory
+   x = poly([0,1])
+   powerTerm = x
+   isUnit = lambda p: p.degree() == 0
+
+   for _ in range(int(polynomial.degree() / 2)):
+      powerTerm = powerTerm.powmod(p, polynomial)
+      gcdOverZmodp = gcd(polynomial, powerTerm - x)
+      if not isUnit(gcdOverZmodp):
+         return False
+
+   return True
+
+
+# generateIrreduciblePolynomial: int, int -> Polynomial
+# generate a random irreducible polynomial of a given degree over Z/p, where p
+# is given by the integer 'modulus'. This algorithm is expected to terminate
+# after 'degree' many irreducilibity tests. By Chernoff bounds the probability
+# it deviates from this by very much is exponentially small.
+def generateIrreduciblePolynomial(modulus, degree):
+   Zp = IntegersModP(modulus)
+   Polynomial = polynomialsOver(Zp)
+
+   while True:
+      coefficients = [Zp(random.randint(0, modulus-1)) for _ in range(degree)]
+      randomMonicPolynomial = Polynomial(coefficients + [Zp(1)])
+      print(randomMonicPolynomial)
+
+      if isIrreducible(randomMonicPolynomial, modulus):
+         return randomMonicPolynomial
+
+
+# create a type constructor for the finite field of order p^m for p prime, m >= 1
+@memoize
+def FiniteField(p, m=1, polynomialModulus=None):
+   Zp = IntegersModP(p)
+   if m == 1:
+      return Zp
+
+   Polynomial = polynomialsOver(Zp)
+   if polynomialModulus is None:
+      polynomialModulus = generateIrreduciblePolynomial(modulus=p, degree=m)
+
+   class Fq(FieldElement):
+      fieldSize = int(p ** m)
+      primeSubfield = Zp
+      idealGenerator = polynomialModulus
+      operatorPrecedence = 3
+
+      def __init__(self, poly):
+         if type(poly) is Fq:
+            self.poly = poly.poly
+         elif type(poly) is int or type(poly) is Zp:
+            self.poly = Polynomial([Zp(poly)])
+         elif isinstance(poly, Polynomial):
+            self.poly = poly % polynomialModulus
+         else:
+            self.poly = Polynomial([Zp(x) for x in poly]) % polynomialModulus
+
+         self.field = Fq
+
+      @typecheck
+      def __add__(self, other): return Fq(self.poly + other.poly)
+      @typecheck
+      def __sub__(self, other): return Fq(self.poly - other.poly)
+      @typecheck
+      def __mul__(self, other): return Fq(self.poly * other.poly)
+      @typecheck
+      def __eq__(self, other): return isinstance(other, Fq) and self.poly == other.poly
+
+      def __pow__(self, n): return Fq(pow(self.poly, n))
+      def __neg__(self): return Fq(-self.poly)
+      def __abs__(self): return abs(self.poly)
+      def __repr__(self): return repr(self.poly) + ' \u2208 ' + self.__class__.__name__
+
+      @typecheck
+      def __divmod__(self, divisor):
+         q,r = divmod(self.poly, divisor.poly)
+         return (Fq(q), Fq(r))
+
+
+      def inverse(self):
+         if self == Fq(0):
+            raise ZeroDivisionError
+
+         x,y,d = extendedEuclideanAlgorithm(self.poly, self.idealGenerator)
+         if d.degree() != 0:
+            raise Exception('Somehow, this element has no inverse! Maybe intialized with a non-prime?')
+
+         return Fq(x) * Fq(d.coefficients[0].inverse())
+
+
+   Fq.__name__ = 'F_{%d^%d}' % (p,m)
+   return Fq
+
+
+if __name__ == "__main__":
+   F23 = FiniteField(2,3)
+   x = F23([1,1])
+
+   F35 = FiniteField(3,5)
+   y = F35([1,1,2])

finitefield/modp.py

@@ -0,0 +1,80 @@
+
+from .euclidean import *
+from .numbertype import *
+
+# so all IntegersModP are instances of the same base class
+class _Modular(FieldElement):
+   pass
+
+
+@memoize
+def IntegersModP(p):
+   # assume p is prime
+
+   class IntegerModP(_Modular):
+      def __init__(self, n):
+         try:
+            self.n = int(n) % IntegerModP.p
+         except:
+            raise TypeError("Can't cast type %s to %s in __init__" % (type(n).__name__, type(self).__name__))
+
+         self.field = IntegerModP
+
+      @typecheck
+      def __add__(self, other):
+         return IntegerModP(self.n + other.n)
+
+      @typecheck
+      def __sub__(self, other):
+         return IntegerModP(self.n - other.n)
+
+      @typecheck
+      def __mul__(self, other):
+         return IntegerModP(self.n * other.n)
+
+      def __neg__(self):
+         return IntegerModP(-self.n)
+
+      @typecheck
+      def __eq__(self, other):
+         return isinstance(other, IntegerModP) and self.n == other.n
+
+      @typecheck
+      def __ne__(self, other):
+         return isinstance(other, IntegerModP) is False or self.n != other.n
+
+      @typecheck
+      def __divmod__(self, divisor):
+         q,r = divmod(self.n, divisor.n)
+         return (IntegerModP(q), IntegerModP(r))
+
+      def inverse(self):
+         # need to use the division algorithm *as integers* because we're
+         # doing it on the modulus itself (which would otherwise be zero)
+         x,y,d = extendedEuclideanAlgorithm(self.n, self.p)
+
+         if d != 1:
+            raise Exception("Error: p is not prime in %s!" % (self.__name__))
+
+         return IntegerModP(x)
+
+      def __abs__(self):
+         return abs(self.n)
+
+      def __str__(self):
+         return str(self.n)
+
+      def __repr__(self):
+         return '%d (mod %d)' % (self.n, self.p)
+
+      def __int__(self):
+         return self.n
+
+   IntegerModP.p = p
+   IntegerModP.__name__ = 'Z/%d' % (p)
+   IntegerModP.englishName = 'IntegersMod%d' % (p)
+   return IntegerModP
+
+
+if __name__ == "__main__":
+   mod7 = IntegersModP(7)