"Fossies" - the Fresh Open Source Software Archive

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

    1 """
    2 Elliptic Curve Method Factorization
    3 
    4 It's using Montgomery's parametrization.
    5 """
    6 
    7 from __future__ import division
    8 import logging
    9 import random
   10 import nzmath.arith1 as _arith1
   11 import nzmath.gcd as _gcd
   12 import nzmath.prime as _prime
   13 
   14 _log = logging.getLogger('nzmath.factor.ecm')
   15 
   16 
   17 # curve types
   18 S = SUYAMA = 1
   19 B = BERNSTEIN = 2
   20 A1 = ASUNCION1 = 5
   21 A2 = ASUNCION2 = 6
   22 A3 = ASUNCION3 = 8
   23 A4 = ASUNCION4 = 9
   24 A5 = ASUNCION5 = 10
   25 N3 = NAKAMURA3 = 13
   26 
   27 class Curve (object):
   28     """
   29     Elliptic curves for factorization.
   30     """
   31 
   32     _CURVE_TYPES = (S, B, A1, A2, A3, A4, A5, N3)
   33     
   34     def __init__(self, c):
   35         """
   36         Initialize a Curve object with Montgomery's parameter c.
   37         """
   38         self.c = c
   39         self.c2 = (c + 2)//4
   40 
   41     @classmethod
   42     def get_random_curve_with_point(cls, curve_type, n, bounds):
   43         """
   44         Return the curve with parameter C and a point Q on the curve,
   45         according to the curve_type, factorization target n and the
   46         bounds for stages.
   47 
   48         curve_type should be one of the module constants corresponding
   49         to parameters:
   50           S: Suyama's parameter selection strategy
   51           B: Bernstein's [2:1], [16,18,4,2]
   52           A1: Asuncion's [2:1], [4,14,1,1]
   53           A2: Asuncion's [2:1], [16,174,4,41]
   54           A3: Asuncion's [3:1], [9,48,1,2]
   55           A4: Asuncion's [3:1], [9,39,1,1]
   56           A5: Asuncion's [4:1], [16,84,1,1]
   57           N3: Nakamura's [2:1], [28,22,7,3]
   58           This is a class method.
   59         """
   60         bound = bounds.first
   61         if curve_type not in cls._CURVE_TYPES:
   62             raise ValueError("Input curve_type is wrong.")
   63         if curve_type == SUYAMA:
   64             t = n
   65             while _gcd.gcd(t, n) != 1:
   66                 sigma = random.randrange(6, bound + 1)
   67                 u, v = (sigma**2 - 5) % n, (4*sigma) % n
   68                 t = 4*(u**3)*v
   69             d = _arith1.inverse(t, n)
   70             curve = cls((((u - v)**3 * (3*u + v)) * d - 2) % n)
   71             start_point = Point(pow(u, 3, n), pow(v, 3, n))
   72         elif curve_type == BERNSTEIN:
   73             d = random.randrange(1, bound + 1)
   74             start_point = Point(2, 1)
   75             curve = cls((4*d + 2) % n)
   76         elif curve_type == ASUNCION1:
   77             d = random.randrange(1, bound + 1)
   78             start_point = Point(2, 1)
   79             curve = cls((d + 1) % n)
   80         elif curve_type == ASUNCION2:
   81             d = random.randrange(1, bound + 1)
   82             start_point = Point(2, 1)
   83             curve = cls((4*d + 41) % n)
   84         elif curve_type == ASUNCION3:
   85             d = random.randrange(1, bound + 1)
   86             start_point = Point(3, 1)
   87             curve = cls((d + 2) % n)
   88         elif curve_type == ASUNCION4:
   89             d = random.randrange(1, bound + 1)
   90             start_point = Point(3, 1)
   91             curve = cls((d + 1) % n)
   92         elif curve_type == ASUNCION5:
   93             d = random.randrange(1, bound + 1)
   94             start_point = Point(4, 1)
   95             curve = cls((d + 1) % n)
   96         elif curve_type == NAKAMURA3:
   97             d = random.randrange(1, bound + 1)
   98             start_point = Point(2, 1)
   99             curve = cls((7*d + 3) % n)
  100         return curve, start_point
  101 
  102 
  103 class Point (tuple):
  104     """
  105     x and z projective coordinates of elliptic curves.
  106     """
  107     def __new__(cls, x, z):
  108         """
  109         Create a new Point object.
  110         """
  111         return super(Point, cls).__new__(Point, (x, z))
  112 
  113     def _get_x(self):
  114         "getter for x"
  115         return self[0]
  116 
  117     def _get_z(self):
  118         "getter for z"
  119         return self[1]
  120 
  121     x = property(_get_x, None, None, "x")
  122     z = property(_get_z, None, None, "z")
  123 
  124 POINT_AT_INFINITY = Point(0, 0)
  125 
  126 
  127 class Bounds (object):
  128     """
  129     Bounds for ECM.
  130 
  131     public attributes:
  132       first: bound for first stage
  133       second: bound for second stage
  134     """
  135     def __init__(self, n):
  136         """
  137         Create a Bounds object from target number n.
  138         """
  139         if _arith1.log(n, 10) < 20:
  140             self.first = 1000
  141         elif _arith1.log(n, 10) < 25:
  142             self.first = 10000
  143         else:
  144             self.first = 100000
  145         self.second = self.first * 50
  146 
  147     def increment(self, scale=10):
  148         """
  149         Multiply both bounds by optional argument scale (default to 10).
  150         """
  151         self.first *= scale
  152         self.second *= scale
  153 
  154     def __getitem__(self, key):
  155         """
  156         Return first by key 1, and second key 2.
  157         Otherwise KeyError is raised.
  158         """
  159         if key == 1:
  160             return self.first
  161         if key == 2:
  162             return self.second
  163         raise KeyError("Key must be 1 or 2.")
  164 
  165     def __hash__(self):
  166         """
  167         Bounds objects are mutable and thus unhashable.
  168         """
  169         raise TypeError("Bounds objects are unhashable")
  170 
  171     def __str__(self):
  172         """
  173         (first, second)
  174         """
  175         return "(%d, %d)" % (self.first, self.second)
  176 
  177 
  178 def ecm(n, curve_type=A1, incs=3, trials=20, **options):
  179     """
  180     Find a factor of n with Elliptic Curve Method.
  181     An unsuccessful factorization returns 1.
  182 
  183     There are a few optional arguments.
  184 
  185     By 'curve_type', the function choose a family of curves.
  186     Please use a module constant to specify the curve_type.
  187 
  188     The second optional argument 'incs' specifies the number of times
  189     for changing bounds.  The function repeats factorization trials
  190     several times changing curves with a fixed bounds.
  191 
  192     The last optional argument 'trials' can control how quickly move
  193     on to the next higher bounds.
  194     """
  195     # verbosity
  196     verbose = options.get('verbose', False)
  197     if verbose:
  198         _log.setLevel(logging.DEBUG)
  199         _log.debug("verbose")
  200     else:
  201         _log.setLevel(logging.NOTSET)
  202 
  203     # trivial checks
  204     if _prime.primeq(n):
  205         _log.info("%d is prime!" % n)
  206         return n
  207     if _gcd.gcd(n, 6) != 1: 
  208         _log.info("%d is not coprime to 6!" % n)
  209         if n % 2 == 0:
  210             return 2
  211         if n % 3 == 0:
  212             return 3
  213 
  214     # main loop
  215     bounds = Bounds(n)
  216     for inc in range(incs):
  217         _log.info("bounds increment %d times" % inc)
  218         _log.debug("Bounds B1, B2 = %s" % (bounds,))
  219         for trial in range(trials):
  220             _log.info("Trial #: %d" % (trial + 1))
  221             curve, point = Curve.get_random_curve_with_point(curve_type, n, bounds)
  222             _log.debug("Curve param: %d" % curve.c)
  223             g = stage1(n, bounds, curve, point)
  224             if 1 < g < n:
  225                 return g
  226             g = stage2(n, bounds, curve, point)
  227             if 1 < g < n:
  228                 return g
  229         bounds.increment()
  230     _log.info("ECM 2 step test failed!")
  231 
  232     if not verbose:
  233         _log.setLevel(logging.DEBUG)
  234     return 1
  235 
  236 def stage1(n, bounds, C, Q):
  237     """
  238     ECM stage 1 for factoring n.
  239     The upper bound for primes to be tested is bounds.first.
  240     It uses curve C and starting point Q.
  241     """
  242     for p in _prime.generator_eratosthenes(bounds.first):
  243         q = p
  244         while q < bounds.first:
  245             Q = mul(Q, p, C, n)
  246             q *= p
  247     g = _gcd.gcd(Q.z, n)
  248     _log.debug("Stage 1: %d" % g)
  249     return g
  250 
  251 def stage2(n, bounds, C, Q):
  252     """
  253     ECM stage 2 for factoring n.
  254     The upper bounds for primes to be tested are stored in bounds.
  255     It uses curve C and starting point Q.
  256     """
  257     d1 = bounds.first + random.randrange(1, 16)
  258     d2 = 0
  259     while _gcd.gcd(d1, d2) != 1:
  260         d2 = random.randrange(2, d1//5 + 1) # We want to keep d2 small
  261     for i in range(1, bounds.second//d1):
  262         if _gcd.gcd(i, d2) != 1:
  263             continue
  264         for j in range(1, d1//2):
  265             if _gcd.gcd(j, d1) == 1:
  266                 Q = mul(Q, i*d1 + j*d2, C, n)
  267                 if i*d1 - j*d2 > bounds.first:
  268                     Q = mul(Q, i*d1 - j*d2, C, n)
  269         g = _gcd.gcd(Q.z, n)
  270         if 1 < g < n:
  271             _log.debug("Stage 2: %d" % g)
  272             return g
  273     _log.debug("Stage 2: %d" % g)
  274     return 1
  275 
  276 def mul(Q, x, C, n):
  277     """
  278     Return x*Q on the curve C mod n.
  279 
  280     m*Q and (m+1)*Q are being tracked in the main loop.
  281     """
  282     if x == 0:
  283         return POINT_AT_INFINITY
  284     if x == 1:
  285         return Q
  286     if x == 2:
  287         return double(Q, C, n)
  288     minor = Q
  289     major = double(Q, C, n)
  290     binary = _arith1.expand(x, 2)
  291     lastbit, binary = binary[0], binary[1:]
  292     while binary:
  293         if binary.pop() == 1:
  294             minor = add(major, minor, Q, n)
  295             major = double(major, C, n)
  296         else:
  297             major = add(minor, major, Q, n)
  298             minor = double(minor, C, n)
  299     if lastbit:
  300         return add(minor, major, Q, n)
  301     return double(minor, C, n)
  302 
  303 # Montgomery's Arithmetic
  304 
  305 def double(Q, C, n):
  306     """
  307     Return the doubled Q on the curve C mod n.
  308     """
  309     u = (Q.x + Q.z)**2
  310     v = (Q.x - Q.z)**2
  311     t = C.c2 * (u - v) + v
  312     return Point(u * v % n, (u - v) * t % n)
  313 
  314 def add(P, Q, M, n):
  315     """
  316     Return the sum of P and Q mod n with a help of auxiliary point M.
  317     M must be P - Q. The curve is not explicitly necessary.
  318     """
  319     u = (P.x + P.z) * (Q.x - Q.z)
  320     v = (P.x - P.z) * (Q.x + Q.z)
  321     w = (u + v)**2
  322     t = (u - v)**2
  323     return Point(M.z * w % n, M.x * t % n)