#include <stdio.h>
#include <limits.h>

// http://mpir.org/ 

#include "mpirxx.h"


unsigned primes[] = { 
  3,5,7, // up to 2^32, skiped...
};

size_t primesSize = sizeof(primes)/sizeof(*primes);

struct pellSolution { int n, y, x; }; //solution of x^2-n*y^2=1
pellSolution pellSolutions[] = {
 {2,2,3},
 {3,1,2},
 {5,4,9},
 {6,2,5},
 {7,3,8},
 {8,1,3},
};

size_t pellSolutionsSize = sizeof(pellSolutions)/sizeof(pellSolution);

// вход - решение ax^2 - n * ay^2  == 1 по модулю mod, "возводится в степень" mul,
// выход - новое решение bx^2 - n * by == 1 по модулю mod.
void pair_mul(mpz_class &bx, mpz_class &by, mpz_class ax, mpz_class ay, 
   unsigned long long n, unsigned long long mul, const mpz_class & mod)
{
  mpz_class t;
  bx = 1;
  by = 0;

  for(;;) 
  {
    if (mul & 1)
    {
      t = (ax*bx + ay*by*n) % mod;
      by = (ax*by + ay*bx) % mod;
      bx = t;
    }
    if (0 == (mul >>= 1)) break;
    t = (ax*ax + ay*ay*n) % mod;
    ay = ax*ay*2 % mod;
    ax = t;
  }
}


int main(int argc, char* argv[])
{
 mpz_class p;
 mpz_class a, b, a1, b1, t;
 unsigned long long i, j, k, n;

 if (argc<3) 
 {
	 return 1;
	 //p = "59649589127497219";
	 //k = 3; //тест
 }
 else
 {
	 p = argv[1];
	 k = strtoul(argv[2],NULL,10);
 }

 for (j = 0; j < pellSolutionsSize; ++j)
 {
	if (pellSolutions[j].n == k)
	{
		a = pellSolutions[j].x;
		b = pellSolutions[j].y;
		break;
	}
 }
 if (j >= pellSolutionsSize) return 1;
 //first stage
 //2^64
 pair_mul(a1,b1,a,b,k,0x100000000ui64,p);
 pair_mul(a,b,a1,b1,k,0x100000000ui64,p);

 mpz_gcd(t.get_mpz_t(), b.get_mpz_t(), p.get_mpz_t());
 if (t != 1) {
	mpz_out_str (NULL, 10, t.get_mpz_t());
	getchar();
	return 0;
 }

 for (i=0; i<primesSize; ++i) {
	n = primes[i];
	j = ULLONG_MAX/primes[i];
	while (n < j) n *= primes[i];
	pair_mul(a1,b1,a,b,k,n,p);
	if (b1 == 0) { 
		for(;;)
		{
		   pair_mul(a1,b1,a,b,k,primes[i],p);
		   if (b1 == 0) { puts("Strange"); break; }
		   mpz_gcd(t.get_mpz_t(), b1.get_mpz_t(), p.get_mpz_t());
		   if (t != 1) {
			printf("On %u step\n", i);
			mpz_out_str (NULL, 10, t.get_mpz_t());
			getchar();
  			return 0;
		   }
		   a = a1;
		   b = b1;
		}
		break;
	}
	mpz_gcd(t.get_mpz_t(), b1.get_mpz_t(), p.get_mpz_t());
    if (t != 1) {
	  printf("On %u step\n", i);
	  mpz_out_str (NULL, 10, t.get_mpz_t());
      getchar();
  	  return 0;
    }
	a = a1;
	b = b1;
 }

   puts("Second stage...");
   b = b*b % p;
   n = 4*k;
   b1 = (n*b+4)*b%p;

   i = 0;
   j = 1;
   t = abs(b-b1);
   while (true)
   {
     if (i == j ){
       mpz_gcd(a.get_mpz_t(), t.get_mpz_t(), p.get_mpz_t());
       if (a != 1) {
         printf("On %u step\n",i);
         mpz_out_str (NULL, 10, a.get_mpz_t());
         break;
       }
       b = b1;
       j = j*2; 
     }
     b1 = (n*b1+4)*b1%p;
	 a = t*abs(b-b1)%p;
	 if (a==0) {
       mpz_gcd(a.get_mpz_t(), t.get_mpz_t(), p.get_mpz_t());
       if (a != 1) {
         printf("On %u step\n",i);
         mpz_out_str (NULL, 10, a.get_mpz_t());
         break;
       } else {
	     printf("Fail on %u step\n",i);
		 break;
	   }
	 }
	 t = a;
     i = i + 1;
   }
 
   getchar();
}