plonk circuit setup routines

scripts/finite_fields/euclidean.py

@@ -2,14 +2,14 @@
 # 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)
+    if abs(a) < abs(b):
+        return gcd(b, a)
 
-   while abs(b) > 0:
-      _,r = divmod(a,b)
-      a,b = b,r
+    while abs(b) > 0:
+        _,r = divmod(a,b)
+        a,b = b,r
 
-   return a
+    return a
 
 
 # extendedEuclideanAlgorithm: int, int -> int, int, int
@@ -17,19 +17,18 @@ def 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) > abs(a):
+        (x,y,d) = extendedEuclideanAlgorithm(b, a)
+        return (y,x,d)
 
-   if abs(b) == 0:
-      return (1, 0, a)
+    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)
+    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

@@ -4,9 +4,9 @@ from polynomial import *
 from modp import *
 
 def p(L, q):
-   f = IntegersModP(q)
-   Polynomial = polynomialsOver(f).factory
-   return Polynomial(L)
+    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))

scripts/finite_fields/finitefield.py

@@ -1,6 +1,6 @@
 import random
-from polynomial import polynomialsOver
-from modp import *
+from .polynomial import polynomialsOver
+from .modp import *
 
 
 
@@ -9,23 +9,23 @@ from modp import *
 # 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__))
+    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
+    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
+    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
+    return True
 
 
 # generateIrreduciblePolynomial: int, int -> Polynomial
@@ -34,95 +34,95 @@ def isIrreducible(polynomial, p):
 # 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)
+    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)
+    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
+        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
+    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])
+    F23 = FiniteField(2,3)
+    x = F23([1,1])
 
-   F35 = FiniteField(3,5)
-   y = F35([1,1,2])
+    F35 = FiniteField(3,5)
+    y = F35([1,1,2])

scripts/finite_fields/modp-test.py

@@ -5,9 +5,8 @@ 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(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

@@ -4,77 +4,81 @@ from .numbertype import *
 
 # so all IntegersModP are instances of the same base class
 class _Modular(FieldElement):
-   pass
+    pass
 
 
 @memoize
 def IntegersModP(p):
-   # assume p is prime
+    # 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__))
+    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
+            self.field = IntegerModP
 
-      @typecheck
-      def __add__(self, other):
-         return IntegerModP(self.n + other.n)
+        @typecheck
+        def __add__(self, other):
+            return IntegerModP(self.n + other.n)
 
-      @typecheck
-      def __sub__(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)
+        @typecheck
+        def __mul__(self, other):
+            return IntegerModP(self.n * other.n)
 
-      def __neg__(self):
-         return IntegerModP(-self.n)
+        def __neg__(self):
+            return IntegerModP(-self.n)
 
-      @typecheck
-      def __eq__(self, other):
-         return isinstance(other, IntegerModP) and self.n == other.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 __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))
+        @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)
+        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__))
+            if d != 1:
+                raise Exception("Error: p is not prime in %s!" % (self.__name__))
 
-         return IntegerModP(x)
+            return IntegerModP(x)
 
-      def __abs__(self):
-         return abs(self.n)
+        def __abs__(self):
+            return abs(self.n)
 
-      def __str__(self):
-         return str(self.n)
+        def __str__(self):
+            return str(self.n)
 
-      def __repr__(self):
-         return '%d (mod %d)' % (self.n, self.p)
+        def __repr__(self):
+            return '%d (mod %d)' % (self.n, self.p)
 
-      def __int__(self):
-         return self.n
+        def __int__(self):
+            return self.n
 
-   IntegerModP.p = p
-   IntegerModP.__name__ = 'Z/%d' % (p)
-   IntegerModP.englishName = 'IntegersMod%d' % (p)
-   return IntegerModP
+        def __hash__(self):
+            return hash((self.n, self.p))
+
+    IntegerModP.p = p
+    IntegerModP.__name__ = 'Z/%d' % (p)
+    IntegerModP.englishName = 'IntegersMod%d' % (p)
+    return IntegerModP
 
 
 if __name__ == "__main__":
-   mod7 = IntegersModP(7)
+    mod7 = IntegersModP(7)

scripts/finite_fields/numbertype.py

@@ -2,38 +2,38 @@
 # this helps typechecking by never creating two separate
 # instances of a number class.
 def memoize(f):
-   cache = {}
+    cache = {}
 
