bls scalar field

scripts/finite_fields/README.md

@@ -0,0 +1,4 @@
+finite-fields
+=============
+
+Python code and tests for the post ["Programming with Finite Fields"](http://jeremykun.com/2014/03/13/programming-with-finite-fields/)

scripts/finite_fields/euclidean-test.py

@@ -0,0 +1,45 @@
+from test import test
+from euclidean import *
+
+test(1, gcd(7, 9))
+test(2, gcd(8, 18))
+test(-12, gcd(-12, 24))
+test(12, gcd(12, -24)) # gcd is only unique up to multiplication by a unit, and so sometimes we'll get negatives.
+test(38, gcd(4864, 3458))
+
+test((32, -45, 38), extendedEuclideanAlgorithm(4864, 3458))
+test((-45, 32, 38), extendedEuclideanAlgorithm(3458, 4864))
+
+from modp import *
+
+Mod2 = IntegersModP(2)
+test(Mod2(1), gcd(Mod2(1), Mod2(0)))
+test(Mod2(1), gcd(Mod2(1), Mod2(1)))
+test(Mod2(0), gcd(Mod2(2), Mod2(2)))
+
+Mod7 = IntegersModP(7)
+test(Mod7(6), gcd(Mod7(6), Mod7(14)))
+test(Mod7(2), gcd(Mod7(6), Mod7(9)))
+
+ModHuge = IntegersModP(9923)
+test(ModHuge(38), gcd(ModHuge(4864), ModHuge(3458)))
+test((ModHuge(32), ModHuge(-45), ModHuge(38)),
+     extendedEuclideanAlgorithm(ModHuge(4864), ModHuge(3458)))
+
+from polynomial import *
+
+p = polynomialsOver(Mod7).factory
+test(p([-1, 1]), gcd(p([-1,0,1]), p([-1,0,0,1])))
+f = p([-1,0,1])
+g = p([-1,0,0,1])
+test((p([0,-1]), p([1]), p([-1, 1])), extendedEuclideanAlgorithm(f, g))
+test(p([-1,1]), f * p([0,-1]) + g * p([1]))
+
+p = polynomialsOver(Mod2).factory
+f = p([1,0,0,0,1,1,1,0,1,1,1]) # x^10 + x^9 + x^8 + x^6 + x^5 + x^4 + 1
+g = p([1,0,1,1,0,1,1,0,0,1])   # x^9 + x^6 + x^5 + x^3 + x^1 + 1
+theGcd = p([1,1,0,1]) # x^3 + x + 1
+x = p([0,0,0,0,1]) # x^4
+y = p([1,1,1,1,1,1]) # x^5 + x^4 + x^3 + x^2 + x + 1
+
+test((x, y, theGcd), extendedEuclideanAlgorithm(f, g))

scripts/finite_fields/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)
+

scripts/finite_fields/finitefield-test.py

@@ -0,0 +1,28 @@
+from test import test
+from finitefield import *
+from polynomial import *
+from modp import *
+
+def p(L, q):
+   f = IntegersModP(q)
+   Polynomial = polynomialsOver(f).factory
+   return Polynomial(L)
+
+test(True, isIrreducible(p([0,1], 2), 2))
+test(False, isIrreducible(p([1,0,1], 2), 2))
+test(True, isIrreducible(p([1,0,1], 3), 3))
+
+test(False, isIrreducible(p([1,0,0,1], 5), 5))
+test(False, isIrreducible(p([1,0,0,1], 7), 7))
+test(False, isIrreducible(p([1,0,0,1], 11), 11))
+
+
+test(True, isIrreducible(p([-2, 0, 1], 13), 13))
+
+
+Z5 = IntegersModP(5)
+Poly = polynomialsOver(Z5).factory
+f = Poly([3,0,1])
+F25 = FiniteField(5, 2, polynomialModulus=f)
+x = F25([2,1])
+test(Poly([1,2]), x.inverse())

scripts/finite_fields/finitefield.py

