NZMATH  1.2.0
About: NZMATH is a Python based number theory oriented calculation system.
  Fossies Dox: NZMATH-1.2.0.tar.gz  ("inofficial" and yet experimental doxygen-generated source code documentation)  

elliptic.py
Go to the documentation of this file.
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
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 
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 
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)
nzmath.elliptic.ECoverGF.BSGS
def BSGS(self, P)
Definition: elliptic.py:1326
nzmath.elliptic.ECGeneric.b
b
Definition: elliptic.py:133
nzmath.bigrandom
Definition: bigrandom.py:1
nzmath.elliptic.ECoverGF._sub1
def _sub1(self, k, t, lth_div)
Definition: elliptic.py:871
nzmath.bigrange.range
def range(start, stop=None, step=None)
Definition: bigrange.py:19
nzmath.ring
Definition: ring.py:1
nzmath.elliptic.ECoverGF.issupersingular
def issupersingular(self)
Definition: elliptic.py:1446
nzmath.elliptic.ECoverGF
Definition: elliptic.py:540
nzmath.finitefield.FinitePrimeField
Definition: finitefield.py:205
nzmath.elliptic.ECoverGF.Miller
def Miller(self, P, m, Q, R)
Definition: elliptic.py:1221
nzmath.elliptic.ECGeneric.changeCurve
def changeCurve(self, V)
Definition: elliptic.py:200
nzmath.elliptic._PolyMod
def _PolyMod(f, g)
Definition: elliptic.py:48
nzmath.poly.multiutil
Definition: multiutil.py:1
nzmath.elliptic.ECGeneric.ch
ch
Definition: elliptic.py:113
nzmath.elliptic.ECoverGF.WeilPairing
def WeilPairing(self, m, P, Q)
Definition: elliptic.py:1302
nzmath.factor.misc
Definition: misc.py:1
nzmath.elliptic._PolyMulRed
def _PolyMulRed(multipliees, poly)
Definition: elliptic.py:64
nzmath.elliptic._UniVarPolynomial
def _UniVarPolynomial(dict, coeffring=None)
Definition: elliptic.py:29
nzmath.elliptic.ECoverGF._trace_to_order
def _trace_to_order(self, trace)
Definition: elliptic.py:1095
nzmath.elliptic.ECGeneric.a1
a1
Definition: elliptic.py:119
nzmath.elliptic.ECGeneric.c6
c6
Definition: elliptic.py:129
nzmath.elliptic.ECoverGF.abelian
abelian
Definition: elliptic.py:695
nzmath.elliptic.ECGeneric.sub
def sub(self, P, Q)
Definition: elliptic.py:352
nzmath.elliptic.ECGeneric.a6
a6
Definition: elliptic.py:123
nzmath.poly.termorder
Definition: termorder.py:1
nzmath.elliptic.EC
def EC(coefficient, basefield=None)
Definition: elliptic.py:1455
nzmath.elliptic.ECoverGF.TatePairing_Extend
def TatePairing_Extend(self, m, P, Q)
Definition: elliptic.py:1296
nzmath.elliptic.MultiVarPolynomial
MultiVarPolynomial
Definition: elliptic.py:26
nzmath.gcd
Definition: gcd.py:1
nzmath.elliptic.ECoverGF._order_3
def _order_3(self)
Definition: elliptic.py:1058
nzmath.elliptic.ECoverGF.division_polynomials
division_polynomials
Definition: elliptic.py:745
nzmath.elliptic.ECoverGF._sub2
def _sub2(self, k, g, t, lth_div)
Definition: elliptic.py:894
nzmath.elliptic.ECGeneric.j
j
Definition: elliptic.py:150
nzmath.rational
Definition: rational.py:1
nzmath.rational.Rational
Definition: rational.py:10
nzmath.elliptic.ECoverQ.__repr__
def __repr__(self)
Definition: elliptic.py:517
nzmath.elliptic.ECGeneric.a4
a4
Definition: elliptic.py:122
nzmath.elliptic.ECoverGF.point
def point(self)
Definition: elliptic.py:699
nzmath.elliptic.ECGeneric.disc
disc
Definition: elliptic.py:130
nzmath.elliptic.ECoverGF.Shanks_Mestre
def Shanks_Mestre(self)
Definition: elliptic.py:959
nzmath.elliptic.ECGeneric.coefficient
coefficient
Definition: elliptic.py:116
nzmath.elliptic.ECoverGF.DLP_BSGS
def DLP_BSGS(self, n, P, Q)
Definition: elliptic.py:1366
nzmath.elliptic.ECGeneric.simple
def simple(self)
Definition: elliptic.py:185
nzmath.elliptic.ECoverGF._step
def _step(self, P, W)
Definition: elliptic.py:936
nzmath.elliptic.ECGeneric.a
a
Definition: elliptic.py:132
nzmath.elliptic.ECoverGF.ord
ord
Definition: elliptic.py:694
nzmath.elliptic.ECGeneric.infpoint
infpoint
Definition: elliptic.py:114
nzmath.elliptic.ECGeneric.__repr__
def __repr__(self)
Definition: elliptic.py:161
nzmath.elliptic.ECGeneric.whetherOn
def whetherOn(self, P)
Definition: elliptic.py:280
nzmath.elliptic.ECoverGF._Schoof_mod2
def _Schoof_mod2(self)
Definition: elliptic.py:761
nzmath.elliptic.ECGeneric.cubic
cubic
Definition: elliptic.py:151
nzmath.elliptic._strMultiPoly
def _strMultiPoly(poly, symbol=["x","y"], asc=False)
Definition: elliptic.py:37
nzmath.elliptic.ECoverGF.naive
def naive(self)
Definition: elliptic.py:1017
nzmath.elliptic.ECoverGF._Z_y
def _Z_y(self, t, D, g0, bfsize, lth_div)
Definition: elliptic.py:926
nzmath.elliptic.ECGeneric.divPoly
def divPoly(self, Number=None)
Definition: elliptic.py:395
nzmath.elliptic.ECoverGF.Schoof
def Schoof(self)
Definition: elliptic.py:733
nzmath.elliptic.ECoverQ
Definition: elliptic.py:498
nzmath.elliptic.ECoverQ.__init__
def __init__(self, coefficient)
Definition: elliptic.py:502
nzmath.elliptic.ECoverGF._Schoof_mod_l
def _Schoof_mod_l(self, l)
Definition: elliptic.py:781
nzmath.compatibility.cmp
cmp
Definition: compatibility.py:20
nzmath.arith1
Definition: arith1.py:1
nzmath.poly.uniutil
Definition: uniutil.py:1
nzmath.arith1.inverse
def inverse(x, n)
Definition: arith1.py:183
nzmath.elliptic.ECGeneric.changePoint
def changePoint(self, P, V)
Definition: elliptic.py:227
nzmath.elliptic.ECGeneric.b4
b4
Definition: elliptic.py:125
nzmath.elliptic._PolyGCD
def _PolyGCD(f, g)
Definition: elliptic.py:54
nzmath.prime
Definition: prime.py:1
nzmath.elliptic.ECGeneric.b8
b8
Definition: elliptic.py:127
nzmath.compatibility
Definition: compatibility.py:1
nzmath.elliptic.ECoverGF._order_to_trace
def _order_to_trace(self, order)
Definition: elliptic.py:1088
nzmath.elliptic.ECoverGF.TatePairing
def TatePairing(self, m, P, Q)
Definition: elliptic.py:1274
nzmath.elliptic.ECGeneric.a3
a3
Definition: elliptic.py:121
nzmath.elliptic.ECoverGF.order
def order(self, index=None)
Definition: elliptic.py:1145
nzmath.elliptic.ECoverGF._Z_x
def _Z_x(self, t, D, e_q, bfsize, lth_div)
Definition: elliptic.py:917
nzmath.factor.methods
Definition: methods.py:1
nzmath.rational.RationalField
Definition: rational.py:517
nzmath.elliptic.ECoverGF.structure
def structure(self)
Definition: elliptic.py:1398
nzmath.elliptic.ECGeneric.c4
c4
Definition: elliptic.py:128
nzmath.elliptic.ECoverGF.line
def line(self, P, Q, R)
Definition: elliptic.py:1181
nzmath.elliptic.ECGeneric.mul
def mul(self, k, P)
Definition: elliptic.py:374
nzmath.elliptic.ECGeneric
Definition: elliptic.py:85
nzmath.elliptic.ECGeneric.__len__
def __len__(self)
Definition: elliptic.py:158
nzmath.elliptic.ECoverGF.trace
def trace(self, index=None)
Definition: elliptic.py:1102
nzmath.elliptic.ECGeneric.__init__
def __init__(self, coefficient, basefield=None)
Definition: elliptic.py:90
nzmath.elliptic.ECGeneric.__str__
def __str__(self)
Definition: elliptic.py:168
nzmath.elliptic._strUniPoly
def _strUniPoly(poly, symbol="x", asc=True)
Definition: elliptic.py:33
nzmath.elliptic._PolyPow
def _PolyPow(f, d, g)
Definition: elliptic.py:58
nzmath.elliptic.ECoverGF.__init__
def __init__(self, coefficient, basefield=None)
Definition: elliptic.py:544
nzmath.elliptic.ECGeneric.add
def add(self, P, Q)
Definition: elliptic.py:303
nzmath.elliptic.ECGeneric.b2
b2
Definition: elliptic.py:124
nzmath.elliptic.ECGeneric.basefield
basefield
Definition: elliptic.py:99
nzmath.elliptic.ECoverQ.point
def point(self, limit=1000)
Definition: elliptic.py:523
nzmath.elliptic.ECGeneric.b6
b6
Definition: elliptic.py:126
nzmath.elliptic.ECoverGF._order_2
def _order_2(self)
Definition: elliptic.py:1039
nzmath.elliptic.ECoverGF.pointorder
def pointorder(self, P, ord_factor=None)
Definition: elliptic.py:1203
nzmath.compatibility.card
def card(virtualset)
Definition: compatibility.py:28
nzmath.elliptic.ECGeneric.a2
a2
Definition: elliptic.py:120
nzmath.finitefield
Definition: finitefield.py:1
nzmath.elliptic.ECGeneric.coordinateY
def coordinateY(self, x)
Definition: elliptic.py:245