Coverage for sympy/ntheory/factor_.py : 26%

""" 

Integer factorization 

""" 

from __future__ import print_function, division 

import random 

import math 

from .primetest import isprime 

from .generate import sieve, primerange, nextprime 

from sympy.core import sympify 

from sympy.core.evalf import bitcount 

from sympy.core.logic import fuzzy_and 

from sympy.core.numbers import igcd, ilcm, Rational 

from sympy.core.power import integer_nthroot, Pow 

from sympy.core.mul import Mul 

from sympy.core.compatibility import as_int, SYMPY_INTS, range 

from sympy.core.singleton import S 

from sympy.core.function import Function 

small_trailing = [i and max(int(not i % 2**j) and j for j in range(1, 8)) 

for i in range(256)] 

def smoothness(n): 

""" 

Return the B-smooth and B-power smooth values of n. 

The smoothness of n is the largest prime factor of n; the power- 

smoothness is the largest divisor raised to its multiplicity. 

>>> from sympy.ntheory.factor_ import smoothness 

>>> smoothness(2**7*3**2) 

(3, 128) 

>>> smoothness(2**4*13) 

(13, 16) 

>>> smoothness(2) 

(2, 2) 

See Also 

======== 

factorint, smoothness_p 

""" 

if n == 1: 

return (1, 1) # not prime, but otherwise this causes headaches 

facs = factorint(n) 

return max(facs), max(m**facs[m] for m in facs) 

def smoothness_p(n, m=-1, power=0, visual=None): 

""" 

Return a list of [m, (p, (M, sm(p + m), psm(p + m)))...] 

where: 

1. p**M is the base-p divisor of n 

2. sm(p + m) is the smoothness of p + m (m = -1 by default) 

3. psm(p + m) is the power smoothness of p + m 

The list is sorted according to smoothness (default) or by power smoothness 

if power=1. 

The smoothness of the numbers to the left (m = -1) or right (m = 1) of a 

factor govern the results that are obtained from the p +/- 1 type factoring 

methods. 

>>> from sympy.ntheory.factor_ import smoothness_p, factorint 

>>> smoothness_p(10431, m=1) 

(1, [(3, (2, 2, 4)), (19, (1, 5, 5)), (61, (1, 31, 31))]) 

>>> smoothness_p(10431) 

(-1, [(3, (2, 2, 2)), (19, (1, 3, 9)), (61, (1, 5, 5))]) 

>>> smoothness_p(10431, power=1) 

(-1, [(3, (2, 2, 2)), (61, (1, 5, 5)), (19, (1, 3, 9))]) 

If visual=True then an annotated string will be returned: 

>>> print(smoothness_p(21477639576571, visual=1)) 

p**i=4410317**1 has p-1 B=1787, B-pow=1787 

p**i=4869863**1 has p-1 B=2434931, B-pow=2434931 

This string can also be generated directly from a factorization dictionary 

and vice versa: 

>>> factorint(17*9) 

{3: 2, 17: 1} 

>>> smoothness_p(_) 

'p**i=3**2 has p-1 B=2, B-pow=2\\np**i=17**1 has p-1 B=2, B-pow=16' 

>>> smoothness_p(_) 

{3: 2, 17: 1} 

The table of the output logic is: 

====== ====== ======= ======= 

| Visual 

------ ---------------------- 

Input True False other 

====== ====== ======= ======= 

dict str tuple str 

str str tuple dict 

tuple str tuple str 

n str tuple tuple 

mul str tuple tuple 

====== ====== ======= ======= 

See Also 

======== 

factorint, smoothness 

""" 

from sympy.utilities import flatten 

# visual must be True, False or other (stored as None) 

if visual in (1, 0): 

visual = bool(visual) 

elif visual not in (True, False): 

visual = None 

if type(n) is str: 

if visual: 

return n 

d = {} 

for li in n.splitlines(): 

k, v = [int(i) for i in 

li.split('has')[0].split('=')[1].split('**')] 

d[k] = v 

if visual is not True and visual is not False: 

return d 

return smoothness_p(d, visual=False) 