@@ -0,0 +1,128 @@
+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, 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
+      @typecheck
+      def __ne__(self, other): return not self == other
+      
+      def __pow__(self, n):
+         if n==0: return Fq([1])
+         if n==1: return self
+         if n%2==0:
+            sqrut = self**(n//2)
+            return sqrut*sqrut
+         if n%2==1: return (self**(n-1))*self
+      
+      #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])

scripts/finite_fields/modp-test.py

@@ -0,0 +1,13 @@
+from modp import *
+from test import test
+
+mod7 = IntegersModP(7)
+
+test(mod7(5), mod7(5)) # Sanity check
+test(mod7(5), 1 / mod7(3))
+test(mod7(1), mod7(3) * mod7(5)) 
+test(mod7(3), mod7(3) * 1)
+test(mod7(2), mod7(5) + mod7(4))
+
+test(True, mod7(0) == mod7(3) + mod7(4))
+

scripts/finite_fields/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)

scripts/finite_fields/numbertype.py

@@ -0,0 +1,98 @@
+# memoize calls to the class constructors for fields
+# this helps typechecking by never creating two separate
+# instances of a number class.
+def memoize(f):
+   cache = {}
+
+   def memoizedFunction(*args, **kwargs):
+      argTuple = args + tuple(kwargs)
+      if argTuple not in cache:
+         cache[argTuple] = f(*args, **kwargs)
+      return cache[argTuple]
+
+   memoizedFunction.cache = cache
+   return memoizedFunction
+
+
+# type check a binary operation, and silently typecast 0 or 1
+def typecheck(f):
+   def newF(self, other):
+      if (hasattr(other.__class__, 'operatorPrecedence') and
+            other.__class__.operatorPrecedence > self.__class__.operatorPrecedence):
+         return NotImplemented
+
+      if type(self) is not type(other):
+         try:
+            other = self.__class__(other)
+         except TypeError:
+            message = 'Not able to typecast %s of type %s to type %s in function %s'
+            raise TypeError(message % (other, type(other).__name__, type(self).__name__, f.__name__))
+         except Exception as e:
+            message = 'Type error on arguments %r, %r for functon %s. Reason:%s'
+            raise TypeError(message % (self, other, f.__name__, type(other).__name__, type(self).__name__, e))
+
+      return f(self, other)
+
+   return newF
+
+
+
+# require a subclass to implement +-* neg and to perform typechecks on all of
+# the binary operations finally, the __init__ must operate when given a single
+# argument, provided that argument is the int zero or one
+class DomainElement(object):
+   operatorPrecedence = 1
+
+   # the 'r'-operators are only used when typecasting ints
+   def __radd__(self, other): return self + other
+   def __rsub__(self, other): return -self + other
+   def __rmul__(self, other): return self * other
+
+   # square-and-multiply algorithm for fast exponentiation
+   def __pow__(self, n):
+      if type(n) is not int:
+         raise TypeError
+
+      Q = self
+      R = self if n & 1 else self.__class__(1)
+
+      i = 2
+      while i <= n:
+         Q = (Q * Q)
+
+         if n & i == i:
+            R = (Q * R)
+
+         i = i << 1
+
+      return R
+
+
+   # requires the additional % operator (i.e. a Euclidean Domain)
+   def powmod(self, n, modulus):
+      if type(n) is not int:
+         raise TypeError
+
+      Q = self
+      R = self if n & 1 else self.__class__(1)
+
+      i = 2
+      while i <= n:
+         Q = (Q * Q) % modulus
+
+         if n & i == i:
+            R = (Q * R) % modulus
+
+         i = i << 1
+
+      return R
+
+
+
+# additionally require inverse() on subclasses
+class FieldElement(DomainElement):
+   def __truediv__(self, other): return self * other.inverse()
+   def __rtruediv__(self, other): return self.inverse() * other
+   def __div__(self, other): return self.__truediv__(other)
+   def __rdiv__(self, other): return self.__rtruediv__(other)
+

scripts/finite_fields/polynomial-test.py

