Brent's variation of a factorization algorithm

[n00m]

Maybe someone'll make use of it:


def gcd(x, y):
if y == 0:
return x
return gcd(y, x % y)

def brent(n):
c = 11
y, r, q, m = 1, 1, 1, 137
while 1:
x = y
for i in range(1, r + 1):
y = (y * y + c) % n
k = 0
while 1:
ys = y
for i in range(1, min(m, r - k) + 1):
y = (y * y + c) % n
q = (q * abs(x - y)) % n
g = gcd(q, n)
k += m
if k >=r or g > 1:
break
r *= 2
if g > 1:
break
if g == n:
while 1:
ys = (ys * ys + c) % n
g = gcd(abs(x - ys), n)
if g > 1:
break
return g


while 1:
n = eval(raw_input())
g = brent(n)
print '==', g, '*', n / g



IDLE 1.2 ==== No Subprocess ====

1170999422783 * 10001
== 73 * 160426920921271
1170999422783 * 254885996264477
== 1170999422783 * 254885996264477
1170999422783 * 415841978209842084233328691123
== 1170999422783 * 415841978209842084233328691123
51539607551 * 80630964769
== 51539607551 * 80630964769
304250263527209 * 792606555396977
== 304250263527209 * 792606555396977

 

[Gabriel Genellina]

Maybe someone'll make use of it:


def gcd(x, y):
if y == 0:
return x
return gcd(y, x % y)

def brent(n): ...

A better place to publish this code would be the Python Cookbook:
http://code.activestate.com

 

[MRAB]

A better place to publish this code would be the Python Cookbook:
http://code.activestate.com

An iterative alternative is:

def gcd(x, y):
while y != 0:
x, y = y, x % y
return x

 

[Gabriel Genellina]

An iterative alternative is:

def gcd(x, y):
while y != 0:
x, y = y, x % y
return x

(note that the interesting part was the snipped code, the brent()
function...)

 

[n00m]

Being an absolute dummy in Theory of Number
for me ***c'est fantastique*** that brent() works =)

PS
1.
Values of magic parameters c = 11 and m = 137
almost don't matter. Usually they choose c = 2
(what about to run brent() in parallel with different
values of "c" waiting for "n" is cracked?)

2.
Before calling brent() "n" should be tested for its
primality. If it is a prime brent(n) may freeze for good.

3.

A better place to publish this code would be the Python Cookbook:

It requires a tedious registration etc.
Gabriel, don't you mind to publish the code there by yourself?
In the long run, it is an invention by Richard Brent (b.1946) =)
I just rewrote it to Python from a pseudo-code once available in
Wiki but which for some vague reason was later on removed from there.

 

[n00m]

PPS
The code was successfully tested e.g. here:
http://www.spoj.pl/ranks/FACT1/ (see my 2nd and 4th places).
They confused versions: the 2nd is in Python 2.5, not 2.6.2.
PPPS
Funnilly... almost only Python on the 1st page =)

 

[Irmen de Jong]

Maybe someone'll make use of it:


def gcd(x, y):
if y == 0:
return x
return gcd(y, x % y)

def brent(n):

[...]

[D:\Projects]python brentfactor.py
999999999
== 27 * 37037037

What gives? Isn't this thing supposed to factor numbers into the product
of two primes?

-irmen

 

[n00m]

999999999
== 27 * 37037037

What gives? Isn't this thing supposed to factor numbers into the product
of two primes?

-irmen

Only if you yield to it a SEMIprime =)

27 * 37037037

Now you can apply brent() to these numbers, and so on

 

[Irmen de Jong]

Only if you yield to it a SEMIprime =)

A 'semiprime' being a product of 2 prime numbers, I suppose.

Now you can apply brent() to these numbers, and so on

Right I more or less expected it to do that by itself (didn't look
at the algorithm)

But you wrote that it might run indefinately if you feed it a prime
number. There's no safeguard then against getting into an endless loop
if you keep applying brent() to the factors it produces? Because in the
end one or both factors will be prime...

-irmen

 

[n00m]

A 'semiprime' being a product of 2 prime numbers, I suppose.


Right   I more or less expected it to do that by itself (didn't look
at the algorithm)

But you wrote that it might run indefinately if you feed it a prime
number. There's no safeguard then against getting into an endless loop
if you keep applying brent() to the factors it produces? Because in the
end one or both factors will be prime...

-irmen

Just to use e.g. Rabin-Miller's before supplying a number to brent()

A 'semiprime' being a product of 2 prime numbers, I suppose.

Aha, exactly.
==============================
Also it's worthy to test a num is it a perfect square?
But all this is obvious trifles...