/*
 *Pollard's rho Algorithm
 *Richard Brent's variant
 *
 * n = integer to be factored
 * 0 <= x0 <= n
 * m > 0
 *
 * y = x0, r = 1, q = 1
 * Do: 
 *    x = y
 *    for i = 1 to r
 *       y = f(y)
 *    k = 0
 *    Do:
 *       ys = y
 *       for i = 1 to min(m, r-k)
 *          y = f(y)
 *          q = ( q x |x - y| ) mod n
 *       g = GCD(q, n)
 *       k = k + m
 *    until (k >= r or g > 1)
 *    r = 2r
 * Until g > 1
 * if g == n then
 *    Do:
 *       ys = f(ys)
 *       g = GCD(| x - ys| ,n)
 *    until g > 1
 * if g == n
 *    return failure
 * else
 *    return g
 *
 */

#include <stdio.h>
#include "gmp.h"

char buff[100];
mpz_t n;

int i,r,k,m;

int min(int a, int b){
   if (a < b)
      return a;
   return b;
}

//f(x) = (x^2 + 123) mod n
void f(mpz_t a){
   mpz_t result;
   mpz_init(result);
   mpz_powm_ui(result, a, 2, n);
   mpz_add_ui(result, result, 123);
   mpz_mod(a, result, n);
}

int pollardrho_brent(mpz_t res, mpz_t n){
   long long k, r = 1, m = 43;
   mpz_t x, y, ys, q, gcd, tmp;
   mpz_init(gcd);
   mpz_init(tmp);
   mpz_init(ys);
   mpz_init(x);
   mpz_init_set(y,n);
   mpz_init_set_ui(q,1);
   do {
      mpz_set(x, y);
      for (int i=0; i<r; i++)
         f(y);
      k = 0;
      do {
         mpz_set(ys, y);
	 for (int i=0; i< min(m, r-k); i++){
	    f(y);
	    mpz_sub(tmp, x, y);
	    mpz_abs(tmp, tmp);
	    mpz_mul(q, q, tmp);
	    mpz_mod(q, q, n);
	 }
	 mpz_gcd(gcd, q, n);
	 k += m;
      }while (k<r && mpz_cmp_ui(gcd,1)<=1);
      r *= 2;
   }while(mpz_cmp_ui(gcd,1)<=0);
   if (mpz_cmp(gcd, n)==0){
      do {
         f(ys);
	 mpz_sub(tmp, x, ys);
	 mpz_abs(tmp, tmp);
	 mpz_gcd(gcd, tmp, n);
      }while (mpz_cmp_ui(gcd,1)>0);
   }
   if (mpz_cmp(gcd,n)==0){
      return 0;
   }else{
      mpz_set(res, gcd);
      return 1;
   }
}

int main(){
   mpz_init(n);
   printf("Input a number to be factored: ");
   scanf("%s",buff);
   mpz_set_str(n,buff,10);
   if (pollardrho_brent(n,n)==0){
      printf("FAILURE\n");
   }else{
      mpz_get_str(buff, 10, n);
      printf("Factor: %s\n",buff);
   }
   return 0;
}