/*** This program finds a prime number by using Euler's criterion
     10 times  ***/
int a,i,n,x,y,w;
int primality;
int count=0;
/** This computes x^e **/
exp(int x, int e){
  int y; int i,m; int b[100];
  i=0;
  while(e>0){ // This part gets binary expansion of e into array b
    if((e/2)*2==e) b[i]=0; else b[i]=1;
    e=e/2; i++;
  }
  /** The following computes x^e by repeated squaring **/
  m=i-1; y=1; 
  for(i=m;i>=0;i--){y=y*y; y=y%n; if(b[i]==1){y=y*x; y=y%n;}}
  return y;
}
int J(int a, int n){ // Jaqcobi function
  if(a==0)return 0;
  else if(a==1)return 1;
  else if((a/2)*2==a)return J(a/2, n)*exp(-1, (n*n-1)/8);
  else return J(n % a, a)*exp(-1, (a-1)*(n-1)/4);
}
main(){
  n=random()%10000;  // Choose an initial candidate for prime
  if((n/2)*2==n)n++; // If n is even, add 1
  printf("n %d",n); getchar();
  do{
    primality=1; //Suppose n is prime tentatively
    count++;
    for(i=1; i<=10; i++){
      a=random()%n;
      x=J(a,n);
      if(x!=0){
        y=exp(a, (n-1)/2);
//        printf("a %d n %d jacobi x %d exp y %d", a, n, x, y); getchar();
//        printf("(n-1)/2  %d ", (n-1)/2); getchar();
        if(y==n-1)y=-1;
        if(x!=y){primality=0; break;} // n is not prime
      } else {primality=0; break;}
    } 
    /** take the next odd number **/
   n=n+2;
  printf("n, count %d %d \n",n-2,count); getchar();
  }
  while(primality==0);
  printf("n, count %d %d \n",n-2,count);
}