## "Fossies" - the Fresh Open Source Software Archive

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

```    1 """
2 All methods defined here return one of a factor of given integer.
3 When 1 is returned, the method has failed to factor,
4 but 1 is a factor anyway.
5
6 'verbose' boolean flag can be specified for verbose reports.
7 """
8
9 import logging
10 import nzmath.arith1 as arith1
11 import nzmath.bigrandom as bigrandom
12 import nzmath.gcd as gcd
13 import nzmath.prime as prime
14
15 _log = logging.getLogger('nzmath.factor.find')
16
17 # Pollard's rho method
18 def rhomethod(n, **options):
19     """
20     Find a non-trivial factor of n using Pollard's rho algorithm.
21     The implementation refers the explanation in C.Pomerance's book.
22     """
23     # verbosity
24     verbose = options.get('verbose', False)
25     if not verbose:
26         _silence()
27
28     if n <= 3:
29         return 1
30
31     g = n
32     while g == n:
33         # x^2 + a is iterated. Starting value x = u.
34         a = bigrandom.randrange(1, n-2)
35         u = v = bigrandom.randrange(0, n-1)
36         _log.info("%d %d" % (a, u))
37         g = gcd.gcd((v**2 + v + a) % n - u, n)
38         while g == 1:
39             u = (u**2 + a) % n
40             v = ((pow(v, 2, n) + a)**2 + a) % n
41             g = gcd.gcd(v - u, n)
42     if not verbose:
43         _verbose()
44     return g
45
46 # p-1 method
47 def pmom(n, **options):
48     """
49     This function tries to find a non-trivial factor of n using
50     Algorithm 8.8.2 (p-1 first stage) of Cohen's book.
51     In case of N = pow(2,i), this program will not terminate.
52     """
53     # verbosity
54     verbose = options.get('verbose', False)
55     if not verbose:
56         _silence()
57
58     # initialize
59     x = y = 2
60     primes = []
61     if 'B' in options:
62         B = options['B']
63     else:
64         B = 10000
65
66     for q in prime.generator():
67         primes.append(q)
68         if q > B:
69             if gcd.gcd(x-1, n) == 1:
70                 if not verbose:
71                     _verbose()
72                 return 1
73             x = y
74             break
75         q1 = q
76         l = B // q
77         while q1 <= l:
78             q1 *= q
79         x = pow(x, q1, n)
80         if len(primes) >= 20:
81             if gcd.gcd(x-1, n) == 1:
82                 primes, y = [], x
83             else:
84                 x = y
85                 break
86
87     for q in primes:
88         q1 = q
89         while q1 <= B:
90             x = pow(x, q, n)
91             g = gcd.gcd(x-1, n)
92             if g != 1:
93                 if not verbose:
94                     _verbose()
95                 if g == n:
96                     return 1
97                 return g
98             q1 *= q
99
100 def trialDivision(n, **options):
101     """
102     Return a factor of given integer by trial division.
103
104     options can be either:
105     1) 'start' and 'stop' as range parameters.
106     2) 'iterator' as an iterator of primes.
107     If both options are not given, prime factor is searched from 2
108     to the square root of the given integer.
109     """
110     # verbosity
111     verbose = options.get('verbose', False)
112     if not verbose:
113         _silence()
114
115     if 'start' in options and 'stop' in options:
116         if 'step' in options:
117             trials = range(options['start'], options['stop'], options['step'])
118         else:
119             trials = range(options['start'], options['stop'])
120     elif 'iterator' in options:
121         trials = options['iterator']
122     elif n < 1000000:
123         trials = prime.generator_eratosthenes(arith1.floorsqrt(n))
124     else:
125         trials = prime.generator()
126
127     limit = arith1.floorsqrt(n)
128     for p in trials:
129         if limit < p:
130             break
131         if 0 == n % p:
132             if not verbose:
133                 _verbose()
134             return p
135     if not verbose:
136         _verbose()
137     return 1
138
139 def _silence():
140     """
141     Stop verbose outputs.
142     """
143     _log.setLevel(logging.NOTSET)
144
145 def _verbose():
146     """
147     Stop silencing.
148     """
149     _log.setLevel(logging.DEBUG)
```