## "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
304         """
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:
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):
383                 if l[j] == 1:
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):
391                 if l[j] == 1:
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):
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)
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
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:
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)]
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)
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:
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)
```