-   def memoizedFunction(*args, **kwargs):
-      argTuple = args + tuple(kwargs)
-      if argTuple not in cache:
-         cache[argTuple] = f(*args, **kwargs)
-      return cache[argTuple]
+    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
+    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
+    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))
+        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 f(self, other)
 
-   return newF
+    return newF
 
 
 
@@ -41,58 +41,57 @@ def typecheck(f):
 # 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
+    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
+    # 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
+    # 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)
+        Q = self
+        R = self if n & 1 else self.__class__(1)
 
-      i = 2
-      while i <= n:
-         Q = (Q * Q)
+        i = 2
+        while i <= n:
+            Q = (Q * Q)
 
-         if n & i == i:
-            R = (Q * R)
+            if n & i == i:
+                R = (Q * R)
 
-         i = i << 1
+            i = i << 1
 
-      return R
+        return R
 
 
-   # requires the additional % operator (i.e. a Euclidean Domain)
-   def powmod(self, n, modulus):
-      if type(n) is not int:
-         raise TypeError
+    # 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)
+        Q = self
+        R = self if n & 1 else self.__class__(1)
 
-      i = 2
-      while i <= n:
-         Q = (Q * Q) % modulus
+        i = 2
+        while i <= n:
+            Q = (Q * Q) % modulus
 
-         if n & i == i:
-            R = (Q * R) % modulus
+            if n & i == i:
+                R = (Q * R) % modulus
 
-         i = i << 1
+            i = i << 1
 
-      return R
+        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)
-
+    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

@@ -13,37 +13,37 @@ 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]))
+    # 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]))

scripts/finite_fields/polynomial.py

@@ -8,13 +8,13 @@ from .numbertype import *
 
 # strip all copies of elt from the end of the list
 def strip(L, elt):
-   if len(L) == 0: return L
+    if len(L) == 0: return L
 
-   i = len(L) - 1
-   while i >= 0 and L[i] == elt:
-      i -= 1
+    i = len(L) - 1
+    while i >= 0 and L[i] == elt:
+        i -= 1
 
-   return L[:i+1]
+    return L[:i+1]
 
 
 # create a polynomial with coefficients in a field; coefficients are in
@@ -23,111 +23,121 @@ def strip(L, elt):
 @memoize
 def polynomialsOver(field=fractions.Fraction):
 
-   class Polynomial(DomainElement):
-      operatorPrecedence = 2
+    class Polynomial(DomainElement):
+        operatorPrecedence = 2
 
-      @classmethod
-      def factory(cls, L):
-         return Polynomial([cls.field(x) for x in L])
+        @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
+        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))
+            self.coefficients = strip(self.coefficients, field(0))
 
 
-      def isZero(self): return self.coefficients == []
+        def isZero(self): return self.coefficients == []
 
-      def __repr__(self):
-         if self.isZero():
-            return '0'
+        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)])
+            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 __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
+        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 __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 __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 __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()
+        @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)]
+            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
+            for i,a in enumerate(self):
+                for j,b in enumerate(other):
+                    newCoeffs[i+j] += a*b
 
-         return Polynomial(newCoeffs)
+            return Polynomial(newCoeffs)
 
 
-      @typecheck
-      def __divmod__(self, divisor):
-         quotient, remainder = Zero(), self
-         divisorDeg = divisor.degree()
-         divisorLC = divisor.leadingCoefficient()
+        @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])
+            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
+                quotient += monomialDivisor
+                remainder -= monomialDivisor * divisor
 
-         return quotient, remainder
+            return quotient, remainder
 
 
-      @typecheck
-      def __truediv__(self, divisor):
-         if divisor.isZero():
-            raise ZeroDivisionError
-         return divmod(self, divisor)[0]
+        @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]
+        @typecheck
+        def __mod__(self, divisor):
+            if divisor.isZero():
+                raise ZeroDivisionError
+            return divmod(self, divisor)[1]
 
+        def __call__(self, x):
+            if type(x) is int:
+                x = self.field(x)
+            assert type(x) is self.field
+            if self.isZero():
+                return self.field(0)
+            y = self.leadingCoefficient()
+            for coeff in self.coefficients[-2::-1]:
+                y = y * x + coeff
+            return y
 
-   def Zero():
-      return Polynomial([])
 
+    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
 
+    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

@@ -1,9 +1,8 @@
 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)
+    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

@@ -9,4 +9,3 @@ p = Polynomial([1,2])
 
 x+p
 p+x
-

scripts/halo/misc.py

