"Fossies" - the Fresh Open Source Software Archive

Member "NZMATH-1.2.0/nzmath/elliptic.py" (19 Nov 2012, 53620 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "elliptic.py" see the Fossies "Dox" file reference documentation.

    1 """ Elliptic Curves over finite field.
    2 """
    3 
    4 from __future__ import division
    5 import logging
    6 import random
    7 
    8 import nzmath.poly.uniutil as uniutil
    9 import nzmath.poly.multiutil as multiutil
   10 import nzmath.poly.termorder as termorder
   11 import nzmath.arith1 as arith1
   12 import nzmath.bigrandom as bigrandom
   13 import nzmath.gcd as gcd
   14 import nzmath.prime as prime
   15 import nzmath.rational as rational
   16 import nzmath.ring as ring
   17 import nzmath.factor.methods as factor_methods
   18 import nzmath.factor.misc as factor_misc
   19 import nzmath.finitefield as finitefield
   20 import nzmath.compatibility
   21 
   22 _log = logging.getLogger('nzmath.elliptic')
   23 
   24 # symbol aliases
   25 uniutil.special_ring_table[finitefield.FinitePrimeField] = uniutil.FinitePrimeFieldPolynomial
   26 MultiVarPolynomial = multiutil.MultiVariableSparsePolynomial
   27 
   28 # polynomial wrapper
   29 def _UniVarPolynomial(dict,coeffring=None):
   30     return uniutil.polynomial(dict,coeffring)
   31 
   32 # string format wrapper
   33 def _strUniPoly(poly, symbol="x", asc=True):
   34     """return format string of UniutilPolynomial"""
   35     return termorder.ascending_order.format(poly, symbol, asc)
   36 
   37 def _strMultiPoly(poly, symbol=["x","y"], asc=False):
   38     """return format string of MultiVarPolynomial for EC"""
   39     def Weierstrasscmp(left, right):
   40         sum_left, sum_right = sum(left), sum(right)
   41         if sum_left != sum_right:
   42             return -cmp(sum_left, sum_right)
   43         return -cmp(right, left)
   44 
   45     return termorder.MultivarTermOrder(Weierstrasscmp).format(MultiVarPolynomial(poly,symbol), symbol, asc)
   46 
   47 
   48 def _PolyMod(f, g):
   49     """
   50     return f (mod g)
   51     """
   52     return f % g
   53 
   54 def _PolyGCD(f, g):
   55     # other cases
   56     return f.gcd(g)
   57 
   58 def _PolyPow(f, d, g):
   59     """
   60     returns (f^d)%g
   61     """
   62     return g.mod_pow(f, d)
   63 
   64 def _PolyMulRed(multipliees, poly):
   65     """
   66     multipliees[*] is (OneSparsePoly,int,long)
   67     poly is OneSparsePoly
   68     """
   69     if poly.degree() < 1:
   70         return poly.getRing().zero
   71     product = multipliees.pop()
   72     for factor in multipliees:
   73         #if factor.degree() >= poly.degree():
   74         #factor = PolyMod(factor, poly)
   75         #if factor == 0:
   76         #    return 0
   77         product = product * factor
   78         if product.degree() >= poly.degree():
   79             product = _PolyMod(product, poly)
   80             if not product:
   81                 break
   82     return product
   83 
   84 
   85 class ECGeneric:
   86     """
   87     Definition of Elliptic curves over generic field.
   88     this class is fundamental class, normally only called sub class.
   89     """
   90     def __init__(self, coefficient, basefield=None):
   91         """
   92         Initialize an elliptic curve. If coefficient has 5 elements,
   93         it represents E:y**2+a1*x*y+a3*y=x**3+a2*x**2+a4*x+a6 or 2
   94         elements, E:y*2=x*3+a*x+b.
   95         """
   96 
   97         try:
   98             character = basefield.getCharacteristic()
   99             self.basefield = basefield
  100         except:
  101             # backward compatibility support
  102             if isinstance(basefield, rational.RationalField) or (not basefield):
  103                 character = 0
  104                 self.basefield = rational.theRationalField
  105             elif isinstance(basefield, (int,long)):
  106                 character = basefield
  107                 if character == 1 or character < 0:
  108                     raise ValueError("basefield characteristic must be 0 or prime.")
  109                 self.basefield = finitefield.FinitePrimeField.getInstance(character)
  110             else:
  111                 raise ValueError("basefield must be FiniteField.")
  112 
  113         self.ch = character
  114         self.infpoint = [self.basefield.zero]
  115         if isinstance(coefficient, list):
  116             self.coefficient = coefficient
  117             if self.ch == 0:
  118                 if len(self) == 5:
  119                     self.a1 = self.coefficient[0]
  120                     self.a2 = self.coefficient[1]
  121                     self.a3 = self.coefficient[2]
  122                     self.a4 = self.coefficient[3]
  123                     self.a6 = self.coefficient[4]
  124                     self.b2 = self.a1**2+4*self.a2
  125                     self.b4 = self.a1*self.a3+2*self.a4
  126                     self.b6 = self.a3**2+4*self.a6
  127                     self.b8 = self.a1**2*self.a6+4*self.a2*self.a6-self.a1*self.a3*self.a4+self.a2*self.a3**2-self.a4**2
  128                     self.c4 = self.b2**2-24*self.b4
  129                     self.c6 = -self.b2**3+36*self.b2*self.b4-216*self.b6
  130                     self.disc = -self.b2**2*self.b8-8*self.b4**3-27*self.b6**2+9*self.b2*self.b4*self.b6
  131                 elif len(self) == 2:
  132                     self.a = self.coefficient[0]
  133                     self.b = self.coefficient[1]
  134                     self.a1 = 0
  135                     self.a2 = 0
  136                     self.a3 = 0
  137                     self.a4 = self.coefficient[0]
  138                     self.a6 = self.coefficient[1]
  139                     self.b2 = 0
  140                     self.b4 = 2*self.a
  141                     self.b6 = 4*self.b
  142                     self.b8 = -self.a**2
  143                     self.c4 = -48*self.a
  144                     self.c6 = -864*self.b
  145                     self.disc = (self.c4**3-self.c6**2)/1728
  146                 else:
  147                     raise ValueError("coefficient is less or more, can't defined EC.")
  148                 if self.disc == 0:
  149                     raise ValueError("this curve is singular.")
  150                 self.j = (self.c4**3)/self.disc
  151                 self.cubic = _UniVarPolynomial({0:self.a6, 1:self.a4,
  152                                                3:self.basefield.one},
  153                                               self.basefield)
  154             else:
  155                 pass # support for subclass
  156         else:
  157             raise ValueError("parameters must be (coefficient, basefield)")
  158     def __len__(self):
  159         return len(self.coefficient)
  160 
  161     def __repr__(self):
  162         """
  163         return represantation form.
  164         only defined generic representation.
  165         """
  166         return "EC(%s, %s)" % (str(self.coefficient), repr(self.basefield))
  167 
  168     def __str__(self):
  169         dictL = {(0, 2):self.basefield.one}
  170         dictR = {3:self.basefield.one}
  171         if self.a1:
  172             dictL[(1, 1)] = self.a1
  173         if self.a2:
  174             dictR[2] = self.a2
  175         if self.a3:
  176             dictL[(0, 1)] = self.a3
  177         if self.a4:
  178             dictR[1] = self.a4
  179         if self.a6:
  180             dictR[0] = self.a6
  181         return ("%s = %s" %
  182                 (str(_strMultiPoly(dictL)),
  183                  str(_strUniPoly(_UniVarPolynomial(dictR, self.basefield)))))
  184 
  185     def simple(self):
  186         """
  187         Transform
  188           E:y^2+a1*x*y+a3*y = x^3+a2*x^2+a4*x+a6
  189         to
  190           E':Y^2 = X^3+(-27*c4)*X+(-54*c6),
  191         if ch is not 2 or 3
  192         """
  193         if len(self) == 2 or not self.a1 and not self.a2 and not self.a3:
  194             return self
  195         elif self.ch == 2 or self.ch == 3:
  196             return self
  197         else:
  198             return EC([-27*self.c4, -54*self.c6], self.basefield)
  199 
  200     def changeCurve(self, V):
  201         """
  202         this transforms E to E' by V=[u,r,s,t]
  203           x->u^2*x'+r,y->u^3*y'+s*u^2*x'+t
  204         """
  205         # FIXME: is "(not V[0] != 0)" a correct condition?
  206         if (not isinstance(V, list)) or (not len(V) == 4) or (not V[0] != 0):
  207             raise ValueError("([u, r, s, t]) with u must not be 0.")
  208 
  209         # FIXME: this is dependent rational.Rational
  210         if self.ch == 0:
  211             return EC([rational.Rational(self.a1+2*V[2], V[0]),
  212                        rational.Rational(self.a2-V[2]*self.a1+3*V[1]-V[2]**2, V[0]**2),
  213                        rational.Rational(self.a3+V[1]*self.a1+2*V[3], V[0]**3),
  214                        rational.Rational(self.a4-V[2]*self.a3+2*V[1]*self.a2-(V[3]+V[1]*V[2])*self.a1+3*V[1]**2-2*V[2]*V[3], V[0]**4),
  215                        rational.Rational(self.a6+V[1]*self.a4+V[1]**2*self.a2+V[1]**3-V[3]*self.a3-V[3]**2-V[1]*V[3]*self.a1, V[0]**6)])
  216         else:
  217             for v in V:
  218                 if not isinstance(v, (int, long)) and not (v in self.basefield):
  219                     raise ValueError("transform V must be integer sequence.")
  220             v = self.basefield.createElement(V[0]).inverse()
  221             return EC([(self.a1+2*V[2])*v,
  222                        ((self.a2-V[2]*self.a1+3*V[1]-V[2]**2)*v**2),
  223                        ((self.a3+V[1]*self.a1+2*V[3])*v**3),
  224                        ((self.a4-V[2]*self.a3+2*V[1]*self.a2-(V[3]+V[1]*V[2])*self.a1+3*V[1]**2-2*V[2]*V[3])*v**4),
  225                        ((self.a6+V[1]*self.a4+V[1]**2*self.a2+V[1]**3-V[3]*self.a3-V[3]**2-V[1]*V[3]*self.a1)*v**6)], self.basefield)
  226 
  227     def changePoint(self, P, V):
  228         """
  229         this transforms P to P' by V=[u,r,s,t] ; x->u^2*x'+r,y->u^3*y'+s*u^2*x'+t
  230         """
  231         if (not (isinstance(P, list) and isinstance(V, list))) or \
  232                not (len(P) == 2 and len(V) == 4 and V[0] != 0):
  233             raise ValueError("(P,V) must be ([px, py], [u, r, s, t]) with u != 0.")
  234 
  235         if self.ch == 0:
  236             Q0 = rational.IntegerIfIntOrLong(P[0]-V[1])/rational.IntegerIfIntOrLong(V[0]**2)
  237             Q1 = rational.IntegerIfIntOrLong(P[1]-V[2]*(P[0]-V[1])-V[3])/rational.IntegerIfIntOrLong(V[0]**3)
  238         else:
  239             v = self.basefield.createElement(V[0]).inverse()
  240             Q0 = ((P[0]-V[1])*v**2)
  241             Q1 = ((P[1]-V[2]*(P[0]-V[1])-V[3])*v**3)
  242         Q = [Q0, Q1]
  243         return Q
  244 
  245     def coordinateY(self, x):
  246         """
  247         this returns the y(P)>0,x(P)==x
  248         """
  249         if self.ch > 2:
  250             y1 = self.a1*x + self.a3
  251             y2 = x**3 + self.a2*x**2 + self.a4*x + self.a6
  252             if len(self) == 2 or (self.a1 == self.a2 == self.a3 == 0):
  253                 if self.basefield.Legendre(y2) == 1:
  254                     return self.basefield.sqrt(y2)
  255                 else:
  256                     return False
  257             else:
  258                 if y1**2 + 4*y2 >= 0:
  259                     d = y1**2 + 4*y2
  260                     return (-y1 - self.basefield.sqrt(d)) // (2*self.basefield.one)
  261                 else:
  262                     return False
  263         elif self.ch == 2:
  264             raise NotImplementedError("This is not implemented.")
  265         y1 = self.a1*x + self.a3
  266         y2 = x**3 + self.a2*x**2 + self.a4*x + self.a6
  267         Y = y1**2 + 4*y2
  268         if Y >= 0:
  269             if isinstance(Y, rational.Rational):
  270                 yn = arith1.issquare(Y.numerator)
  271                 yd = arith1.issquare(Y.denominator)
  272                 if yn and yd:
  273                     return ((-1)*y1 + rational.Rational(yn, yd)) / 2
  274             else:
  275                 Z = arith1.issquare(Y)
  276                 if Z:
  277                     return rational.Rational((-1)*y1 + Z, 2)
  278         return False
  279 
  280     def whetherOn(self, P):
  281         """
  282         Determine whether P is on curve or not.
  283         Return True if P is on, return False otherwise.
  284         """
  285         if isinstance(P, list):
  286             if len(P) == 2:
  287                 if self.ch == 0:
  288                     if P[1]**2+self.a1*P[0]*P[1]+self.a3*P[1] == P[0]**3+self.a2*P[0]**2+self.a4*P[0]+self.a6:
  289                         return True
  290                     else:
  291                         #_log.debug(str(P[1]**2+self.a1*P[0]*P[1]+self.a3*P[1]))
  292                         #_log.debug(str(P[0]**3+self.a2*P[0]**2+self.a4*P[0]+self.a6))
  293                         return False
  294                 else:
  295                     if P[1]**2+self.a1*P[0]*P[1]+self.a3*P[1] == P[0]**3+self.a2*P[0]**2+self.a4*P[0]+self.a6:
  296                         return True
  297                     else:
  298                         return False
  299             elif P == [self.basefield.zero]:
  300                 return True
  301         raise ValueError("point P must be [px, py] or [0].")
  302 
  303     def add(self, P, Q):
  304         """
  305         return point addition P+Q
  306         """
  307         if not (isinstance(P, list) and isinstance(Q, list)):
  308             raise ValueError("point P (resp. Q) must be [px, py] (resp. [qx, qy])")
  309         #if not (self.whetherOn(P) and self.whetherOn(Q)):
  310         #    raise ValueError("either points must not be point on curve.")
  311 
  312         if (P == self.infpoint) and (Q != self.infpoint):
  313             return Q
  314         elif (P != self.infpoint) and (Q == self.infpoint):
  315             return P
  316         elif (P == self.infpoint) and (Q == self.infpoint):
  317             return self.infpoint
  318 
  319         if self.ch == 0:
  320             # FIXME
  321             if P[0] == Q[0]:
  322                 if P[1]+Q[1]+self.a1*Q[0]+self.a3 == 0:
  323                     return self.infpoint
  324                 else:
  325                     s = (3*P[0]**2+2*self.a2*P[0]+self.a4-self.a1*P[1])/(2*P[1]+self.a1*P[0]+self.a3)
  326                     t = (-P[0]**3+self.a4*P[0]+2*self.a6-self.a3*P[1])/(2*P[1]+self.a1*P[0]+self.a3)
  327             else:
  328                 s = (Q[1]-P[1])/(Q[0]-P[0])
  329                 t = (P[1]*Q[0]-Q[1]*P[0])/(Q[0]-P[0])
  330             x3 = s**2+self.a1*s-self.a2-P[0]-Q[0]
  331             y3 = -(s+self.a1)*x3-t-self.a3
  332             R = [x3, y3]
  333             return R
  334         else:
  335             if not (P[0] - Q[0]):
  336                 # FIXME: the condition is P[0] == Q[0] intuitively,
  337                 #        but sometimes there are int vs FiniteFieldElement
  338                 #        comparisons ...
  339                 if not (P[1]+Q[1]+self.a1*Q[0]+self.a3):
  340                     return self.infpoint
  341                 else:
  342                     s = (3*P[0]**2+2*self.a2*P[0]+self.a4-self.a1*P[1])/(2*P[1]+self.a1*P[0]+self.a3)
  343                     t = (-P[0]**3+self.a4*P[0]+2*self.a6-self.a3*P[1])/(2*P[1]+self.a1*P[0]+self.a3)
  344             else:
  345                 s = (Q[1] - P[1]*self.basefield.one) / (Q[0] - P[0])
  346                 t = (P[1]*Q[0] - Q[1]*P[0]*self.basefield.one)/ (Q[0] - P[0])
  347             x3 = s**2+self.a1*s-self.a2-P[0]-Q[0]
  348             y3 = -(s+self.a1)*x3-t-self.a3
  349             R = [x3, y3]
  350             return R
  351 
  352     def sub(self, P, Q):
  353         """
  354         return point subtraction P-Q
  355         """
  356         if not (isinstance(P, list) and isinstance(Q, list)):
  357             raise ValueError("point P (resp. Q) must be [px, py] (resp. [qx, qy])")
  358         #if not (self.whetherOn(P) and self.whetherOn(Q)):
  359         #    raise ValueError("either points must not be point on curve.")
  360 
  361         if (P != self.infpoint) and (Q == self.infpoint):
  362             return P
  363         elif (P == self.infpoint) and (Q == self.infpoint):
  364             return self.infpoint
  365 
  366         x = Q[0]
  367         y = -Q[1]-self.a1*Q[0]-self.a3
  368         R = [x, y]
  369         if (P == self.infpoint) and (Q != self.infpoint):
  370             return R
  371         else:
  372             return self.add(P, R)
  373 
  374     def mul(self, k, P):
  375         """
  376         this returns [k]*P
  377         """
  378         if k >= 0:
  379             l = arith1.expand(k, 2)
  380             Q = self.infpoint
  381             for j in range(len(l)-1, -1, -1):
  382                 Q = self.add(Q, Q)
  383                 if l[j] == 1:
  384                     Q = self.add(Q, P)
  385             return Q
  386         else:
  387             l = arith1.expand(-k, 2)
  388             Q = self.infpoint
  389             for j in range(len(l)-1, -1, -1):
  390                 Q = self.add(Q, Q)
  391                 if l[j] == 1:
  392                     Q = self.add(Q, P)
  393             return self.sub(self.infpoint, Q)
  394 
  395     def divPoly(self, Number=None):
  396         """ Return division polynomial
  397         """
  398         def heart(q):
  399             """
  400             this is for Schoof's, internal function
  401             """
  402             l = []
  403             j = 1
  404             bound = 4 * arith1.floorsqrt(q)
  405             for p in prime.generator():
  406                 if p != q:
  407                     l.append(p)
  408                     j *= p
  409                     if j > bound:
  410                         break
  411             return l
  412 
  413         # divPoly mainblock
  414         one = self.basefield.one
  415         Kx = ring.getRing(self.cubic)
  416         f = {-1: -Kx.one, 0: Kx.zero, 1: Kx.one, 2:Kx.one}
  417 
  418         # initialize f[3], f[4] and e
  419         if not Number:
  420             if self.ch <= 3:
  421                 raise ValueError("You must input (Number)")
  422             H = heart(card(self.basefield))
  423             loopbound = H[-1] + 2
  424             E = self.simple()
  425             e = 4 * E.cubic
  426             f[3] = _UniVarPolynomial({4:3*one,
  427                                      2:6*E.a,
  428                                      1:12*E.b,
  429                                      0:-E.a**2},
  430                                     self.basefield)
  431             f[4] = _UniVarPolynomial({6:2*one,
  432                                      4:10*E.a,
  433                                      3:40*E.b,
  434                                      2:-10*E.a**2,
  435                                      1:-8*E.a*E.b,
  436                                      0:-2*(E.a**3 + 8*E.b**2)},
  437                                     self.basefield)
  438         else:
  439             loopbound = Number + 1
  440             if self.ch == 0:
  441                 E = self.simple()
  442                 e = E.cubic
  443                 f[3] = _UniVarPolynomial({4:3*one,
  444                                          2:6*E.a,
  445                                          1:12*E.b,
  446                                          0:-E.a**2},
  447                                         self.basefield)
  448                 f[4] = _UniVarPolynomial({6:2*one,
  449                                          4:10*E.a,
  450                                          3:40*E.b,
  451                                          2:-10*E.a**2,
  452                                          1:-8*E.a*E.b,
  453                                          0:-2*(E.a**3 + 8*E.b**2)},
  454                                         self.basefield)
  455             else:
  456                 e = _UniVarPolynomial({3:4*one,
  457                                       2:self.b2,
  458                                       1:2*self.b4,
  459                                       0:self.b6},
  460                                      self.basefield)
  461                 f[3] = _UniVarPolynomial({4:3*one,
  462                                          3:self.b2,
  463                                          2:3*self.b4,
  464                                          1:3*self.b6,
  465                                          0:self.b8},
  466                                         self.basefield)
  467                 f[4] = _UniVarPolynomial({6:2*one,
  468                                          5:self.b2,
  469                                          4:5*self.b4,
  470                                          3:10*self.b6,
  471                                          2:10*self.b8,
  472                                          1:self.b2*self.b8 - self.b4*self.b6,
  473                                          0:self.b4*self.b8 - self.b6**2},
  474                                         self.basefield)
  475 
  476         # recursive calculation
  477         i = 5
  478         while i < loopbound:
  479             if i & 1:
  480                 j = (i - 1) >> 1
  481                 if j & 1:
  482                     f[i] = f[j+2]*f[j]**3 - e**2*f[j-1]*f[j+1]**3
  483                 else:
  484                     f[i] = e**2*f[j+2]*f[j]**3 - f[j-1]*f[j+1]**3
  485             else:
  486                 j = i >> 1
  487                 f[i] = (f[j+2]*f[j-1]**2 - f[j-2]*f[j+1]**2) * f[j]
  488             i += 1
  489 
  490         # final result
  491         if not Number:
  492             for i in range(2, loopbound, 2):
  493                 f[i] = 2*f[i]
  494             return (f, H)
  495         else:
  496             return f[Number]
  497 
  498 class ECoverQ(ECGeneric):
  499     """
  500     Elliptic curves over Q.
  501     """
  502     def __init__(self, coefficient):
  503         field = rational.RationalField()
  504         coeffs_list = []
  505         if isinstance(coefficient, list):
  506             for c in coefficient:
  507                 if isinstance(c, (int, long)):
  508                     coeff = field.createElement(c)
  509                 elif c in field:
  510                     coeff = c
  511                 else:
  512                     raise ValueError("coefficient not in basefield.")
  513                 coeffs_list.append(coeff)
  514 
  515         ECGeneric.__init__(self, coeffs_list, field)
  516 
  517     def __repr__(self):
  518         if len(self) == 2 or self.a1 == self.a2 == self.a3 == 0:
  519             return "ECoverQ(["+repr(self.a4)+","+repr(self.a6)+"])"
  520         else:
  521             return "ECoverQ(["+repr(self.a1)+","+repr(self.a2)+","+repr(self.a3)+","+repr(self.a4)+","+repr(self.a6)+"])"
  522 
  523     def point(self, limit=1000):
  524         """
  525         this returns a random point on eliiptic curve over Q.
  526         limit set maximal find time.
  527         """
  528         i = 9
  529         while i < limit:
  530             s = random.randrange(1, i)
  531             t = random.randrange(1, i)
  532             x = rational.Rational(s, t);
  533             y = self.coordinateY(x)
  534             if y != False:
  535                 return [x, y]
  536             i = i+10
  537         raise ValueError("Times exceeded for limit.")
  538 
  539 
  540 class ECoverGF(ECGeneric):
  541     """
  542     Elliptic curves over Galois Field.
  543     """
  544     def __init__(self, coefficient, basefield=None):
  545         """ create ECoverGF object.
  546         coefficient must be length 5 or 2 list(represent as Weierstrass form),
  547         any coefficient is in basefield.
  548         basefield must be FiniteField subclass object
  549         (i.e. FinitePrimeField or FiniteExtendedField object.)
  550         """
  551 
  552         # parameter parse
  553         try:
  554             character = basefield.getCharacteristic()
  555             field = basefield
  556         except AttributeError:
  557             # backward compatibility
  558             if isinstance(basefield, (int, long)):
  559                 field = finitefield.FinitePrimeField.getInstance(basefield)
  560                 character = basefield
  561             else:
  562                 raise ValueError("basefield must be FiniteField object.")
  563 
  564         coeffs_list = []
  565         if isinstance(coefficient, list):
  566             for c in coefficient:
  567                 if isinstance(c, (int, long)):
  568                     coeff = field.createElement(c)
  569                 elif c in field:
  570                     coeff = c
  571                 else:
  572                     raise ValueError("coefficient not in basefield.")
  573                 coeffs_list.append(coeff)
  574 
  575         # general initialize
  576         ECGeneric.__init__(self, coeffs_list, field)
  577 
  578         zero = self.basefield.zero
  579         one = self.basefield.one
  580 
  581         # format attribute
  582         if self.ch == 2:
  583             if len(self) == 5:
  584                 # FIXME
  585                 if coeffs_list[0] & 1 == one and coeffs_list[2] & 1 == coeffs_list[3] & 1 == zero and coeffs_list[4]:
  586                     self.a1 = one
  587                     self.a2 = coeffs_list[1]
  588                     self.a3 = zero
  589                     self.a4 = zero
  590                     self.a6 = coeffs_list[4]
  591                     self.b2 = one
  592                     self.b4 = zero
  593                     self.b6 = zero
  594                     self.b8 = self.a6
  595                     self.c4 = one
  596                     self.c6 = one
  597                     self.disc = self.a6
  598                     self.j = self.disc.inverse()
  599                 elif coeffs_list[0] & 1 == coeffs_list[1] & 1 == zero and coeffs_list[2]:
  600                     self.a1 = zero
  601                     self.a2 = zero
  602                     self.a3 = coeffs_list[2]
  603                     self.a4 = coeffs_list[3]
  604                     self.a6 = coeffs_list[4]
  605                     self.b2 = zero
  606                     self.b4 = zero
  607                     self.b6 = self.a3**2
  608                     self.b8 = self.a4**2
  609                     self.c4 = zero
  610                     self.c6 = zero
  611                     self.disc = self.a3**4
  612                     self.j = zero
  613                 else:
  614                     raise ValueError("coefficient may be not representation of EC.")
  615             else:
  616                 raise ValueError("coefficient may only use full Weierstrass form for characteristic 2.")
  617         elif self.ch == 3: # y^2=x^3+a2*x^2+a6 or y^2=x^3+a4*x+a6
  618             # FIXME
  619             if len(self) == 5:
  620                 if coeffs_list[0] % 3 == coeffs_list[2] % 3 == coeffs_list[3] % 3 == 0 and coeffs_list[1] and coeffs_list[4]:
  621                     self.a1 = zero
  622                     self.a2 = coeffs_list[1]
  623                     self.a3 = zero
  624                     self.a4 = zero
  625                     self.a6 = coeffs_list[4]
  626                     self.b2 = self.a2
  627                     self.b4 = zero
  628                     self.b6 = self.a6
  629                     self.b8 = self.a2*self.a6
  630                     self.c4 = self.b2**2
  631                     self.c6 = 2*self.b2**3
  632                     self.disc = -self.a2**3*self.a6
  633                     self.j = (-self.a2**3)*self.a6.inverse()
  634                 elif coeffs_list[0] == coeffs_list[1] == coeffs_list[2] == 0 and coeffs_list[3]:
  635                     self.a1 = zero
  636                     self.a2 = zero
  637                     self.a3 = zero
  638                     self.a4 = coeffs_list[3]
  639                     self.a6 = coeffs_list[4]
  640                     self.b2 = zero
  641                     self.b4 = 2*self.a4
  642                     self.b6 = self.a6
  643                     self.b8 = 2*self.a4**2
  644                     self.c4 = zero
  645                     self.c6 = zero
  646                     self.disc = -self.a4**3
  647                     self.j = zero
  648                 else:
  649                     raise ValueError("can't defined EC.")
  650                 if not self.disc:
  651                     raise ValueError("this curve is singular.")
  652             else:
  653                 raise ValueError("coefficient is less or more, can't defined EC.")
  654         else:
  655             if len(self) == 5:
  656                 self.a1 = coeffs_list[0]
  657                 self.a2 = coeffs_list[1]
  658                 self.a3 = coeffs_list[2]
  659                 self.a4 = coeffs_list[3]
  660                 self.a6 = coeffs_list[4]
  661                 self.b2 = self.a1**2+4*self.a2
  662                 self.b4 = self.a1*self.a3+2*self.a4
  663                 self.b6 = self.a3**2+4*self.a6
  664                 self.b8 = self.a1**2*self.a6+4*self.a2*self.a6-self.a1*self.a3*self.a4+self.a2*self.a3**2-self.a4**2
  665                 self.c4 = self.b2**2-24*self.b4
  666                 self.c6 = -self.b2**3+36*self.b2*self.b4-216*self.b6
  667                 self.disc = -self.b2**2*self.b8-8*self.b4**3-27*self.b6**2+9*self.b2*self.b4*self.b6
  668                 if self.disc:
  669                     self.j = self.c4**3*self.disc.inverse()
  670                 else:
  671                     raise ValueError("coefficients creates singular curve.")
  672             elif len(self) == 2:
  673                 self.a = coeffs_list[0]
  674                 self.b = coeffs_list[1]
  675                 self.a1 = zero
  676                 self.a2 = zero
  677                 self.a3 = zero
  678                 self.a4 = self.a
  679                 self.a6 = self.b
  680                 self.b2 = zero
  681                 self.b4 = 2*self.a
  682                 self.b6 = 4*self.b
  683                 self.b8 = -(self.a**2)
  684                 self.c4 = -48*self.a
  685                 self.c6 = -864*self.b
  686                 self.disc = -self.b2**2*self.b8-8*self.b4**3-27*self.b6**2+9*self.b2*self.b4*self.b6
  687                 if self.disc:
  688                     self.j = self.c4**3*self.disc.inverse()
  689                 else:
  690                     raise ValueError("coefficients creates singular curve.")
  691             else:
  692                 raise ValueError("coefficient is less or more, can't defined EC.")
  693 
  694         self.ord = None
  695         self.abelian = None
  696         self.cubic = _UniVarPolynomial({0:self.a6, 1:self.a4, 2:self.a2, 3:one},
  697                                       self.basefield)
  698 
  699     def point(self):
  700         """
  701         Return a random point on eliiptic curve over ch(field)>3
  702         """
  703         bfsize = card(self.basefield)
  704         one = self.basefield.one
  705         t = self.basefield.zero
  706         if len(self) == 2 or (self.a1 == self.a2 == self.a3 == self.basefield.zero):
  707             while self.basefield.Legendre(t) != 1:
  708                 s = self.basefield.createElement(bigrandom.randrange(bfsize))
  709                 t = self.cubic(s)
  710                 if not t:
  711                     return [s, t]
  712             t = self.basefield.sqrt(t)
  713             r = bigrandom.randrange(2)
  714             if r:
  715                 return [s, -t]
  716             return [s, t]
  717         elif self.ch != 2 and self.ch != 3:
  718             sform = self.simple()
  719             while sform.basefield.Legendre(t) != 1:
  720                 s = sform.basefield.createElement(bigrandom.randrange(bfsize))
  721                 t = (s**3+sform.a*s+sform.b)
  722             x = (s-3*self.b2) // (36*one)
  723             y = (sform.basefield.sqrt(t) // (108*one)-self.a1*x-self.a3)//(2*one)
  724             return [x, y]
  725         elif self.ch == 3:
  726             while sform.basefield.Legendre(t) != 1:
  727                 s = self.basefield.createElement(bigrandom.randrange(bfsize))
  728                 t = (s**3+self.a2*s**2+self.a4*s+self.a6)
  729             return [s, self.basefield.sqrt(t)]
  730         else:
  731             raise NotImplementedError("This is not implemented.")
  732 
  733     def Schoof(self):
  734         """
  735         Return t = q + 1 - #E(F_q). q is basefield order >>1.
  736         """
  737         if self.ch in (2, 3):
  738             raise NotImplementedError("characteristic should be >> 1")
  739 
  740         if len(self) != 2:
  741             return self.simple().Schoof()
  742 
  743         traces = []
  744         D, prime_moduli = self.divPoly()
  745         self.division_polynomials = D # temporary attribute
  746 
  747         # the main loop
  748         for l in prime_moduli:
  749             _log.debug("l = %d" % l)
  750             traces.append(self._Schoof_mod_l(l))
  751 
  752         del self.division_polynomials # temporary attribute deleted
  753 
  754         trace = arith1.CRT(traces)
  755         modulus = arith1.product(prime_moduli)
  756         if trace > (modulus >> 1):
  757             trace -= modulus
  758         assert abs(trace) <= 2*arith1.floorsqrt(card(self.basefield)), str(self)
  759         return trace
  760 
  761     def _Schoof_mod2(self):
  762         """
  763         Return (#E(F_q) mod 2, 2).  char(F_q) > 3 is required.
  764 
  765         For odd characteristic > 3, t = #E(F_q) mod 2 is determined by
  766         gcd(self.cubic, X^q - X) == 1 <=> t = 1 mod 2.
  767         """
  768         if not self.b:
  769             result = 0
  770             _log.debug("(%d, 2) #" % result)
  771         else:
  772             linearfactors = _UniVarPolynomial({card(self.basefield):self.basefield.one, 1:-self.basefield.one}, self.basefield)
  773             if _PolyGCD(self.cubic, linearfactors).degree() == 0:
  774                 result = 1
  775                 _log.debug("(%d, 2) ##" % result)
  776             else:
  777                 result = 0
  778                 _log.debug("(%d, 2) ###" % result)
  779         return (result, 2)
  780 
  781     def _Schoof_mod_l(self, l):
  782         """
  783         Return q + 1 - #E(F_q) mod l.
  784         """
  785         if l == 2:
  786             return self._Schoof_mod2()
  787         E = self.cubic
  788         D = self.division_polynomials
  789         lth_div = self.division_polynomials[l]
  790         field = self.basefield
  791         bfsize = card(field)
  792         x = _UniVarPolynomial({1:field.one}, field)
  793         k = bfsize % l
  794         x_frob = _PolyPow(x, bfsize, lth_div) #x_frob=x^q
  795         x_frobfrob = _PolyPow(x_frob, bfsize, lth_div) #x_frobfrob=x^{q^2}
  796 
  797         # test for x^{q^2} - x
  798         f, P = self._sub1(k, x_frobfrob - x, lth_div)
  799         f0, f3 = f[0], f[3]
  800 
  801         if _PolyGCD(lth_div, P).degree() > 0:
  802             if arith1.legendre(k, l) == -1:
  803                 _log.debug("%s $" % str((0, l)))
  804                 return (0, l)
  805 
  806             # arith1.legendre(k, l) == 1 <=> k is QR
  807             w = arith1.modsqrt(k, l)
  808             f, P = self._sub1(w, x_frob - x, lth_div)
  809 
  810             if _PolyGCD(lth_div, P).degree() == 0: # coprime
  811                 _log.debug("%s $$$$" % str((0, l)))
  812                 return (0, l)
  813 
  814             # there exist non trivial common divisors
  815             g0 = _PolyPow(E, (bfsize - 1) >> 1, lth_div) #y^(q-1)
  816             P = self._sub2(w, g0, f[3], lth_div)
  817 
  818             if _PolyGCD(lth_div, P).degree() > 0:
  819                 _log.debug("%s $$" % str((2*w % l, l)))
  820                 return (2*w % l, l)
  821             else:
  822                 _log.debug("%s $$$" % str((-2*w % l, l)))
  823                 return (-2*w % l, l)
  824 
  825         else: # coprime (GCD(P, lth_div).degree() == 0)
  826             Y = x - x_frobfrob
  827             g0 = _PolyPow(E, (bfsize - 1) >> 1, lth_div) #y^(q-1)
  828             g1 = _PolyPow(g0, bfsize + 1, lth_div) #y^(q^2-1)
  829             f = -self._sub2(k, g1, f3, lth_div)
  830             h1 = _PolyMulRed([f, f], lth_div)
  831             if k & 1 == 0:
  832                 g = (_PolyMulRed([Y, E, f3], lth_div) - f0) * 4
  833                 h0 = _PolyMulRed([g, g], lth_div)
  834                 aux1 = _PolyMulRed([f0, h0], lth_div) + h1
  835                 X_d = _PolyMulRed([E, f3, h0], lth_div)
  836             else:
  837                 g = (_PolyMulRed([Y, f3], lth_div) - _PolyMulRed([E, f0], lth_div)) * 4
  838                 h0 = _PolyMulRed([g, g], lth_div)
  839                 aux1 = _PolyMulRed([E, _PolyMulRed([f0, h0], lth_div) + h1], lth_div)
  840                 X_d = _PolyMulRed([f3, h0], lth_div)
  841             X_n = _PolyMulRed([X_d, x_frobfrob + x_frob + x], lth_div) - aux1
  842 
  843             # loop of t
  844             e_q = _PolyPow(self.cubic, bfsize, lth_div)
  845             for t in range(1, ((l - 1) >> 1) + 1):
  846                 Z_d_x, Z_n_x = self._Z_x(t, D, e_q, bfsize, lth_div)
  847                 # X_n * Z_d_x == X_d * Z_n_x (mod lth_div)?
  848                 if not _PolyMod(X_n * Z_d_x - X_d * Z_n_x, lth_div):
  849                     break
  850             else: # loop of t exhausted
  851                 _log.debug("%s @@@" % str((0, l)))
  852                 return (0, l)
  853 
  854             # found: X_n * Z_d_x == X_d * Z_n_x (mod lth_div)
  855             y0 = _PolyMulRed([-2*x_frobfrob - x, X_d], lth_div) + aux1
  856             if k & 1 == 0:
  857                 Y_d = _PolyMulRed([E, D[k], g, X_d], lth_div)
  858             else:
  859                 Y_d = _PolyMulRed([D[k], g, X_d], lth_div)
  860             Y_n = -_PolyMulRed([g1, Y_d], lth_div) - _PolyMulRed([f, y0], lth_div)
  861             Z_d_y, Z_n_y = self._Z_y(t, D, g0, bfsize, lth_div)
  862 
  863             # Y_n * Z_d_y == Y_d * Z_n_y (mod lth_div)?
  864             if _PolyMod(Y_n * Z_d_y - Y_d * Z_n_y, lth_div):
  865                 _log.debug("%s @@" % str((l-t, l)))
  866                 return (l-t, l)
  867             else:
  868                 _log.debug("%s @" % str((t, l)))
  869                 return (t, l)
  870 
  871     def _sub1(self, k, t, lth_div):
  872         """
  873         Compute in F_q[X]/(D_l):
  874         P = t * D_{k}^2 * f_E + D_{k-1} * D_{k+1}  (k is even)
  875         P = t * D_{k}^2 + D{k-1} * D{k+1} * f_E    (k is odd)
  876         and
  877         return a dict f and P, where f includes
  878         f[0] = D{k-1} * D{k+1} and f[3] = D_{k}^2
  879         for later computation.
  880 
  881         The aim is to compute x^q - x mod D_{l} or x^{q^2} - x mod D_{l}.
  882         """
  883         D = self.division_polynomials
  884         f = {}
  885         f[0] = _PolyMulRed([D[k-1], D[k+1]], lth_div)
  886         f[3] = _PolyMulRed([D[k], D[k]], lth_div)
  887         f[4] = _PolyMulRed([t, f[3]], lth_div)
  888         if k & 1 == 0:
  889             P = _PolyMulRed([f[4], self.cubic], lth_div) + f[0]
  890         else:
  891             P = f[4] + _PolyMulRed([f[0], self.cubic], lth_div)
  892         return f, P
  893 
  894     def _sub2(self, k, g, t, lth_div):
  895         """
  896         Compute in F_q[X]/(D[l]):
  897         P = 4 * g * t * D_{k} * f_E^2 - f_1 + f_2  (k is even)
  898         P = 4 * g * t * D_{k} - f_1 + f_2    (k is odd)
  899         where
  900         f_1 = D_{k-1}^2 * D_{k+2}
  901         f_2 = D_{k-2} * D_{k+1}^2
  902         and return P.
  903 
  904         The aim is to compute y^q - y mod D_{l} or y^{q^2} - y mod D_{l}.
  905         """
  906         D = self.division_polynomials
  907         fk = _PolyMulRed([4 * g, t, D[k]], lth_div)
  908         f1 = _PolyMulRed([D[k-1], D[k-1], D[k+2]], lth_div)
  909         f2 = _PolyMulRed([D[k-2], D[k+1], D[k+1]], lth_div)
  910         if k & 1 == 0:
  911             cubic = self.cubic
  912             P = _PolyMulRed([fk, cubic, cubic], lth_div) - f1 + f2
  913         else:
  914             P = fk - f1 + f2
  915         return P
  916 
  917     def _Z_x(self, t, D, e_q, bfsize, lth_div):
  918         d = _PolyPow(D[t], 2 * bfsize, lth_div)
  919         n = _PolyPow(_PolyMulRed([D[t-1], D[t+1]], lth_div), bfsize, lth_div)
  920         if t & 1 == 0:
  921             d = _PolyMulRed([e_q, d], lth_div)
  922         else:
  923             n = _PolyMulRed([e_q, n], lth_div)
  924         return d, n
  925 
  926     def _Z_y(self, t, D, g0, bfsize, lth_div):
  927         E = self.cubic
  928         if t & 1 == 0:
  929             d = _PolyPow(_PolyMulRed([4 * E, E, D[t], D[t], D[t]], lth_div), bfsize, lth_div)
  930         else:
  931             d = _PolyPow(_PolyMulRed([4 * D[t], D[t], D[t]], lth_div), bfsize, lth_div)
  932         z0 = _PolyPow(_PolyMulRed([D[t-1], D[t-1], D[t+2]], lth_div) - _PolyMulRed([D[t-2], D[t+1], D[t+1]], lth_div), bfsize, lth_div)
  933         n = _PolyMulRed([g0, z0], lth_div)
  934         return d, n
  935 
  936     def _step(self, P, W):
  937         """
  938         Return three component list [A, B, C] used in Shanks_Mestre,
  939         where A is a list of x-coordinates of [q+1]P to [q+1+W]P,
  940         B is a list of x-coordinates of [0]P, [W]P to [W*W]P and
  941         C is the intersection set of A and B.
  942         """
  943         Qa = self.mul(self.ch + 1, P)
  944         if len(Qa) == 1:
  945             _log.debug("[%d]P is zero" % (self.ch + 1))
  946         Qb = R = self.mul(W, P)
  947         A = [Qa[0]]
  948         B = [0, Qb[0]] # 0 = [0][0] ((infinity)_x)
  949         for i in range(1, W):
  950             Qa = self.add(Qa, P)
  951             Qb = self.add(Qb, R)
  952             A.append(Qa[0])
  953             B.append(Qb[0])
  954             if len(Qa) == 1:
  955                 _log.debug("[%d]P is zero" % (self.ch + 1 + i))
  956                 break
  957         return [A, B, set(A).intersection(set(B))]
  958 
  959     def Shanks_Mestre(self):
  960         """
  961         Return t = p + 1 - #E(F_p) for prime field F_p.
  962 
  963         Use in the range 229 <= self.ch <=10**5+o(1).
  964 
  965         This method is using
  966         Algorithm 7.5.3(Shanks-Mestre assessment of curve order)
  967         Crandall & Pomerance, PRIME NUMBERS
  968         """
  969         if type(self.basefield) != finitefield.FinitePrimeField:
  970             raise NotImplementedError("FIXME: this is only implemented for FinitePrimeField.")
  971 
  972         if self.ch <= 229:
  973             return self.naive()
  974         else: #E.ch>229
  975             if len(self) != 2:
  976                 return self.simple().Shanks_Mestre()
  977             g = 0
  978             while arith1.legendre(g, self.ch) != -1:
  979                 g = bigrandom.randrange(2, self.ch)
  980             W = arith1.floorsqrt(2 * (arith1.floorsqrt(self.ch) + 1)) + 1
  981 
  982             c, d = g**2*self.a, g**3*self.b
  983             f = self.cubic
  984             used = set()
  985             while True:
  986                 x = bigrandom.randrange(self.ch)
  987                 while x in used:
  988                     x = bigrandom.randrange(self.ch)
  989                 used.add(x)
  990                 y2 = f(x)
  991                 if not y2:
  992                     continue
  993                 if arith1.legendre(y2.n, self.ch) == 1:
  994                     E = self
  995                     cg = 1
  996                 else: #arith1.legendre(f(x),self.ch)==-1
  997                     E = EC([c, d], self.ch)
  998                     cg = -1
  999                     x = g*x % self.ch
 1000                 x = self.basefield.createElement(x)
 1001                 P = [x, E.coordinateY(x)]
 1002                 A, B, S = E._step(P, W)
 1003                 _log.debug("#S = %d" % len(S))
 1004                 if len(S) == 1:
 1005                     s = S.pop()
 1006                     if B.count(s) <= 2:
 1007                         break
 1008             aa = A.index(s)
 1009             bb = B.index(s)
 1010             t = aa - bb*W
 1011             if E.mul(self.ch + 1 + t, P) == [0]:
 1012                 return -cg*t
 1013             else:
 1014                 t = aa + bb*W
 1015                 return -cg*t
 1016 
 1017     def naive(self):
 1018         """
 1019         Return t = p + 1 - #E(Fp).
 1020 
 1021         This method computes t by counting up Legendre symbols and thus
 1022         needs exponential time.
 1023 
 1024         RESTRICTIONS:
 1025         The field cannot be of characteristic two nor three.
 1026         """
 1027         if self.ch in (2, 3):
 1028             raise TypeError("cannot be defined over characteristic two/three.")
 1029 
 1030         sform = self.simple()
 1031 
 1032         k = 0
 1033         f = sform.cubic
 1034         for i in range(card(sform.basefield)):
 1035             x = sform.basefield.createElement(i)
 1036             k += sform.basefield.Legendre(f(x))
 1037         return -k
 1038 
 1039     def _order_2(self):
 1040         """
 1041         Return #E(F_2).
 1042 
 1043         E is in long Weierstrass form:
 1044         Y^2 + a1 XY + a3 Y = X^3 + a2 X^2 + a4 X + a6
 1045         """
 1046         assert self.ch == 2
 1047         result = 1 # for infinite place
 1048         if not self.a6: # (0, 0)
 1049             result += 1
 1050         if self.a3 != self.a6: # (0, 1)
 1051             result += 1
 1052         if self.a2 + self.a4 + self.a6: # (1, 0)
 1053             result += 1
 1054         if self.a1 + self.a3 == self.a2 + self.a4 + self.a6: # (1, 1)
 1055             result += 1
 1056         return result
 1057 
 1058     def _order_3(self):
 1059         """
 1060         Return #E(F_3).
 1061 
 1062         E is in long Weierstrass form:
 1063         Y^2 + a1 XY + a3 Y = X^3 + a2 X^2 + a4 X + a6
 1064         """
 1065         assert self.ch == 3
 1066         one = self.basefield.one
 1067         result = 1 # for infinite place
 1068         if not self.a6: # (0, 0)
 1069             result += 1
 1070         if one + self.a3 == self.a6: # (0, 1)
 1071             result += 1
 1072         if one - self.a3 == self.a6: # (0, 2)
 1073             result += 1
 1074         if not (one + self.a2 + self.a4 + self.a6): # (1, 0)
 1075             result += 1
 1076         if self.a1 + self.a3 == self.a2 + self.a4 + self.a6: # (1, 1)
 1077             result += 1
 1078         if -(self.a1 + self.a3) == self.a2 + self.a4 + self.a6: # (1, 2)
 1079             result += 1
 1080         if one == self.a2 - self.a4 + self.a6: # (2, 0)
 1081             result += 1
 1082         if -self.a1 + self.a3 == one + self.a2 - self.a4 + self.a6: # (2, 1)
 1083             result += 1
 1084         if self.a1 - self.a3 == one + self.a2 - self.a4 + self.a6: # (2, 2)
 1085             result += 1
 1086         return result
 1087 
 1088     def _order_to_trace(self, order):
 1089         """
 1090         Return the trace calculated from the given order:
 1091           t = q + 1 - #E(F_q)
 1092         """
 1093         return card(self.basefield) + 1 - order
 1094 
 1095     def _trace_to_order(self, trace):
 1096         """
 1097         Return the order calculated from the given trace:
 1098           #(E_q) = q + 1 - t
 1099         """
 1100         return card(self.basefield) + 1 - trace
 1101 
 1102     def trace(self, index=None):
 1103         """
 1104         Return the Frobenius trace t = q + 1 - #E(F_q),
 1105         where q is basefield cardinality.
 1106 
 1107         If index is an integer greater than 1, then return the trace
 1108         t = q^r + 1 - #E(F_q^r) for a subfield curve defined over F_q.
 1109         """
 1110         bfsize = card(self.basefield)
 1111 
 1112         if not self.ord:
 1113             if bfsize == self.ch: # prime field
 1114                 # special cases
 1115                 if bfsize == 2 or bfsize == 3:
 1116                     trace = self._order_to_trace(self.order())
 1117                 # trace main block
 1118                 elif bfsize < 10**4:
 1119                     trace = self.naive()
 1120                 elif bfsize < 10**30:
 1121                     trace = self.Shanks_Mestre()
 1122                 else: # self.ch>=10**30
 1123                     trace = self.Schoof()
 1124             else:
 1125                 if self.ch in (2, 3):
 1126                     error_message = "no E/F_{%d} trace" % bfsize
 1127                     raise NotImplementedError(error_message)
 1128                 else:
 1129                     trace = self.Schoof()
 1130 
 1131             self.ord = self._trace_to_order(trace) # cached
 1132         else:
 1133             trace = self._order_to_trace(self.ord)
 1134 
 1135         # final result
 1136         if index is not None:
 1137             # for subfield curve
 1138             basetrace = trace
 1139             trace, oldtrace = basetrace, 2
 1140             for i in range(2, index + 1):
 1141                 trace, oldtrace = basetrace*trace - bfsize*oldtrace, trace
 1142 
 1143         return trace
 1144 
 1145     def order(self, index=None):
 1146         """
 1147         Return #E(F_q) or #E(F_{q^r}).
 1148 
 1149         E is defined over F_q.
 1150         If the method is called as E.order(), the result is #E(F_q).
 1151         If the method is called as E.order(r), the result is #E(F_{q^r}).
 1152 
 1153         RESTRICTIONS:
 1154         F_q cannot be q = 2^k or q = 3^k for k > 1.
 1155         """
 1156         bfsize = card(self.basefield)
 1157 
 1158         if not self.ord:
 1159             if self.ch in (2, 3):
 1160                 if bfsize == self.ch == 2:
 1161                     self.ord = self._order_2()
 1162                 elif bfsize == self.ch == 3:
 1163                     self.ord = self._order_3()
 1164                 else:
 1165                     error_message = "no E/F_{%d} order" % bfsize
 1166                     raise NotImplementedError(error_message)
 1167             else:
 1168                 self.ord = self._trace_to_order(self.trace())
 1169 
 1170         # final result
 1171         if index:
 1172             # for subfield curve
 1173             basetrace = self._order_to_trace(self.ord)
 1174             trace, oldtrace = basetrace, 2
 1175             for i in range(2, index + 1):
 1176                 trace, oldtrace = basetrace*trace - bfsize*oldtrace, trace
 1177             return bfsize ** index + 1 - trace
 1178 
 1179         return self.ord
 1180 
 1181     def line(self, P, Q, R):
 1182         """
 1183         this use to compute weil pairing.
 1184         return line_{P,Q}(R).
 1185         """
 1186         if not R or R == self.infpoint:
 1187             raise ValueError("R must not be zero")
 1188         if P == Q == self.infpoint:
 1189             return self.basefield.one
 1190         if P == self.infpoint:
 1191             return R[0] - Q[0]
 1192         if Q == self.infpoint:
 1193             return R[0] - P[0]
 1194         if P[0] != Q[0]:
 1195             return (Q[0] - P[0]) * R[1] - (Q[1] - P[1]) * \
 1196                      R[0]- Q[0] * P[1] + P[0] * Q[1]
 1197         if P[1] != Q[1]:
 1198             return R[0] - P[0]
 1199         return (3 * P[0] ** 2 + 2 * self.a2 * P[0] + self.a4 - self.a1 * P[1] ) * \
 1200                 R[0] - (2 * P[1] + self.a1 * P[0] + self.a3 ) * R[1] - \
 1201                 (P[0] ** 3) + self.a4 * P[0] + 2 * self.a6 - self.a3 * P[1]
 1202 
 1203     def pointorder(self, P, ord_factor=None):
 1204         """
 1205         find point order of P and return order.
 1206         """
 1207         # parameter ord_factor is extension for structre.
 1208         if ord_factor is not None:
 1209             N = int(ord_factor)
 1210         else:
 1211             N = self.order()
 1212             ord_factor = factor_misc.FactoredInteger(N)
 1213         o = 1
 1214         for p, e in ord_factor:
 1215             B = self.mul(N//(p**e), P)
 1216             while B != self.infpoint:
 1217                 o = o*p
 1218                 B = self.mul(p, B)
 1219         return o
 1220 
 1221     def Miller(self, P, m, Q, R):
 1222         """
 1223         return f_P(D_Q)
 1224         this is used for Tate/Weil pairing
 1225 
 1226         suppose that P be an m-torsion point .
 1227         support point R must be in neither groups generated by P nor Q.
 1228 
 1229         return False if parameters lack any conditions above.
 1230 
 1231         NOTE: algorithm limitations forced that assume R is not P-Q .
 1232         """
 1233         # check order
 1234         if m < 2 or not (m % self.ch):
 1235             raise ValueError("order more than 1 and not be divisible by characteristic")
 1236 
 1237         O = self.infpoint
 1238 
 1239         # support point must not be P-Q
 1240         S = self.add(R, Q)
 1241         if S == O:
 1242             return False
 1243 
 1244         # j = 1
 1245         jP = P
 1246         v = self.basefield.one
 1247         for k in arith1.expand(m, 2)[-2::-1]:
 1248             j2P = self.mul(2, jP)
 1249             denominator = self.line(jP, jP, R) * self.line(j2P, O, S)
 1250             if not denominator:
 1251                 return False
 1252             numerator = self.line(jP, jP, S) * self.line(j2P, O, R)
 1253             if not numerator:
 1254                 return False
 1255             f = numerator / denominator
 1256             v = v**2 * f
 1257             # j *= 2
 1258             jP = j2P
 1259             if k:
 1260                 kjP = self.add(P, jP)
 1261                 denominator = self.line(P, jP, R) * self.line(kjP, O, S)
 1262                 if not denominator:
 1263                     return False
 1264                 numerator = self.line(P, jP, S) * self.line(kjP, O, R)
 1265                 if not numerator:
 1266                     return False
 1267                 f = numerator / denominator
 1268                 v = v * f
 1269                 # j += 1
 1270                 jP = kjP
 1271         # now j == m
 1272         return v
 1273 
 1274     def TatePairing(self, m, P, Q):
 1275         """
 1276         computing the Tate-Lichetenbaum pairing with Miller's algorithm.
 1277         parameters satisfies that mul(m,P)==[0].
 1278         """
 1279         O = self.infpoint
 1280         if self.mul(m, P) != O:
 1281             raise ValueError("sorry, not mP=[0].")
 1282 
 1283         if P == O or Q == O:
 1284             return self.basefield.one
 1285 
 1286         forbidden = [O, P, self.mul(-1, Q), self.sub(P, Q)]
 1287         R = self.add(P, Q)
 1288         T = False
 1289         while (not T):
 1290             while R in forbidden:
 1291                 R = self.point()
 1292             forbidden.append(R)
 1293             T = self.Miller(P, m, Q, R)
 1294         return T
 1295 
 1296     def TatePairing_Extend(self, m, P, Q):
 1297         """
 1298         computing the Tate-Lichtenbaum pairing extended original Tate Pairing.
 1299         """
 1300         return self.TatePairing(m, P, Q)**((card(self.basefield)-1)//m)
 1301 
 1302     def WeilPairing(self, m, P, Q):
 1303         """
 1304         computing the Weil pairing with Miller's algorithm.
 1305         we assume point P and Q that be in different m-tortion group .
 1306         """
 1307         O = self.infpoint
 1308         if self.mul(m, P) != O or self.mul(m, Q) != O:
 1309             raise ValueError("sorry, not mP=[0] or mQ=[0].")
 1310 
 1311         if P == O or Q == O or P == Q:
 1312             return self.basefield.one
 1313 
 1314         T = U = False
 1315         forbidden = [O, P, Q, self.sub(Q, P)]
 1316         R = self.sub(P,Q) # assume Q not in group generated P
 1317         while (not T) or (not U):
 1318             while R in forbidden:
 1319                 R = self.point()
 1320             T = self.Miller(Q, m, P, R)
 1321             U = self.Miller(P, m, Q, self.mul(-1, R))
 1322             forbidden.append(R)
 1323         F = U/T
 1324         return F
 1325 
 1326     def BSGS(self, P):
 1327         """
 1328         returns order of P such that kP=[0]
 1329         refered to Washington 4.3.4.
 1330         """
 1331         if P == self.infpoint:
 1332             return 1
 1333 
 1334         bfsize = card(self.basefield)
 1335 
 1336         Q = self.mul(bfsize + 1, P)
 1337         m = arith1.floorpowerroot(bfsize, 4) + 1
 1338         Plist = [self.infpoint]
 1339         R = P
 1340         j = 1
 1341         while j <= m:
 1342             Plist.append(R)
 1343             R = self.add(R, P)
 1344             j += 1
 1345         R = self.mul(2*m, P)
 1346         k = -m
 1347         Plist_rev = map(self.mul, [-1]*(m+1), Plist) # make reverse point mapping
 1348         j = 0
 1349         while k <= m:
 1350             S = self.add(Q, self.mul(k, R))
 1351             if S in Plist:
 1352                 j = Plist.index(S)
 1353                 break
 1354             elif S in Plist_rev:
 1355                 j = -Plist_rev.index(S)
 1356                 break
 1357             k += 1
 1358         M = self.ch + 1 + 2*m*k - j
 1359         Flist = factor_methods.factor(M)
 1360         for p, e in Flist:
 1361             for i in range(e):
 1362                 if self.mul(M//p, P) == self.infpoint:
 1363                     M = M//p
 1364         return M
 1365 
 1366     def DLP_BSGS(self, n, P, Q):
 1367         """
 1368         returns k such that Q=kP
 1369         P,Q is the elements of the same group
 1370         """
 1371         B = []
 1372         m = arith1.floorsqrt(n) + 1
 1373         R = Q
 1374         B.append(Q)
 1375         i = 1
 1376         while i <= m:
 1377             R = self.sub(R, P)
 1378             if R == [0]:
 1379                 return i
 1380             B.append(R)
 1381             i += 1
 1382         R = self.mul(m, P)
 1383         P = R
 1384         if R in B:
 1385             return m+B.index(R)
 1386         else:
 1387             i = 2
 1388             while i <= m:
 1389                 R = self.add(R, P)
 1390                 if R in B:
 1391                     return i*m+B.index(R)
 1392                 i = i+1
 1393             if self.mul(n, P) == Q:
 1394                 return n
 1395             else:
 1396                 return False
 1397 
 1398     def structure(self):
 1399         """
 1400         returns group structure E(K)=Z_d x Z_m with d|m|#E(K)
 1401         """
 1402         if self.abelian:
 1403             return self.abelian
 1404 
 1405         # step 1. find order E/F_p.
 1406         simplified = self.simple()
 1407         N = simplified.order()
 1408         if prime.primeq(N):
 1409             return (1, N)
 1410 
 1411         # step 2. decompose N.
 1412         r = gcd.gcd(card(simplified.basefield) - 1, N)
 1413         _log.debug("r = %d, N = %d" % (r, N))
 1414         r_factor = factor_misc.FactoredInteger(r)
 1415         N1, N2 = 1, N
 1416         for p in r_factor.prime_divisors():
 1417             k, N2 = arith1.vp(N2, p=p)
 1418             N1 *= p**k
 1419 
 1420         _log.debug("loop")
 1421         while 1:
 1422             P1 = self.infpoint
 1423             while P1 == self.infpoint:
 1424                 P1 = simplified.point()
 1425             P2 = self.infpoint
 1426             while P2 == self.infpoint:
 1427                 P2 = simplified.point()
 1428             P1, P2 = simplified.mul(N2, P1), simplified.mul(N2, P2)
 1429             s = simplified.pointorder(P1, r_factor)
 1430             t = simplified.pointorder(P2, r_factor)
 1431             m = gcd.lcm(s, t)
 1432             if m > 1:
 1433                 e = simplified.WeilPairing(m, P1, P2)
 1434                 if e != self.basefield.one:
 1435                     d = e.order()
 1436                 else:
 1437                     d = 1
 1438                 if m*d == N1:
 1439                     _log.debug("N1 = %d" % N1)
 1440                     _log.debug("P1 = %s (pointorder=%d)" % (P1, s))
 1441                     _log.debug("P2 = %s (pointorder=%d)" % (P2, t))
 1442                     assert (not (N//d) % d), d
 1443                     self.abelian = (d, N//d)
 1444                     return self.abelian
 1445 
 1446     def issupersingular(self):
 1447         """
 1448         returns supersingularities.
 1449         """
 1450         if self.order() % self.ch == 1:
 1451             return True
 1452         else:
 1453             return False
 1454 
 1455 def EC(coefficient, basefield=None):
 1456     """
 1457     generate new elliptic curve class.
 1458     """
 1459     try:
 1460         character = basefield.getCharacteristic()
 1461         field = basefield
 1462     except:
 1463         # backward compatiblity
 1464         if isinstance(basefield, (int, long)):
 1465             field = finitefield.FinitePrimeField(basefield)
 1466             character = basefield
 1467         elif isinstance(basefield, rational.RationalField) or not basefield:
 1468             character = 0 # necessary to be defined
 1469         else:
 1470             raise ValueError("basefield must be RationalFieid or FiniteField.")
 1471 
 1472     if isinstance(coefficient, list):
 1473         if not character:
 1474             return ECoverQ(coefficient)
 1475         else:
 1476             return ECoverGF(coefficient, field)