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

/*****************************************************************************
 *                      Shank's square form factorization                    *
 * 									     *
 * 				Christer Karlsson                            *
 *                                                                           *								
 * This is an improvement of Fermat's factorization method. In order to      * 
 * succeed Fermat's method depends on finding integers s and t such that     * 
 * n = s^2 - t^2, a possible improvement is to search for integers s and t   *
 * such that $s^2 \equiv t^2$ (mod $n$). Finding such a pair $(s,y)$ does    * 
 * not necessarily guarantee a factorization of n, but it do imply that n    *
 * is a factor of s^2-t^2 =(s-t)(s+t, and that therefore there is a really   *
 * good chance that the prime divisors of n are distributed between these    *
 * two factors, in such a way that calculation of the highest common factor  *
 * of n and s-t will yield a nontrivial factor of n.                         *
 *                                                                           *
 *****************************************************************************/
public class Shank_Square {
   public static void main( String[] args)
   {
	BigInteger one      = new BigInteger("1");		// Constants		
	SecureRandom random = new SecureRandom();
    	// Variables
	boolean solved = false;
	boolean pfsq   = false;
	int bitlength = 32;					// Size of our prime numbers
	BigInteger i = one;					// 1st Counter
	BigInteger j = one;					// 2nd Counter
	BigInteger L = one;					// Upper limit		
       	BigInteger prime1 = new BigInteger("4088422847");	// 1st Prime number
       	BigInteger prime2 = new BigInteger("4088422847");	// 2nd Prime number
    	BigInteger n      = new BigInteger("5157437");		// Our random composite number      	
    	BigInteger sn;						// Square-root of n
    	BigInteger p0 = one;					// p[0]
       	BigInteger p1 = one;					// p[1]       	
       	BigInteger q0 = one;					// q[0]
       	BigInteger q1 = one;					// q[1]
       	BigInteger q2 = one;					// q[1]       	
       	BigInteger b0 = one;					// b[0]
       	BigInteger r  = one;					// r      	
       	
        // generate a test value for p1
        // 32 bits long, probably prime.
       	prime1 = new BigInteger( "16" );  // Set p1 to composite number
        while(!prime1.isProbablePrime(100))
        {
        	prime1 = new BigInteger( bitlength,	// size in bits       		//
                    			 100,		// 1 in 100 chance not prime 	//
                    			 random );
        }
        if(prime1.isProbablePrime(100))
        	System.out.print("Prime           p1: " );
        System.out.println( prime1 +" " + prime1.bitLength() +" bits long");
        // generate a test value for p2
        // 32 bits long, probably prime.
        prime2 = new BigInteger( "16" );  // Set p2 to composite number
        while(!prime2.isProbablePrime(100))
        {
        	prime2 = new BigInteger( bitlength,	// size in bits       		//
       			 		 100,		// 1 in 100 chance not prime //
       			 		 random );
        }
        if( prime2.isProbablePrime(100))
        	System.out.print("Prime           p2: " );
        System.out.println( prime2 +" " + prime2.bitLength() +" bits long");

        // Create a composite number
        n = prime1.multiply(prime2);
        System.out.print("Composite Number n: " );
        System.out.println( n +" " + n.bitLength() +" bits long");
        System.out.println("...");
        // Let p[0] = sqrt(n);
        sn = sqrt(n);
        // Is n a perfect square?
        if( (sn.multiply(sn)).compareTo(n) == 0)
        {
    		solved = true; pfsq = true; 
        }        

        // Initialize p[0], q[0] and q[1]
        p0 = sn;
        q1 = n.subtract(p0.pow(2));
        // Start solving
        // Step one a q that is a perfect square
        while (!solved && !pfsq)
        {
      		b0 = (sn.add(p0)).divide(q1);
      		p1 = (b0.multiply(q1)).subtract(p0);
      		q2 = q0.add(b0.multiply(p0.subtract(p1)));
      		// Is q2 a perfect square?
      		r = sqrt(q2);
            if( (r.multiply(r)).compareTo(q2) == 0)        		
        		solved = true;        	
            if(!solved)
            {
            	p0 = p1;
            	q0 = q1;
            	q1 = q2;
            	i = i.add(one);
            }
        }
        System.out.println("Found a perfect square after "
        		             +i +" iterations");        
        // Initialize b[0], p[0], q[0] and q[1]
        // Unless n was a perfect square
        if(!pfsq)
        {
        	b0 = (sn.subtract(p1)).divide(sqrt(q2));
        	p0 = b0.multiply(sqrt(q2)).add(p1);
        	q0 = sqrt(q2);
        	q1 = (n.subtract(p0.pow(2))).divide(q0);
            solved = false;
        }
        // Set upper limit of iterations
        L = sqrt(n);
        // Step two find p[i+1] = p[i]
        while (!solved && !pfsq && j.compareTo(L) < 0)
        {       	
        	b0 = (sn.add(p0)).divide(q1);
        	p1 = (b0.multiply(q1)).subtract(p0);
        	q2 = q0.add(b0.multiply(p0.subtract(p1)));
        	// Is p[i+1] = p[i]?
        	if(p1.compareTo(p0) == 0)
        		solved = true;
        	if(!solved)
        	{
            	p0 = p1;
            	q0 = q1;
            	q1 = q2;
            	j = j.add(one);
        	}
        }       
        // We found no factor
        if(!solved || (n.gcd(p0)).compareTo(one) == 0)
        {
        	if((n.gcd(p0)).compareTo(one) == 0)
        		System.out.println("Found only trivial solution after " 
        				              +i.add(j) +" iterations");
        	else
        	    System.out.println("Unable to find a factorization with " +L +" iterations");
        }
        else
        {
        	System.out.println( "We found a solution in " +i  +" + " +j 
        			           +" = " +i.add(j) +" iterations");
        	if(pfsq)
        		r = p0;
        	else
        		r = n.gcd(p0);
        	System.out.println( n +" = " + r + " * " +n.divide(r));        	
        }
	}

	// Simple Square root function
	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);
	}	
}