import java.math.BigInteger;
import java.security.SecureRandom;

/******************************************************************************
 *                            Pollard Rho Factorization                       *
 * 									      *
 *  				Christer Karlsson                             *
 *                                                                            *								
 *  Let n be a large composite number, and p its smallest prime divisor.      *
 *  Choose integers x[0],x[1],...x[s] s.t that they have distinct least       * 
 *  non-negative residues mod n, but their least non-negative rsidues mod p   *
 *  is not all distinct.  Once integers x[i] and x[j] 0 <= i <= j <= s        *
 *  such that x[i] ~ x[j] (mod p) but x[i] != x[j] (mod n) are found          *
 *  it follows that (x[i]-x[j],n) is a non trivial divisor of n as p divides  *
 *  x[i] - x[j], but n does not. To find (x[i]-x[j],n) can be done quicly     *
 *  using Euclidian algorith. To find x[i] and x[j] we first seed x[0] which  * 
 *  is chosen randomly and a polynomial function f(x) with integer coef. > 1  * 
 *  Computer x[k], k = 1,2,3, recursivle using:                               *
 *                           x[k+1] ~ f(x[k]) (mod n) 0 <= x[k+1] < n         *
 *  The polynomial f(x) should be picked such that the probability is high    *
 *  that a large number of integers x[i] are generated before they repeat.    *
 *  Tests have shown that f(x) = x^2 +1 performs well.                        *
 *                                                                            *
 ******************************************************************************/

public class PollardRho {

   public static void main( String[] args)
   {
	BigInteger one      = new BigInteger("1");
	BigInteger two	    = new BigInteger("2");		
	SecureRandom random = new SecureRandom();

    	// Variables
	boolean solved = false;		
	int bitlength = 32;				// Size of our prime numbers
	BigInteger i  = one;				// 1st Counter
    	BigInteger k  = two;				// 2nd Counter
    	BigInteger L  = two;				// Upper limit
    	BigInteger p1;					// 1st Prime number
       	BigInteger p2;					// 2nd Prime number    	
       	BigInteger n  = new BigInteger("5157437");	// p*q The composite number
       	BigInteger n1;					// n-1											
       	BigInteger d;					// Remainder
       	BigInteger x1;
    	BigInteger y;					// Part two       	
       	BigInteger temp;

       	
        // generate a test value for p1
        // 32 bits long, probably prime.
        p1 = new BigInteger( "16" );  // Set p1 to composite number
        while(!p1.isProbablePrime(100))
        {
        	p1 = new BigInteger( bitlength,		// size in bits       		//
                                       100,		// 1 in 100 chance not prime 	//
                                       random );
        }
        if(p1.isProbablePrime(100))
        	System.out.print("Prime           p1: " );
        System.out.println( p1 +" " + p1.bitLength() +" bits long");
        // generate a test value for p2
        // 32 bits long, probably prime.
        p2 = new BigInteger( "16" );  // Set p2 to composite number
        while(!p2.isProbablePrime(100))
        {
        	p2 = new BigInteger( bitlength,		// size in bits       		//
        			       100,		// 1 in 100 chance not prime	//
                                       random );
        }
        if(p2.isProbablePrime(100))
        	System.out.print("Prime           p2: " );
        System.out.println( p2 +" " + p2.bitLength() +" bits long");
        // Create a composite number
        n = p1.multiply(p2);
        System.out.print("Composite Number n: " );
        System.out.println( n +" " + n.bitLength() +" bits long");
              
        // Let x1 equal a random number between 0 and n-1
        n1 = n.subtract(one);
        temp = n1;
        // Set the bitlength to the size of n1
        bitlength = n1.bitLength();
        // temp >= n-1
        while (temp.compareTo(n1) >=0)
        {	
        	temp = new BigInteger(bitlength, random);
        }
        x1 = temp;
        System.out.println("...");
        // y = x[i]
        y = x1;
        
        // Set upper limit to n^(1/2)
        L = sqrt(n);
        // Start solving using generating polynomial f(x) = x^2+1
        while (!solved && i.compareTo(L) < 0)
        {
        	// i = i + 1
        	i = i.add(one);
        	// x[i] = (x[i-1]^2+1) mod n
        	x1 = ((x1.pow(2)).add(one)).mod(n);
        	// d = gcd (y - x[i], n)
        	d = (y.subtract(x1)).gcd(n);
        	// if ( d != 1 and d != n)
        	if (d.compareTo(one) != 0 && d.compareTo(n) !=0)
        	{
        		System.out.println( "We found a solution in " +i +" iterations!");
        		System.out.print( n +" = " + d + " * " +n.divide(d));
        		solved = true;
        	}
        	// if ( i == k)
        	if ( i.compareTo(k) == 0)
        	{
        		// y = x[i]
        		y = x1;
        		// k = k*2
        		k = k.multiply(two);
        	}
        }
        // We found no factor
        if(!solved)
        	System.out.println("Unable to find a factorization with " +L +" iterations");        
	}

	public static BigInteger sqrt(BigInteger x)
	{
	    BigInteger p, s, one, two;	    
		one = new BigInteger("1");
		two	= new BigInteger("2");	
		
	    p = x.divide(two);
	    s = (x.divide(p.pow(1))).add(p.multiply(two.subtract(one)));
	    s = s.divide(two);

	    while (p.compareTo(s) != 0 && one.compareTo(p.subtract(s)) !=0 )
	    {
	    	p = s;
	    	s = (x.divide(p.pow(1))).add(p.multiply(two.subtract(one)));
	    	s = s.divide(two);
	    }
	    // Was x a perfect square?
	    if( (s.multiply(s)).compareTo(x) == 0)
	    	return (s);
	    // Make sure that s is the floor of sqrt(x)
	    while( (s.multiply(s)).compareTo(x) < 0)
	    		s = s.add(one);
        s = s.subtract(one);
        return (s);
   }		
}