/* test_mpz.c   basic functions                                 */
/*           compile:  gcc -o test_mpz -I. test_mpz.c -L. -lgmp */
/*           assumes gmp.h and libgmp.a in local directory      */
#include <stdio.h>
#include "gmp.h"

int miller_rabin(mpz_t n, int s); /* 0 may be prime, 1 is composite */
int witness(mpz_t a, mpz_t n);    /* 0 for false, 1 for true */

int main()
{
  mpz_t a, b, c, d;
  char my_str[1000];

  mpz_init(a); /* must init every variable */
  mpz_init(b);
  mpz_init(c);
  mpz_init(d);

  mpz_set_str(a, "100000000000000000000", 10); /* radix 10 */
  mpz_set_str(b,"1", 10);

  mpz_sub(c, a, b);
  mpz_get_str(my_str, 10, c);
  printf("my_str=%s \n", my_str);

  mpz_mul(b, c, c);
  mpz_get_str(my_str, 10, b);
  printf("b=c*c=%s \n", my_str);

  mpz_add(a, b, b);
  mpz_div(d, a, c);
  mpz_get_str(my_str, 10, d);
  printf("d=(a=b+b)/c=%s \n", my_str);

  mpz_gcd(d, a, c);
  mpz_get_str(my_str, 10, d);
  printf("d=gcd(a,c)=%s \n", my_str);

  mpz_get_str(my_str, 2, d);
  printf("d_binary=%s \n", my_str);



  return 0;
}

int miller_rabin(mpz_t n, int s) /* 0 may be prime, 1 is composite */
{
  int i;
  mpz_t nm1;
  mpz_t one;
  mpz_t a;

  mpz_init(nm1);
  mpz_init(one);
  mpz_init(a);
  mpz_set_str(one, "1", 10);
  mpz_sub(nm1, n, one);
  for(i=0; i<s; i++)
  {
    mpz_random(a, 10);
    if( witness(a, n) ) return 1; /* composite number */
  }
  return 0; /* probably a prime, error estimate 2^(-s) */
}

int witness(mpz_t a, mpz_t n)    /* 0 for false, 1 for true */
{
  mpz_t d;
  mpz_t dd;
  mpz_t x;
  char b[100000];
  int i;
  mpz_t nm1;
  mpz_t one;

  mpz_init(nm1);
  mpz_init(one);
  mpz_set_str(one, "1", 10);
  mpz_sub(nm1, n, one);
  mpz_get_str(b, 2, nm1); /* n-1 in binary characters */

  mpz_init(d);
  mpz_init(dd);
  mpz_init(x);
  mpz_set_str(d, "1", 10);

  i=0;
  while(b[i] != '\0')
  {
    mpz_set(x, d);
    mpz_mul(dd, d, d);
    mpz_mod(d, dd, n);


    if( b[i] == '1' )
    {
      mpz_mul(dd, d, a);
      mpz_mod(d, dd, n);
    }
  }

  if( mpz_cmp(d, one) ) return 1;
  return 0;
}