@@ -0,0 +1,51 @@
+from __future__ import division
+from test import test
+from fractions import Fraction
+from polynomial import *
+
+from modp import *
+
+Mod5 = IntegersModP(5)
+Mod11 = IntegersModP(11)
+
+polysOverQ = polynomialsOver(Fraction).factory
+polysMod5 = polynomialsOver(Mod5).factory
+polysMod11 = polynomialsOver(Mod11).factory
+
+for p in [polysOverQ, polysMod5, polysMod11]:
+   # equality
+   test(True, p([]) == p([]))
+   test(True, p([1,2]) == p([1,2]))
+   test(True, p([1,2,0]) == p([1,2,0,0]))
+
+   # addition
+   test(p([1,2,3]), p([1,0,3]) + p([0,2]))
+   test(p([1,2,3]), p([1,2,3]) + p([]))
+   test(p([5,2,3]), p([4]) + p([1,2,3]))
+   test(p([1,2]), p([1,2,3]) + p([0,0,-3]))
+
+   # subtraction
+   test(p([1,-2,3]), p([1,0,3]) - p([0,2]))
+   test(p([1,2,3]), p([1,2,3]) - p([]))
+   test(p([-1,-2,-3]), p([]) - p([1,2,3]))
+
+   # multiplication
+   test(p([1,2,1]), p([1,1]) * p([1,1]))
+   test(p([2,5,5,3]), p([2,3]) * p([1,1,1]))
+   test(p([0,7,49]), p([0,1,7]) * p([7]))
+
+   # division
+   test(p([1,1,1,1,1,1]), p([-1,0,0,0,0,0,1]) / p([-1,1]))
+   test(p([-1,1,-1,1,-1,1]), p([1,0,0,0,0,0,1]) / p([1,1]))
+   test(p([]), p([]) / p([1,1]))
+   test(p([1,1]), p([1,1]) / p([1]))
+   test(p([1,1]), p([2,2]) / p([2]))
+
+   # modulus
+   test(p([]), p([1,7,49]) % p([7]))
+   test(p([-7]), p([-3,10,-5,3]) % p([1,3]))
+
+
+test(polysOverQ([Fraction(1,7), 1, 7]), polysOverQ([1,7,49]) / polysOverQ([7]))
+test(polysMod5([1 / Mod5(7), 1, 7]), polysMod5([1,7,49]) / polysMod5([7]))
+test(polysMod11([1 / Mod11(7), 1, 7]), polysMod11([1,7,49]) / polysMod11([7]))

scripts/finite_fields/polynomial.py

