/******************************************************************************/ /* */ /* 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); }