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

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

#define K_BOUND 1000000
#define A_BOUND 10
#define GCD_LOOP_RATE 1 


int main(int argc, char **argv)
{
	mpz_t composite, k_lcm, a, gcd, a_minus_1;
	int a_plug = 2, prime_count = 0;
	long k = 2;
	int gcd_loop = 0;

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

	mpz_init_set_str(composite, argv[1], 10);
	mpz_init(a);
	mpz_init(a_minus_1);
	mpz_init_set_ui(k_lcm, 1);
	mpz_init(gcd);

	printf("---\nPerforming P-1 Algorithm\n---\n");
	mpz_out_str(stdout, 10, composite);
	printf(" = \n");

	while (mpz_mod_ui(a, composite, 2) == 0)
	{
		prime_count ++;
		printf("\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: ");
		mpz_out_str(stdout, 10, composite);
		printf("\n");
		exit(prime_count);
	}

	for (a_plug = 2; a_plug < A_BOUND; a_plug ++) 
	{
		printf("---\na = %d\n---\n", a_plug);
		mpz_set_ui(a, a_plug);
		for (k = 1; k < K_BOUND; k ++) 
		{
			mpz_gcd_ui(gcd, k_lcm, k); /* find gcd of k and (k-1)!*/
			mpz_divexact(k_lcm, k_lcm, gcd); /* make it lcm */
			mpz_mul_ui(k_lcm, k_lcm, k); /* put in new k */ 
			mpz_powm(a, a, k_lcm, composite);
			mpz_sub_ui(a_minus_1, a, 1);
			gcd_loop ++;
			if (gcd_loop == GCD_LOOP_RATE) 
			{
				gcd_loop = 0;
				mpz_gcd(gcd, a_minus_1, composite);
				if ((mpz_cmp(gcd, composite) < 0) && (mpz_cmp_ui(gcd, 1) > 0))
				{
					prime_count ++;
					printf("\tfactor: ");
					mpz_out_str(stdout, 10, gcd);
					printf("\n");
					mpz_div(composite, composite, gcd);
					if (mpz_probab_prime_p(composite, 25) == 1)
					{
						prime_count ++;
						printf("\tfactor: ");
						mpz_out_str(stdout, 10, composite);
						printf("\n");
						exit(prime_count);
					}
				}
			}
		}
	}
	printf("Remaining unfactored component:\n");
	mpz_out_str(stdout, 10, composite);
	exit(prime_count);
}