classes for CocksPinch curves: k=6 D=3, k=8 D=4, other k,D

CocksPinch_k.py

+# New Cocks-Pinch curves with rho = 2 but optimal ate pairing available.
+# k=5 or k=7
+# any j possible
+import sage
+
+from exceptions import ValueError
+from sage.functions.log import log
+from sage.functions.other import ceil
+from sage.functions.other import sqrt
+from sage.arith.misc import XGCD, xgcd
+from sage.rings.integer import Integer
+from sage.rings.finite_rings.finite_field_constructor import FiniteField
+from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+from sage.schemes.elliptic_curves.ell_finite_field import EllipticCurve_finite_field
+from sage.schemes.elliptic_curves.constructor import EllipticCurve
+from sage.misc.functional import cyclotomic_polynomial
+
+import PairingFriendlyCurve
+from PairingFriendlyCurve import print_parameters, print_parameters_for_RELIC
+from PairingFriendlyCurve import get_curve_generator_order_r
+
+class CocksPinch_k(EllipticCurve_finite_field):
+    """
+    A Cocks-Pinch curve of embedding degree k (=5,7), any given D, rho>=2 and optimal ate pairing available.
+    """
+
+    def __init__(self, k, T, D, exp_tr, p, r, tr, y, a, b):
+        """
+        :param k     : embedding degree
+        :param T     : seed
+        :param exp_tr: integer such that tr = u^exp_tr + 1 mod r
+        :param a     : curve parameter in E: y^2 = x^3 + a*x + b
+        :param b     : curve parameter in E: y^2 = x^3 + a*x + b
+        :param p     : prime, field characteristic
+        :param r     : prime, divides the curve order
+        :param tr    : trace of the curve
+        :param y     : integer s.t. p = (tr^2 + D*y^2)/4
+        """
+        #T = Integer(u)
+        self._k = int(k) # embedding degree
+        self._D = Integer(D) # curve discriminant
+        self._a = Integer(a)
+        self._b = Integer(b)
+        self._u = Integer(T)
+        self._i = int(exp_tr)
+
+        self._r = Integer(r)
+        self._p = Integer(p)
+        self._tr= Integer(tr)
+        self._y = Integer(y)
+        
+        #self.polynomial_r = [1,1,1,1,1] # 1+u+u^2+u^3+u^4
+        self.polynomial_r = cyclotomic_polynomial(self._k).list()
+        r_ = sum([Integer(self.polynomial_r[i])*self._u**i for i in range(len(self.polynomial_r))])
+        if self._r != r_:
+            raise ValueError("Error: r does not match the polynomial r(T)\nr = {:= #x} #(input)\nr = {:= #x} # from polynomial r(T) = {}".format(r_,self._r,self.polynomial_r))
+        else:
+            print("valid r")
+        # can check that tr mod r = T^i mod r...
+
+        if ((self._p+1-self._tr) % self._r) != 0:
+            raise ValueError("Error: p+1-tr mod r != 0,\np+1-tr={}\np     ={}\nr     ={}\ntr    ={}\np+1-tr % r = {}\n".format((self._p+1-self._tr), self._p, self._r, self._tr, ((self._p+1-self._tr) % self._r) ))
+            
+        self._c = (self._p+1-self._tr) // self._r
+        
+        self._pbits = self._p.nbits()
+        
+        try:
+            self._Fp = FiniteField(self._p)
+        except ValueError as err:
+            print("ValueError creating Fp: {}".format(err))
+            raise
+        except:
+            print("Error creating Fp")
+            raise
+        if not self._r.is_prime():
+            raise ValueError("Error r is not prime")
+
+        self._Fpz = PolynomialRing(self._Fp, names=('z',))
+        (self._z,) = self._Fpz._first_ngens(1)
+        
+        self._ap = self._Fp(a)
+        self._bp = self._Fp(b)
+        try:
+            # this init function of super inherits from class EllipticCurve_generic defined in ell_generic.py
+            # __init__ method inherited from ell_generic
+            EllipticCurve_finite_field.__init__(self, self._Fp, [0,0,0,self._ap,self._bp])
+        except ValueError as err:
+            print("ValueError at EllipticCurve_finite_field.__init__: {}".format(err))
+            raise
+        except:
+            print("An error occured when initialising the elliptic curve")
+            raise
+        if self._pbits <= 384:
+            print("check curve order")
+            self.order_checked = super(CocksPinch_k,self).order()
+            if self.order_checked != (self._p+1-self._tr):
+                #if (self.order_checked % self._r) != 0:
+                print("Error, wrong order")
+                raise ValueError("Wrong curve order: this one is a twist")
+        else:
+            print("check order (randomized because pnbits = {})".format(self._pbits))
+            self.curve_order = self._p + Integer(1) - self._tr
+            self.twist_order = self._p + Integer(1) + self._tr
+            for i in range(10):
+                P = self.random_element()
+                if self.curve_order*P != self(0):
+                    if self.twist_order*P == self(0):
+                        raise ValueError("Wrong curve order: this one is a twist: (p+1+tr)*P = 0\ntr={}\np+1-tr={}\np+1+tr={}\n".format(self._tr,self.curve_order,self.twist_order))
+                    else:
+                        print("computing actual curve order to compare")
+                        self.order_checked = super(CocksPinch_k,self).order()
+                        raise ValueError("Wrong curve order:\np+1-tr        = {}\np+1+tr        = {}\nchecked order = {}\np             = {}".format(self.curve_order,self.twist_order,self.order_checked,self._p))
+        print("ok")
+
+        # computes a generator
+        #print("compute a generator")
+        self._G = get_curve_generator_order_r(self)
+        self._Gx = self._G[0]
+        self._Gy = self._G[1]
+        
+    def _repr_(self):
+        return "CP"+str(self._k)+" p"+str(self._pbits)+" (modified Cocks-Pinch k="+str(self._k)+") curve with seed "+str(self._u)+"\n"+super(CocksPinch_k,self)._repr_()
+
+    def u(self):
+        return self._u
+    def p(self):
+        return self._p
+    def r(self):
+        return self._r
+    def c(self):
+        return self._c
+    def tr(self):
+        return self._tr
+    def y(self):
+        return self._y
+    def a(self):
+        return self._a # 0
+    def ap(self):
+        return self._ap # 0
+    def b(self):
+        return self._b # Integer
+    def bp(self):
+        return self._bp # in Fp (finite field element)
+    def beta(self):
+        return None
+    def lamb(self):
+        return None
+    
+    def k(self):
+        return self._k
+    def Fp(self):
+        return self._Fp
+    def Fpz(self):
+        return self._Fpz, self._z
+    def G(self):
+        return self._G
+
+    def print_parameters(self):
+        PairingFriendlyCurve.print_parameters(self)
+        
+    def print_parameters_for_RELIC(self):
+        PairingFriendlyCurve.print_parameters_for_RELIC(self)
+