@@ -0,0 +1,33 @@
+import random
+
+def sample_random(fp, seed):
+    rnd = random.Random(seed)
+    # Range of the field is 0 ... p - 1
+    return fp(rnd.randint(0, fp.p - 1))
+
+def is_power_of_two(n):
+    # Power of two number is represented by a single digit
+    # followed by zeroes.
+    return n & (n - 1) == 0
+
+#| ## Choosing roots of unity
+def get_omega(fp, n, seed=None):
+    """
+    Given a field, this method returns an n^th root of unity.
+    If the seed is not None then this method will return the
+    same n'th root of unity for every run with the same seed
+
+    This only makes sense if n is a power of 2.
+    """
+    assert is_power_of_two(n)
+    # https://crypto.stackexchange.com/questions/63614/finding-the-n-th-root-of-unity-in-a-finite-field
+    while True:
+        # Sample random x != 0
+        x = sample_random(fp, seed)
+        # Compute g = x^{(q - 1)/n}
+        y = pow(x, (fp.p - 1) // n)
+        # If g^{n/2} != 1 then g is a primitive root
+        if y != 1 and pow(y, n // 2) != 1:
+            assert pow(y, n) == 1, "omega must be 2nd root of unity"
+            return y
+

scripts/halo/pasta.py

@@ -0,0 +1,5 @@
+from finite_fields import finitefield
+
+p = 0x40000000000000000000000000000000224698fc094cf91b992d30ed00000001
+fp = finitefield.IntegersModP(p)
+

scripts/halo/polynomial_evalrep.py

@@ -0,0 +1,301 @@
+#| # Evaluation Representation of Polynomials and FFT optimizations
+#| In addition to the coefficient-based representation of polynomials used
+#| in babysnark.py, for performance we will also use an alternative
+#| representation where the polynomial is evaluated at a fixed set of points.
+#| Some operations, like multiplication and division, are significantly more
+#| efficient in this form.
+#| We can use FFT-based tools for efficiently converting
+#| between coefficient and evaluation representation.
+#|
+#| This library provides:
+#|  - Fast fourier transform for finite fields
+#|  - Interpolation and evaluation using FFT
+
+from finite_fields.finitefield import FiniteField
+from finite_fields.polynomial import polynomialsOver
+from finite_fields.euclidean import extendedEuclideanAlgorithm
+import random
+from finite_fields.numbertype import typecheck, memoize, DomainElement
+from functools import reduce
+import numpy as np
+
+#| ## Fast Fourier Transform on Finite Fields
+def fft_helper(a, omega, field):
+    """
+    Given coefficients A of polynomial this method does FFT and returns
+    the evaluation of the polynomial at [omega^0, omega^(n-1)]
+
+    If the polynomial is a0*x^0 + a1*x^1 + ... + an*x^n then the coefficients
+    list is of the form [a0, a1, ... , an].
+    """
+    n = len(a)
+    assert not (n & (n - 1)), "n must be a power of 2"
+
+    if n == 1:
+        return a
+
+    b, c = a[0::2], a[1::2]
+    b_bar = fft_helper(b, pow(omega, 2), field)
+    c_bar = fft_helper(c, pow(omega, 2), field)
+    a_bar = [field(1)] * (n)
+    for j in range(n):
+        k = j % (n // 2)
+        a_bar[j] = b_bar[k] + pow(omega, j) * c_bar[k]
+    return a_bar
+
+
+#| ## Representing a polynomial by evaluation at fixed points
+@memoize
+def make_polynomial_evalrep(field, omega, n):
+    assert n & n - 1 == 0, "n must be a power of 2"
+
+    # Check that omega is an n'th primitive root of unity
+    assert type(omega) is field
+    omega = field(omega)
+    assert omega**(n) == 1
+    _powers = [omega**i for i in range(n)]
+    assert len(set(_powers)) == n
+
+    _poly_coeff = polynomialsOver(field)
+
+    class PolynomialEvalRep(object):
+
+        def __init__(self, xs, ys):
+            # Each element of xs must be a power of omega.
+            # There must be a corresponding y for every x.
+            if type(xs) is not tuple:
+                xs = tuple(xs)
+            if type(ys) is not tuple:
+                ys = tuple(ys)
+
+            assert len(xs) <= n+1
+            assert len(xs) == len(ys)
+            for x in xs:
+                assert x in _powers
+            for y in ys:
+                assert type(y) is field
+
+            self.evalmap = dict(zip(xs, ys))
+
+        @classmethod
+        def from_coeffs(cls, poly):
+            assert type(poly) is _poly_coeff
+            assert poly.degree() <= n
+            padded_coeffs = poly.coefficients + [field(0)] * (n - len(poly.coefficients))
+            ys = fft_helper(padded_coeffs, omega, field)
+            xs = [omega**i for i in range(n) if ys[i] != 0]
+            ys = [y for y in ys if y != 0]
+            return cls(xs, ys)
+
+        def to_coeffs(self):
+            # To convert back to the coefficient form, we use polynomial interpolation.
+            # The non-zero elements stored in self.evalmap, so we fill in the zero values
+            # here.
+            ys = [self.evalmap[x] if x in self.evalmap else field(0) for x in _powers]
+            coeffs = [b / field(n) for b in fft_helper(ys, 1 / omega, field)]
+            return _poly_coeff(coeffs)
+
+        _lagrange_cache = {}
+        def __call__(self, x):
+            if type(x) is int:
+                x = field(x)
+            assert type(x) is field
+            xs = _powers
+
+            def lagrange(x, xi):
+                # Let's cache lagrange values
+                if (x,xi) in PolynomialEvalRep._lagrange_cache:
+                    return PolynomialEvalRep._lagrange_cache[(x,xi)]
+
+                mul = lambda a,b: a*b
+                num = reduce(mul, [x  - xj for xj in xs if xj != xi], field(1))
+                den = reduce(mul, [xi - xj for xj in xs if xj != xi], field(1))
+                PolynomialEvalRep._lagrange_cache[(x,xi)] = num / den
+                return PolynomialEvalRep._lagrange_cache[(x,xi)]
+
+            y = field(0)
+            for xi, yi in self.evalmap.items():
+                y += yi * lagrange(x, xi)
+            return y
+
+        def __mul__(self, other):
+            # Scale by integer
+            if type(other) is int:
+                other = field(other)
+            if type(other) is field:
+                return PolynomialEvalRep(self.evalmap.keys(),
+                                         [other * y for y in self.evalmap.values()])
+
+            # Multiply another polynomial in the same representation
+            if type(other) is type(self):
+                xs = []
+                ys = []
+                for x, y in self.evalmap.items():
+                    if x in other.evalmap:
+                        xs.append(x)
+                        ys.append(y * other.evalmap[x])
+                return PolynomialEvalRep(xs, ys)
+
+        @typecheck
+        def __iadd__(self, other):
+            # Add another polynomial to this one in place.
+            # This is especially efficient when the other polynomial is sparse,
+            # since we only need to add the non-zero elements.
+            for x, y in other.evalmap.items():
+                if x not in self.evalmap:
+                    self.evalmap[x] = y
+                else:
+                    self.evalmap[x] += y
+            return self
+
+        @typecheck
+        def __add__(self, other):
+            res = PolynomialEvalRep(self.evalmap.keys(), self.evalmap.values())
+            res += other
+            return res
+
+        def __sub__(self, other): return self + (-other)
+        def __neg__(self): return PolynomialEvalRep(self.evalmap.keys(),
+                                                    [-y for y in self.evalmap.values()])
+
+        def __truediv__(self, divisor):
+            # Scale by integer
+            if type(divisor) is int:
+                other = field(divisor)
+            if type(divisor) is field:
+                return self * (1/divisor)
+            if type(divisor) is type(self):
+                res = PolynomialEvalRep((),())
+                for x, y in self.evalmap.items():
+                    assert x in divisor.evalmap
+                    res.evalmap[x] = y / divisor.evalmap[x]
+                return res
+            return NotImplemented
+
+        def __copy__(self):
+            return PolynomialEvalRep(self.evalmap.keys(), self.evalmap.values())
+
+        def __repr__(self):
+            return f'PolyEvalRep[{hex(omega.n)[:15]}...,{n}]({len(self.evalmap)} elements)'
+
+        @classmethod
+        def divideWithCoset(cls, p, t, c=field(3)):
+            """
+              This assumes that p and t are polynomials in coefficient representation,
+            and that p is divisible by t.
+               This function is useful when t has roots at some or all of the powers of omega,
+            in which case we cannot just convert to evalrep and use division above
+            (since it would cause a divide by zero.
+               Instead, we evaluate p(X) at powers of (c*omega) for some constant cofactor c.
+            To do this efficiently, we create new polynomials, pc(X) = p(cX), tc(X) = t(cX),
+            and evaluate these at powers of omega. This conversion can be done efficiently
+            on the coefficient representation.
+               See also: cosetFFT in libsnark / libfqfft.
+               https://github.com/scipr-lab/libfqfft/blob/master/libfqfft/evaluation_domain/domains/extended_radix2_domain.tcc
+            """
+            assert type(p) is _poly_coeff
+            assert type(t) is _poly_coeff
+            # Compute p(cX), t(cX) by multiplying coefficients
+            c_acc = field(1)
+            pc = _poly_coeff(list(p.coefficients))  # make a copy
+            for i in range(p.degree() + 1):
+                pc.coefficients[-i-1] *= c_acc
+                c_acc *= c
+            c_acc = field(1)
+            tc = _poly_coeff(list(t.coefficients))  # make a copy
+            for i in range(t.degree() + 1):
+                tc.coefficients[-i-1] *= c_acc
+                c_acc *= c
+
+            # Divide using evalrep
+            pc_rep = cls.from_coeffs(pc)
+            tc_rep = cls.from_coeffs(tc)
+            hc_rep = pc_rep / tc_rep
+            hc = hc_rep.to_coeffs()
+
+            # Compute h(X) from h(cX) by dividing coefficients
+            c_acc = field(1)
+            h = _poly_coeff(list(hc.coefficients))  # make a copy
+            for i in range(hc.degree() + 1):
+                h.coefficients[-i-1] /= c_acc
+                c_acc *= c
+
+            # Correctness checks
+            # assert pc == tc * hc
+            # assert p == t * h
+            return h
+
+
+    return PolynomialEvalRep
+
+#| ## Sparse Matrix
+#| In our setting, we have O(m*m) elements in the matrix, and expect the number of
+#| elements to be O(m).
+#| In this setting, it's appropriate to use a rowdict representation - a dense
+#| array of dictionaries, one for each row, where the keys of each dictionary
+#| are column indices.
+
+class RowDictSparseMatrix():
+    # Only a few necessary methods are included here.
+    # This could be replaced with a generic sparse matrix class, such as scipy.sparse,
+    # but this does not work as well with custom value types like Fp
+
+    def __init__(self, m, n, zero=None):
+        self.m = m
+        self.n = n
+        self.shape = (m,n)
+        self.zero = zero
+        self.rowdicts = [dict() for _ in range(m)]
+
+    def __setitem__(self, key, v):
+        i, j = key
+        self.rowdicts[i][j] = v
+
+    def __getitem__(self, key):
+        i, j = key
+        return self.rowdicts[i][j] if j in self.rowdicts[i] else self.zero
+
+    def items(self):
+        for i in range(self.m):
+            for j, v in self.rowdicts[i].items():
+                yield (i,j), v
+
+    def dot(self, other):
+        if isinstance(other, np.ndarray):
+            assert other.dtype == 'O'
+            assert other.shape in ((self.n,),(self.n,1))
+            ret = np.empty((self.m,), dtype='O')
+            ret.fill(self.zero)
+            for i in range(self.m):
+                for j, v in self.rowdicts[i].items():
+                    ret[i] += other[j] * v
+            return ret
+
+    def to_dense(self):
+        mat = np.empty((self.m, self.n), dtype='O')
+        mat.fill(self.zero)
+        for (i,j), val in self.items():
+            mat[i,j] = val
+        return mat
+
+    def __repr__(self): return repr(self.rowdicts)
+
+#-
+# Examples
+if __name__ == '__main__':
+    import misc
+
+    Fp = FiniteField(52435875175126190479447740508185965837690552500527637822603658699938581184513,1)  # (# noqa: E501)
+    Poly = polynomialsOver(Fp)
+
+    n = 8
+    omega = misc.get_omega(Fp, n)
+    PolyEvalRep = make_polynomial_evalrep(Fp, omega, n)
+
+    f = Poly([1,2,3,4,5])
+    xs = tuple([omega**i for i in range(n)])
+    ys = tuple(map(f, xs))
+    # print('xs:', xs)
+    # print('ys:', ys)
+
+    assert f == PolyEvalRep(xs, ys).to_coeffs()

scripts/halo/test.py

@@ -0,0 +1,119 @@
+import sys
+sys.path.append("..")
+
+import random
+import misc
+import pasta
+from polynomial_evalrep import make_polynomial_evalrep
+
+n = 8
+omega_base = misc.get_omega(pasta.fp, 2**32, seed=0)
+assert misc.is_power_of_two(8)
+omega = omega_base ** (2 ** 32 // n)
+# Order of omega is n
+assert omega ** n == 1
+# Compute complete roots of this group
+ROOTS = [omega ** i for i in range(n)]
+
+PolyEvalRep = make_polynomial_evalrep(pasta.fp, omega, n)
+
+import numpy as np
+from tabulate import tabulate
+
+# Wires
+a = ["x",  "v1", "v2", "1", "1",  "v3", "e1", "e2"]
+b = ["x",  "x",  "x",  "5", "35", "5",  "e3", "e4"]
+c = ["v1", "v2", "v3", "5", "35", "35", "e5", "e6"]
+
+wires = a + b + c
+
+# Gates
+# La + Rb + Oc + Mab + C = 0
+add =           np.array([1, 1, 0, -1,  0])
+mul =           np.array([0, 0, 1, -1,  0])
+const5 =        np.array([0, 1, 0,  0, -5])
+public_input =  np.array([0, 1, 0,  0,  0])
+empty =         np.array([0, 0, 0,  0,  0])
+
+gates_matrix = np.array(
+    [mul, mul, add, const5, public_input, add, empty, empty])
+print("Wires:")
+print(tabulate([["a ="] + a, ["b ="] + b, ["c ="] + c]))
+print()
+print("Gates:")
+print(gates_matrix)
+print()
+
+# The index of the public input in the gates_matrix
+# We specify its position and its value
+public_input_values = [(4, 35)]
+
+def permute_indices(wires):
+    size = len(wires)
+    permutation = [i + 1 for i in range(size)]
+    for i in range(size):
+        for j in range(i + 1, size):
+            if wires[i] == wires[j]:
+                permutation[i], permutation[j] = permutation[j], permutation[i]
+                break
+    return permutation
+
+permutation = permute_indices(wires)
+
+table = [
+    ["Wires"] + wires,
+    ["Indices"] + list(i + 1 for i in range(len(wires))),
+    ["Permutations"] + permutation
+]
+print(tabulate(table))
+print()
+
+import misc
+from pasta import fp
+
+def setup(wires, gates_matrix):
+    # Section 8.1
+    # The selector polynomials that define the circuit's arithmetisation
+    gates_matrix = gates_matrix.transpose()
+    ql = PolyEvalRep(ROOTS, [fp(i) for i in gates_matrix[0]])
+    qr = PolyEvalRep(ROOTS, [fp(i) for i in gates_matrix[1]])
+    qm = PolyEvalRep(ROOTS, [fp(i) for i in gates_matrix[2]])
+    qo = PolyEvalRep(ROOTS, [fp(i) for i in gates_matrix[3]])
+    qc = PolyEvalRep(ROOTS, [fp(i) for i in gates_matrix[4]])
+    selector_polys = [ql, qr, qm, qo, qc]
+
+    public_input = [fp(0) for i in range(len(ROOTS))]
+    for (index, value) in public_input_values:
+        # This is negative because the value is added to
+        # the output of the const selector poly:
+        # La + Rb + Oc + Mab + (C + PI) = 0
+        public_input[index] = fp(-value)
+    public_input_poly = PolyEvalRep(ROOTS, public_input)
+
+    # Identity permutations applied to a, b, c
+    # Ideally H, k_1 H, k_2 H are distinct cosets of H
+    # Here we just sample k and assume it's high-order
+    # Random high order k to form distinct cosets
+    k = misc.sample_random(fp)
+    id_domain_a = ROOTS
+    id_domain_b = [k * root for root in ROOTS]
+    id_domain_c = [k**2 * root for root in ROOTS]
+    id_domain = id_domain_a + id_domain_b + id_domain_c
+
+    # Intermediate step where we permute the positions of the domain
+    # generated above
+    permuted_domain = [id_domain[i - 1] for i in permutation]
+    permuted_domain_a = permuted_domain[:n]
+    permuted_domain_b = permuted_domain[n:2 * n]
+    permuted_domain_c = permuted_domain[2*n:3 * n]
+
+    # The copy permuation applied to a, b, c
+    # Returns the permuted index value (corresponding root of unity coset)
+    # when evaluated on the domain.
+    ssigma_1 = PolyEvalRep(ROOT, permuted_domain_a)
+    ssigma_2 = PolyEvalRep(ROOT, permuted_domain_b)
+    ssigma_3 = PolyEvalRep(ROOT, permuted_domain_c)
+    copy_permutes = [ssigma_1, ssigma_2, ssigma_3]
+
+setup(wires, gates_matrix)
+