Reputation: 117801
Fast prime factorization module
I am looking for an implementation or clear algorithm for getting the prime factors of N in either python, pseudocode or anything else well-readable. There are a few requirements/constraints:
- N is between 1 and ~20 digits
- No pre-calculated lookup table, memoization is fine though
- Need not to be mathematically proven (e.g. could rely on the Goldbach conjecture if needed)
- Need not to be precise, is allowed to be probabilistic/deterministic if needed
I need a fast prime factorization algorithm, not only for itself, but for usage in many other algorithms like calculating the Euler phi(n).
I have tried other algorithms from Wikipedia and such but either I couldn't understand them (ECM) or I couldn't create a working implementation from the algorithm (Pollard-Brent).
I am really interested in the Pollard-Brent algorithm, so any more information/implementations on it would be really nice.
Thanks!
EDIT
After messing around a little I have created a pretty fast prime/factorization module. It combines an optimized trial division algorithm, the Pollard-Brent algorithm, a miller-rabin primality test and the fastest primesieve I found on the internet. gcd is a regular Euclid's GCD implementation (binary Euclid's GCD is much slower then the regular one).
Bounty
Oh joy, a bounty can be acquired! But how can I win it?
- Find an optimization or bug in my module.
- Provide alternative/better algorithms/implementations.
The answer which is the most complete/constructive gets the bounty.
And finally the module itself:
import random
def primesbelow(N):
# http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
#""" Input N>=6, Returns a list of primes, 2 <= p < N """
correction = N % 6 > 1
N = {0:N, 1:N-1, 2:N+4, 3:N+3, 4:N+2, 5:N+1}[N%6]
sieve = [True] * (N // 3)
sieve[0] = False
for i in range(int(N ** .5) // 3 + 1):
if sieve[i]:
k = (3 * i + 1) | 1
sieve[k*k // 3::2*k] = [False] * ((N//6 - (k*k)//6 - 1)//k + 1)
sieve[(k*k + 4*k - 2*k*(i%2)) // 3::2*k] = [False] * ((N // 6 - (k*k + 4*k - 2*k*(i%2))//6 - 1) // k + 1)
return [2, 3] + [(3 * i + 1) | 1 for i in range(1, N//3 - correction) if sieve[i]]
smallprimeset = set(primesbelow(100000))
_smallprimeset = 100000
def isprime(n, precision=7):
# http://en.wikipedia.org/wiki/Miller-Rabin_primality_test#Algorithm_and_running_time
if n < 1:
raise ValueError("Out of bounds, first argument must be > 0")
elif n <= 3:
return n >= 2
elif n % 2 == 0:
return False
elif n < _smallprimeset:
return n in smallprimeset
d = n - 1
s = 0
while d % 2 == 0:
d //= 2
s += 1
for repeat in range(precision):
a = random.randrange(2, n - 2)
x = pow(a, d, n)
if x == 1 or x == n - 1: continue
for r in range(s - 1):
x = pow(x, 2, n)
if x == 1: return False
if x == n - 1: break
else: return False
return True
# https://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/
def pollard_brent(n):
if n % 2 == 0: return 2
if n % 3 == 0: return 3
y, c, m = random.randint(1, n-1), random.randint(1, n-1), random.randint(1, n-1)
g, r, q = 1, 1, 1
while g == 1:
x = y
for i in range(r):
y = (pow(y, 2, n) + c) % n
k = 0
while k < r and g==1:
ys = y
for i in range(min(m, r-k)):
y = (pow(y, 2, n) + c) % n
q = q * abs(x-y) % n
g = gcd(q, n)
k += m
r *= 2
if g == n:
while True:
ys = (pow(ys, 2, n) + c) % n
g = gcd(abs(x - ys), n)
if g > 1:
break
return g
smallprimes = primesbelow(1000) # might seem low, but 1000*1000 = 1000000, so this will fully factor every composite < 1000000
def primefactors(n, sort=False):
factors = []
for checker in smallprimes:
while n % checker == 0:
factors.append(checker)
n //= checker
if checker > n: break
if n < 2: return factors
while n > 1:
if isprime(n):
factors.append(n)
break
factor = pollard_brent(n) # trial division did not fully factor, switch to pollard-brent
factors.extend(primefactors(factor)) # recurse to factor the not necessarily prime factor returned by pollard-brent
n //= factor
if sort: factors.sort()
return factors
def factorization(n):
factors = {}
for p1 in primefactors(n):
try:
factors[p1] += 1
except KeyError:
factors[p1] = 1
return factors
totients = {}
def totient(n):
if n == 0: return 1
try: return totients[n]
except KeyError: pass
tot = 1
for p, exp in factorization(n).items():
tot *= (p - 1) * p ** (exp - 1)
totients[n] = tot
return tot
def gcd(a, b):
if a == b: return a
while b > 0: a, b = b, a % b
return a
def lcm(a, b):
return abs((a // gcd(a, b)) * b)
Upvotes: 107
Views: 81454
Answers (9)
Reputation: 11
At most 1 of the factors is bigger than sqrt(n), so all methods that test numbers up to n are too slow. Test up to sqrt(n), and check the last factor separately.
def prime_factors(n: int) -> dict:
factors = {}
i = 2
while i * i <= n:
if not n % i:
factors[i] = factors.get(i, 0) + 1
n //= i
else:
i += 1
factors[n] = factors.get(n, 0) + 1
return factors
Upvotes: 0
Reputation: 16767
Very interesting question, thanks!
Today I decided to implement for you from scratch Elliptic Curve ECM factorization method in Python.
I tried to make clean and easy understandable code. You can just copy-paste whole function FactorECM() from my code into your code and it will work without any other tweaking or dependencies.
Due to some simplifications in the code algorithm is NOT highly optimized, e.g. it doesn't use multiprocessing module to use processes and all CPU cores. Basically my algorithm is single core for single factored number.
I used following sub-algorithms in my code: Trial Division Factorization, Fermat Probability Test, Sieve of Eratosthenes (prime numbers generator), Euclidean Algorithm (computing Greatest Common Divisor, GCD), Extended Euclidean Algorithm (GCD with Bezu coefficients), Modular Multiplicative Inverse, Right-to-Left Binary Exponentation (for elliptic point multiplication), Elliptic Curve Arithmetics (point addition and multiplication), Lenstra Elliptic Curve Factorization.
NOTE. You can speedup my code 2.5x times if you install gmpy2 module through python -m pip install gmpy2 command line. But don't forget to uncomment line #import gmpy2 in my code, if uncommented then 2.5x boost will be applied. I used gmpy2 for its considerably faster implementation of gmpy2.gcdext() functions, which implements Extended Euclidean Algorithm needed for Modular Multiplicative Inverse.
ECM algorithm in general, if well optimized (for example written in C++), is capable of factoring quite big numbers, even hardest 100-bit number (30 digits), consisting of two 50-bit primes, can be factored within several seconds.
My algorithm is written according to ECM in Wikipedia and is as following:
Check if number is smaller than 2^16, then factor it through Trial Division method. Return result.
Check if number is probably prime with high condifence, for that I use Fermat Test with 32 trials. If number is primes return it as only factor.
Generate curve parameters
A, X, Yrandomly and deriveBfrom curve equationY^2 = X^3 + AX + B (mod N). Check if curve is OK, value4 * A ** 3 - 27 * B ** 2should be non-zero.Generate small primes through Sieve of Eratosthenes, primes below our
Bound. Each prime raise to some small power, this raised prime would be calledK. I do raising into power while it is smaller than someBound2, which isSqrt(Bound)in my case.Compute elliptic point multiplication
P = k * P, where K taken from previous step and P is (X, Y). Compute according to Elliptic Curve Arithmetics Wiki.Point multiplication uses Modular Multiplicative Inverse, which computes
GCD(SomeValue, N)according to Extended Euclidean Algorithm Wiki. If this GCD is not 1, then it gives non-trivial factor of N. Collect this factor, remove it from number and re-run ECM algorithm (steps1.-6.above) for remaining number.If all primes till Bound were multiplied and gave no factor then re-run ECM factorization algorithm (steps
1.-6.above) again with another random curve and bigger Bound. In my code I take new bound by adding 512 to old bound.
See example of usage of my factoring function inside Test() function. After code below is located example console output for factoring hardest 72-bit random number (consisting of two 36-bit primes). If you want different size of number than modify bits = 72 in my code to your desired bit size of input random number.
def FactorECM(N0, *, verbose = False):
# https://en.wikipedia.org/wiki/Lenstra_elliptic-curve_factorization
import math, random, time; gmpy2 = None
#import gmpy2
def Factor_TrialDiv(x):
# https://en.wikipedia.org/wiki/Trial_division
fs = []
while (x & 1) == 0:
fs.append(2)
x >>= 1
for d in range(3, x + 1, 2):
if d * d > x:
break
while x % d == 0:
fs.append(d)
x //= d
if x > 1:
fs.append(x)
return sorted(fs)
def IsProbablyPrime_Fermat(n, trials = 32):
# https://en.wikipedia.org/wiki/Fermat_primality_test
for i in range(trials):
if pow(random.randint(2, n - 2), n - 1, n) != 1:
return False
return True
def GenPrimes_SieveOfEratosthenes(end):
# https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
composite = [False] * end
for p in range(2, int(math.sqrt(end) + 3)):
if composite[p]:
continue
for i in range(p * p, end, p):
composite[i] = True
for p in range(2, end):
if not composite[p]:
yield p
def GCD(a, b):
# https://en.wikipedia.org/wiki/Euclidean_algorithm
while b != 0:
a, b = b, a % b
return a
def EGCD(a, b):
# https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
if gmpy2 is None:
ro, r, so, s = a, b, 1, 0
while r != 0:
ro, (q, r) = r, divmod(ro, r)
so, s = s, so - q * s
return ro, so, (ro - so * a) // b
else:
return tuple(map(int, gmpy2.gcdext(a, b)))
def ModularInverse(a, n):
# https://en.wikipedia.org/wiki/Modular_multiplicative_inverse
g, s, t = EGCD(a, n)
if g != 1:
raise ValueError(a)
return s % n
def EllipticCurveAdd(N, A, B, X0, Y0, X1, Y1):
# https://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication
if X0 == X1 and Y0 == Y1:
# Double
l = ((3 * X0 ** 2 + A) * ModularInverse(2 * Y0, N)) % N
x = (l ** 2 - 2 * X0) % N
y = (l * (X0 - x) - Y0) % N
else:
# Add
l = ((Y1 - Y0) * ModularInverse(X1 - X0, N)) % N
x = (l ** 2 - X0 - X1) % N
y = (l * (X0 - x) - Y0) % N
return x, y
def EllipticCurveMul(N, A, B, X, Y, k):
# https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method
assert k >= 2, k
k -= 1
BX, BY = X, Y
while k != 0:
if k & 1:
X, Y = EllipticCurveAdd(N, A, B, X, Y, BX, BY)
BX, BY = EllipticCurveAdd(N, A, B, BX, BY, BX, BY)
k >>= 1
return X, Y
bound_start = 1 << 9
def Main(N, *, bound = bound_start, icurve = 0):
def NextFactorECM(x):
return Main(x, bound = bound + bound_start, icurve = icurve + 1)
def PrimePow(p, *, bound2 = int(math.sqrt(bound) + 1.01)):
mp = p
while True:
mp *= p
if mp >= bound2:
return mp // p
if N < (1 << 16):
fs = Factor_TrialDiv(N)
if verbose and len(fs) >= 2:
print('Factors from TrialDiv:', fs, flush = True)
return fs
if IsProbablyPrime_Fermat(N):
return [N]
if verbose:
print(f'Curve {icurve:>4}, bound 2^{math.log2(bound):>7.3f}', flush = True)
while True:
X, Y, A = [random.randrange(N) for i in range(3)]
B = (Y ** 2 - X ** 3 - A * X) % N
if 4 * A ** 3 - 27 * B ** 2 == 0:
continue
break
for p in GenPrimes_SieveOfEratosthenes(bound):
k = PrimePow(p)
try:
X, Y = EllipticCurveMul(N, A, B, X, Y, k)
except ValueError as ex:
g = GCD(ex.args[0], N)
assert g > 1, g
if g != N:
if verbose:
print('Factor from ECM:', g, flush = True)
return sorted(NextFactorECM(g) + NextFactorECM(N // g))
else:
return NextFactorECM(N)
return NextFactorECM(N)
if verbose:
print(f'ECM factoring N: {N0} (2^{math.log2(max(1, N0)):.2f})', flush = True)
tb = time.time()
fs = Main(N0)
tb = time.time() - tb
print(f'Factored N {N0}: {fs}')
print(f'Time {tb:.3f} sec.', flush = True)
return fs
else:
return Main(N0)
def Test():
import random
def IsProbablyPrime_Fermat(n, trials = 32):
for i in range(trials):
if pow(random.randint(2, n - 2), n - 1, n) != 1:
return False
return True
def RandPrime(bits):
while True:
p = random.randrange(1 << (bits - 1), 1 << bits) | 1
if IsProbablyPrime_Fermat(p):
return p
bits = 72
N = RandPrime((bits + 1) // 2) * RandPrime(bits // 2) * random.randrange(1 << 16)
FactorECM(N, verbose = True)
if __name__ == '__main__':
Test()
Example console output:
ECM factoring N: 25005272280974861996134424 (2^84.37)
Curve 0, bound 2^ 9.000
Factor from ECM: 2
Curve 1, bound 2^ 10.000
Factor from ECM: 4
Factors from TrialDiv: [2, 2]
Curve 2, bound 2^ 10.585
Factor from ECM: 1117
Curve 3, bound 2^ 11.000
Curve 4, bound 2^ 11.322
Curve 5, bound 2^ 11.585
Curve 6, bound 2^ 11.807
Factor from ECM: 54629318837
Factored N 25005272280974861996134424: [2, 2, 2, 1117, 51222720707, 54629318837]
Time 0.281 sec.
Upvotes: 1
Reputation: 137752
If you don't want to reinvent the wheel, use the library sympy
pip install sympy
Use the function sympy.ntheory.factorint
Given a positive integer
n,factorint(n)returns a dict containing the prime factors ofnas keys and their respective multiplicities as values. For example:
Example:
>>> from sympy.ntheory import factorint
>>> factorint(10**20+1)
{73: 1, 5964848081: 1, 1676321: 1, 137: 1}
You can factor some very large numbers:
>>> factorint(10**100+1)
{401: 1, 5964848081: 1, 1676321: 1, 1601: 1, 1201: 1, 137: 1, 73: 1, 129694419029057750551385771184564274499075700947656757821537291527196801: 1}
Upvotes: 113
Reputation: 763
You could factorize up to a limit then use brent to get higher factors
from fractions import gcd
from random import randint
def brent(N):
if N%2==0: return 2
y,c,m = randint(1, N-1),randint(1, N-1),randint(1, N-1)
g,r,q = 1,1,1
while g==1:
x = y
for i in range(r):
y = ((y*y)%N+c)%N
k = 0
while (k<r and g==1):
ys = y
for i in range(min(m,r-k)):
y = ((y*y)%N+c)%N
q = q*(abs(x-y))%N
g = gcd(q,N)
k = k + m
r = r*2
if g==N:
while True:
ys = ((ys*ys)%N+c)%N
g = gcd(abs(x-ys),N)
if g>1: break
return g
def factorize(n1):
if n1==0: return []
if n1==1: return [1]
n=n1
b=[]
p=0
mx=1000000
while n % 2 ==0 : b.append(2);n//=2
while n % 3 ==0 : b.append(3);n//=3
i=5
inc=2
while i <=mx:
while n % i ==0 : b.append(i); n//=i
i+=inc
inc=6-inc
while n>mx:
p1=n
while p1!=p:
p=p1
p1=brent(p)
b.append(p1);n//=p1
if n!=1:b.append(n)
return sorted(b)
from functools import reduce
#n= 2**1427 * 31 #
n= 67898771 * 492574361 * 10000223 *305175781* 722222227*880949 *908909
li=factorize(n)
print (li)
print (n - reduce(lambda x,y :x*y ,li))
Upvotes: 2
Reputation: 27340
I just ran into a bug in this code when factoring the number 2**1427 * 31.
File "buckets.py", line 48, in prettyprime
factors = primefactors.primefactors(n, sort=True)
File "/private/tmp/primefactors.py", line 83, in primefactors
limit = int(n ** .5) + 1
OverflowError: long int too large to convert to float
This code snippet:
limit = int(n ** .5) + 1
for checker in smallprimes:
if checker > limit: break
while n % checker == 0:
factors.append(checker)
n //= checker
limit = int(n ** .5) + 1
if checker > limit: break
should be rewritten into
for checker in smallprimes:
while n % checker == 0:
factors.append(checker)
n //= checker
if checker > n: break
which will likely perform faster on realistic inputs anyway. Square root is slow — basically the equivalent of many multiplications —, smallprimes only has a few dozen members, and this way we remove the computation of n ** .5 from the tight inner loop, which is certainly helpful when factoring numbers like 2**1427. There's simply no reason to compute sqrt(2**1427), sqrt(2**1426), sqrt(2**1425), etc. etc., when all we care about is "does the [square of the] checker exceed n".
The rewritten code doesn't throw exceptions when presented with big numbers, and is roughly twice as fast according to timeit (on sample inputs 2 and 2**718 * 31).
Also notice that isprime(2) returns the wrong result, but this is okay as long as we don't rely on it. IMHO you should rewrite the intro of that function as
if n <= 3:
return n >= 2
...
Upvotes: 0
Reputation: 7790
There's a python library with a collection of primality tests (including incorrect ones for what not to do). It's called pyprimes. Figured it's worth mentioning for posterity's purpose. I don't think it includes the algorithms you mentioned.
Upvotes: 4
Reputation: 4165
There is no need to calculate smallprimes using primesbelow, use smallprimeset for that.
smallprimes = (2,) + tuple(n for n in xrange(3,1000,2) if n in smallprimeset)
Divide your primefactors into two functions for handling smallprimes and other for pollard_brent, this can save a couple of iterations as all the powers of smallprimes will be divided from n.
def primefactors(n, sort=False):
factors = []
limit = int(n ** .5) + 1
for checker in smallprimes:
print smallprimes[-1]
if checker > limit: break
while n % checker == 0:
factors.append(checker)
n //= checker
if n < 2: return factors
else :
factors.extend(bigfactors(n,sort))
return factors
def bigfactors(n, sort = False):
factors = []
while n > 1:
if isprime(n):
factors.append(n)
break
factor = pollard_brent(n)
factors.extend(bigfactors(factor,sort)) # recurse to factor the not necessarily prime factor returned by pollard-brent
n //= factor
if sort: factors.sort()
return factors
By considering verified results of Pomerance, Selfridge and Wagstaff and Jaeschke, you can reduce the repetitions in isprime which uses Miller-Rabin primality test. From Wiki.
- if n < 1,373,653, it is enough to test a = 2 and 3;
- if n < 9,080,191, it is enough to test a = 31 and 73;
- if n < 4,759,123,141, it is enough to test a = 2, 7, and 61;
- if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11;
- if n < 3,474,749,660,383, it is enough to test a = 2, 3, 5, 7, 11, and 13;
- if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17.
Edit 1: Corrected return call of if-else to append bigfactors to factors in primefactors.
Upvotes: 19
Reputation: 10663
Even on the current one, there are several spots to be noticed.
- Don't do
checker*checkerevery loop, uses=ceil(sqrt(num))andchecher < s - checher should plus 2 each time, ignore all even numbers except 2
- Use
divmodinstead of%and//
Upvotes: 7
Reputation: 3043
You should probably do some prime detection which you could look here, Fast algorithm for finding prime numbers?
You should read that entire blog though, there is a few algorithms that he lists for testing primality.
Upvotes: 5
Related Questions
- Prime Factorisation in Python
- Fast Algorithm to Factorize All Numbers Up to a Given Number
- Prime factorization python
- Integer factorization using Prime list (trying to find fastest algorithm)
- Algorithm recommendation for calculating prime factors (Python Beginner)
- Prime Factorization in Python
- Prime Factorization algorithm
- python prime factorization performance
- Optimizing a Prime Number Factorization algorithm
- Speeding up my Fermat Factorization function (Python)