PairingFriendlyCurve.py

@@ -3,6 +3,8 @@ import sage
 from exceptions import ValueError
 from sage.functions.log import log
 from sage.functions.other import ceil
+from sage.arith.functions import lcm
+from sage.arith.misc import GCD, gcd
 from sage.arith.misc import XGCD, xgcd
 from sage.rings.integer import Integer
 from sage.rings.integer_ring import Z, ZZ
@@ -210,74 +212,168 @@ def print_parameters_for_RELIC(E):
     print("Gy ={:= #x}".format(Integer(E._Gy)))
 
 
-def compute_a_b(D,p,tr):
-    if ((-D) % 4) == 2 or ((-D) % 4) == 3:
-        D = 4*D
+
+def check_order(Fp,p,tr,a,b):
+    curve_order = p+1-tr
+    twist_order = p+1+tr
+    E = EllipticCurve([0,0,0,Fp(a),Fp(b)])
+    E0 = E(0)
+    ok = True
+    twist_ok = False
+    for i in range(10):
+        P = E.random_element()
+        if curve_order*P != E0:
+            ok = False
+            if twist_order*P == E0:
+                twist_ok = True
+            break
+    return ok, twist_ok
+
+def compute_a_b(D,p,tr,y=None):
+    if D < 0:
+        D = -D
     p = Integer(p)
     tr = Integer(tr)