elif type(n) is not tuple: 

facs = factorint(n, visual=False) 

if power: 

k = -1 

else: 

k = 1 

if type(n) is not tuple: 

rv = (m, sorted([(f, 

tuple([M] + list(smoothness(f + m)))) 

for f, M in [i for i in facs.items()]], 

key=lambda x: (x[1][k], x[0]))) 

else: 

rv = n 

if visual is False or (visual is not True) and (type(n) in [int, Mul]): 

return rv 

lines = [] 

for dat in rv[1]: 

dat = flatten(dat) 

dat.insert(2, m) 

lines.append('p**i=%i**%i has p%+i B=%i, B-pow=%i' % tuple(dat)) 

return '\n'.join(lines) 

def trailing(n): 

"""Count the number of trailing zero digits in the binary 

representation of n, i.e. determine the largest power of 2 

that divides n. 

Examples 

======== 

>>> from sympy import trailing 

>>> trailing(128) 

7 

>>> trailing(63) 

0 

""" 

n = int(n) 

if not n: 

return 0 

low_byte = n & 0xff 

if low_byte: 

return small_trailing[low_byte] 

# 2**m is quick for z up through 2**30 

z = bitcount(n) - 1 

if isinstance(z, SYMPY_INTS): 

if n == 1 << z: 

return z 

t = 0 

p = 8 

while not n & 1: 

while not n & ((1 << p) - 1): 

n >>= p 

t += p 

p *= 2 

p //= 2 

return t 

def multiplicity(p, n): 

""" 

Find the greatest integer m such that p**m divides n. 

Examples 

======== 

>>> from sympy.ntheory import multiplicity 

>>> from sympy.core.numbers import Rational as R 

>>> [multiplicity(5, n) for n in [8, 5, 25, 125, 250]] 

[0, 1, 2, 3, 3] 

>>> multiplicity(3, R(1, 9)) 

-2 

""" 

try: 

p, n = as_int(p), as_int(n) 

except ValueError: 

if all(isinstance(i, (SYMPY_INTS, Rational)) for i in (p, n)): 

try: 

p = Rational(p) 

n = Rational(n) 

if p.q == 1: 

if n.p == 1: 

return -multiplicity(p.p, n.q) 

return S.Zero 

elif p.p == 1: 

return multiplicity(p.q, n.q) 

else: 

like = min( 

multiplicity(p.p, n.p), 

multiplicity(p.q, n.q)) 

cross = min( 

multiplicity(p.q, n.p), 

multiplicity(p.p, n.q)) 

return like - cross 

except AttributeError: 

pass 

raise ValueError('expecting ints or fractions, got %s and %s' % (p, n)) 

if n == 0: 

raise ValueError('no such integer exists: multiplicity of %s is not-defined' %(n)) 

if p == 2: 

return trailing(n) 

if p < 2: 

raise ValueError('p must be an integer, 2 or larger, but got %s' % p) 

if p == n: 

return 1 

m = 0 

n, rem = divmod(n, p) 

while not rem: 

m += 1 

if m > 5: 

# The multiplicity could be very large. Better 

# to increment in powers of two 

e = 2 

while 1: 

ppow = p**e 

if ppow < n: 

nnew, rem = divmod(n, ppow) 

if not rem: 

m += e 

e *= 2 

n = nnew 

continue 

return m + multiplicity(p, n) 

n, rem = divmod(n, p) 

return m 

def perfect_power(n, candidates=None, big=True, factor=True): 

""" 

Return ``(b, e)`` such that ``n`` == ``b**e`` if ``n`` is a 

perfect power; otherwise return ``False``. 

By default, the base is recursively decomposed and the exponents 

collected so the largest possible ``e`` is sought. If ``big=False`` 

then the smallest possible ``e`` (thus prime) will be chosen. 

If ``candidates`` for exponents are given, they are assumed to be sorted 

and the first one that is larger than the computed maximum will signal 

failure for the routine. 

If ``factor=True`` then simultaneous factorization of n is attempted 

since finding a factor indicates the only possible root for n. This 

is True by default since only a few small factors will be tested in 

the course of searching for the perfect power. 

Examples 

======== 

>>> from sympy import perfect_power 

>>> perfect_power(16) 

(2, 4) 

>>> perfect_power(16, big = False) 

(4, 2) 

""" 

