/*
 * rho.c
 * by jeff hamblin (hamblin@cs.wisc.edu)
 *
 * pollard's rho method of factoring 
 *
 * be sure to link in libgmp
 */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <gmp.h>

#define DEBUG 0

void f(mpz_t rop, mpz_t x);

int main(int argc, char **argv)
{
	mpz_t composite, x_2k, x_k, diff, gcd, a, bound;
	int x_0 = 2, prime_count = 0, k;
	FILE *log;
	char filename[80];

	if (argc != 2)
	{
		printf("usage: %s composite\n", argv[0]);
		exit(-1);
	}

	mpz_init_set_str(composite, argv[1], 10);
	mpz_init(x_k);
	mpz_init(x_2k);
	mpz_init(diff);
	mpz_init(gcd);
	mpz_init(a);
	mpz_init(bound);
	
	sprintf(filename, "rho.%d.log", getpid());
	log = fopen(filename, "w");

	printf("---\nPerforming rho Algorithm\n---\n");
	fprintf(log, "---\nPerforming rho Algorithm\n---\n");
	mpz_out_str(stdout, 10, composite);
	mpz_out_str(log, 10, composite);
	printf(" = \n");
	fprintf(log, " = \n");
	
	while (mpz_mod_ui(a, composite, 2) == 0)
	{
		prime_count ++;
		printf("\tfactor: 2\n");
		fprintf(log, "\tfactor: 2\n");
		mpz_div_ui(composite, composite, 2);
	}

	if (mpz_cmp_ui(composite, 1) == 0)
	{
		exit(prime_count);
	}

	if (mpz_probab_prime_p(composite, 25) == 1) 
	{
		prime_count ++;
		printf("\tfactor: ");
		fprintf(log, "\tfactor: ");
		mpz_out_str(stdout, 10, composite);
		mpz_out_str(log, 10, composite);
		printf("\n");
		fprintf(log, "\n");
		exit(prime_count);
	}

	while (1)	
	{
		mpz_sqrt(bound, composite);
		printf("---\nusing x_0 = %d\n---\n", x_0);
		fprintf(log, "---\nusing x_0 = %d\n---\n", x_0);
		x_0 ++;
		mpz_set_ui(x_k, x_0);
		mpz_set_ui(x_2k, x_0);
		fflush(log);
		while (mpz_cmp_ui(bound, 0) > 0)	
		{
			mpz_sub_ui(bound, bound, 1);
			f(x_k, x_k);
			mpz_mod(x_k, x_k, composite);
			f(x_2k, x_2k);
			mpz_mod(x_2k, x_2k, composite);
			f(x_2k, x_2k);
			mpz_mod(x_2k, x_2k, composite);
			mpz_sub(diff, x_2k, x_k);
			mpz_gcd(gcd, diff, composite);
			if (mpz_cmp_ui(gcd, 1) > 0)
			{
				prime_count ++;
				printf("\tfactor: ");
				fprintf(log, "\tfactor: ");
				mpz_out_str(stdout, 10, gcd);
				mpz_out_str(log, 10, gcd);
				printf("\n");
				fprintf(log, "\n");
				mpz_div(composite, composite, gcd);
	
				if (mpz_cmp_ui(composite, 1) == 0)
				{
					exit(prime_count);
				}
	
				if (mpz_probab_prime_p(composite, 25) == 1)
				{
					prime_count ++;
					printf("\tfactor: ");
					fprintf(log, "\tfactor: ");
					mpz_out_str(stdout, 10, composite);
					mpz_out_str(log, 10, composite);
					printf("\n");
					fprintf(log, "\n");
					exit(prime_count);
				}
				fflush(log);
				mpz_sqrt(bound, composite);
			}
		}
	}
}

void f(mpz_t rop, mpz_t x)
{
	mpz_mul(rop, x, x);
	mpz_add_ui(rop, rop, 1);
}