+    Fp = GF(p, proof=False)
+    if D == 1 or D == 4 or D == 3:
+        if y == None:
+            assert ((4*p-tr**2) % D) == 0
+            y = ZZ(sqrt((4*p-tr**2)//D))
+        else:
+            y = Integer(y)
+        assert tr**2+D*y**2 == 4*p
+    if D==1 or D == 4:
+        # beta s.t. beta^2 = -1 mod p -> beta = tr/y mod p
+        # invert y mod p
+        #g,u,v = xgcd(y,p) # such that g = u*y + v*p
+        #beta = (u*tr) % p
+        beta = Fp(tr)/Fp(y)
+        a, ap = get_curve_parameter_a_j1728(tr, y, beta, p, Fp)
+        b=0
+        return [(a,b)]
+    if D == 3:
+        # beta: s.t. beta^2+beta+1 = 0 mod p: (-1+sqrt(-3))/2
+        beta = (Fp(tr)/Fp(y) -Fp(1))/Fp(2)
+        b, bp = get_curve_parameter_b_j0(tr, y, beta, p, Fp)
+        a=0
+        return [(a,b)]
+    if ((-D) % 4) == 2 or ((-D) % 4) == 3:
+        D = 4*D
     H=hilbert_class_polynomial(-D)
     print("#    H(-D).degree() is {}".format(H.degree()))
-    Fp = GF(p, proof=False)
-    #Fp = GF(p)
     rH=H.change_ring(Fp).roots()
     min_bnbits = p.nbits()+1
-    min_a = None
-    min_b = None
-    for i,j_expj in enumerate(rH):
-        # y^2 + xy = x^3 - (36/(j - 1728))*x - 1/(j - 1728)
-        # (y+1/2*x)^2 - 1/4*x^2 = x^3  - (36/(j - 1728))*x - 1/(j - 1728)
-        # y' = y+1/2*x
-        # y'^2 = x^3 + 1/4*x^2  - (36/(j - 1728))*x - 1/(j - 1728)
-        # x' = (x-1/12)
-        # y'^2 = x^3 - 1/48*j/(j - 1728)*x + 1/864*j/(j - 1728)
-        # (a,b) ~ (a/u^4, b/u^6)
-        j = Fp(j_expj[0])
-        if j == 0:
-            print("bizarre j=0")
-            continue
-        elif j == 1728:
-            print("bizarre j=1728")
-            continue
-        fac_j = Fp(j)/Fp(j - 1728)
-        a = -fac_j/Fp(48)
-        b = fac_j/Fp(864)
-        u4 = a/Fp(-3)
-        if not u4.is_square():
-            continue
-        u2=u4.sqrt()
-        if not u2.is_square():
-            continue
-        u=u2.sqrt()
-        #print("j[%d]: ok" % i)
-        #print("j = {}".format(j))
-        a = -3
-        #print("a = {}".format(a))
-        #print("b = {} in {}".format(b, b.parent()))
-        b = b/u**6
+    # can we have a=-3?
+    ab = []
+    a_wanted = -3
+    while len(ab) == 0 and a_wanted < 10:
+        min_a = None
+        min_b = None
+        print("#looking for a={}".format(a_wanted))
+        for i,j_expj in enumerate(rH):
+            # y^2 + xy = x^3 - (36/(j - 1728))*x - 1/(j - 1728)
+            # (y+1/2*x)^2 - 1/4*x^2 = x^3  - (36/(j - 1728))*x - 1/(j - 1728)
+            # y' = y+1/2*x
+            # y'^2 = x^3 + 1/4*x^2  - (36/(j - 1728))*x - 1/(j - 1728)
+            # x' = (x-1/12)
+            # y'^2 = x^3 - 1/48*j/(j - 1728)*x + 1/864*j/(j - 1728)
+            # (a,b) ~ (a/u^4, b/u^6)
+            j = Fp(j_expj[0])
+            if j == 0:
+                print("unexpected j=0, a=0")
+                continue
+            elif j == 1728:
+                print("unexpected j=1728, b=0")
+                continue
+            #fac_j = Fp(j)/Fp(j - 1728)
+            #a = -fac_j/Fp(48)
+            #b = fac_j/Fp(864)
+            a = -3*j*(j - 1728)
+            b = 2*j*(j - 1728)**2
+            u4 = a/Fp(a_wanted)
+            if not u4.is_square():
+                continue
+            u2=u4.sqrt()
+            if (p % 4) == 3 and not u2.is_square():
+                u2 = -u2
+            if not u2.is_square():
+                continue
+            u=u2.sqrt()
+            a = a_wanted
+            b = b/u**6
+            # check that this is the curve and not the twist
+            ok, twist_ok = check_order(Fp,p,tr,a,b)
+            if not ok:
+                continue
+            b_ZZ = ZZ(b)
+            if abs(b_ZZ-p) < abs(b_ZZ):
+                b_ZZ -= p
+            bnbits = b_ZZ.nbits()
+            if bnbits < min_bnbits:
+                min_bnbits = bnbits
+                min_a = a
+                min_b = b_ZZ
+        if min_a != None and min_b != None:
+            ab.append((min_a, min_b))
+        if a_wanted == -3:
+            a_wanted = 1
+        else:
+            a_wanted += 1
+    if len(ab) == 0:
+        j, je = rH[0]
+        j = ZZ(j)
+        if abs(j-p) < abs(j):
+            j = j-p
+        print("j={}".format(j))
+        a = ZZ(-3*j*(j - 1728))
+        b = ZZ(2*j*(j - 1728)**2)
+        print("a,b = {},{}".format(a,b))
+        # simplify
+        g = gcd(a,b)
+        gf = g.factor()
+        u = 1
+        for fac,ei in gf:
+            if fac > 1:
+                while (a % (fac**2)) == 0 and (b%(fac**3)) == 0:
+                    a = a // fac**2 # 4
+                    b = b // fac**3 # 6
+                    u *= fac
+        print("a={}, b= {}, simplified by u = {}".format(a,b,u))
+        if abs((a % p)) < abs(a):
+            a = a % p
+            if abs(a-p) < abs(a):
+                a = a-p
+        if abs((b % p)) < abs(b):
+            b = b % p
+            if abs(b-p) < abs(b):
+                b = b-p
         # check that this is the curve and no the twist
-        E = EllipticCurve([0,0,0,Fp(a),Fp(b)])
-        E0 = E(0)
-        curve_order = p+1-tr
-        twist_order = p+1+tr
-        ok = True
-        for i in range(10):
-            P = E.random_element()
-            if curve_order*P != E0:
-                ok = False
-                break
-        if not ok:
-            continue
-        b_ZZ = ZZ(b)
-        if abs(b_ZZ-p) < abs(b_ZZ):
-            b_ZZ -= p
-        bnbits = b_ZZ.nbits()
-        if bnbits < min_bnbits:
-            min_bnbits = bnbits
-            min_a = a
-            min_b = b_ZZ
-    return min_a, min_b
+        ok, twist_ok = check_order(Fp,p,tr,a,b)
+        if not ok and not twist_ok:
+            print("pb with a,b")
+            return []
+        ns = -2
+        nsp= Fp(-2)
+        while not ok and twist_ok: # this is the twist
+            # find a non-square
+            ns += 1
+            nsp += 1
+            while nsp.is_square():
+                ns += 1
+                nsp += 1
+            a = a*ns**2
+            b = b*ns**3
+            ok, twist_ok = check_order(Fp,p,tr,a,b)
+            print("ns = {}, try with (a,b) = ({},{})".format(ns,a,b))
+        ab = [(a,b)]
+    return ab
 
-def convert_CocksPinch57(k,T,r,i,tr,y,D,p):
+def convert_CocksPinch567(k,T,r,i,tr,y,D,p):
+    p = Integer(p)
+    T = Integer(T)
+    r = Integer(r)
+    y = Integer(y)
+    tr = Integer(tr)
     if k == 5:
         # Phi_5(T) = T^4+T^3+T^2+T+1
         if i >= 1 and i <= 3:
@@ -286,6 +382,13 @@ def convert_CocksPinch57(k,T,r,i,tr,y,D,p):
             t0 = -(T**3+T**2+T)
         else:
             raise ValueError("k=5, i s.t. tr=T^i+1 mod r should be 1 <= i <= 4 but got i = {}".format(i))
+    elif k == 6:# i=1, or i=5
+        if i == 1:
+            t0 = T+1
+        elif i ==5:
+            t0 = -T+2
+        else:
+            raise ValueError("k=6, i s.t. tr=T^i+1 mod r should be 1 or 5 but got i = {}".format(i))        
     elif k == 7:
         # Phi_7(T) = T^6+T^5+T^4+T^3+T^2+T+1
         if i >= 1 and i <= 5:
@@ -294,6 +397,15 @@ def convert_CocksPinch57(k,T,r,i,tr,y,D,p):
             t0 = -(T**5+T**4+T**3+T**2+T)
         else:
             raise ValueError("k=7, i s.t. tr=T^i+1 mod r should be 1 <= i <= 6 but got i = {}".format(i))
+    elif k == 8:# i=1,3,5,7
+        if i == 1 or i == 3:
+            t0 = T**i+1
+        elif i ==5:
+            t0 = -T+1
+        elif i == 7:
+            t0 = -T**3 + 1
+        else:
+            raise ValueError("k=8, i s.t. tr=T^i+1 mod r should be 1,3,5, or 7 but got i = {}".format(i))
     
     assert (tr-t0) % r == 0
     ht = (tr-t0)//r
@@ -302,9 +414,34 @@ def convert_CocksPinch57(k,T,r,i,tr,y,D,p):
         y0 = y0-r
     assert (y-y0) % r == 0
     hy = (y-y0)//r
+    if k == 6 and D == 3:
+        if (T % 3) == 2:
+            raise ValueError("Error T = 2 mod 3, should be 1 or 2")
+        if  i == 1:
+            if (T % 3) == 0:
+                expected_y0 = -(T**2-2*T)//3
+            elif (T % 3) == 1:
+                expected_y0 = (2+T**2)//3
+        elif i == 5:
+            if (T % 3) == 0:
+                expected_y0 = (-3+2*T-T**2)//3
+            elif (T % 3) == 1:
+                expected_y0 = (T**2-1)//3
+        if y0 == expected_y0:
+            pass
+        elif y0 == -expected_y0: # y = y0 + hy*r, -y = -y0-hy*r
+            y = -y
+            hy = -hy
+            y0 = -y0
+            print("changed y,y0,hy -> -y,-y0,-hy")
+        else:
+            raise ValueError("Error not able to find compatible y0, k={}, y % r = {}, -y mod r = {}, expected y0 = {}, (y-y0)%r = {}, (y+y0)%r = {}".format(k, (y%r), (-y%r), expected_y0, (y-expected_y0)%r, (y+expected_y0)%r))
+        
     assert ((t0 + ht*r)**2 + D*(y0+hy*r)**2) % 4 == 0
     assert p == ((t0 + ht*r)**2 + D*(y0+hy*r)**2)//4
 
+    p = Integer(p)
+    r = Integer(r)
     print("    {{\'T\':{:#x},".format(Integer(T)))
     print("     \'k\':{},".format(k))
     print("     \'ht\':{:#x},".format(Integer(ht)))
@@ -318,13 +455,15 @@ def convert_CocksPinch57(k,T,r,i,tr,y,D,p):
     print("     \'y0\':{:#x},".format(Integer(y0)))
     print("     \'y\':{:#x},".format(Integer(y)))
     print("     \'pnbits\':{},".format(p.nbits()))
-    a, b = compute_a_b(D,p)
-    if a != None and b != None:
-        print("     \'a\':{},".format(Integer(a)))
-        print("     \'b\':{:#x},".format(Integer(b)))
-    else:
-        print("     #\'a\':,")
-        print("     #\'b\':,")
+    print("     \'rnbits\':{},".format(r.nbits()))
+    pairs_ab = compute_a_b(D,p,tr)
+    for a,b in pairs_ab:
+        if a != None and b != None:
+            print("     \'a\':{:#x},".format(Integer(a)))
+            print("     \'b\':{:#x},".format(Integer(b)))
+        else:
+            print("     #\'a\':,")
+            print("     #\'b\':,")
     print("    },")
 
 def print_JL_polys(C,Y):