n = int(n) 

if n < 3: 

return False 

logn = math.log(n, 2) 

max_possible = int(logn) + 2 # only check values less than this 

not_square = n % 10 in [2, 3, 7, 8] # squares cannot end in 2, 3, 7, 8 

if not candidates: 

candidates = primerange(2 + not_square, max_possible) 

afactor = 2 + n % 2 

for e in candidates: 

if e < 3: 

if e == 1 or e == 2 and not_square: 

continue 

if e > max_possible: 

return False 

# see if there is a factor present 

if factor: 

if n % afactor == 0: 

# find what the potential power is 

if afactor == 2: 

e = trailing(n) 

else: 

e = multiplicity(afactor, n) 

# if it's a trivial power we are done 

if e == 1: 

return False 

# maybe the bth root of n is exact 

r, exact = integer_nthroot(n, e) 

if not exact: 

# then remove this factor and check to see if 

# any of e's factors are a common exponent; if 

# not then it's not a perfect power 

n //= afactor**e 

m = perfect_power(n, candidates=primefactors(e), big=big) 

if m is False: 

return False 

else: 

r, m = m 

# adjust the two exponents so the bases can 

# be combined 

g = igcd(m, e) 

if g == 1: 

return False 

m //= g 

e //= g 

r, e = r**m*afactor**e, g 

if not big: 

e0 = primefactors(e) 

if len(e0) > 1 or e0[0] != e: 

e0 = e0[0] 

