fix: some bugs in numpy implementation of RS lib (still not enough, doesn't pass rstest.py).

Signed-off-by: Stephen L. <lrq3000@gmail.com>

Author: lrq3000

lib/brownanrs/nff.py

@@ -42,27 +42,46 @@
         137, 180, 124, 184, 38, 119, 153, 227, 165, 103, 74, 237, 222, 197,
         49, 254, 24, 13, 99, 140, 128, 192, 247, 112, 7]
 
-class nGF256int(int):
+class GF256int(int):
     """Instances of this object are elements of the field GF(2^8)
     Instances are integers in the range 0 to 255
     This field is defined using the irreducable polynomial
     x^8 + x^4 + x^3 + x + 1
     and using 3 as the generator for the exponent table and log table.
     """
+    
+    def __add__(a, b):
+        "Addition in GF(2^8) is the xor of the two"
+        return GF256int(a ^ b)
+    __sub__ = __add__
+    __radd__ = __add__
+    __rsub__ = __add__
+    def __neg__(self):
+        return self
+
+    def __mul__(a, b):
+        "Multiplication in GF(2^8)"
+        if a == 0 or b == 0:
+            return GF256int(0)
+        x = nGF256int_logtable[a]
+        y = nGF256int_logtable[b]
+        z = (x + y) % 255
+        return GF256int(nGF256int_exptable[z])
+    __rmul__ = __mul__
 
     def __pow__(self, power, modulo=None):
-        if isinstance(power, nGF256int):
-            raise TypeError("Raising a Field element to another Field element is not defined. power must be a regular integer")
+        #if isinstance(power, GF256int):
+            #raise TypeError("Raising a Field element to another Field element is not defined. power must be a regular integer")
         x = nGF256int_logtable[self]
         z = (x * power) % 255
-        return nGF256int(nGF256int_exptable[z])
+        return GF256int(nGF256int_exptable[z])
 
     def inverse(self):
         e = nGF256int_logtable[self]
-        return nGF256int(nGF256int_exptable[255 - e])
+        return GF256int(nGF256int_exptable[255 - e])
 
     def __div__(self, other):
-        return self * nGF256int(other).inverse()
+        return self * GF256int(other).inverse()
     def __rdiv__(self, other):
         return self.inverse() * other
 

lib/brownanrs/npolynomial.py

@@ -47,14 +47,42 @@
         137, 180, 124, 184, 38, 119, 153, 227, 165, 103, 74, 237, 222, 197,
         49, 254, 24, 13, 99, 140, 128, 192, 247, 112, 7], dtype=int )
 
-class nPolynomial(object):
+def GF256int_add(a, b):
+    "Addition in GF(2^8) is the xor of the two"
+    return a ^ b
+
+def GF256int_neg(a):
+    return a
+
+def GF256int_mul(a, b):
+    "Multiplication in GF(2^8)"
+    if a == 0 or b == 0:
+        return 0
+    x = GF256int_logtable[a]
+    y = GF256int_logtable[b]
+    z = (x + y) % 255
+    return GF256int_exptable[z]
+
+def GF256int_pow(a, power):
+    x = GF256int_logtable[a]
+    z = (x * power) % 255
+    return GF256int_exptable[z]
+
+def GF256int_inverse(a):
+    e = GF256int_logtable[a]
+    return GF256int_exptable[255 - e]
+
+def GF256int_div(a, b):
+    return GF256int_mul(a, GF256int_inverse(b))
+
+class Polynomial(object):
     """Polynomial for GF256int class with numpy implementation.
     This class implements vectorized operations on GF256int, thus the type is "defined" implitly in the operations involving GF256int_logtable and GF256int_exptable.
 
     Polynomial objects are immutable.
 
     Implementation note: this class is NOT type agnostic, contrariwise to polynomial.py and cpolynomial.py."""
