Pollard's rho algorithm
Pollard's rho algorithm is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
In the area of number decomposition algorithms, Pollard's rho algorithm is a fast, yet fairly simple method to find factors of a number, using a pseudorandom number generator modulo the number.
Except for trivial factors like 2, 3, 5 etc., it is much faster than trial division, finding a factor quite often, but may terminate without a factor. Thus it cannot prove that a number is prime. If it is known that a number is not prime, the algorithm can be re-iterated with changed parameters to finally find a factor (not necessarily a prime factor) with quite high probability.
- Task
Write a function to demonstrate the algorithm (see pseudocode in Wikipedia). The intent is to show the principle as clear as possible in the respective language, but avoid an endless loop. Enhanced versions may be presented as extra examples.
- Warning
The variables 'x' and 'y' may be equal during an iteration, calling 'gcd()' with first argument 0, while the GCD is only defined for postitive numbers, so a library gdc() function might fail. If in this case, the other argument is returned, the algorithm reports failure, which is acceptable. The gcd() functions supplied for AWK and C have this property.
Example numbers, if supported:
- 4294967213 = 75167 * 57139 (32 bits)
- 9759463979 = 98779 * 98801 (34 bits)
- 34225158206557151 = 130051349 * 263166499 (55 bits)
- 763218146048580636353 = 25998278833 * 29356487441 (70 bits)
- Related tasks
- prime decomposition
- AKS test for primes
- factors of an integer
- factors of a Mersenne number
- trial factoring of a Mersenne number
- partition an integer X into N primes
- sequence of primes by Trial Division
- References
Includes additional test cases as in some of the other samples, also includes 4 which (possibly surprisingly) can't be factorised by the algorithm.
This slightly modifies the algorithm by checking for even numbers first and returning 2 in that case (as in some of the other samples). The number itself is returned instead of 0 if a factor couldn't be found, e.g. if the number is a prime.
If a factor can't found on the first attempt, the number added to the square in the g function is incremented and the algorithm tried again, repeating up to 20 times if necessary. However after the second attempt and no factor found, if the number is a square or a prime, the search stops (avoiding an up-front primality test which could be expensive for larger numbers).
Following the task test cases, the algorithm is applied to n up to 100 000 and where a factor couldn't be found the number is checked to be prime: 1 is the only number that isn't prime that can't be (partially) factorised.
In ALGOL 68 Genie versions 2 and 3, LONG INT is large enough for the task test cases but squaring a large LONG INT will need to have a LONG LONG INT result, hence the use of the LENG and SHORTEN operators to ensure the result will not overflow. Alternatively, LONG LONG INT could be used for everything and a suitable precision specified in a pragmatic comment.
This should also work with other Algol 68 implementations that have suitable sizes of INT.
Note, the RC ALGOL 68 primes library (primes.incl.a68) is on a separate page on Rosetta Code - see the above link,
BEGIN # Pollard's rho algorithm #
PR read "primes.incl.a68" PR
# select a suitable size INT mode for the number of digits required #
MODE INTEGER = LONG INT;
PROC gcd = ( INTEGER x, y )INTEGER: # iterative gcd #
BEGIN
INTEGER a := x, b := y;
WHILE b /= 0 DO
INTEGER next a = b;
b := a MOD b;
a := next a
OD;
ABS a
END # gcd # ;
# returns a prime factor (if possible) of n via Pollard's rho algorithm #
PROC pollard rho = ( INTEGER n )INTEGER:
IF n < 0
THEN pollard rho( - n )
ELIF n < 2
THEN n
ELIF NOT ODD n
THEN 2
ELSE
INTEGER d := 1;
INT number of attempts := 1;
BOOL square or prime := FALSE;
PROC g = ( INTEGER gx )INTEGER: SHORTEN ( ( ( LENG gx * LENG gx ) + number of attempts ) MOD n );
WHILE
INTEGER x := 2, y := 2;
d := 1;
WHILE d = 1 DO
x := g( x );
y := g( g( y ) );
d := gcd( ABS( x - y ), n )
OD;
d = n AND number of attempts < 20 AND number of attempts < n
AND NOT square or prime
DO
IF ( number of attempts +:= 1 ) = 3 THEN
# second try failed - before we try again, check for a #
# square or prime #
LONG REAL root n = long sqrt( n );
IF ENTIER root n = root n THEN
# n is a square #
square or prime := TRUE;
d := ENTIER root n
ELIF is probably prime( n ) THEN
# n is a prime #
square or prime := TRUE
FI
FI
OD;
d
FI # pollard rho # ;
BEGIN # test cases #
[]INTEGER test = ( 4, 13, 125, 4294967213, 9759463979, 34225158206557151, 763218146048580636353 );
FOR test pos FROM LWB test TO UPB test DO
INTEGER n = test[ test pos ];
INTEGER d1 = pollard rho( n );
print( ( whole( n, -24 ), " = ", whole( d1, -12 ), " * ", whole( n OVER d1, 0 ), newline ) )
OD;
INT max prime check = 100 000;
print( ( "Checking for primes up to ", whole( max prime check, 0 ), newline ) );
FOR i TO max prime check DO
IF pollard rho( i ) = i
THEN # pollard's rho couldn't find a factor or i is prime #
IF NOT is probably prime( i )
THEN print( ( " ?", whole( i, 0 ) ) )
ELSE IF i < 100 THEN print( ( " ", whole( i, 0 ) ) ) FI
FI
FI
OD;
print( ( newline ) )
END
END- Output:
4 = 2 * 2
13 = 13 * 1
125 = 5 * 25
4294967213 = 57139 * 75167
9759463979 = 98779 * 98801
34225158206557151 = 130051349 * 263166499
763218146048580636353 = 25998278833 * 29356487441
Checking for primes up to 100000
?1 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
Note that AWK uses internally floating point numbers, thus the programme fails for larger numbers.
However, Gawk currently supports arbitrary precision arithmetic, enabled by specifying -M
or --bignum on the command line.
Obtaining the numbers from terminal input is normal for AWK.
# the modulo generator function
function g(x, m) {
return (x*x + 1) % m
}
# the rho iteration
function rho(n, x, y, d, s, lim) {
x = 1; y = x; d = 1; s = 0
lim =sqrt(n)
while (d == 1 && s < lim) {
s += 1
x = g(x,n)
y = g( g(y,n), n)
d = gcd( abs(x-y), n)
if (d == n)
return 0 # failure
}
return d # factor found
}
{
i = $0
p = rho(i)
if( p == 0)
print i ": no factor found"
else
printf "%d = %d * %d\n", i, p, i/p
}
# common helper functions
function gcd(x, y) {
if (x < y) return gcd(y, x)
while (y > 0) {
t = x % y; x = y; y = t
}
return x
}
function abs(x) {
if (x < 0)
return -x
return x
}
- Output:
Using Gawk with -M:
234571407 = 3 * 78190469 4294967213 = 57139 * 75167 9759463979 = 98779 * 98801 34225158206557151 = 130051349 * 263166499 763218146048580636353 = 25998278833 * 29356487441
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
typedef long long int int64;
int64 gcd(int64 x, int64 y); // helper function at end
// modulo generator function
int64 g(int64 x, int64 m) {
return (x*x + 1) % m;
}
// the rho iteration
int64 rho(int64 n) {
int64 x = 1;
int64 y = x;
int64 d = 1;
int64 s = 0;
int64 lim = sqrt(n);
while (d == 1 && s < lim) {
s += 1;
x = g(x,n);
y = g( g(y,n), n);
d = gcd( abs(x-y), n);
if (d == n)
return 0; // failure
}
return d; //factor found
}
int64 main(int64 argc, char** argv) {
int64 i = 9759463979;
int64 p = rho(i);
if( p == 0)
printf("%lld: no factor found", i);
else
printf("%lld = %lld * %lld\n", i, p, i/p);
}
// well-known helper function
int64 gcd(int64 x, int64 y) {
if (x < y)
return gcd(y, x);
int64 t;
while (y > 0) {
t = x % y;
x = y;
y = t;
}
return x;
}
#include <bit>
#include <cmath>
#include <iostream>
#include <random>
#include <vector>
const uint32_t WORD_SIZE = 64;
std::random_device random;
std::mt19937 generator(random());
uint64_t gcd(uint64_t a, uint64_t b) {
if ( a < b ) {
return gcd(b, a);
}
while ( b > 0 ) {
const uint64_t temp = a % b;
a = b;
b = temp;
}
return a;
}
uint64_t pollards_rho(const uint64_t& number) {
if ( number % 2 == 0 ) {
return 2;
}
const uint32_t bit_length = WORD_SIZE - std::countl_zero(number);
std::uniform_int_distribution<uint64_t> distribution(0, std::pow(2, bit_length - 1));
const uint64_t constant = distribution(generator);
uint64_t x = distribution(generator);
uint64_t y = x;
uint64_t divisor = 1;
do {
x = ( x * x + constant ) % number;
y = ( y * y + constant ) % number;
y = ( y * y + constant ) % number;
divisor = gcd(x - y, number);
} while ( divisor == 1 );
return divisor;
}
int main() {
const std::vector<uint64_t> tests = { 4294967213, 9759463979, 34225158206557151, 13 };
for ( const uint64_t& test : tests ) {
const uint64_t divisor_one = pollards_rho(test);
const uint64_t divisor_two = test / divisor_one;
std::cout << test << " = " << divisor_one << " * " << divisor_two << std::endl;
}
}
- Output:
4294967213 = 75167 * 57139 9759463979 = 98779 * 98801 34225158206557151 = 130051349 * 263166499 13 = 13 * 1
import 'dart:math';
int gcd(int n, int d) {
if (n < d) return gcd(d, n);
while (d > 0) {
final temp = n % d;
n = d;
d = temp;
}
return n;
'return d == 0 ? n : gcd(d, n % d); MAS LENTO
}
int bitLength(int n) {
var count = 0;
while (n > 0) {
n >>= 1;
count++;
}
return count;
}
int pollardsRho(int number) {
if (number % 2 == 0) return 2;
final random = Random();
final bitLen = number.toRadixString(2).length;
final constant = random.nextInt(bitLen - 1);
var x = random.nextInt(bitLen - 1);
var y = x;
var divisor = 1;
do {
x = (x * x + constant) % number;
y = (y * y + constant) % number;
y = (y * y + constant) % number;
divisor = gcd((x - y).abs(), number);
} while (divisor == 1);
return divisor;
}
void main() {
final tests = [4294967213, 9759463979, 34225158206557151, 13];
final stopwatch = Stopwatch()..start();
for (final test in tests) {
final divisor1 = pollardsRho(test);
final divisor2 = test ~/ divisor1;
final bits = bitLength(test);
print('$test = $divisor1 * $divisor2 ($bits bits)');
}
print('\nTook ${stopwatch.elapsedMilliseconds / 1000} seconds');
}
Const As Integer WORD_SIZE = 64
Function GCD(Byval n As Longint, Byval d As Longint) As Longint
If n < d Then Return GCD(d, n)
Return Iif(d = 0, n, GCD(d, n Mod d))
End Function
Function countlZero(Byval n As Longint) As Integer
Dim As Integer count = 0
For i As Integer = WORD_SIZE-1 To 0 Step -1
If (n And (1LL Shl i)) = 0 Then
count += 1
Else
Exit For
End If
Next
Return count
End Function
Function randomRange(Byval min As Longint, Byval max As Longint) As Longint
Return min + (Rnd * (max - min + 1))
End Function
Function pollardsRho(Byval number As Longint) As Longint
If (number Mod 2) = 0 Then Return 2
Dim As Integer bit_length = WORD_SIZE - countlZero(number)
Dim As Longint constant = randomRange(0, bit_length - 1)
Dim As Longint x = randomRange(0, bit_length - 1)
Dim As Longint y = x
Dim As Longint divisor = 1
Do
x = ((x * x + constant) Mod number)
y = ((y * y + constant) Mod number)
y = ((y * y + constant) Mod number)
divisor = Iif(x > y, GCD(x - y, number), GCD(y - x, number))
Loop While divisor = 1
Return divisor
End Function
Function bitLength(Byval n As Longint) As Byte
Dim As Byte count = 0
While n > 0
n Shr= 1
count += 1
Wend
Return count
End Function
Randomize Timer
' Test values
Dim tests(3) As Longint = {4294967213LL, 9759463979LL, 34225158206557151LL, 13LL}
For i As Integer = 0 To Ubound(tests)
Dim As Longint test = tests(i)
Dim As Longint divisor1 = pollardsRho(test)
Dim As Longint divisor2 = test \ divisor1
Dim As Byte bits = bitLength(test)
Print Using "& = & * & (& bits)"; test; divisor1; divisor2; bits
Next
Sleep
- Output:
4294967213 = 75167 * 57139 (32 bits) 9759463979 = 98779 * 98801 (34 bits) 34225158206557151 = 130051349 * 263166499 (55 bits) 13 = 13 * 1 (4 bits)
Uses GMP/MPFR framework: https://www.brilorsoftware.com/fb/utilities/gmp-mpfr.zip
FB 7.0.36, macOS 14.7.2 on M2 Max MacStudio
include "gmp.incl"
CFMutableStringRef mutStr
mutStr = fn MutableStringNew
void local fn Compute_GCD( result as mpz_t, a as mpz_t, b as mpz_t )
mpz_gcd( result, a, b )
end fn
void local fn Pollards_Rho( factor as mpz_t, n as mpz_t )
if ( fn mpz_even_p( n ) )
mpz_set_ui( factor, 2 )
return
end if
gmp_randstate_t state
gmp_randinit_default( state )
gmp_randseed_ui( state, 254893255633 )
mpz_t x, y, c, d, abs_diff
mpz_inits( x, y, c, d, abs_diff, NULL )
mpz_urandomm( x, state, n )
mpz_set( y, x )
mpz_urandomm( c, state, n )
mpz_set_ui( d, 1 )
while ( fn mpz_cmp_ui( d, 1 ) == 0 )
mpz_mul( x, x, x )
mpz_add( x, x, c )
mpz_mod( x, x, n )
for int i = 0 to 1
mpz_mul( y, y, y )
mpz_add( y, y, c )
mpz_mod( y, y, n )
next
mpz_sub( abs_diff, x, y )
mpz_abs( abs_diff, abs_diff )
fn Compute_GCD( d, abs_diff, n )
wend
if ( fn mpz_cmp( d, n ) == 0 )
mpz_set_ui( factor, 0 )
else
mpz_set( factor, d )
end if
mpz_clears( x, y, c, d, abs_diff, NULL )
gmp_randclear( state )
end fn
void local fn Factor_Recursive( n as mpz_t )
if ( fn mpz_cmp_ui( n, 1 ) == 0 ) then return
if ( fn mpz_probab_prime_p( n, 25 ) )
MutableStringAppendFormat( mutStr, @"%@ * \b", fn mpz_cf2( n ) )
exit fn
end if
mpz_t factor
mpz_init( factor )
while ( 1 )
fn Pollards_Rho( factor, n )
if ( fn mpz_cmp_ui( factor, 0 ) != 0 ) then break
wend
mpz_t cofactor
mpz_init( cofactor )
mpz_div( cofactor, n, factor )
fn Factor_Recursive( factor )
fn Factor_Recursive( cofactor )
mpz_clear( factor )
mpz_clear( cofactor )
end fn
local fn GetFactors( numberToFactor as CFStringRef )
mpz_t number
mpz_init( number )
fn mpz_init_set_cf( number, numberToFactor, 10 )
fn Factor_Recursive( number )
printf @"%@ = %@", numberToFactor, fn StringSubstringToIndex( mutStr, len(mutStr) -3 )
MutableStringSetString( mutStr, @"" )
mpz_clear( number )
end fn
CFTimeInterval t : t = fn CACurrentMediaTime
fn GetFactors( @"13" )
fn GetFactors( @"4294967213" )
fn GetFactors( @"9759463979" )
fn GetFactors( @"34225158206557151" )
fn GetFactors( @"5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497" )
printf @"\nCompute time: %.3f seconds",(fn CACurrentMediaTime-t)
HandleEvents- Output:
13 = 13 4294967213 = 75167 * 57139 9759463979 = 98779 * 98801 34225158206557151 = 130051349 * 263166499 5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497 = 165901 * 10424087 * 18830281 * 95174413 * 91931221 * 59663291 * 56402249 * 53204737 * 790130065009 * 444161842339 * 305293727939 Compute time: 0.279 seconds
package main
import (
"crypto/rand"
"fmt"
"math/big"
"time"
)
func main() {
// Create test numbers
tests := []*big.Int{}
n1, _ := new(big.Int).SetString("4294967213", 10)
n2, _ := new(big.Int).SetString("9759463979", 10)
n3, _ := new(big.Int).SetString("34225158206557151", 10)
n4, _ := new(big.Int).SetString("763218146048580636353", 10)
n5, _ := new(big.Int).SetString("13", 10)
tests = append(tests, n1, n2, n3, n4, n5)
// Test each number
for _, test := range tests {
divisorOne := pollardsRho(test)
divisorTwo := new(big.Int).Div(test, divisorOne)
fmt.Printf("%s = %s * %s\n", test, divisorOne, divisorTwo)
}
}
func pollardsRho(aNumber *big.Int) *big.Int {
// Check if even
if new(big.Int).Mod(aNumber, big.NewInt(2)).Cmp(big.NewInt(0)) == 0 {
return big.NewInt(2)
}
// Generate random numbers in the appropriate range
bitLength := aNumber.BitLen()
constant := randomBigInt(bitLength)
constant.Mod(constant, aNumber)
x := randomBigInt(bitLength)
x.Mod(x, aNumber)
y := new(big.Int).Set(x)
divisor := big.NewInt(1)
for divisor.Cmp(big.NewInt(1)) == 0 {
// x = (x*x + constant) % aNumber
x.Mul(x, x)
x.Add(x, constant)
x.Mod(x, aNumber)
// y = (y*y + constant) % aNumber, applied twice
y.Mul(y, y)
y.Add(y, constant)
y.Mod(y, aNumber)
y.Mul(y, y)
y.Add(y, constant)
y.Mod(y, aNumber)
// Calculate difference and find GCD
diff := new(big.Int).Sub(x, y)
if diff.Sign() < 0 {
diff.Neg(diff)
}
divisor.GCD(nil, nil, diff, aNumber)
}
return divisor
}
func randomBigInt(bitLength int) *big.Int {
// Create a new big.Int with random bits
max := new(big.Int).Lsh(big.NewInt(1), uint(bitLength))
n, err := rand.Int(rand.Reader, max)
if err != nil {
// Fallback to a simpler method if crypto/rand fails
r := new(big.Int)
// This is not cryptographically secure but works as a fallback
r.SetInt64(time.Now().UnixNano())
return r.Mod(r, max)
}
return n
}
- Output:
4294967213 = 57139 * 75167 9759463979 = 98779 * 98801 34225158206557151 = 263166499 * 130051349 763218146048580636353 = 25998278833 * 29356487441 13 = 13 * 1
The following works in both languages.
link random
procedure main(A)
a := (*A > 0,A) | [4294967213,9759463979,34225158206557151,763218146048580636353]
every n := !a do write(n,": ",x := prho(n)," x ",n/x)
end
procedure prho(n)
if n == 1 then return n
if n%2 == 0 then return 2
y := x := rand_num()%(n-2) + 2
c := rand_num()%(n-1) + 1
d := 1
while d = 1 do {
x := (pow(x,2,n)+c+n)%n
every (y|y) := (pow(y,2,n)+c+n)%n
d := gcd(+(x-y),n)
}
return d
end
procedure pow(b,e,m)
r := 1
while e > 0 do
if e%2 = 0 then (e /:= 2, b := (b*b)%m)
else (e -:= 1, r := (r*b)%m, e /:= 2, b := (b*b)%m)
return r
end
Sample run:
->prho 4294967213: 57139 x 75167 9759463979: 98801 x 98779 34225158206557151: 130051349 x 263166499 763218146048580636353: 25998278833 x 29356487441 ->
import java.math.BigInteger;
import java.util.List;
import java.util.concurrent.ThreadLocalRandom;
public final class PollardsRhoAlgorithm {
public static void main(String[] args) {
List<BigInteger> tests = List.of( new BigInteger("4294967213"),
new BigInteger("9759463979"),
new BigInteger("34225158206557151"),
new BigInteger("763218146048580636353"),
new BigInteger("13")
);
tests.forEach( test -> {
final BigInteger divisorOne = pollardsRho(test);
final BigInteger divisorTwo = test.divide(divisorOne);
System.out.println(test + " = " + divisorOne + " * " + divisorTwo);
} );
}
private static BigInteger pollardsRho(BigInteger aNumber) {
if ( ! aNumber.testBit(0) ) {
return BigInteger.TWO;
}
final BigInteger constant = new BigInteger(aNumber.bitLength(), RANDOM);
BigInteger x = new BigInteger(aNumber.bitLength(), RANDOM);
BigInteger y = x;
BigInteger divisor = BigInteger.ONE;
do {
x = x.multiply(x).add(constant).mod(aNumber);
y = y.multiply(y).add(constant).mod(aNumber);
y = y.multiply(y).add(constant).mod(aNumber);
divisor = x.subtract(y).gcd(aNumber);
} while ( divisor.compareTo(BigInteger.ONE) == 0 );
return divisor;
}
private static final ThreadLocalRandom RANDOM = ThreadLocalRandom.current();
}
- Output:
4294967213 = 57139 * 75167 9759463979 = 98801 * 98779 34225158206557151 = 130051349 * 263166499 763218146048580636353 = 25998278833 * 29356487441 13 = 13 * 1
function pollards_rho(number::Integer)
iseven(number) && return 2
bit_length = ndigits(number, base=2)
x, constant = rand(0:2^(bit_length-1), 2)
y = x
divisor = 1
while divisor == 1
x = ( x * x + constant ) % number
y = ( y * y + constant ) % number
y = ( y * y + constant ) % number
divisor = gcd(x - y, number)
end
return divisor;
end
const tests = [ 4294967213, 9759463979, 34225158206557151, 13 ]
for test in tests
divisor_one = pollards_rho(test)
divisor_two = div(test, divisor_one)
@assert rem(test, divisor_one) == 0
println(test, " = ", divisor_one, " * ", divisor_two)
end
- Output:
Same as C++ example.
g[x_Integer, m_Integer] := Mod[x^2 + 1, m];
rho[n_Integer] := Module[{x = 1, y = 1, d = 1, s = 0, lim},
lim = Floor[Sqrt[Abs[n]]];
While[d == 1 && s < lim,
s++;
x = g[x, n];
y = g[g[y, n], n];
d = GCD[Abs[x - y], n];
If[d == n, Return[0]];
];
Return[d];
];
tests = {4294967213, 9759463979, 34225158206557151, 13};
For[i = 1, i <= Length[tests], i++,
p = rho[tests[[i]]];
If[p == 0,
Print[tests[[i]], " no factor found"],
Print[tests[[i]], " = ", p, " * ", tests[[i]]/p]
];
]
- Output:
4294967213 = 57139 * 75167 9759463979 = 98779 * 98801 34225158206557151 = 130051349 * 263166499 13 = 1 * 13
with javascript_semantics
include mpfr.e
for t in {"13","125","217","4294967213","9759463979","34225158206557151","763218146048580636353",
"5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497"} do
integer b = mpz_sizeinbase(mpz_init(t),2)
printf(1,"%s = %s (%d bits)\n",{t,join(mpz_pollard_rho(t,true)," * "),b})
end for
- Output:
13 = 13 (4 bits) 125 = 5 * 5 * 5 (7 bits) 217 = 7 * 31 (8 bits) 4294967213 = 57139 * 75167 (32 bits) 9759463979 = 98779 * 98801 (34 bits) 34225158206557151 = 130051349 * 263166499 (55 bits) 763218146048580636353 = 25998278833 * 29356487441 (70 bits) 5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497 = 165901 * 10424087 * 18830281 * 53204737 * 56402249 * 59663291 * 91931221 * 95174413 * 305293727939 * 444161842339 * 790130065009 (312 bits)
[be warned the last test gets you a blank screen for 16s under p2js/in a web browser, otherwise 0.5s, on my box]
The 'bigint.factors' method uses Pollard's Rho algorithm under the hood and runs this task in around 31 seconds on my machine (core i7).
If 'bigint.factors' is replaced by 'mpz.factors' then the runtime is only 0.17 seconds. This invokes the 'factor' utility which is written in C using GMP and is also said to be based on Pollard's Rho.
require "bignum"
local fmt = require "fmt"
local tests = {13, 4294967213, 9759463979, 34225158206557151, "763218146048580636353", "5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497"}
for tests as test do
local bi = bigint.new(test)
local factors = bigint.factors(bi)
local factstr = factors:map(|f| -> f:tostring()):concat(" * ")
fmt.print("%s = %s (%d bits)", bi, factstr, bi:bitlength())
end
- Output:
13 = 13 (4 bits) 4294967213 = 57139 * 75167 (32 bits) 9759463979 = 98779 * 98801 (34 bits) 34225158206557151 = 130051349 * 263166499 (55 bits) 763218146048580636353 = 25998278833 * 29356487441 (70 bits) 5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497 = 165901 * 10424087 * 18830281 * 53204737 * 56402249 * 59663291 * 91931221 * 95174413 * 305293727939 * 444161842339 * 790130065009 (312 bits)
This solution was based on the C solution.
import math
def g(x: int, m: int):
return (x ** 2 + 1) % m
def rho(n: int):
x = 1
y = x
d = 1
s = 0
lim = int(math.sqrt(abs(n)))
while d == 1 and s < lim:
s += 1
x = g(x, n)
y = g(g(y, n), n)
d = math.gcd(abs(x - y), n)
if d == n:
return 0
return d
def main():
tests = { 4294967213, 9759463979, 34225158206557151, 13 }
for i in tests:
p = rho(i)
if p == 0:
print(f"{i} no factor found")
else:
print(f"{i} = {p} * {i / p:.0f}")
if __name__ == '__main__':
main()
- Output:
34225158206557151 = 130051349 * 263166499 9759463979 = 98779 * 98801 13 = 13 * 1 4294967213 = 57139 * 75167
rho takes two arguments; 2OS is the number to be factored, TOS is the seed value for the PRNG g, typically 2. This allows rho to be run again with a different seed if the algorithm fails, as described in the wikipedia article. This is not illustrated here.
rho returns two results. TOS is a boolean indicating success (true = 1), or failure (false = 0). 2OS is the found factor, if the algorithm succeeded, and should be discarded if the algorithm failed.
The PRNG, g, requires a number on the top of the ancillary stack n.
[ [ dup while
tuck mod again ]
drop abs ] is gcd ( n n --> n )
[ stack ] is n ( --> s )
[ stack ] is d ( --> s )
[ dup * 1+ n share mod ] is g ( n --> n )
[ swap n put 1 d put
dup
[ d share 1 = while
dip g g g
2dup - abs
n share gcd
d replace again ]
2drop
d take n take over != ] is rho ( n n --> n b )
' [ 1234571407 4294967213 9759463979
34225158206557151 763218146048580636353 ]
witheach
[ cr dup echo
dup 2 rho not iff
[ 2drop say " might be prime." ]
done
say " = " dup echo say " * " / echo ]- Output:
1234571407 might be prime. 4294967213 = 57139 * 75167 9759463979 = 98779 * 98801 34225158206557151 = 130051349 * 263166499 763218146048580636353 = 25998278833 * 29356487441
pollards_rho <- function(number) {
# Check if even
if (number %% 2 == 0) return(2)
# Get bit length
bit_length <- floor(log2(number)) + 1
# Generate random values
x <- sample(0:(2^(bit_length-1)), 1)
constant <- sample(0:(2^(bit_length-1)), 1)
y <- x
divisor <- 1
while (divisor == 1) {
x <- (x * x + constant) %% number
y <- (y * y + constant) %% number
y <- (y * y + constant) %% number
divisor <- gcd(x - y, number)
}
return(divisor)
}
# GCD function (not in base R)
gcd <- function(a, b) {
while (b != 0) {
temp <- b
b <- a %% b
a <- temp
}
return(abs(a))
}
# Disable scientific notation
options(scipen = 999)
# Test cases
tests <- c(4294967213, 9759463979, 34225158206557151, 13)
for (test in tests) {
divisor_one <- pollards_rho(test)
divisor_two <- test %/% divisor_one
stopifnot(test %% divisor_one == 0)
cat(test, " = ", divisor_one, " * ", divisor_two, "\n")
}
- Output:
4294967213 = 4294967213 * 1 9759463979 = 98779 * 98801 34225158206557152 = 2 * 17112579103278576 13 = 13 * 1
Copied verbatim (with different test values) from the Prime decomposition page. This finds all prime factors of an integer of any size. Very large integers, especially those with two or more factors larger than 40 bit may take a while, but it will find them eventually.
Raku handles big integers seamlessly and transparently. No special handling required.
sub prime-factors ( Int $n where * > 0 ) {
return $n if $n.is-prime;
return () if $n == 1;
my $factor = find-factor( $n );
sort flat ( $factor, $n div $factor ).map: &prime-factors;
}
sub find-factor ( Int $n, $constant = 1 ) {
return 2 unless $n +& 1;
if (my $gcd = $n gcd 6541380665835015) > 1 { # magic number: [*] primes 3 .. 43
return $gcd if $gcd != $n
}
my $x = 2;
my $rho = 1;
my $factor = 1;
while $factor == 1 {
$rho = $rho +< 1;
my $fixed = $x;
my int $i = 0;
while $i < $rho {
$x = ( $x * $x + $constant ) % $n;
$factor = ( $x - $fixed ) gcd $n;
last if 1 < $factor;
$i = $i + 1;
}
}
$factor = find-factor( $n, $constant + 1 ) if $n == $factor;
$factor;
}
.put for (125, 217, 4294967213, 9759463979, 34225158206557151, 763218146048580636353,
5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497)\
.hyper(:1batch).map: -> $n {
my $start = now;
"factors of $n: ({1+$n.msb} bits) ",
prime-factors($n).join(' × '), " \t in ", (now - $start).fmt("%0.3f"), " sec."
}
- Output:
factors of 125: (7 bits) 5 × 5 × 5 in 0.002 sec. factors of 217: (8 bits) 7 × 31 in 0.001 sec. factors of 4294967213: (32 bits) 57139 × 75167 in 0.001 sec. factors of 9759463979: (34 bits) 98779 × 98801 in 0.001 sec. factors of 34225158206557151: (55 bits) 130051349 × 263166499 in 0.021 sec. factors of 763218146048580636353: (70 bits) 25998278833 × 29356487441 in 0.303 sec. factors of 5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497: (312 bits) 165901 × 10424087 × 18830281 × 53204737 × 56402249 × 59663291 × 91931221 × 95174413 × 305293727939 × 444161842339 × 790130065009 in 5.524 sec.
Include: How to use
Include: Source code
Pollard-Rho tends to find low factors first, but you cannot rely on that. To avoid an endless loop, my implementation (Pollardrho in Numbers) has an outer loop of 3 attempts with different random parameters and an inner loop of 1 million times. If no factor is found within these borders, the procedure fails. On a typical machine, this will be after a few minutes.
-- 15 Nov 2025
include Setting
numeric digits 200
SAY 'POLLARD''S RHO ALGORITHM'
say version
say
call TestNumbers
call Selected
call Randomized
exit
TestNumbers:
procedure expose Test.
t='31',
'4294967213',
'9759463979',
'34225158206557151',
'576460752303423487',
'147573952589676412927',
'763218146048580636353',
'5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497'
do i=1 to Words(t)
Test.i=Word(t,i)
end
Test.0=i-1
return
Selected:
procedure expose Mult. Test. Memo.
call Time('r')
say 'Find a factor for 8 selected numbers...'
do t=1 for Test.0
n=Test.t
call Time('r')
f=Pollardrho(n)
if f=0 then
u='failed'
else
u=f 'x' n/f
say Format(Time('e'),4,3)'s' n '('Xpon(n)+1 'digits) =' u
end
say
return
Randomized:
procedure expose Mult. Memo.
say 'Find a factor for 30 Random numbers...'
x=0
do until x=30
n=''
do Random(1,30)
n||=Random(9)
end
n+=0
if Pos(Right(n,1),024568) > 0 | Digitsum(n)//3 = 0 then
iterate
call Time('r')
f=Pollardrho(n)
if f=0 then
u='failed'
else
u=f 'x' n/f
say Format(Time('e'),4,3)'s' n '('Xpon(n)+1 'digits) =' u
x+=1
end
say
return
-- Xpon (get exponent)
include Basic
-- Pollardrho (find a factor)
include Ntheory
-- Digitsum (calculate digital sum)
include Special
- Output:
POLLARD'S RHO ALGORITHM REXX-ooRexx_5.2.0(MT)_64-bit 6.05 25 Oct 2025 Find a factor for 8 selected numbers... 0.000s 31 (2 digits) = 31 x 1 0.005s 4294967213 (10 digits) = 75167 x 57139 0.000s 9759463979 (10 digits) = 98801 x 98779 0.376s 34225158206557151 (17 digits) = 130051349 x 263166499 0.027s 576460752303423487 (18 digits) = 179951 x 3203431780337 0.619s 147573952589676412927 (21 digits) = 193707721 x 761838257287 3.211s 763218146048580636353 (21 digits) = 25998278833 x 29356487441 0.228s 5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497 (94 digits) = 165901 x 32945014744179408009060442759193408327096567196437142390968926669880637276718920200636197 Find a factor for 30 Random numbers... 0.002s 24259868215004138240819873 (26 digits) = 373 x 65039861166230933621501 0.000s 59557 (5 digits) = 59557 x 1 0.000s 869 (3 digits) = 11 x 79 0.001s 709731792904816829693 (21 digits) = 7 x 101390256129259547099 0.001s 5966934138985286315110860533 (28 digits) = 13 x 458994933768098947316220041 0.000s 353210170615905187739 (21 digits) = 7 x 50458595802272169677 0.001s 6462829 (7 digits) = 47 x 137507 0.000s 796416043693508693 (18 digits) = 7 x 113773720527644099 0.000s 59 (2 digits) = 59 x 1 0.000s 51761 (5 digits) = 271 x 191 0.005s 832366642689786692938445707 (27 digits) = 661 x 1259253619803005586896287 0.001s 5208674690037789311 (19 digits) = 11 x 473515880912526301 0.001s 237215059152773003651255867 (27 digits) = 23 x 10313698224033608854402429 0.001s 8052474889482857215709096749 (28 digits) = 7 x 1150353555640408173672728107 8.049s 978417738904120867499291 (24 digits) = 46126308661 x 21211706882831 0.003s 283405191016889805549881357 (27 digits) = 317 x 894022684595866894479121 0.000s 621889 (6 digits) = 19 x 32731 0.000s 45107 (5 digits) = 43 x 1049 0.000s 1 (1 digits) = 1 x 1 0.000s 852031 (6 digits) = 852031 x 1 0.001s 412888731886589 (15 digits) = 7 x 58984104555227 0.000s 4274020223 (10 digits) = 11 x 388547293 0.000s 768352124663 (12 digits) = 191 x 4022785993 0.001s 390637859057255220859194877 (27 digits) = 13 x 30049066081327324681476529 0.000s 46411 (5 digits) = 46411 x 1 0.006s 4288367478478251313519427 (25 digits) = 5209 x 823261178436984318203 0.000s 634817738029 (12 digits) = 137 x 4633706117 0.000s 605104291 (9 digits) = 11 x 55009481 0.001s 79995062083846819 (17 digits) = 13 x 6153466314142063 0.001s 697145365593908803910591 (24 digits) = 11 x 63376851417628073082781
« SQ 1 + ←n MOD » 'G' STO « → ←n « 2 2 1 WHILE DUP 1 == REPEAT ROT G ROT G G ROT DROP DUP2 - ←n GCD END ←n MOD NIP NIP » » 'RHO' STO « DUP RHO IF DUP THEN SWAP OVER / 2 →LIST ELSE DROP 1 →LIST END » 'TASK' STO
{ 3082021 4294967213 9759463979 34225158206557151 } 1 « TASK » DOLIST
- Output:
1: { { 3082021 } { 57139 75167 } { 98779 98801 } { 130051349 263166499 } }
Breaking 763218146048580636353 in a reasonable time seems out of range.
use num_integer::gcd;
use rand::Rng;
fn pollards_rho(number: i128) -> i128 {
if number & 1 == 0 {
return 2;
}
let bit_length = number.ilog2() as u32;
let mut x = rand::thread_rng().gen_range(0..2_i64.pow(bit_length)) as i128;
let constant = rand::thread_rng().gen_range(0..2_i64.pow(bit_length)) as i128;
let mut y = x;
let mut divisor = 1;
while divisor == 1 {
x = ( x * x + constant ) % number;
y = ( y * y + constant ) % number;
y = ( y * y + constant ) % number;
divisor = gcd(x - y, number);
}
return divisor;
}
fn main() {
let tests = [ 4294967213, 9759463979, 34225158206557151, 13 ];
for test in tests {
let divisor_one = pollards_rho(test);
let divisor_two = test / divisor_one;
println!("{} = {} * {}", test, divisor_one, divisor_two);
}
}
- Output:
4294967213 = 57139 * 75167 9759463979 = 98801 * 98779 34225158206557151 = 263166499 * 130051349 13 = 13 * 1
Demonstration version; only limit number of iterations, no different modulo function, no repeat. C mode is used (with keywords instead of symbols) intended for publication.
\C C mode on
pollards rho (n):
x = random next integer %% n \ starting value, for calling again
y = x
d = 1
s = 0
lim = math int sqrt n \ reverse step counter
while d = 1 && s < lim
s += 1
\ debug print s, x, y, d
x = g x mod n
y = g (g y mod n) mod n
d = math gcd (x-y).abs, n
if d = n
debug print "not found after ", s, "rounds"
return () \ failure
debug print "Found", d, "after", s, "interations"
return d \ factor found
\ the generator function
g (x) mod (m):
y = x^2 + 1
return y %% m
\c C mode off
\ testing
main (parm):+
\ debug enable
i =: string parm[1] as integer else 9759463979
p =: pollards rho i
? p = ()
print i, "no factor found"
|
print i, '=', p, '*', i//p
- Output:
9759463979 = 98779 * 98801
Wren can only deal accurately with integers less than 2^53 so we need to use BigInt to find factors of the 55 bit and 70 bit numbers in the set examples.
The BigInt class already contains a version of the Pollard's Rho algorithm which is normally used as part of a routine to find the prime factors of a given number. However, it can be called directly so that's what we do here.
As it's written entirely in Wren, it's not very fast but does find a factor for the 55 bit number in under 5 seconds though the 70 bit number takes around 155 seconds on my Core i7 machine.
import "./big" for BigInt
import "./fmt" for Fmt
var tests = [13, 4294967213, 9759463979, "34225158206557151", "763218146048580636353"]
for (test in tests) {
var bi = BigInt.new(test)
var fac = BigInt.pollardRho(bi)
if (fac > BigInt.zero) {
Fmt.print("$i = $i * $i ($i bits)", bi, fac, bi/fac, bi.bitLength)
} else {
Fmt.print("$i : $s", bi, "Pollard's Rho failed to find a factor.")
}
}
- Output:
13 = 13 * 1 (4 bits) 4294967213 = 57139 * 75167 (32 bits) 9759463979 = 98779 * 98801 (34 bits) 34225158206557151 = 130051349 * 263166499 (55 bits) 763218146048580636353 = 25998278833 * 29356487441 (70 bits)
Alternatively, we can speed things up around 240 times by using GMP which produces identical output in under 0.7 seconds.
import "./gmp" for Mpz
import "./fmt" for Fmt
var tests = ["13", "4294967213", "9759463979", "34225158206557151", "763218146048580636353"]
for (test in tests) {
var bi = Mpz.fromStr(test)
var fac = Mpz.pollardRho(bi)
if (fac > Mpz.zero) {
Fmt.print("$i = $i * $i ($i bits)", bi, fac, bi/fac, bi.sizeInBits)
} else {
Fmt.print("$i : $s", bi, "Pollard's Rho failed to find a factor.")
}
}
const std = @import("std");
fn gcd(a: i128, b: i128) i128 {
var x = @abs(a);
var y = @abs(b);
while (y != 0) {
const temp = y;
y = x % y;
x = temp;
}
return @as(i128, @intCast(x));
}
fn pollards_rho(number: i128) i128 {
if (number & 1 == 0) {
return 2;
}
var prng = std.Random.DefaultPrng.init(@as(u64, @truncate(@as(u128, @bitCast(std.time.nanoTimestamp())))));
const rand = prng.random();
const bit_length: u6 = @intCast(std.math.log2_int(u128, @intCast(number)));
var x: i128 = @intCast(rand.intRangeAtMost(i64, 0, std.math.pow(i64, 2, bit_length)));
const constant: i128 = @intCast(rand.intRangeAtMost(i64, 0, std.math.pow(i64, 2, bit_length)));
var y: i128 = x;
var divisor: i128 = 1;
while (divisor == 1) {
x = @rem( (x * x + constant) , number);
y = @rem( (y * y + constant) , number);
y = @rem( (y * y + constant) , number);
divisor = gcd(x - y, number);
}
return divisor;
}
pub fn main() !void {
const tests = [_]i128{ 4294967213, 9759463979, 34225158206557151, 13 };
for (tests) |my_test| {
const divisor_one = pollards_rho(my_test);
const divisor_two = @divTrunc(my_test, divisor_one);
std.debug.print("{} = {} * {}\n", .{ my_test, divisor_one, divisor_two });
}
}
- Output:
4294967213 = 57139 * 75167 9759463979 = 98801 * 98779 34225158206557151 = 263166499 * 130051349 13 = 13 * 1