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)  

ecpp.py
Go to the documentation of this file.
1 from __future__ import division
2 import os
3 import csv
4 import logging
5 import nzmath.arith1 as arith1
6 import nzmath.bigrandom as bigrandom
7 import nzmath.bigrange as bigrange
8 import nzmath.equation as equation
9 import nzmath.factor as factor
10 import nzmath.intresidue as intresidue
11 import nzmath.prime as prime
12 import nzmath.quad as quad
13 import nzmath.squarefree as squarefree
15 from nzmath.config import DATADIR, HAVE_MPMATH, HAVE_NET
16 
17 if HAVE_MPMATH:
18  import mpmath
19 
20  log = mpmath.log
21 
23  """
24  Round x to nearest integer.
25  """
26  return mpmath.floor(x + 0.5)
27 
28 else:
29  import math
30 
31  log = math.log
32 
33  def nearest_integer(x):
34  """
35  Round x to nearest integer.
36  """
37  return math.floor(x + 0.5)
38 
39 
40 _log = logging.getLogger('nzmath.ecpp')
41 
42 # mapping from disc to list of order functions.
43 default_orders = [lambda n, u, v: n + 1 + u,
44  lambda n, u, v: n + 1 - u]
45 orders = {-3: default_orders + [lambda n, u, v: n + 1 - ((u + 3*v) >> 1),
46  lambda n, u, v: n + 1 + ((u + 3*v) >> 1),
47  lambda n, u, v: n + 1 - ((u - 3*v) >> 1),
48  lambda n, u, v: n + 1 + ((u - 3*v) >> 1)],
49  -4: default_orders + [lambda n, u, v: n + 1 + (v << 1),
50  lambda n, u, v: n + 1 - (v << 1)]}
51 # default absolute value bound of discriminant
52 DEFAULT_ABS_BOUND = 200000
53 
54 def dedekind(tau, floatpre):
55  """
56  Algorithm 22 (Dedekind eta)
57  Input : tau in the upper half-plane, k in N
58  Output : eta(tau)
59  """
60  a = 2 * mpmath.pi / mpmath.mpf(24)
61  b = mpmath.exp(mpmath.mpc(0, a))
62  p = 1
63  m = 0
64  while m <= 0.999:
65  n = nearest_integer(tau.real)
66  if n != 0:
67  tau -= n
68  p *= b**n
69  m = tau.real*tau.real + tau.imag*tau.imag
70  if m <= 0.999:
71  ro = mpmath.sqrt(mpmath.power(tau, -1)*1j)
72  if ro.real < 0:
73  ro = -ro
74  p = p*ro
75  tau = (-p.real + p.imag*1j) / m
76  q1 = mpmath.exp(a*tau*1j)
77  q = q1**24
78  s = 1
79  qs = mpmath.mpc(1, 0)
80  qn = 1
81  des = mpmath.mpf(10)**(-floatpre)
82  while abs(qs) > des:
83  t = -q*qn*qn*qs
84  qn = qn*q
85  qs = qn*t
86  s += t + qs
87  return p*q1*s
88 
90  """
91  Algorithm 23&24 (Floating-point precision & Hilbert class polynomial)
92  Input : a negative fundamental discriminant D
93  Output : the class number h(D) and the Hilbert class polynomial T in Z[X]
94  """
95  if D >= 0:
96  raise ValueError("Input number must be a negative integer.")
97  T = [complex(1), 0]
98  mpmath.mp.default()
99  h, reduced_forms = quad.class_group(D)
100  harmonic_sum = sum(1/mpmath.mpf(form[0]) for form in reduced_forms)
101  floatpre = nearest_integer(mpmath.pi*mpmath.sqrt(-D)*harmonic_sum / mpmath.log(10)) + 10
102  mpmath.mp.dps = floatpre
103  sqrtD = mpmath.sqrt(D)
104  for a, b, c in reduced_forms:
105  if b < 0:continue
106  tau = (-b + sqrtD)/(2*a)
107  f = mpmath.power(dedekind(2*tau, floatpre) / dedekind(tau, floatpre), 24)
108  j = ((256*f + 1)**3) / f
109  TC = tuple(T)
110  if b == a or c == a or b == 0:
111  T[0] = TC[0]*(-j)
112  if len(T) > 2:
113  for i in range(1, len(T)):
114  T[i] = TC[i - 1] + (TC[i]*(-j))
115  T.append(1)
116  elif T[1] != 0:
117  T[1] = TC[0] + (TC[1]*(-j))
118  T.append(1)
119  else:
120  T[1] = 1
121  else:
122  rej = -2*j.real
123  abj = j.real**2 + j.imag**2
124  T[0] = TC[0]*abj
125  T[1] = TC[0]*rej + TC[1]*abj
126  if len(T) > 2:
127  for i in range(2, len(T)):
128  T[i] = TC[i-2] + TC[i-1]*rej + TC[i]*abj
129  T.extend([TC[-2] + TC[-1]*rej, 1])
130  T = [(int(nearest_integer(c.real)) if hasattr(c, 'real') else c) for c in T]
131  return h, T
132 
133 def hilbert(D):
134  """
135  wrapper to split mpmath part and others.
136  """
137  if HAVE_NET:
138  try:
139  import urllib2
140  url = 'http://hilbert-class-polynomial.appspot.com/%d/' % D
141  result = urllib2.urlopen(url).read()
142  if result == 'null':
143  pass
144  else:
145  hcp = eval(result)
146  return len(hcp) - 1, hcp
147  except Exception:
148  pass
149 
150  if HAVE_MPMATH:
151  return hilbert_mpmath(D)
152 
153  raise prime.ImplementLimit("impliment limit")
154 
155 
156 def outputHP(absbound=DEFAULT_ABS_BOUND):
157  """
158  Output Hilbert class polynomials to file 'Hilbert_Poly.txt'.
159  """
160  f = open('Hilbert_Poly.txt', 'a')
161  for disc in disc_gen(absbound):
162  h, T = hilbert(disc)
163  f.write("%s %s %s\n" % (str(disc), str(h), str(T)))
164  f.close()
165 
166 
167 class Elliptic(object):
168  """
169  The class is for operations of elliptic curve points.
170  """
171  def __init__(self, coefficient, modulus):
172  self.modulus = modulus
173  self.a, self.b = coefficient
174  self.f = lambda x: ((x*x + self.a)*x + self.b) % self.modulus
175 
176  def add(self, P, Q):
177  """
178  Return P + Q
179  """
180  infinity = [0]
181  if P == infinity:
182  return Q
183  elif Q == infinity:
184  return P
185  else:
186  if P[0] == Q[0]:
187  if not (P[1] + Q[1]):
188  return infinity
189  s = 3*P[0]*P[0] + self.a
190  t = 2*P[1]
191  else:
192  s = Q[1] - P[1]
193  t = Q[0] - P[0]
194  m = s/t
195  x_3 = m*m - P[0] - Q[0]
196  return [x_3, m*(P[0] - x_3) - P[1]]
197 
198  def sub(self, P, Q):
199  """
200  return P - Q
201  """
202  infinity = [0]
203  if Q == infinity:
204  return P
205  R = [Q[0], -Q[1]]
206  if P == infinity:
207  return R
208  else:
209  return self.add(P, R)
210 
211  def mul(self, k, P):
212  """
213  this returns [k]*P
214  """
215  if k >= 0:
216  l = arith1.expand(k, 2)
217  Q = [0]
218  for j in range(len(l) - 1, -1, -1):
219  Q = self.add(Q, Q)
220  if l[j] == 1:
221  Q = self.add(Q, P)
222  return Q
223  else:
224  return self.sub([0], self.mul(-k, P))
225 
226  def choose_point(self):
227  """
228  Choose point on E_{a,b}(Z_n)
229  Algorithm 27 (Atkin-morain ECPP) Step5
230  """
231  n, f = self.modulus, self.f
232  x = bigrandom.randrange(n)
233  Q = f(x)
234  while arith1.legendre(Q, n) == -1:
235  x = bigrandom.randrange(n)
236  Q = f(x)
237  y = arith1.modsqrt(Q, n)
238  return [intresidue.IntegerResidueClass(t, n) for t in (x, y)]
239 
240 
241 def cmm(p):
242  """
243  Algorithm 25 (CM methods)
244  Input : a prime number p (>=3)
245  Output : curve parameters for CM curves
246  """
247  disc, g = _cmm_find_disc(p)[:2]
248  return [param for param in param_gen(disc, p, g)]
249 
250 def cmm_order(p):
251  """
252  Algorithm 25 (CM methods)
253  Input : a prime number p (>=3)
254  Output : curve parameters for CM curves and orders
255  """
256  disc, g, u, v = _cmm_find_disc(p)
257  m = [f(p, u, v) for f in orders.get(disc, default_orders)]
258  params = []
259  for param in param_gen(disc, p, g):
260  E = Elliptic(param, p)
261  base_point = E.choose_point()
262  for mi in m:
263  Q = E.mul(mi, base_point)
264  if Q == [0]:
265  params.append(param + (mi,))
266  m.remove(mi)
267  break
268  return params
269 
270 def _cmm_find_disc(p, absbound=DEFAULT_ABS_BOUND):
271  """
272  Input : a prime number p (>=3)
273  Output : (disc, g, u, v) where:
274  - disc: discriminant appropriate for CM method
275  - g: quasi primitive element of p
276  - u, v: a solution of u**2 + disc * v**2 = 4*p
277  """
278  for disc in disc_gen(absbound):
279  if arith1.legendre(disc, p) == 1:
280  try:
281  u, v = cornacchiamodify(disc, p)
282  g = quasi_primitive(p, disc == -3)
283  return disc, g, u, v
284  except ValueError:
285  continue
286  raise ValueError("no discriminant found")
287 
289  """
290  Algorithm 26 (Modified cornacchia)
291  Input : p be a prime and d be an integer such that d < 0 and d > -4p with
292  d = 0, 1 (mod 4)
293  Output : the solution of u^2 -d * v^2 = 4p.
294  """
295  q = 4 * p
296  if (d >= 0) or (d <= -q):
297  raise ValueError("invalid input")
298  if p == 2:
299  b = arith1.issquare(d + 8)
300  if b:
301  return (b, 1)
302  else:
303  raise ValueError("no solution")
304  if arith1.legendre(d, p) == -1:
305  raise ValueError("no solution")
306  x0 = arith1.modsqrt(d, p)
307  if (x0 - d) & 1:
308  x0 = p - x0
309  a = 2 * p
310  b = x0
311  l = arith1.floorsqrt(q)
312  while b > l:
313  a, b = b, a % b
314  c, r = divmod(q - b * b, -d)
315  if r:
316  raise ValueError("no solution")
317  t = arith1.issquare(c)
318  if t:
319  return (b, t)
320  else:
321  raise ValueError("no solution")
322 
323 def quasi_primitive(p, disc_is_minus3):
324  """
325  Return g s.t.
326  0) 1 < g < p
327  1) quadratic nonresidue modulo p
328  and
329  2) cubic nonresidue modulo p if p == 1 mod 3
330  with p >= 3, a prime.
331 
332  For p, not a prime, condition 1 means the Jacobi symbol
333  (g/p) != -1, and condition 2 means g doesn't have order
334  dividing (p - 1)//3.
335 
336  if disc_is_minus3 is True, cubic nonresidue check becomes
337  stricter so that g**((p-1)//3) must match with one of the
338  primitive third roots of unity.
339  """
340  third_order = (p - 1) // 3
341  for g in bigrange.range(2, p):
342  # if g is not quadratic nonresidue, then skipped.
343  legendre_jacobi = arith1.legendre(g, p)
344  if legendre_jacobi != -1:
345  continue
346  # if p != 1 mod 3 or g is cubic nonresidue, it's what we need.
347  if p % 3 != 1:
348  return g
349  cubic_symbol = pow(g, third_order, p)
350  if cubic_symbol != 1:
351  # g is cubic nonresidue.
352  if disc_is_minus3 and pow(cubic_symbol, 3, p) != 1:
353  # stricter check
354  continue
355  return g
356  else:
357  raise ValueError("p is not prime.")
358 
359 def param_gen(disc, n, g):
360  """
361  Return generator to generate curve params.
362  """
363  if disc == -3:
364  return _generate_params_for_minus3(n, g)
365  elif disc == -4:
366  return _generate_params_for_minus4(n, g)
367  else:
368  return _generate_params_for_general_disc(disc, n, g)
369 
371  """
372  Generate curve params for discriminant -3.
373  """
374  yield (0, n - 1)
375  gk = -1
376  for i in range(5):
377  gk = (gk*g)%n
378  yield (0, gk)
379 
381  """
382  Generate curve params for discriminant -4.
383  """
384  yield (n - 1, 0)
385  gk = -1
386  for i in range(3):
387  gk = (gk*g)%n
388  yield (gk, 0)
389 
391  """
392  Generate curve params for discriminant other than -3 or -4.
393  """
394  jr = equation.root_Fp([c % n for c in hilbert(disc)[-1]], n)
395  c = (jr * arith1.inverse(jr - 1728, n)) % n
396  r = (-3 * c) % n
397  s = (2 * c) % n
398  yield (r, s)
399  g2 = (g * g) % n
400  yield ((r * g2) % n, (s * g2 * g) % n)
401 
402 
403 def ecpp(n, era=None):
404  """
405  Algorithm 27 (Atkin-Morain ECPP)
406  Input : a natural number n
407  Output : If n is prime, retrun True. Else, return False
408  """
409  _log.info('n = %d' % n)
410  disc_and_order_info = _ecpp_find_disc(n, era)
411  if not disc_and_order_info:
412  return False
413  # disc, m = k*q, era
414  disc, m, k, q, era = disc_and_order_info
415  # non(quadratic|cubic) residue
416  g = quasi_primitive(n, disc == -3)
417 
418  try:
419  # find param corresponding to order m
420  for param in param_gen(disc, n, g):
421  E = Elliptic(param, n)
422  P = E.choose_point()
423  if E.mul(m, P) == [0]:
424  break
425  # find a point [k]P is not the unit.
426  U = E.mul(k, P)
427  while U == [0]:
428  P = E.choose_point()
429  U = E.mul(k, P)
430  # compute [q]([k]P)
431  V = E.mul(q, U)
432  except ZeroDivisionError:
433  return False
434  if V == [0]:
435  if q > 1000000000:
436  return ecpp(q, era)
437  elif prime.trialDivision(q):
438  return True
439  else:
440  return False
441  else:
442  return False
443 
444 def _ecpp_find_disc(n, era=None):
445  """
446  Return (disc, m, k, q, era) where:
447  - disc: appropriate discriminant for ecpp
448  - m: one of possible orders for disc
449  - k: smooth part of m
450  - q: rough part of m which is greater than (n^(1/4) + 1)^2
451  - era: list of small primes
452  """
453  up = (arith1.floorpowerroot(n, 4) + 1)**2
454  for disc in disc_gen(n):
455  legendre = arith1.legendre(disc, n)
456  if legendre == 0:
457  return False
458  elif legendre == -1:
459  continue
460  # legendre == 1:
461  try:
462  u, v = cornacchiamodify(disc, n)
463  except ValueError:
464  continue
465  for m in (f(n, u, v) for f in orders.get(disc, default_orders)):
466  k, q, era = _factor_order(m, up, era)
467  if k:
468  return disc, m, k, q, era
469  else:
470  # no usable discriminat found
471  return False
472 
473 def _factor_order(m, u, era=None):
474  """
475  Return triple (k, q, primes) if m is factored as m = kq,
476  where k > 1 and q is a probable prime > u.
477  u is expected to be (n^(1/4)+1)^2 for Atkin-Morain ECPP.
478  If m is not successfully factored into desired form,
479  return (False, False, primes).
480  The third component of both cases is a list of primes
481  to do trial division.
482 
483  Algorithm 27 (Atkin-morain ECPP) Step3
484  """
485  if era is None:
486  v = min(500000, int(log(m)))
487  era = factor.mpqs.eratosthenes(v)
488  k = 1
489  q = m
490  for p in era:
491  if p > u:
492  break
493  e, q = arith1.vp(q, p)
494  k *= p**e
495  if q < u or k == 1:
496  return False, False, era
497  if not prime.millerRabin(q):
498  return False, False, era
499  return k, q, era
500 
501 def next_disc(d, absbound):
502  """
503  Return fundamental discriminant D, such that D < d.
504  If there is no D with |d| < |D| < absbound, return False
505  """
506  # -disc % 16
507  negdisc_mod16 = (3, 4, 7, 8, 11, 15)
508  for negdisc in bigrange.range(-d + 1, absbound):
509  if negdisc & 15 not in negdisc_mod16:
510  continue
511  if negdisc & 1 and not squarefree.trial_division(negdisc):
512  continue
513  if arith1.issquare(negdisc):
514  continue
515  return -negdisc
516  return False
517 
518 def disc_gen(absbound=DEFAULT_ABS_BOUND):
519  """
520  Generate negative discriminants.
521 
522  The order of discriminants depends on how to produce them.
523  By default, discriminants are in ascending order of class
524  numbers while reading from precomputed data file, then
525  in descending order of discriminant itself.
526 
527  If discriminant reach the absbound, then StopIteration will be
528  raised. The default value of absbound is DEFAULT_ABS_BOUND.
529  """
530  csvfile = open(os.path.join(DATADIR, "discriminant.csv"))
531  for disc_str in csv.reader(csvfile):
532  disc = int(disc_str[0])
533  if -disc >= absbound:
534  raise StopIteration("absbound reached")
535  yield disc
536  disc = next_disc(disc, absbound)
537  while disc:
538  yield disc
539  disc = next_disc(disc, absbound)
540  raise StopIteration("absbound reached")
nzmath.factor
Definition: __init__.py:1
nzmath.ecpp.Elliptic.add
def add(self, P, Q)
Definition: ecpp.py:176
nzmath.bigrandom
Definition: bigrandom.py:1
nzmath.ecpp.Elliptic.__init__
def __init__(self, coefficient, modulus)
Definition: ecpp.py:171
nzmath.bigrange.range
def range(start, stop=None, step=None)
Definition: bigrange.py:19
nzmath.ecpp.Elliptic.f
f
Definition: ecpp.py:174
nzmath.ecpp.quasi_primitive
def quasi_primitive(p, disc_is_minus3)
Definition: ecpp.py:323
nzmath.ecpp.log
log
Definition: ecpp.py:20
nzmath.ecpp._generate_params_for_minus4
def _generate_params_for_minus4(n, g)
Definition: ecpp.py:380
nzmath.ecpp.Elliptic.modulus
modulus
Definition: ecpp.py:172
nzmath.ecpp.hilbert_mpmath
def hilbert_mpmath(D)
Definition: ecpp.py:89
nzmath.ecpp.cmm
def cmm(p)
Definition: ecpp.py:241
nzmath.ecpp._cmm_find_disc
def _cmm_find_disc(p, absbound=DEFAULT_ABS_BOUND)
Definition: ecpp.py:270
nzmath.squarefree
Definition: squarefree.py:1
nzmath.ecpp.cmm_order
def cmm_order(p)
Definition: ecpp.py:250
nzmath.ecpp.disc_gen
def disc_gen(absbound=DEFAULT_ABS_BOUND)
Definition: ecpp.py:518
nzmath.ecpp.dedekind
def dedekind(tau, floatpre)
Definition: ecpp.py:54
nzmath.ecpp.hilbert
def hilbert(D)
Definition: ecpp.py:133
nzmath.ecpp.nearest_integer
def nearest_integer(x)
Definition: ecpp.py:22
nzmath.intresidue.IntegerResidueClass
Definition: intresidue.py:12
nzmath.ecpp.Elliptic.sub
def sub(self, P, Q)
Definition: ecpp.py:198
nzmath.prime.ImplementLimit
Definition: prime.py:1207
nzmath.ecpp._generate_params_for_general_disc
def _generate_params_for_general_disc(disc, n, g)
Definition: ecpp.py:390
nzmath.ecpp.outputHP
def outputHP(absbound=DEFAULT_ABS_BOUND)
Definition: ecpp.py:156
nzmath.ecpp.Elliptic.mul
def mul(self, k, P)
Definition: ecpp.py:211
nzmath.ecpp.Elliptic.b
b
Definition: ecpp.py:173
nzmath.ecpp._generate_params_for_minus3
def _generate_params_for_minus3(n, g)
Definition: ecpp.py:370
nzmath.intresidue
Definition: intresidue.py:1
nzmath.quad
Definition: quad.py:1
nzmath.arith1
Definition: arith1.py:1
nzmath.equation
Definition: equation.py:1
nzmath.ecpp.Elliptic.choose_point
def choose_point(self)
Definition: ecpp.py:226
nzmath.ecpp._factor_order
def _factor_order(m, u, era=None)
Definition: ecpp.py:473
nzmath.ecpp.cornacchiamodify
def cornacchiamodify(d, p)
Definition: ecpp.py:288
nzmath.prime
Definition: prime.py:1
nzmath.compatibility
Definition: compatibility.py:1
nzmath.ecpp.ecpp
def ecpp(n, era=None)
Definition: ecpp.py:403
nzmath.ecpp.param_gen
def param_gen(disc, n, g)
Definition: ecpp.py:359
nzmath.bigrange
Definition: bigrange.py:1
nzmath.ecpp.next_disc
def next_disc(d, absbound)
Definition: ecpp.py:501
nzmath.ecpp.Elliptic
Definition: ecpp.py:167
nzmath.ecpp._ecpp_find_disc
def _ecpp_find_disc(n, era=None)
Definition: ecpp.py:444
nzmath.config
Definition: config.py:1