r, e = r**(e//e0), e0 

return r, e 

else: 

# get the next factor ready for the next pass through the loop 

afactor = nextprime(afactor) 

# Weed out downright impossible candidates 

if logn/e < 40: 

b = 2.0**(logn/e) 

if abs(int(b + 0.5) - b) > 0.01: 

continue 

# now see if the plausible e makes a perfect power 

r, exact = integer_nthroot(n, e) 

if exact: 

if big: 

m = perfect_power(r, big=big, factor=factor) 

if m is not False: 

r, e = m[0], e*m[1] 

return int(r), e 

else: 

return False 

def pollard_rho(n, s=2, a=1, retries=5, seed=1234, max_steps=None, F=None): 

r""" 

Use Pollard's rho method to try to extract a nontrivial factor 

of ``n``. The returned factor may be a composite number. If no 

factor is found, ``None`` is returned. 

The algorithm generates pseudo-random values of x with a generator 

function, replacing x with F(x). If F is not supplied then the 

function x**2 + ``a`` is used. The first value supplied to F(x) is ``s``. 

Upon failure (if ``retries`` is > 0) a new ``a`` and ``s`` will be 

supplied; the ``a`` will be ignored if F was supplied. 

The sequence of numbers generated by such functions generally have a 

a lead-up to some number and then loop around back to that number and 

begin to repeat the sequence, e.g. 1, 2, 3, 4, 5, 3, 4, 5 -- this leader 

and loop look a bit like the Greek letter rho, and thus the name, 'rho'. 

For a given function, very different leader-loop values can be obtained 

so it is a good idea to allow for retries: 

>>> from sympy.ntheory.generate import cycle_length 

>>> n = 16843009 

>>> F = lambda x:(2048*pow(x, 2, n) + 32767) % n 

>>> for s in range(5): 

... print('loop length = %4i; leader length = %3i' % next(cycle_length(F, s))) 

... 

loop length = 2489; leader length = 42 

loop length = 78; leader length = 120 

loop length = 1482; leader length = 99 

loop length = 1482; leader length = 285 

loop length = 1482; leader length = 100 

Here is an explicit example where there is a two element leadup to 

a sequence of 3 numbers (11, 14, 4) that then repeat: 

>>> x=2 

>>> for i in range(9): 

... x=(x**2+12)%17 

... print(x) 

... 

16 

13 

11 

14 

4 

11 

14 

4 

11 

>>> next(cycle_length(lambda x: (x**2+12)%17, 2)) 

(3, 2) 

>>> list(cycle_length(lambda x: (x**2+12)%17, 2, values=True)) 

[16, 13, 11, 14, 4] 

Instead of checking the differences of all generated values for a gcd 

with n, only the kth and 2*kth numbers are checked, e.g. 1st and 2nd, 

2nd and 4th, 3rd and 6th until it has been detected that the loop has been 

traversed. Loops may be many thousands of steps long before rho finds a 

factor or reports failure. If ``max_steps`` is specified, the iteration 

is cancelled with a failure after the specified number of steps. 

Examples 

======== 

>>> from sympy import pollard_rho 

>>> n=16843009 

>>> F=lambda x:(2048*pow(x,2,n) + 32767) % n 

>>> pollard_rho(n, F=F) 

257 

Use the default setting with a bad value of ``a`` and no retries: 

>>> pollard_rho(n, a=n-2, retries=0) 

If retries is > 0 then perhaps the problem will correct itself when 

new values are generated for a: 

>>> pollard_rho(n, a=n-2, retries=1) 

257 

References 

========== 

- Richard Crandall & Carl Pomerance (2005), "Prime Numbers: 

A Computational Perspective", Springer, 2nd edition, 229-231 

""" 

n = int(n) 

if n < 5: 

raise ValueError('pollard_rho should receive n > 4') 

prng = random.Random(seed + retries) 

V = s 

for i in range(retries + 1): 

U = V 

if not F: 

F = lambda x: (pow(x, 2, n) + a) % n 

j = 0 

while 1: 

if max_steps and (j > max_steps): 

break 

j += 1 

U = F(U) 

V = F(F(V)) # V is 2x further along than U 

g = igcd(U - V, n) 

if g == 1: 

continue 

if g == n: 

break 

return int(g) 

V = prng.randint(0, n - 1) 

a = prng.randint(1, n - 3) # for x**2 + a, a%n should not be 0 or -2 

F = None 

return None 

def pollard_pm1(n, B=10, a=2, retries=0, seed=1234): 

""" 

Use Pollard's p-1 method to try to extract a nontrivial factor 

of ``n``. Either a divisor (perhaps composite) or ``None`` is returned. 

The value of ``a`` is the base that is used in the test gcd(a**M - 1, n). 

The default is 2. If ``retries`` > 0 then if no factor is found after the 

first attempt, a new ``a`` will be generated randomly (using the ``seed``) 

and the process repeated. 

Note: the value of M is lcm(1..B) = reduce(ilcm, range(2, B + 1)). 

A search is made for factors next to even numbers having a power smoothness 

less than ``B``. Choosing a larger B increases the likelihood of finding a 

larger factor but takes longer. Whether a factor of n is found or not 

depends on ``a`` and the power smoothness of the even mumber just less than 

the factor p (hence the name p - 1). 

Although some discussion of what constitutes a good ``a`` some 

descriptions are hard to interpret. At the modular.math site referenced 

below it is stated that if gcd(a**M - 1, n) = N then a**M % q**r is 1 

for every prime power divisor of N. But consider the following: 

>>> from sympy.ntheory.factor_ import smoothness_p, pollard_pm1 

>>> n=257*1009 

>>> smoothness_p(n) 

(-1, [(257, (1, 2, 256)), (1009, (1, 7, 16))]) 

So we should (and can) find a root with B=16: 

>>> pollard_pm1(n, B=16, a=3) 

1009 

If we attempt to increase B to 256 we find that it doesn't work: 

>>> pollard_pm1(n, B=256) 

>>> 

But if the value of ``a`` is changed we find that only multiples of 

257 work, e.g.: 

>>> pollard_pm1(n, B=256, a=257) 

1009 

Checking different ``a`` values shows that all the ones that didn't 

work had a gcd value not equal to ``n`` but equal to one of the 

factors: 

>>> from sympy.core.numbers import ilcm, igcd 

>>> from sympy import factorint, Pow 

>>> M = 1 

>>> for i in range(2, 256): 

... M = ilcm(M, i) 

... 

>>> set([igcd(pow(a, M, n) - 1, n) for a in range(2, 256) if 

... igcd(pow(a, M, n) - 1, n) != n]) 

set([1009]) 

But does aM % d for every divisor of n give 1? 

>>> aM = pow(255, M, n) 

>>> [(d, aM%Pow(*d.args)) for d in factorint(n, visual=True).args] 

[(257**1, 1), (1009**1, 1)] 

No, only one of them. So perhaps the principle is that a root will 

be found for a given value of B provided that: 

1) the power smoothness of the p - 1 value next to the root 

does not exceed B 

2) a**M % p != 1 for any of the divisors of n. 

By trying more than one ``a`` it is possible that one of them 

will yield a factor. 

Examples 

======== 

With the default smoothness bound, this number can't be cracked: 

>>> from sympy.ntheory import pollard_pm1, primefactors 

>>> pollard_pm1(21477639576571) 

Increasing the smoothness bound helps: 

>>> pollard_pm1(21477639576571, B=2000) 

4410317 

Looking at the smoothness of the factors of this number we find: 

>>> from sympy.utilities import flatten 

>>> from sympy.ntheory.factor_ import smoothness_p, factorint 

>>> print(smoothness_p(21477639576571, visual=1)) 

p**i=4410317**1 has p-1 B=1787, B-pow=1787 

p**i=4869863**1 has p-1 B=2434931, B-pow=2434931 

The B and B-pow are the same for the p - 1 factorizations of the divisors 

because those factorizations had a very large prime factor: 

>>> factorint(4410317 - 1) 

{2: 2, 617: 1, 1787: 1} 

>>> factorint(4869863-1) 

{2: 1, 2434931: 1} 

Note that until B reaches the B-pow value of 1787, the number is not cracked; 

>>> pollard_pm1(21477639576571, B=1786) 

>>> pollard_pm1(21477639576571, B=1787) 

4410317 

The B value has to do with the factors of the number next to the divisor, 

not the divisors themselves. A worst case scenario is that the number next 

to the factor p has a large prime divisisor or is a perfect power. If these 

conditions apply then the power-smoothness will be about p/2 or p. The more 

realistic is that there will be a large prime factor next to p requiring 

a B value on the order of p/2. Although primes may have been searched for 

up to this level, the p/2 is a factor of p - 1, something that we don't 

know. The modular.math reference below states that 15% of numbers in the 

range of 10**15 to 15**15 + 10**4 are 10**6 power smooth so a B of 10**6 

will fail 85% of the time in that range. From 10**8 to 10**8 + 10**3 the 

percentages are nearly reversed...but in that range the simple trial 

division is quite fast. 

References 

========== 

- Richard Crandall & Carl Pomerance (2005), "Prime Numbers: 

A Computational Perspective", Springer, 2nd edition, 236-238 

- http://modular.math.washington.edu/edu/2007/spring/ent/ent-html/node81.html 

- http://www.cs.toronto.edu/~yuvalf/Factorization.pdf 

""" 

n = int(n) 

if n < 4 or B < 3: 

raise ValueError('pollard_pm1 should receive n > 3 and B > 2') 

prng = random.Random(seed + B) 

# computing a**lcm(1,2,3,..B) % n for B > 2 

# it looks weird, but it's right: primes run [2, B] 

# and the answer's not right until the loop is done. 

for i in range(retries + 1): 

aM = a 

for p in sieve.primerange(2, B + 1): 

e = int(math.log(B, p)) 

aM = pow(aM, pow(p, e), n) 

g = igcd(aM - 1, n) 

if 1 < g < n: 

return int(g) 

# get a new a: 

# since the exponent, lcm(1..B), is even, if we allow 'a' to be 'n-1' 

# then (n - 1)**even % n will be 1 which will give a g of 0 and 1 will 

# give a zero, too, so we set the range as [2, n-2]. Some references 

# say 'a' should be coprime to n, but either will detect factors. 

a = prng.randint(2, n - 2) 

def _trial(factors, n, candidates, verbose=False): 

""" 

Helper function for integer factorization. Trial factors ``n` 

against all integers given in the sequence ``candidates`` 

and updates the dict ``factors`` in-place. Returns the reduced 

value of ``n`` and a flag indicating whether any factors were found. 

""" 

if verbose: 

factors0 = list(factors.keys()) 

nfactors = len(factors) 

for d in candidates: 

if n % d == 0: 

m = multiplicity(d, n) 

n //= d**m 

factors[d] = m 

if verbose: 

for k in sorted(set(factors).difference(set(factors0))): 

print(factor_msg % (k, factors[k])) 

return int(n), len(factors) != nfactors 

def _check_termination(factors, n, limitp1, use_trial, use_rho, use_pm1, 

verbose): 

""" 

Helper function for integer factorization. Checks if ``n`` 

is a prime or a perfect power, and in those cases updates 

the factorization and raises ``StopIteration``. 

""" 

if verbose: 

print('Check for termination') 

# since we've already been factoring there is no need to do 

# simultaneous factoring with the power check 

p = perfect_power(n, factor=False) 

if p is not False: 

base, exp = p 

if limitp1: 

limit = limitp1 - 1 

else: 

limit = limitp1 

facs = factorint(base, limit, use_trial, use_rho, use_pm1, 

verbose=False) 

for b, e in facs.items(): 

if verbose: 

print(factor_msg % (b, e)) 

factors[b] = exp*e 

raise StopIteration 

if isprime(n): 

factors[int(n)] = 1 

raise StopIteration 

if n == 1: 

raise StopIteration 

trial_int_msg = "Trial division with ints [%i ... %i] and fail_max=%i" 

trial_msg = "Trial division with primes [%i ... %i]" 

rho_msg = "Pollard's rho with retries %i, max_steps %i and seed %i" 

pm1_msg = "Pollard's p-1 with smoothness bound %i and seed %i" 

factor_msg = '\t%i ** %i' 

fermat_msg = 'Close factors satisying Fermat condition found.' 

complete_msg = 'Factorization is complete.' 

def _factorint_small(factors, n, limit, fail_max): 

""" 

Return the value of n and either a 0 (indicating that factorization up 

to the limit was complete) or else the next near-prime that would have 

been tested. 

Factoring stops if there are fail_max unsuccessful tests in a row. 

If factors of n were found they will be in the factors dictionary as 

{factor: multiplicity} and the returned value of n will have had those 

factors removed. The factors dictionary is modified in-place. 

""" 

def done(n, d): 

"""return n, d if the sqrt(n) wasn't reached yet, else 

n, 0 indicating that factoring is done. 

""" 

if d*d <= n: 

return n, d 

return n, 0 

d = 2 

m = trailing(n) 

if m: 

factors[d] = m 

n >>= m 

d = 3 

if limit < d: 

if n > 1: 

factors[n] = 1 

return done(n, d) 

# reduce 

m = 0 

while n % d == 0: 

n //= d 

m += 1 

if m == 20: 

mm = multiplicity(d, n) 

m += mm 

n //= d**mm 

break 

if m: 

factors[d] = m 

# when d*d exceeds maxx or n we are done; if limit**2 is greater 

# than n then maxx is set to zero so the value of n will flag the finish 

if limit*limit > n: 

maxx = 0 

else: 

maxx = limit*limit 

dd = maxx or n 

d = 5 

fails = 0 

while fails < fail_max: 

if d*d > dd: 

break 

# d = 6*i - 1 

# reduce 

m = 0 

while n % d == 0: 

n //= d 

m += 1 

if m == 20: 

mm = multiplicity(d, n) 

m += mm 

n //= d**mm 

break 

if m: 

factors[d] = m 

dd = maxx or n 

fails = 0 

else: 

fails += 1 

d += 2 

if d*d > dd: 

break 

# d = 6*i - 1 

# reduce 

m = 0 

while n % d == 0: 

n //= d 

m += 1 

if m == 20: 

mm = multiplicity(d, n) 

m += mm 

n //= d**mm 

break 

if m: 

factors[d] = m