# Brent's variation of a factorization algorithm

Discussion in 'Python' started by n00m, Nov 27, 2009.

1. ### n00mGuest

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

n00m, Nov 27, 2009

2. ### Gabriel GenellinaGuest

En Fri, 27 Nov 2009 12:36:29 -0300, n00m <> escribió:

> 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

--
Gabriel Genellina

Gabriel Genellina, Dec 8, 2009

3. ### MRABGuest

Gabriel Genellina wrote:
> En Fri, 27 Nov 2009 12:36:29 -0300, n00m <> escribió:
>
>> 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
>

An iterative alternative is:

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

MRAB, Dec 8, 2009
4. ### Gabriel GenellinaGuest

En Tue, 08 Dec 2009 15:51:29 -0300, MRAB <>
escribió:
> Gabriel Genellina wrote:
>> En Fri, 27 Nov 2009 12:36:29 -0300, n00m <> escribió:
>>
>>> 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
>>

> 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...)

--
Gabriel Genellina

Gabriel Genellina, Dec 9, 2009
5. ### n00mGuest

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, Dec 9, 2009
6. ### n00mGuest

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 =)

n00m, Dec 9, 2009
7. ### Irmen de JongGuest

On 27-11-2009 16:36, n00m wrote:
> 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

Irmen de Jong, Dec 9, 2009
8. ### n00mGuest

On Dec 10, 1:11 am, Irmen de Jong <> wrote:
> 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

n00m, Dec 9, 2009
9. ### Irmen de JongGuest

On 12/10/09 12:52 AM, n00m wrote:
> On Dec 10, 1:11 am, Irmen de Jong<> wrote:
>> 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 =)

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

>> 27 * 37037037

> 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

Irmen de Jong, Dec 10, 2009
10. ### n00mGuest

On Dec 10, 2:59 am, Irmen de Jong <> wrote:
> On 12/10/09 12:52 AM, n00m wrote:
>
> > On Dec 10, 1:11 am, Irmen de Jong<>  wrote:
> >> 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 =)

>
> A 'semiprime' being a product of 2 prime numbers, I suppose.
>
> >> 27 * 37037037

> > 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

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...

n00m, Dec 10, 2009