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

#define base 32768
#define base2 1073741824

/* Basic modular arithmetic in Z/pZ */

int residue(int x, int p)
{ int temp;
 temp = x%p;
 if (temp < 0) temp += p;
 return temp;
       }

int prod(int x,int y, int p) 
     /* This calculates product of x and y in Z/pZ
	only valid for p < 2^23 approx 8 million  */
{
  unsigned long int x1,x2,y1,y2,z1,z2, z3;
  unsigned long int x0,y0, b2;
  if ((abs(x) < base) && (abs(y) < base)) return (x*y)%p;
  /* else */
  x0 = residue(x,p);
  y0 = residue(y,p);
  x2 = x0 >> 15;
  y2 = y0 >> 15;
  x1 = x0 - (x2 << 15);
  y1 = y0 - (y2 << 15);
  z1 = x1*y1%p;
  z2 = ((x1<< 15)%p*y2)%p;
  z2 += ((y1<< 15)%p*x2)%p;
  b2 = base2%p;
  z3 = (b2*x2%p)*y2%p;
  return (z1+z2+z3)%p;
 }

int square(int x, int p)
{return prod(x,x, p);}


int power(int x, int N, int p)
{int pow = x;
 int i;
 if (x == 0) return 0;
 if (N == 1) return x;
 if (N == 2) return prod(x,x,p);
 if (N == 3) return prod(x, prod(x,x,p),p);
 if (N == 4) return square(square(x,p),p);
 if (N == 5) return prod(square(square(x,p),p), x,p);
 for(i=1; i < N; i++)
   {pow = prod(x,pow,p);}
 return pow;
}


int num_roots_of_unity(int N, int p )
{int count =1;
 int x;
 for (x=2; x<p ; x++)
   if (power(x,N,p)== 1) 
     {count++;
     if (count == N) break;}
 return count;
}

int inverse(int x, int p) 
     /* calculate inverse using euclidean algorithm */
{int b,r;
 int val=1;
 r =x;
 while(r != 1)
   {if(r == 0) 
     return -1; /* No inverse */
   b = p/r;
   val *= -b;
   r = p%r;
}
 val = residue(x,p);
 return val;}

int order(int x, int p)
{int val, pow;
 pow = residue(x,p);
 val = 1;
 while(pow != 1)
   {
     if ( pow == 0) return -1; /* error */
     pow = prod(pow, x,p);
     val++;}
 return val;
}