-    def __init__(self, coefficients=[], **sparse):
+    def __init__(self, coefficients=None, **sparse):
         """
         There are three ways to initialize a Polynomial object.
         1) With a list, tuple, or other iterable, creates a polynomial using
@@ -80,16 +108,15 @@ def __init__(self, coefficients=[], **sparse):
                     " both")
         if coefficients is not None:
             # Polynomial( [1, 2, 3, ...] )
-            #c = coefficients
-            #if isinstance(coefficients, tuple): c = list(coefficients)
+            #if isinstance(coefficients, tuple): coefficients = list(coefficients)
             # Expunge any leading 0 coefficients
             while len(coefficients) > 0 and coefficients[0] == 0:
                 if isinstance(coefficients, np.ndarray):
                     coefficients = np.delete(coefficients, 0)
                 else:
                     coefficients.pop(0)
             if len(coefficients) == 0:
-                coefficients.append(0)
+                coefficients = [0]
 
             self.coefficients = np.array(coefficients, dtype=int)
         elif sparse:
@@ -109,16 +136,19 @@ def __init__(self, coefficients=[], **sparse):
         else:
             # Polynomial()
             self.coefficients = np.zeros(1, dtype=int)
-        self.degree = len(self.coefficients)
+        # In any case, compute the degree of the polynomial (=the maximum degree)
+        self.degree = len(self.coefficients)-1
 
     def __len__(self):
         """Returns the number of terms in the polynomial"""
-        return self.degree
+        return self.degree+1
+        # return len(self.coefficients)
 
     def degree(self, poly=None):
         """Returns the degree of the polynomial"""
         if not poly:
-            return self.degree - 1
+            return self.degree
+            #return len(self.coefficients) - 1
         elif poly and hasattr("coefficients", poly):
             return len(poly.coefficients) - 1
         else:
@@ -127,7 +157,7 @@ def degree(self, poly=None):
             return len(poly)-1
 
     def __add__(self, other):
-        diff = self.degree - other.degree
+        diff = len(self) - len(other)
 
         # Warning: value must be absolute, there's no negative value accepted here! (in fact they will be accepted but the result will be meaningless! and divmod will be in infinite loop)
         t1 = np.pad(np.absolute(self.coefficients), (np.max([0, -diff]),0), mode='constant')
@@ -218,21 +248,26 @@ def __divmod__(dividend, divisor):
             remainder = dividend
             remainder_power = dividend_power
             remainder_coefficient = dividend_coefficient
-            while remainder_power >= divisor_power: # Until there's no remainder left (or the remainder cannot be divided anymore by the divisor)
-                quotient_power = remainder_power - divisor_power
-                quotient_coefficient = remainder_coefficient / divisor_coefficient
+            quotient_power = remainder_power - divisor_power # need to set at least 1 just to start the loop. Warning if set to remainder_power - divisor_power: because it may skip the loop altogether (and we want to at least do one iteration to set the quotient)
+
+            # Compute how many times the highest order term in the divisor goes into the dividend
+            while quotient_power >= 0 and np.count_nonzero(remainder.coefficients) > 0: # Until there's no remainder left (or the remainder cannot be divided anymore by the divisor)
+                quotient_coefficient = GF256int_div(remainder_coefficient, divisor_coefficient)
                 q = class_( np.pad([quotient_coefficient], (0,quotient_power), mode='constant') ) # construct an array with only the quotient major coefficient (we divide the remainder only with the major coeff)
                 quotient = quotient + q # add the coeff to the full quotient
                 remainder = remainder - q * divisor # divide the remainder with the major coeff quotient multiplied by the divisor, this gives us the new remainder
                 remainder_power = remainder.degree # compute the new remainder degree
                 remainder_coefficient = remainder[0] # Compute the new remainder coefficient
-                #print "quotient: %s remainder: %s" % (quotient, remainder)
+                quotient_power = remainder_power - divisor_power
         return quotient, remainder
 
     def __eq__(self, other):
-        return self.coefficients == other.coefficients
+        if isinstance(self.coefficients, np.ndarray):
+            return self.coefficients.tolist() == other.coefficients.tolist()
+        else:
+            return self.coefficients == other.coefficients
     def __ne__(self, other):
-        return self.coefficients != other.coefficients
+        return not self == other
     def __hash__(self):
         return hash(self.coefficients)
 
@@ -258,7 +293,7 @@ def __str__(self):
         return buf.getvalue()[:-3]
 
     def evaluate(self, x):
-        "Evaluate this polynomial at value x, returning the result."
+        '''Evaluate this polynomial at value x, returning the result.'''
         # Holds the sum over each term in the polynomial
         c = 0
 
@@ -273,7 +308,7 @@ def evaluate(self, x):
         return c
 
     def get_coefficient(self, degree):
-        """Returns the coefficient of the specified term"""
+        '''Returns the coefficient of the specified term'''
         if degree > self.degree:
             return 0
         else:

lib/brownanrs/npolynomialtest.py

@@ -0,0 +1,155 @@
+import unittest
+
+from npolynomial import Polynomial
+from nff import GF256int
+
+class TestGFPoly(unittest.TestCase):
+    """Tests that the Polynomial class works when given GF256int objects
+    instead of regular integers
+    """
+    def test_add(self):
+        one = Polynomial(map(GF256int,     [8,3,5,1]))
+        two = Polynomial(map(GF256int, [5,3,1,1,6,8]))
+
+        r = one + two
+
+        self.assertEqual(list(r.coefficients), [5,3,9,2,3,9])
+
+    def test_sub(self):
+        one = Polynomial(map(GF256int,     [8,3,5,1]))
+        two = Polynomial(map(GF256int, [5,3,1,1,6,8]))
+        r = one - two
+        self.assertEqual(list(r.coefficients), [5,3,9,2,3,9])
+
+    def test_mul(self):
+        one = Polynomial(map(GF256int,     [8,3,5,1]))
+        two = Polynomial(map(GF256int, [5,3,1,1,6,8]))
+        r = one * two
+        self.assertEqual(list(r.coefficients), [40,23,28,1,53,78,7,46,8])
+
+    def test_div(self):
+        one = Polynomial(map(GF256int,     [8,3,5,1]))
+        two = Polynomial(map(GF256int, [5,3,1,1,6,8]))
+        q, r = divmod(two,one)
+        self.assertEqual(list(q.coefficients), [101, 152, 11])
+        self.assertEqual(list(r.coefficients), [183, 185, 3])
+
+        # Make sure they multiply back out okay
+        self.assertEqual(q*one + r, two)
+
+    def test_div_scalar(self):
+        """Tests division by a scalar"""
+        numbers = map(GF256int, [5,20,50,100,134,158,0,148,233,254,4,5,2])
+        scalar = GF256int(17)
+
+        poly = Polynomial(numbers)
+        scalarpoly = Polynomial(x0=scalar)
+
+        self.assertEqual(
+                list((poly // scalarpoly).coefficients),
+                map(lambda x: x / scalar, numbers)
+                )
+
+    def test_div_scalar2(self):
+        """Test that dividing by a scalar is the same as multiplying by the
+        scalar's inverse"""
+        a = Polynomial(map(GF256int, [5,3,1,1,6,8]))
+
+        scalar = GF256int(50)
+
+        self.assertEqual(
+                a * Polynomial(x0=scalar),
+                a // Polynomial(x0=scalar.inverse())
+                )
+
+
+
+# class TestPolynomial(unittest.TestCase):
+    # def test_add_1(self):
+        # one = Polynomial([2,4,7,3])
+        # two = Polynomial([5,2,4,2])
+
+        # r = one + two
+
+        # self.assertEqual(list(r.coefficients), [7, 6, 11, 5])
+
+    # def test_add_2(self):
+        # one = Polynomial([2,4,7,3,5,2])
+        # two = Polynomial([5,2,4,2])
+
+        # r = one + two
+
+        # self.assertEqual(list(r.coefficients), [2,4,12,5,9,4])
+
+    # def test_add_3(self):
+        # one = Polynomial([7,3,5,2])
+        # two = Polynomial([6,8,5,2,4,2])
+
+        # r = one + two
+
+        # self.assertEqual(list(r.coefficients), [6,8,12,5,9,4])
+
+    # def test_mul_1(self):
+        # one = Polynomial([2,4,7,3])
+        # two = Polynomial([5,2,4,2])
+
+        # r = one._mul_standard(two)
+
+        # self.assertEqual(list(r.coefficients),
+                # [10,24,51,49,42,26,6])
+
+    # def test_div_1(self):
+        # one = Polynomial([1,4,0,3])
+        # two = Polynomial([1,0,1])
+
+        # q, r = divmod(one, two)
+        # self.assertEqual(q, one // two)
+        # self.assertEqual(r, one % two)
+
+        # self.assertEqual(list(q.coefficients), [1,4])
+        # self.assertEqual(list(r.coefficients), [-1,-1])
+
+    # def test_div_2(self):
+        # one = Polynomial([1,0,0,2,2,0,1,2,1])
+        # two = Polynomial([1,0,-1])
+
+        # q, r = divmod(one, two)
+        # self.assertEqual(q, one // two)
+        # self.assertEqual(r, one % two)
+
+        # self.assertEqual(list(q.coefficients), [1,0,1,2,3,2,4])
+        # self.assertEqual(list(r.coefficients), [4,5])
+
+    # def test_div_3(self):
+        # # 0 quotient
+        # one = Polynomial([1,0,-1])
+        # two = Polynomial([1,1,0,0,-1])
+
+        # q, r = divmod(one, two)
+        # self.assertEqual(q, one // two)
+        # self.assertEqual(r, one % two)
+
+        # self.assertEqual(list(q.coefficients), [0])
+        # self.assertEqual(list(r.coefficients), [1,0,-1])
+
+    # def test_div_4(self):
+        # # no remander
+        # one = Polynomial([1,0,0,2,2,0,1,-2,-4])
+        # two = Polynomial([1,0,-1])
+
+        # q, r = divmod(one, two)
+        # self.assertEqual(q, one // two)
+        # self.assertEqual(r, one % two)
+
+        # self.assertEqual(list(q.coefficients), [1,0,1,2,3,2,4])
+        # self.assertEqual(list(r.coefficients), [0])
+
+    # def test_getcoeff(self):
+        # p = Polynomial([9,3,3,2,2,3,1,-2,-4])
+        # self.assertEqual(p.get_coefficient(0), -4)
+        # self.assertEqual(p.get_coefficient(2), 1)
+        # self.assertEqual(p.get_coefficient(8), 9)
+        # self.assertEqual(p.get_coefficient(9), 0)
+
+if __name__ == "__main__":
+    unittest.main()

lib/brownanrs/rs.py

@@ -3,9 +3,9 @@
 # See LICENSE.txt for license terms
 
 # try: # Numpy implementation import. NOTE: it's working correctly but it's actually slower than the pure-python implementation! (although the code is vectorized and __mul__ is faster, but it seems it's not enough to speed things up and the numpy call overhead is too much for any gain here).
-    # from npolynomial import nPolynomial as Polynomial
-    # from nff import nGF256int as GF256int
-#except ImportError:
+    # from npolynomial import Polynomial
+    # from nff import GF256int
+# except ImportError:
 try: # Cython implementation import. This should be a bit faster than using PyPy with the pure-python implementation.
     from cff import GF256int
     from cpolynomial import Polynomial