@@ -0,0 +1,133 @@
+try:
+    from itertools import zip_longest
+except ImportError:
+    from itertools import izip_longest as zip_longest
+import fractions
+
+from numbertype import *
+
+# strip all copies of elt from the end of the list
+def strip(L, elt):
+   if len(L) == 0: return L
+
+   i = len(L) - 1
+   while i >= 0 and L[i] == elt:
+      i -= 1
+
+   return L[:i+1]
+
+
+# create a polynomial with coefficients in a field; coefficients are in
+# increasing order of monomial degree so that, for example, [1,2,3]
+# corresponds to 1 + 2x + 3x^2
+@memoize
+def polynomialsOver(field=fractions.Fraction):
+
+   class Polynomial(DomainElement):
+      operatorPrecedence = 2
+
+      @classmethod
+      def factory(cls, L):
+         return Polynomial([cls.field(x) for x in L])
+
+      def __init__(self, c):
+         if type(c) is Polynomial:
+            self.coefficients = c.coefficients
+         elif isinstance(c, field):
+            self.coefficients = [c]
+         elif not hasattr(c, '__iter__') and not hasattr(c, 'iter'):
+            self.coefficients = [field(c)]
+         else:
+            self.coefficients = c
+
+         self.coefficients = strip(self.coefficients, field(0))
+
+
+      def isZero(self): return self.coefficients == []
+
+      def __repr__(self):
+         if self.isZero():
+            return '0'
+
+         return ' + '.join(['%s x^%d' % (a,i) if i > 0 else '%s'%a
+                              for i,a in enumerate(self.coefficients)])
+
+
+      def __abs__(self): return len(self.coefficients) # the valuation only gives 0 to the zero polynomial, i.e. 1+degree
+      def __len__(self): return len(self.coefficients)
+      def __sub__(self, other): return self + (-other)
+      def __iter__(self): return iter(self.coefficients)
+      def __neg__(self): return Polynomial([-a for a in self])
+
+      def iter(self): return self.__iter__()
+      def leadingCoefficient(self): return self.coefficients[-1]
+      def degree(self): return abs(self) - 1
+
+      @typecheck
+      def __eq__(self, other):
+         return self.degree() == other.degree() and all([x==y for (x,y) in zip(self, other)])
+
+      @typecheck
+      def __ne__(self, other):
+          return self.degree() != other.degree() or any([x!=y for (x,y) in zip(self, other)])
+
+      @typecheck
+      def __add__(self, other):
+         newCoefficients = [sum(x) for x in zip_longest(self, other, fillvalue=self.field(0))]
+         return Polynomial(newCoefficients)
+
+
+      @typecheck
+      def __mul__(self, other):
+         if self.isZero() or other.isZero():
+            return Zero()
+
+         newCoeffs = [self.field(0) for _ in range(len(self) + len(other) - 1)]
+
+         for i,a in enumerate(self):
+            for j,b in enumerate(other):
+               newCoeffs[i+j] += a*b
+
+         return Polynomial(newCoeffs)
+
+
+      @typecheck
+      def __divmod__(self, divisor):
+         quotient, remainder = Zero(), self
+         divisorDeg = divisor.degree()
+         divisorLC = divisor.leadingCoefficient()
+
+         while remainder.degree() >= divisorDeg:
+            monomialExponent = remainder.degree() - divisorDeg
+            monomialZeros = [self.field(0) for _ in range(monomialExponent)]
+            monomialDivisor = Polynomial(monomialZeros + [remainder.leadingCoefficient() / divisorLC])
+
+            quotient += monomialDivisor
+            remainder -= monomialDivisor * divisor
+
+         return quotient, remainder
+
+
+      @typecheck
+      def __truediv__(self, divisor):
+         if divisor.isZero():
+            raise ZeroDivisionError
+         return divmod(self, divisor)[0]
+
+
+      @typecheck
+      def __mod__(self, divisor):
+         if divisor.isZero():
+            raise ZeroDivisionError
+         return divmod(self, divisor)[1]
+
+
+   def Zero():
+      return Polynomial([])
+
+
+   Polynomial.field = field
+   Polynomial.__name__ = '(%s)[x]' % field.__name__
+   Polynomial.englishName = 'Polynomials in one variable over %s' % field.__name__
+   return Polynomial
+

scripts/finite_fields/test.py

@@ -0,0 +1,9 @@
+def test(expected, actual):
+   if expected != actual:
+      import sys, traceback
+      (filename, lineno, container, code) = traceback.extract_stack()[-2]
+      print("Test: %r failed on line %d in file %r.\nExpected %r but got %r\n" %
+         (code, lineno, filename, expected, actual))
+
+      sys.exit(1)
+

scripts/finite_fields/typecast-test.py

@@ -0,0 +1,12 @@
+from modp import *
+from polynomial import *
+
+
+mod3 = IntegersModP(3)
+Polynomial = polynomialsOver(mod3)
+x = mod3(1)
+p = Polynomial([1,2])
+
+x+p
+p+x
+

scripts/modp.py

@@ -0,0 +1,10 @@
+from finite_fields.modp import IntegersModP
+
+q = 0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001
+modq = IntegersModP(q)
+
+print("0x%x" % modq(2).inverse().n)
+inv2 = 0x39f6d3a994cebea4199cec0404d0ec02a9ded2017fff2dff7fffffff80000001
+assert modq(2).inverse().n == inv2
+print((2 * inv2) % q)
+

src/eq.rs

@@ -100,6 +100,9 @@ impl Circuit<bls12_381::Scalar> for MyCircuit {
 }
 
 fn main() {
+    let x = Scalar::from(2);
+    println!("{:?}", x.invert().unwrap());
+
     use std::time::Instant;
 
     let start = Instant::now();