# New Cocks-Pinch curves with rho = 2 but optimal ate pairing available. # j = 1728 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 import PairingFriendlyCurve from PairingFriendlyCurve import get_curve_parameter_a_j1728, get_curve_generator_order_r_j1728 from PairingFriendlyCurve import print_parameters, print_parameters_for_RELIC class CocksPinch8(EllipticCurve_finite_field): """ A Cocks-Pinch curve of embedding degree k=8, D=4, rho=2 and optimal ate pairing available. """ def __init__(self, u, hy, ht, exp_tr=None, p=None, r=None, c=None, tr=None, y=None, a=None, beta=None, lamb=None): """ :param u : seed :param hy : integer such that y = y0 + hy*r :param ht : integer such that tr = tr0 + ht*r :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 :param p : prime, field characteristic :param r : prime, divides the curve order :param c : cofactor such that r*c = (p+1-tr), the curve order :param tr : trace of the curve :param y : integer s.t. p = (tr^2 + D*y^2)/4 :param beta : integer s.t. beta^2+1 = 0 mod p :param lamb : integer s.t. lamb^2+1 = 0 mod r """ hy = Integer(hy) ht = Integer(ht) self._k = 8 # embedding degree self._D = 4 # curve discriminant self._b = 0 self._u = Integer(u) self._hy = hy self._ht = ht if exp_tr == None: self._i = 1 else: self._i = exp_tr # see bottom of file for formulas self.polynomial_r = [1, 0, 0, 0, 1] # u^4 + 1 self.polynomial_tr_denom = Integer(1) self.polynomial_r_denom = Integer(1) self.polynomial_p_denom = Integer(4) self.polynomial_y_denom = Integer(2) self.polynomial_c_denom = Integer(4) self.polynomial_lamb = [0,0,1] # T^2 mod r self.polynomial_lamb_denom = 1 if self._i == 1: self.polynomial_p = [ht**2+4*hy**2+2*ht+1, 2*ht+2, 4*hy+1, -4*hy, 2*ht**2+8*hy**2+2*ht+1, 2*ht-2, 4*hy+1, -4*hy, ht**2+4*hy**2] self.polynomial_y = [2*hy, 0, 1, -1, 2*hy] # valid only when T is even self.polynomial_tr = [ht+1, 1, 0, 0, ht] self.polynomial_c = [ht**2+4*hy**2-2*ht+1, 2*ht-2, 4*hy+1, -4*hy, ht**2+4*hy**2] self.polynomial_beta = [ \ ht**5 - 2*ht**4*hy - 8*ht**3*hy**2 - 16*ht**2*hy**3 - 48*ht*hy**4 - 32*hy**5 + 2*ht**4 - 4*ht**3*hy + 32*ht**2*hy**2 - 16*ht*hy**3 - 32*hy**4 + 2*ht**3 - 14*ht**2*hy + 40*ht*hy**2 + 8*hy**3 + 2*ht**2 - 12*ht*hy + ht, \ ht**5 - 6*ht**4*hy + 8*ht**3*hy**2 - 16*ht**2*hy**3 + 16*ht*hy**4 + 32*hy**5 + 4*ht**4 - 12*ht**3*hy - 8*ht**2*hy**2 + 16*ht*hy**3 + 32*hy**4 + 4*ht**3 + 2*ht**2*hy - 16*ht*hy**2 + 8*hy**3 + 2*ht**2 - 4*ht*hy + ht, \ -ht**5 + 6*ht**4*hy - 8*ht**3*hy**2 + 16*ht**2*hy**3 - 16*ht*hy**4 - 32*hy**5 - 2*ht**4 + 12*ht**3*hy - 16*ht**2*hy**2 - 16*ht*hy**3 - 32*hy**4 - 3*ht**3 + 20*ht**2*hy + 4*ht*hy**2 - 16*hy**3 - 4*ht**2 + 12*ht*hy - 2*ht, \ -ht**5 + 2*ht**4*hy + 8*ht**3*hy**2 + 16*ht**2*hy**3 + 48*ht*hy**4 + 32*hy**5 - 2*ht**4 + 32*hy**4 + ht**3 - 16*ht**2*hy + 20*ht*hy**2 + 32*hy**3 + 6*ht**2 - 28*ht*hy + 4*ht, \ ht**5 - 2*ht**4*hy - 8*ht**3*hy**2 - 16*ht**2*hy**3 - 48*ht*hy**4 - 32*hy**5 - 4*ht**3*hy + 80*ht**2*hy**2 - 16*ht*hy**3 - 64*hy**4 + 2*ht**3 - 30*ht**2*hy + 24*ht*hy**2 - 24*hy**3 + 16*ht*hy - 3*ht, \ ht**5 - 6*ht**4*hy + 8*ht**3*hy**2 - 16*ht**2*hy**3 + 16*ht*hy**4 + 32*hy**5 + 2*ht**4 - 4*ht**3*hy - 8*ht**2*hy**2 + 48*ht*hy**3 + 64*hy**4 + 10*ht**2*hy - 48*ht*hy**2 + 8*hy**3 + ht, \ -ht**5 + 6*ht**4*hy - 8*ht**3*hy**2 + 16*ht**2*hy**3 - 16*ht*hy**4 - 32*hy**5 - 2*ht**4 + 12*ht**3*hy - 16*ht**2*hy**2 - 16*ht*hy**3 - 32*hy**4 - ht**3 + 28*ht*hy**2 - 4*ht*hy, \ -ht**5 + 2*ht**4*hy + 8*ht**3*hy**2 + 16*ht**2*hy**3 + 48*ht*hy**4 + 32*hy**5 - 8*ht**3*hy - 32*ht*hy**3 + ht**3 + 4*ht*hy**2] self.polynomial_beta_denom = 2*ht**4 - 16*ht**3*hy + 16*ht**2*hy**2 + 64*ht*hy**3 + 32*hy**4 + 6*ht**3 - 28*ht**2*hy - 8*ht*hy**2 + 16*hy**3 + 6*ht**2 - 12*ht*hy + 2*ht elif self._i == 5: self.polynomial_p = [ht**2+4*hy**2+2*ht+1, -2*ht-2, 4*hy+1, 4*hy, 2*ht**2+8*hy**2+2*ht+1, -2*ht+2, 4*hy+1, 4*hy, ht**2+4*hy**2] self.polynomial_y = [2*hy, 0, 1, 1, 2*hy] self.polynomial_tr = [ht+1, -1, 0, 0, ht] self.polynomial_c = [ht**2+4*hy**2-2*ht+1, -2*ht+2, 4*hy+1, 4*hy, ht**2+4*hy**2] self.polynomial_beta = [\ -ht**5 + 2*ht**4*hy + 8*ht**3*hy**2 + 16*ht**2*hy**3 + 48*ht*hy**4 + 32*hy**5 - 2*ht**4 + 4*ht**3*hy - 32*ht**2*hy**2 + 16*ht*hy**3 + 32*hy**4 - 2*ht**3 + 14*ht**2*hy - 40*ht*hy**2 - 8*hy**3 - 2*ht**2 + 12*ht*hy - ht, \ ht**5 - 6*ht**4*hy + 8*ht**3*hy**2 - 16*ht**2*hy**3 + 16*ht*hy**4 + 32*hy**5 + 4*ht**4 - 12*ht**3*hy - 8*ht**2*hy**2 + 16*ht*hy**3 + 32*hy**4 + 4*ht**3 + 2*ht**2*hy - 16*ht*hy**2 + 8*hy**3 + 2*ht**2 - 4*ht*hy + ht, \ ht**5 - 6*ht**4*hy + 8*ht**3*hy**2 - 16*ht**2*hy**3 + 16*ht*hy**4 + 32*hy**5 + 2*ht**4 - 12*ht**3*hy + 16*ht**2*hy**2 + 16*ht*hy**3 + 32*hy**4 + 3*ht**3 - 20*ht**2*hy - 4*ht*hy**2 + 16*hy**3 + 4*ht**2 - 12*ht*hy + 2*ht, \ -ht**5 + 2*ht**4*hy + 8*ht**3*hy**2 + 16*ht**2*hy**3 + 48*ht*hy**4 + 32*hy**5 - 2*ht**4 + 32*hy**4 + ht**3 - 16*ht**2*hy + 20*ht*hy**2 + 32*hy**3 + 6*ht**2 - 28*ht*hy + 4*ht, \ -ht**5 + 2*ht**4*hy + 8*ht**3*hy**2 + 16*ht**2*hy**3 + 48*ht*hy**4 + 32*hy**5 + 4*ht**3*hy - 80*ht**2*hy**2 + 16*ht*hy**3 + 64*hy**4 - 2*ht**3 + 30*ht**2*hy - 24*ht*hy**2 + 24*hy**3 - 16*ht*hy + 3*ht, \ ht**5 - 6*ht**4*hy + 8*ht**3*hy**2 - 16*ht**2*hy**3 + 16*ht*hy**4 + 32*hy**5 + 2*ht**4 - 4*ht**3*hy - 8*ht**2*hy**2 + 48*ht*hy**3 + 64*hy**4 + 10*ht**2*hy - 48*ht*hy**2 + 8*hy**3 + ht, \ ht**5 - 6*ht**4*hy + 8*ht**3*hy**2 - 16*ht**2*hy**3 + 16*ht*hy**4 + 32*hy**5 + 2*ht**4 - 12*ht**3*hy + 16*ht**2*hy**2 + 16*ht*hy**3 + 32*hy**4 + ht**3 - 28*ht*hy**2 + 4*ht*hy, \ -ht**5 + 2*ht**4*hy + 8*ht**3*hy**2 + 16*ht**2*hy**3 + 48*ht*hy**4 + 32*hy**5 - 8*ht**3*hy - 32*ht*hy**3 + ht**3 + 4*ht*hy**2] self.polynomial_beta_denom = 2*ht**4 - 16*ht**3*hy + 16*ht**2*hy**2 + 64*ht*hy**3 + 32*hy**4 + 6*ht**3 - 28*ht**2*hy - 8*ht*hy**2 + 16*hy**3 + 6*ht**2 - 12*ht*hy + 2*ht else: print("Sorry CocksPinch8 is not implemented for other parameters i except tr = T^i+1 where i=1,5") raise ValueError("Error i={}, sorry but only i=1 and i=5 are implemented, where tr = T^i+1 mod r".format(self._i)) _r = sum([Integer(self.polynomial_r [_i])*self._u**_i for _i in range(len(self.polynomial_r ))])//Integer(self.polynomial_r_denom) _tr= sum([Integer(self.polynomial_tr[_i])*self._u**_i for _i in range(len(self.polynomial_tr))])//Integer(self.polynomial_tr_denom) _y = sum([Integer(self.polynomial_y [_i])*self._u**_i for _i in range(len(self.polynomial_y ))])//Integer(self.polynomial_y_denom) _p = sum([Integer(self.polynomial_p [_i])*self._u**_i for _i in range(len(self.polynomial_p ))])//Integer(self.polynomial_p_denom) _c = sum([Integer(self.polynomial_c [_i])*self._u**_i for _i in range(len(self.polynomial_c ))])//Integer(self.polynomial_c_denom) _beta = sum([Integer(self.polynomial_beta[_i])*self._u**_i for _i in range(len(self.polynomial_beta))])/Integer(self.polynomial_beta_denom) _lamb = sum([Integer(self.polynomial_lamb[_i])*self._u**_i for _i in range(len(self.polynomial_lamb))])#//Integer(self.polynomial_lamb_denom) if r != None: self._r = Integer(r) if _r != self._r: raise ValueError("Error: r does not match the polynomial r(T)\nr = {:= #x} #(input)\nr = {:= #x} # from polynomial r(T,ht,hy) = {}".format(r,_r,self.polynomial_r)) #else: # print("valid r") else: self._r = _r assert (self._r).is_prime() if tr != None: self._tr = Integer(tr) if _tr != self._tr: raise ValueError("Error: tr does not match the polynomial tr(T,ht)\ntr = {:= #x} #(input)\ntr = {:= #x} # from polynomial tr(T,ht,hy) = {}".format(tr,_tr,self.polynomial_tr)) #else: # print("valid tr") else: self._tr = _tr if y != None: self._y = Integer(y) if _y != self._y: raise ValueError("Error: y does not match the polynomial y(T)\ny = {:= #x} #(input)\ny = {:= #x} # from polynomial y(T,ht,hy) = {}".format(y,_y,self.polynomial_y)) #else: # print("valid y") else: self._y = _y if p != None: self._p = Integer(p) if _p != self._p: raise ValueError("Error: p does not match the polynomial P(T,ht,hy)\np = {:= #x} #(input)\np = {:= #x} # from polynomial p(T,ht,hy) = {}/{}".format(p,_p,self.polynomial_p,self.polynomial_p_denom)) else: self._p = _p assert (self._p).is_prime() assert ((self._p+1-self._tr) % self._r) == 0 self._pnbits = self._p.nbits() if c != None: self._c = Integer(c) if _c != self._c: raise ValueError("Error: c does not match the polynomial c(T,ht,hy)\nc = {:= #x} #(input)\nc = {:= #x} # from polynomial c(T,ht,hy)".format(c,_c)) else: self._c = _c assert ((self._p+1-self._tr) == self._c * self._r) if beta != None: self._beta = Integer(beta) if _beta != self._beta: raise ValueError("Error: beta does not match the polynomial beta(T,ht,hy)\nbeta = {:= #x} #(input)\nbeta = {:= #x} # from polynomial beta(T,ht,hy)".format(beta,_beta)) else: self._beta = _beta #print("self._beta = {}".format(self._beta)) if lamb != None: self._lamb = Integer(lamb) if _lamb != self._lamb: raise ValueError("Error: lamb does not match the polynomial lamb(T,ht,hy)\nlamb = {:= #x} #(input)\nlamb = {:= #x} # from polynomial lamb(T,ht,hy)".format(lamb,_lamb)) else: self._lamb = _lamb #print("self._lamb = {:= #x}".format(self._lamb)) try: self._Fp = FiniteField(self._p) except ValueError as err: print("ValueError creating Fp: {}".format(err)) raise except: print("Error creating Fp") raise self._Fpz = PolynomialRing(self._Fp, names=('z',)) (self._z,) = self._Fpz._first_ngens(1) num_b = (self._Fp)(self._beta.numer()) den_b = (self._Fp)(self._beta.denom()) self._beta = Integer(num_b/den_b) if ((self._beta**2 + 1) % self._p) != 0: raise ValueError("Error beta^2 + 1 != 0 mod p") #else: # print("valid beta") if ((self._lamb**2 + 1) % self._r) != 0: raise ValueError("Error lamb^2 + 1 != 0 mod r") #else: # print("valid lamb") if a != None: try: a = Integer(a) except: raise self._a = a self._ap = self._Fp(a) else: self._a, self._ap = get_curve_parameter_a_j1728(self._tr, self._y, self._beta, self._p, self._Fp) self._bp = self._Fp(0) # second curve parameter is 0 because j=1728 # Now self._a is such that E: y^2 = x^3 + a*x has order r 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,0]) 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 self._order_checked = super(CocksPinch8,self).order() expected_order = (self._p+1-self._tr) if self._order_checked != expected_order: #if (self._order_checked % self._r) != 0: print("Error, wrong order") print("order recomputed: {:= #x}".format(self._order_checked)) print("expected order (p+1-tr): {:= #x}".format(expected_order)) print("quadratic twist (p+1+tr): {:= #x}".format(expected_order+2*self._tr)) raise ValueError("Wrong curve order: this one is a twist") # computes a generator self._G = get_curve_generator_order_r_j1728(self) self._Gx = self._G[0] self._Gy = self._G[1] # adjust beta and lamb according to the curve # do we have (-x,y*beta) = lamb*(x,y)? if self([-self._Gx, self._Gy*self._beta]) != self._lamb*self._G: print("adjusting beta, lambda") if self([-self._Gx, self._Gy*(-self._beta)]) == self._lamb*self._G: self._beta = self._p-self._beta print("beta -> -beta") elif self([-self._Gx, self._Gy*self._beta]) == (-self._lamb)*self._G: self._lamb = -self._lamb print("lamb -> -lamb") elif self([-self._Gx, self._Gy*(-self._beta)]) == (-self._lamb)*self._G: self._beta = self._p-self._beta self._lamb = -self._lamb print("lamb -> -lamb") print("beta -> -beta") else: raise ValueError("Error while adjusting beta, lamb: compatibility not found") def _repr_(self): return "CP8 p"+str(self._pnbits)+" (modified Cocks-Pinch k=8) curve with seed "+str(self._u)+"\n"+super(CocksPinch8,self)._repr_() def u(self): return self._u def T(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 D(self): return self._D 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 self._beta def lamb(self): return self._lamb 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) """ R = QQ['T','ht','hy'] T,ht,hy = R.gens() D=4 i=5 r = T^4+1 t0 = (T^i+1) % r print("t0 = {}".format(t0)) t0 = -T+1 tr = t0+ ht*r print("tr = {}".format(tr)) # tr = T^4*ht - T + ht + 1 print tr.polynomial(T).list() # [ht + 1, -1, 0, 0, ht] inv_sqrt_D = (-T^2)/2 assert ((-inv_sqrt_D^2*D) % r) == 1 y0 = ((tr-2)*inv_sqrt_D) % r y = y0+hy*r print("y0 = 1/2*({})".format(2*y0)) print (2*y0).polynomial(T).list(), '/2' print("y = 1/2*({})".format(2*y)) print (2*y).polynomial(T).list(), '/2' #y0 = 1/2*(T^3 + T^2) #[0, 0, 1, 1] /2 #y = 1/2*(2*T^4*hy + T^3 + T^2 + 2*hy) #[2*hy, 0, 1, 1, 2*hy] /2 p = (tr^2 + D*y^2)/4 print("p = 1/4*({})".format(4*p)) print (4*p).polynomial(T).list(), '/4' #[ht^2 + 4*hy^2 + 2*ht + 1, -2*ht - 2, 4*hy + 1, 4*hy, 2*ht^2 + 8*hy^2 + 2*ht + 1, -2*ht + 2, 4*hy + 1, 4*hy, ht^2 + 4*hy^2] /4 print 'p+1-tr =', (p+1-tr).factor() c = (p+1-tr) // r print (4*c).polynomial(T).list(), '/4' # [ht^2 + 4*hy^2 - 2*ht + 1, -2*ht + 2, 4*hy + 1, 4*hy, ht^2 + 4*hy^2] /4 # now see Magma file for computing inverse mod of multivariate polynomials, thanks to Pierre-Jean! den = 2*ht^4 - 16*ht^3*hy + 6*ht^3 + 16*ht^2*hy^2 - 28*ht^2*hy + 6*ht^2 + 64*ht*hy^3 - 8*ht*hy^2 - 12*ht*hy + 2*ht + 32*hy^4 + 16*hy^3 num = (-ht^5 + 2*ht^4*hy + 8*ht^3*hy^2 - 8*ht^3*hy + ht^3 + 16*ht^2*hy^3 + 48*ht*hy^4 - 32*ht*hy^3 + 4*ht*hy^2 + 32*hy^5)*T^7 + (ht^5 - 6*ht^4*hy + 2*ht^4 + 8*ht^3*hy^2 - 12*ht^3*hy + ht^3 - 16*ht^2*hy^3 + 16*ht^2*hy^2 + 16*ht*hy^4 + 16*ht*hy^3 - 28*ht*hy^2 + 4*ht*hy + 32*hy^5 + 32*hy^4)*T^6 + (ht^5 - 6*ht^4*hy + 2*ht^4 + 8*ht^3*hy^2 - 4*ht^3*hy - 16*ht^2*hy^3 - 8*ht^2*hy^2 + 10*ht^2*hy + 16*ht*hy^4 + 48*ht*hy^3 - 48*ht*hy^2 + ht + 32*hy^5 + 64*hy^4 + 8*hy^3)*T^5 + (-ht^5 + 2*ht^4*hy + 8*ht^3*hy^2 + 4*ht^3*hy - 2*ht^3 + 16*ht^2*hy^3 - 80*ht^2*hy^2 + 30*ht^2*hy + 48*ht*hy^4 + 16*ht*hy^3 - 24*ht*hy^2 - 16*ht*hy + 3*ht + 32*hy^5 + 64*hy^4 + 24*hy^3)*T^4 + (-ht^5 + 2*ht^4*hy - 2*ht^4 + 8*ht^3*hy^2 + ht^3 + 16*ht^2*hy^3 - 16*ht^2*hy + 6*ht^2 + 48*ht*hy^4 + 20*ht*hy^2 - 28*ht*hy + 4*ht + 32*hy^5 + 32*hy^4 + 32*hy^3)*T^3 + (ht^5 - 6*ht^4*hy + 2*ht^4 + 8*ht^3*hy^2 - 12*ht^3*hy + 3*ht^3 - 16*ht^2*hy^3 + 16*ht^2*hy^2 - 20*ht^2*hy + 4*ht^2 + 16*ht*hy^4 + 16*ht*hy^3 - 4*ht*hy^2 - 12*ht*hy + 2*ht + 32*hy^5 + 32*hy^4 + 16*hy^3)*T^2 + (ht^5 - 6*ht^4*hy + 4*ht^4 + 8*ht^3*hy^2 - 12*ht^3*hy + 4*ht^3 - 16*ht^2*hy^3 - 8*ht^2*hy^2 + 2*ht^2*hy + 2*ht^2 + 16*ht*hy^4 + 16*ht*hy^3 - 16*ht*hy^2 - 4*ht*hy + ht + 32*hy^5 + 32*hy^4 + 8*hy^3)*T - ht^5 + 2*ht^4*hy - 2*ht^4 + 8*ht^3*hy^2 + 4*ht^3*hy - 2*ht^3 + 16*ht^2*hy^3 - 32*ht^2*hy^2 + 14*ht^2*hy - 2*ht^2 + 48*ht*hy^4 + 16*ht*hy^3 - 40*ht*hy^2 + 12*ht*hy - ht + 32*hy^5 + 32*hy^4 - 8*hy^3 print den.polynomial(T).list() print num.polynomial(T).list() [2*ht^4 - 16*ht^3*hy + 16*ht^2*hy^2 + 64*ht*hy^3 + 32*hy^4 + 6*ht^3 - 28*ht^2*hy - 8*ht*hy^2 + 16*hy^3 + 6*ht^2 - 12*ht*hy + 2*ht] [-ht^5 + 2*ht^4*hy + 8*ht^3*hy^2 + 16*ht^2*hy^3 + 48*ht*hy^4 + 32*hy^5 - 2*ht^4 + 4*ht^3*hy - 32*ht^2*hy^2 + 16*ht*hy^3 + 32*hy^4 - 2*ht^3 + 14*ht^2*hy - 40*ht*hy^2 - 8*hy^3 - 2*ht^2 + 12*ht*hy - ht, \ ht^5 - 6*ht^4*hy + 8*ht^3*hy^2 - 16*ht^2*hy^3 + 16*ht*hy^4 + 32*hy^5 + 4*ht^4 - 12*ht^3*hy - 8*ht^2*hy^2 + 16*ht*hy^3 + 32*hy^4 + 4*ht^3 + 2*ht^2*hy - 16*ht*hy^2 + 8*hy^3 + 2*ht^2 - 4*ht*hy + ht, \ ht^5 - 6*ht^4*hy + 8*ht^3*hy^2 - 16*ht^2*hy^3 + 16*ht*hy^4 + 32*hy^5 + 2*ht^4 - 12*ht^3*hy + 16*ht^2*hy^2 + 16*ht*hy^3 + 32*hy^4 + 3*ht^3 - 20*ht^2*hy - 4*ht*hy^2 + 16*hy^3 + 4*ht^2 - 12*ht*hy + 2*ht, \ -ht^5 + 2*ht^4*hy + 8*ht^3*hy^2 + 16*ht^2*hy^3 + 48*ht*hy^4 + 32*hy^5 - 2*ht^4 + 32*hy^4 + ht^3 - 16*ht^2*hy + 20*ht*hy^2 + 32*hy^3 + 6*ht^2 - 28*ht*hy + 4*ht, \ -ht^5 + 2*ht^4*hy + 8*ht^3*hy^2 + 16*ht^2*hy^3 + 48*ht*hy^4 + 32*hy^5 + 4*ht^3*hy - 80*ht^2*hy^2 + 16*ht*hy^3 + 64*hy^4 - 2*ht^3 + 30*ht^2*hy - 24*ht*hy^2 + 24*hy^3 - 16*ht*hy + 3*ht, \ ht^5 - 6*ht^4*hy + 8*ht^3*hy^2 - 16*ht^2*hy^3 + 16*ht*hy^4 + 32*hy^5 + 2*ht^4 - 4*ht^3*hy - 8*ht^2*hy^2 + 48*ht*hy^3 + 64*hy^4 + 10*ht^2*hy - 48*ht*hy^2 + 8*hy^3 + ht, \ ht^5 - 6*ht^4*hy + 8*ht^3*hy^2 - 16*ht^2*hy^3 + 16*ht*hy^4 + 32*hy^5 + 2*ht^4 - 12*ht^3*hy + 16*ht^2*hy^2 + 16*ht*hy^3 + 32*hy^4 + ht^3 - 28*ht*hy^2 + 4*ht*hy, \ -ht^5 + 2*ht^4*hy + 8*ht^3*hy^2 + 16*ht^2*hy^3 + 48*ht*hy^4 + 32*hy^5 - 8*ht^3*hy - 32*ht*hy^3 + ht^3 + 4*ht*hy^2] i=1 r = T^4+1 t0 = T+1 tr = t0+ ht*r print tr.polynomial(T).list() y0 = -(T-1)*(T^2)/2 print y0.polynomial(T).list() y = y0+hy*r print y.polynomial(T).list() p = (tr^2 + D*y^2)/4 print p.polynomial(T).list() print (4*p).polynomial(T).list() c = (p+1-tr) // r print c.polynomial(T).list() print (4*c).polynomial(T).list() h_t = ht h_y = -hy num_beta = (-h_t^5 - 2*h_t^4*h_y + 8*h_t^3*h_y^2 + 8*h_t^3*h_y + h_t^3 - 16*h_t^2*h_y^3 + 48*h_t*h_y^4 + 32*h_t*h_y^3 + 4*h_t*h_y^2 - 32*h_y^5)*T^7 + (-h_t^5 - 6*h_t^4*h_y - 2*h_t^4 - 8*h_t^3*h_y^2 - 12*h_t^3*h_y - h_t^3 - 16*h_t^2*h_y^3 - 16*h_t^2*h_y^2 - 16*h_t*h_y^4 + 16*h_t*h_y^3 + 28*h_t*h_y^2 + 4*h_t*h_y + 32*h_y^5 - 32*h_y^4)*T^6 + (h_t^5 + 6*h_t^4*h_y + 2*h_t^4 + 8*h_t^3*h_y^2 + 4*h_t^3*h_y + 16*h_t^2*h_y^3 - 8*h_t^2*h_y^2 - 10*h_t^2*h_y + 16*h_t*h_y^4 - 48*h_t*h_y^3 - 48*h_t*h_y^2 + h_t - 32*h_y^5 + 64*h_y^4 - 8*h_y^3)*T^5 + (h_t^5 + 2*h_t^4*h_y - 8*h_t^3*h_y^2 + 4*h_t^3*h_y + 2*h_t^3 + 16*h_t^2*h_y^3 + 80*h_t^2*h_y^2 + 30*h_t^2*h_y - 48*h_t*h_y^4 + 16*h_t*h_y^3 + 24*h_t*h_y^2 - 16*h_t*h_y - 3*h_t + 32*h_y^5 - 64*h_y^4 + 24*h_y^3)*T^4 + (-h_t^5 - 2*h_t^4*h_y - 2*h_t^4 + 8*h_t^3*h_y^2 + h_t^3 - 16*h_t^2*h_y^3 + 16*h_t^2*h_y + 6*h_t^2 + 48*h_t*h_y^4 + 20*h_t*h_y^2 + 28*h_t*h_y + 4*h_t - 32*h_y^5 + 32*h_y^4 - 32*h_y^3)*T^3 + (-h_t^5 - 6*h_t^4*h_y - 2*h_t^4 - 8*h_t^3*h_y^2 - 12*h_t^3*h_y - 3*h_t^3 - 16*h_t^2*h_y^3 - 16*h_t^2*h_y^2 - 20*h_t^2*h_y - 4*h_t^2 - 16*h_t*h_y^4 + 16*h_t*h_y^3 + 4*h_t*h_y^2 - 12*h_t*h_y - 2*h_t + 32*h_y^5 - 32*h_y^4 + 16*h_y^3)*T^2 + (h_t^5 + 6*h_t^4*h_y + 4*h_t^4 + 8*h_t^3*h_y^2 + 12*h_t^3*h_y + 4*h_t^3 + 16*h_t^2*h_y^3 - 8*h_t^2*h_y^2 - 2*h_t^2*h_y + 2*h_t^2 + 16*h_t*h_y^4 - 16*h_t*h_y^3 - 16*h_t*h_y^2 + 4*h_t*h_y + h_t - 32*h_y^5 + 32*h_y^4 - 8*h_y^3)*T + h_t^5 + 2*h_t^4*h_y + 2*h_t^4 - 8*h_t^3*h_y^2 + 4*h_t^3*h_y + 2*h_t^3 + 16*h_t^2*h_y^3 + 32*h_t^2*h_y^2 + 14*h_t^2*h_y + 2*h_t^2 - 48*h_t*h_y^4 + 16*h_t*h_y^3 + 40*h_t*h_y^2 + 12*h_t*h_y + h_t + 32*h_y^5 - 32*h_y^4 - 8*h_y^3 den_beta = 2*(h_t*(h_t+1)^3 + 2*h_t*h_y*((h_t+1)*(4*h_t+3) + 2*h_y*(2*h_t-8*h_y-1)) + 8*h_y^3*(2*h_y-1)) r = T^4 + 1 -> poly [1,0,0,0,1] if i==1: t0 = T+1 -> poly [1,1] tr = t0+ ht*r = 1+ht+T+ht*T^4 -> [1+ht,1,0,0,ht] y0 = (t0-2)/sqrt(-4) - r (why -r ?) y0 = -(T-1)*(T^2)/2 = (-T^3+T^2)/2 -> [0,0,1,-1]/2 y = y0+hy*r = (-T^3+T^2)/2 +hy*(1+T^4) = hy+T^2/2-T^3/2+hy*T^4 -> [hy,0,-1/2,1/2,hy] p = (tr^2 + D*y^2)/4 = (1+ht+T+ht*T^4)^2/4 + (hy+T^2/2-T^3/2+hy*T^4)^2 """