/******************************************************************************/
/*									     */
/*  CHECK RECIPROCITY							     */
/*  06/16/06 (dkc)							     */
/*									     */
/*  This C program checks cubic reciprocity.				     */
/*									     */
/*****************************************************************************/
#include <stdio.h>
int main ()
{
//
//  "Gaussian sum" cubed
//
int input[97*3]={
 7, 1, 3, 
 13, 4, 3, 
 19, -2, 3, 
 31, -5, -6, 
 37, 4, -3, 
 43, 1, -6, 
 61, 4, 9, 
 67, -2, -9, 
 73, 1, 9, 
 79, 10, 3, 
 97, -8, 3, 
 103, -11, -9, 
 109, 7, 12, 
 127, 7, -6, 
 139, 10, -3, 
 151, -14, -9, 
 157, 13, 12, 
 163, -14, -3, 
 181, 4, 15, 
 193, 16, 9, 
 199, -2, -15, 
 211, -14, -15, 
 223, -17, -6, 
 229, -5, 12, 
 241, 16, 15, 
 271, 10, -9, 
 277, 7, -12, 
 283, 13, -6, 
 307, -17, -18, 
 313, 16, -3, 
 331, 10, 21, 
 337, -8, -21, 
 349, -20, -3, 
 367, 22, 9, 
 373, 4, 21, 
 379, 22, 15, 
 397, -23, -12, 
 409, -8, 15, 
 421, -20, -21, 
 433, -11, -24, 
 439, -23, -18, 
 457, -17, -24, 
 463, 22, 21, 
 487, -2, 21, 
 499, 7, -18, 
 523, -26, -9, 
 541, 4, -21, 
 547, -14, -27, 
 571, -5, 21, 
 577, -8, -27, 
 601, 1, -24, 
 607, -23, 3, 
 613, 28, 9, 
 619, 22, 27, 
 631, -29, -15, 
 643, -11, 18, 
 661, -20, 9, 
 673, -29, -21, 
 691, 19, 30, 
 709, 28, 3, 
 727, 13, -18, 
 733, 31, 12, 
 739, 7, 30, 
 751, 31, 21, 
 757, 28, 27, 
 769, -17, 15, 
 787, -2, 27, 
 811, 25, -6, 
 823, -14, -33, 
 829, -20, -33, 
 853, 4, -27, 
 859, 10, 33, 
 877, 28, -3, 
 883, 34, 21, 
 907, -26, -33, 
 919, -35, -18, 
 937, -32, -3, 
 967, 7, -27, 
 991, -26, 9, 
 997, -23, -36, 
 1009, -8, 27, 
 1021, -11, -36, 
 1033, 16, -21, 
 1039, 22, -15, 
 1051, -35, -6, 
 1063, 34, 3, 
 1069, 37, 12, 
 1087, -17, 21, 
 1093, -29, -36, 
 1117, 28, -9, 
 1123, 34, 33, 
 1129, -32, 3, 
 1153, 16, 39, 
 1171, -14, -39, 
 1201, 40, 21, 
 1213, 28, 39, 
 1231, 10, 39};

int count=97;
int A[2],p,q;
int f,g,h,i,j,flag;
int T[2],U[2],temp,t;
FILE *Outfp;
Outfp = fopen("solve3.dat","w");
for (f=0; f<count; f++) {
   p=input[3*f];
   A[0]=input[3*f+1];
   A[1]=input[3*f+2];
   for (g=0; g<count; g++) {
      if (g==f)
	 continue;
      q=input[3*g];
/***************************/
/*  compute p*A modulus q  */
/***************************/
      T[0]=p*A[0];
      T[0]=T[0]-(T[0]/q)*q;
      if (T[0]<0)
	 T[0]+=q;
      T[1]=p*A[1];
      T[1]=T[1]-(T[1]/q)*q;
      if (T[1]<0)
	 T[1]+=q;
/******************/
/*  compute cubes */
/******************/
      for (h=0; h<q; h++) {
	 for (i=0; i<q; i++) {
	    U[0]=h;
	    U[1]=i;
	    for (j=1; j<3; j++) {
	       t=U[1]*i;
	       temp=U[0]*h-t;
	       U[1]=U[0]*i+U[1]*h-t;
	       U[0]=temp;
	       }
	    U[0]=U[0]-(U[0]/q)*q;
	    U[1]=U[1]-(U[1]/q)*q;
	    if (U[0]<0)
	       U[0]+=q;
	    if (U[1]<0)
	       U[1]+=q;
	    flag=0;
	    if ((U[0]==T[0])&&(U[1]==T[1])) {
	       fprintf(Outfp,"solution 1: %d %d \n",p,q);
	       printf("solution 1: %d %d \n",p,q);
	       flag=1;
	       goto askip;
	       }
	    }
	 }
/******************************/
/*  compute q**(p-1)/3 mod p  */
/******************************/
askip:temp=1;
      for (h=0; h<(p-1)/3; h++) {
	 temp=temp*q-((temp*q)/p)*p;
	 }
      if (temp==1) {
	 if (flag==1) {
	    fprintf(Outfp,"solution 2: %d %d \n",p,q);
	    printf("solution 2: %d %d \n",p,q);
	    printf("\n");
	    }
	 else {
	    fprintf(Outfp,"error: %d %d \n",p,q);
	    printf("error: %d %d \n",p,q);
	    }
	 }
      else {
	 if (flag==1) {
	    fprintf(Outfp,"error: %d %d \n",p,q);
	    printf("error: %d %d \n",p,q);
	    }
	 }
      }
   }
fclose(Outfp);
return(0);
}