/* a basic program for testing the primality of a number using trial 
 * division
 *
 * Author: Ron Henry - rlh@ciphermonk.net
 * Ciphermonk Research
 *
 * Notes: Assumes an _odd number_ as input. Program has
 *        a maximum runtime of sqrt(n)/2.
 *
 * Variables: n == the input number to be tested for primality
 *            i == the iterator/counter variable 
 *            root == sqrt of n. used to reduce the length of the test
 *            rem == the remainder captured from modulo operations
 */

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

int main(int argc, char *argv[])
{
   if(argc < 2) {
      fprintf(stderr, "Usage: %s <primality test candidate>\n", *argv);
      return 1;
   }

   /* Create and Initialize variables */
   mpz_t i, n, root, rem;
   
   mpz_init(n);
   mpz_init_set_str(i, "3", 10);
   mpz_init(root);
   mpz_init(rem);

   mpz_set_str(n, argv[1], 10); /* Copy user input to 'n' */
  
   /* Test for even number as input. */

   mpz_cdiv_r_ui(rem, n, (unsigned long) 2);

   if((mpz_cmp_ui(rem, (unsigned long) 0)) == 0) {
      printf("%s expects odd numbers as input!\n", argv[0]);
      return 0;
   }
   /* End even number input check */
   

   /* ============== Tests for Primality Start ============== */

   mpz_init(rem);     /* Zero out the remainder again */
   mpz_sqrt(root, n); /* we shouldn't have to test
                         more than sqrt(n)/2 factors
                         to determine primality
                      */

   mpz_cdiv_q_ui(rem, root, (unsigned long) 2); /* borrowing rem for a quick calculation */
   mpz_add_ui(root,root, (unsigned long) 1);    /* increment root by 1 to compensate for
                                                   rounding */

   gmp_printf("\n- Attempting Trial Division for up to %Zd iterations\n", root);

   while(mpz_cmp(i, root) < 0) {
      mpz_cdiv_r(rem, n, i);
      if((mpz_cmp_ui(rem, (unsigned long) 0)) == 0) {
         mpz_cdiv_q(rem, n, i);
         gmp_printf("Factors of %Zd are %Zd and %Zd\n", n, i, rem);
         return 0;
      }
      else
      mpz_add_ui(i, i, (unsigned long) 2); /*divide by the next odd integer*/
   }
   
   gmp_printf("%Zd is prime!\n", n);

   mpz_clear(i);
   mpz_clear(n);
   mpz_clear(root);
   mpz_clear(rem);

  return 0;
}