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

/* 
 * td.c 
 * jeff hamblin (hamblin@cs.wisc.edu)
 *
 * this is a simple trial division algorithm on arbitrary length integers.
 * given an n to factor, first test it to see if it is prime.  if it is, then
 * quit (since it has no factors other that 1 and itself).  if n fails the
 * primality test, start trial dividing by integers 2 to sqrt(n).  we do
 * trial division by all of these integers (except evens > 2) so that we need 
 * not maintain a list of exactly which ones are prime.  the space savings here 
 * are well worth the extra time spent dividing; this algorithm is really slow 
 * anyhow.
 *
 * don't forget to link in libgmp.a with this (probably just use -lgmp if you 
 * have installed the gmp library in a normal place).
 */

int main(int argc, char **argv)
{
	mpz_t composite, bound, iter, mod;
	int prime_count = 0; /* count the number of primes for fun */
	if (argc != 2)
	{
		printf("usage: %s composite\n", argv[0]);
		exit(-1);
	}

	mpz_init_set_str(composite, argv[1], 10); /* composite is the int to factor */
	mpz_init(bound);
	mpz_sqrt(bound, composite); /* bound is sqrt(composite) + 1*/
	mpz_add_ui(bound, bound, 1);
	mpz_init_set_ui(iter, 3); /* the integers to start trying as factors */
	mpz_init(mod);

	if (mpz_probab_prime_p(composite, 25) == 1) /* make sure !prime */
	{
		printf("number to factor is prime!\n");
		exit(prime_count);
	}
	
	printf("---\nPerforming Trial Division\n---\n");
	mpz_out_str(stdout, 10, composite);
	printf(" = \n");

	while (mpz_mod_ui(mod, composite, 2) == 0)
	{
		prime_count ++;
		printf("\t2\n");
		mpz_div_ui(composite, composite, 2);
	}
		
	for (iter; mpz_cmp(iter, bound) <= 0; mpz_add_ui(iter, iter, 2))
	{
		do {
			mpz_mod(mod, composite, iter);
			if (mpz_cmp_ui(mod, 0) == 0) /* we have a factor */
			{
				prime_count ++;
				printf("\t");
				mpz_out_str(stdout, 10, iter);
				printf("\n");
				mpz_divexact(composite, composite, iter);
			}
		} while (mpz_cmp_ui(mod, 0) == 0);

		if (mpz_probab_prime_p(composite, 25) == 1) /* all that's left is prime */
		{
			prime_count ++;
			printf("\t");
			mpz_out_str(stdout, 10, composite);
			printf("\n---\n");
			printf("Number of primes in composite: %d\n---\n", prime_count);
			exit(prime_count);
		}
	}
}