"Fossies" - the Fresh Open Source Software Archive

Member "NZMATH-1.2.0/nzmath/ecpp.py" (19 Nov 2012, 15612 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 "ecpp.py" see the Fossies "Dox" file reference documentation.

    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
   14 import nzmath.compatibility
   15 from nzmath.config import DATADIR, HAVE_MPMATH, HAVE_NET
   16 
   17 if HAVE_MPMATH:
   18     import mpmath
   19 
   20     log = mpmath.log
   21 
   22     def nearest_integer(x):
   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 
   89 def hilbert_mpmath(D):
   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 
  288 def cornacchiamodify(d, p):
  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 
  370 def _generate_params_for_minus3(n, g):
  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 
  380 def _generate_params_for_minus4(n, g):
  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 
  390 def _generate_params_for_general_disc(disc, n, g):
  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")