Pollard's rho algorithm

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

References

Wikipedia: Pollard's rho algorithm

ALGOL 68

Works with: ALGOL 68 Genie version 2 and 3 - tested with releases 2.8.3.win32 and 3.5.9 Win32

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.

Library: ALGOL 68-primes

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

AWK

Note that AWK uses internally floating point numbers, thus the programme fails for larger numbers.

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
}

C

#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;
}

C++

#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

Dart

Translation of: C++

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');
}

FreeBASIC

Translation of: C++

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)

Java

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

Julia

Translation of: C++

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.

Phix

Translation of: Wren

Library: Phix/mpfr

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]

Quackery

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

Raku

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.

Rust

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

TAV

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
  print system cp

Output:

9759463979 = 98779 * 98801
0.000740373 sec

Wren

Library: Wren-big

Library: Wren-fmt

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)

Library: Wren-gmp

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