reductio.c

calc大数库

				        FREEMPI(TEMP);
					TEMP = MULT_I(C, 2);
					x = CHECK_UVARRAYS(B, TEMP);
					if(x){
					    class_number=class_number+1;
					    z=PERIOD(DD,B,TEMP);
					    if(z % 2){
					       parity_flag = 1;
					    }
					}
				        FREEMPI(TEMP);
				        FREEMPI(C);
				    }
				    FREEMPI(R);
				}
			}
			FREEMPI(H); 
			FREEMPI(B);
		}
		FREEMPI(A);
		FREEMPI(I);
		FREEMPI(J);
	}
	if(e==2){
		d = d / 4;
	}
	if(parity_flag){
             *norm = -1;
        }else{
	     *norm = 1;
        }
        /*printf("h(%lu)=%u, N(eta)=%d\n", d, class_number, *norm);*/
	FREEMPI(ONE);
	FREEMPI(F);
	FREEMPI(G);
	FREEMPI(DD);
	FREEMPIA(GLOBALARRAY_U);
        FREEMPIA(GLOBALARRAY_V);
        return(class_number);
}

USL TABLEPOS(MPI *M, MPI *N){
/* The following gives a table of h(d) for all squarefree d
 * in the range 2<=M<=d<=N<10^6. 
 * The number of squarefree d encountered is returned.
 */
	MPI *DD;
	USL count, m, n, d;
	USI h;
	char buff[20];
	FILE *outfile;
	int norm;

	count = 0;
	m = CONVERTI(M);
	n = CONVERTI(N);

        strcpy(buff, "tablepos.out");
	if(access(buff, R_OK) == 0)
	  unlink(buff);
	outfile = fopen(buff, "a");

	for(d = m; d <= n; d++){
            h = SQUAREFREE_TEST1000(d);
	    if(h == 0){
		continue;
	    }else{
		DD = CHANGE(d);
                h = POS1(DD, &norm);
		FREEMPI(DD);
		count++;
                printf("h(%lu) = %u, N(eta)= %d\n",d, h, norm);
                fprintf(outfile, "h(%lu) = %u, N(eta)= %d\n",d, h, norm);
            }
	}
	if(count==0){
		printf("no squarefree d in the range [%lu,%lu]\n",m,n);
	}
	fclose(outfile);
	return (count);
}

MPI *TABLEPOSX(MPI *M, MPI *N){
	MPI *LIMIT, *COUNT, *ONE;
	USL count;
	int t;

	LIMIT = POWER_I(10, 6);

	if(COMPAREI(M, N) > 0){
		printf("M > N\n");
		FREEMPI(LIMIT);
		return (NULL);
	}
	if(COMPAREI(N, LIMIT) >= 0){
		printf("N >= 10^6\n");
		FREEMPI(LIMIT);
		return (NULL);
	}
	FREEMPI(LIMIT);
	ONE = ONEI();
	t = COMPAREI(M, ONE);
	FREEMPI(ONE);
	if(t <= 0){
		printf("M <= 1\n");
		return (NULL);
	}
	count = TABLEPOS(M, N);
	COUNT = CHANGE(count);
	return(COUNT);
}

USI REDUCE_NEG0(MPI *A, MPI *B, MPI *C, MPI **AA, MPI **BB, MPI **CC, MPI **alpha, MPI **beta, MPI **gamma, MPI **delta, USI print_flag){
/* This is Gauss' algorithm for reducing a positive definite binary
 * quadratic form (A,B,C). 
 * See L.E. Dickson, Introduction to the theory of numbers,
 * page 69. Here d=B^2-4AC <0, A>0,C>0, while the reduced form (AA,BB,CC)
 * satisfies -AA<BB<=AA, CC>=AA, with BB>=0 if CC=AA.
 * The number of steps taken in the reduction is returned.
 */
	MPI *D, *X, *Y, *U, *V, *TEMP1, *TEMP2, *TEMP3;
	MPI *a, *b, *c, *TEMPX, *TEMPY, *TEMPU, *TEMPV, *MINUSD, *DELTA;
	MPI *TEMPB;
	USI i;
	int r, s, t;

	i = 0;
	TEMP1 = MULTI(B, B);
	TEMP2 = MULTI(A, C);
	TEMP3 = MULT_I(TEMP2, 4);
	D = SUBI(TEMP1, TEMP3);
	FREEMPI(TEMP1);
	FREEMPI(TEMP2);
	FREEMPI(TEMP3);
	TEMP1 = MINUSI(A);
	r = COMPAREI(TEMP1, B);
	FREEMPI(TEMP1);
	s = COMPAREI(B, A);
	t = COMPAREI(A, C);
	if(r < 0 && s <= 0 && t <= 0){
		r = COMPAREI(A, C);
		if(t < 0 || (r == 0 && B->S >= 0)){
                   if(print_flag){
			printf("(");
			PRINTI(A);
			printf(",");
			PRINTI(B);
			printf(",");
			PRINTI(C);
			printf(") is reduced\n");
                	printf("Transforming matrix:1,0,0,1\n");
                    }
                    FREEMPI(D);
                    *alpha = ONEI();
                    *beta = ZEROI();
                    *gamma = ZEROI();
                    *delta = ONEI();
                    *AA = COPYI(A);
                    *BB = COPYI(B);
                    *CC = COPYI(C);
                    return (0);
		}
	}
	X = ONEI();
	Y = ZEROI();
	U = ZEROI();
	V = ONEI();
	MINUSD = MINUSI(D);
	a=COPYI(A);
	b=COPYI(B);
	c=COPYI(C);
	while(1){
		TEMPB = b;
		TEMP1 = MINUSI(b);
		TEMP2 = MULT_I(c, 2);
		b = HALFMOD(TEMP1, TEMP2);
		FREEMPI(TEMP1);
		TEMP1 = ADDI(TEMPB, b);
		FREEMPI(TEMPB);
		TEMP3 = MINUSI(TEMP1);
		FREEMPI(TEMP1);
		DELTA = INT(TEMP3, TEMP2);
		FREEMPI(TEMP2);
		FREEMPI(TEMP3);
		TEMPU = U;
		TEMPX = X;
		TEMPY = Y;
		TEMPV = V;
		X = MINUSI(Y);
		Y = MULTABC(TEMPX, Y, DELTA);
		FREEMPI(TEMPX);
		FREEMPI(TEMPY);
		U = MINUSI(V);
		V = MULTABC(TEMPU, V, DELTA);
		FREEMPI(TEMPU);
		FREEMPI(TEMPV);
		FREEMPI(DELTA);
		TEMP1 = a;
		a = COPYI(c);
		FREEMPI(TEMP1);
		TEMP1 = c;
		TEMP2 = MULTABC(MINUSD, b, b);
		TEMP3 = MULT_I(a, 4);
		c = INT(TEMP2, TEMP3);
		FREEMPI(TEMP1);
		FREEMPI(TEMP2);
		FREEMPI(TEMP3);
		i++;
		if(print_flag){
		   printf("(");
		   PRINTI(a);
		   printf(",");
		   PRINTI(b);
		   printf(",");
		   PRINTI(c);
		   printf(")\n");
		}
		if(COMPAREI(c, a) >= 0){
			break;
		}
	}
	if(COMPAREI(a, c) == 0 && b->S < 0){
		TEMP1 = b;
		b = MINUSI(b);
		FREEMPI(TEMP1);
		TEMPX = X;
		TEMPY = Y;
		X = COPYI(Y);
		Y = MINUSI(TEMPX);
		FREEMPI(TEMPX);
		FREEMPI(TEMPY);
		TEMPU = U;
		TEMPV = V;
		U = COPYI(V);
		V = MINUSI(TEMPU);
		FREEMPI(TEMPU);
		FREEMPI(TEMPV);
	}
	*alpha = COPYI(X);
	*beta = COPYI(Y);
	*gamma = COPYI(U);
	*delta = COPYI(V);
	*AA = COPYI(a);
	*BB = COPYI(b);
	*CC = COPYI(c);
	if(print_flag){
           printf("(");
	   PRINTI(a);
	   printf(",");
	   PRINTI(b);
	   printf(",");
	   PRINTI(c);
	   printf(") is reduced\n");
	   printf("Transforming matrix: ");
	   PRINTI(X);
	   printf(",");
	   PRINTI(Y);
	   printf(",");
	   PRINTI(U);
	   printf(",");
	   PRINTI(V);
	   printf("\n");
	}
        FREEMPI(MINUSD);
	FREEMPI(D);
	FREEMPI(X);
	FREEMPI(Y);
	FREEMPI(U);
	FREEMPI(V);
	FREEMPI(a);
	FREEMPI(b);
	FREEMPI(c);
	return(i);
}

USI REP_DEFINITE(MPI *A, MPI *B, MPI *C, MPI *M, USI print_flag){
/* Given a positive definite binary quadratic form Ax^2+Bxy+Cy^2, 
 * we use an algorithm of Gauss to determine if a given positive integer M
 * is expressible as M = AX^2+BXY+CY^2, X and Y integers, gcd(X,Y) = 1.
 * Note: Here D = B^2 - 4AC < 0, A > 0, C > 0.
 * See L.E. Dickson, Introduction to the theory of numbers, pages 74-75.
 */
   USI solution_number;
   MPI *AA1, *BB1, *CC1, *R1, *S1, *T1, *U1, *D;
   MPI *AA2, *BB2, *CC2, *R2, *S2, *T2, *U2;
   MPI *TEMP0, *TEMP1, *TEMP2, *FOURM, *S, *MODULUS, *TWOM;
   MPI *L, *X, *Y, *XX, *YY, *T, *U;
   MPI *ALPHA1, *BETA1, *GAMMA1, *DELTA1, *N, *M0, *TEMP3, *TEMP;
   MPIA SOLUTION, SOL;
   USI l, i, r, s, t, numbr, m0, j, k;
   FILE *outfile;
   char buff[20];

   solution_number = 0;
   REDUCE_NEG0(A, B, C, &AA1, &BB1, &CC1, &R1, &S1, &T1, &U1, 0);
   TEMP0 = MULTI(B, B);
   TEMP1 = MULTI(A, C);
   TEMP2 = MULT_I(TEMP1, 4);
   D = SUBI(TEMP0, TEMP2);
   FREEMPI(TEMP0);
   FREEMPI(TEMP1);
   FREEMPI(TEMP2);
   FOURM = MULT_I(M, 4);
/* First solve x^2=D (mod FOURM). */ 
	S = SQROOT(D, FOURM, &SOLUTION, &MODULUS, &l);
	if(EQMINUSONEI(S)){ /* cleanup of nettbytes */
		FREEMPI(S);
		FREEMPIA(SOLUTION);
		FREEMPI(MODULUS);
		printf("x^2=");
		PRINTI(D);
		FREEMPI(D);
		printf("(mod ");
		PRINTI(FOURM);
		FREEMPI(FOURM);
		printf(") not soluble.\n");
		printf("Hence no primitive solutions.\n");
		return(0);
	}
	FREEMPI(S);
	M0 = INT0(FOURM, MODULUS);
	if(M0->D > 1){
		execerror("M0 >= 2^32", "");
	}
	m0 = CONVERTI(M0); /* valid, as 0 < M0 < 2^32 */
	TEMP = MULT_I(M0, l);
	FREEMPI(M0);
	TEMP1 = INT0_(TEMP, 100);
	FREEMPI(TEMP);
	t = TEMP1->D;
	FREEMPI(TEMP1);
	if(t){
           execerror("number of positive solutions of X^2=N(mod FOURM) > 100", "");
	}
	numbr = 0;
        SOL = BUILDMPIA(); 
	ADD_TO_MPIA(SOL, NULL, 0);
	for(j=0;j<l;j++){
	   for(k=0;k<m0;k++){
		   TEMP1=MULT_I(MODULUS, k);
		   TEMP2=ADD0I(SOLUTION->A[j],TEMP1);
		   TEMP3 = MULT_I(TEMP2, 2);
		   if(COMPAREI(TEMP3, FOURM) <= 0){
                      ADD_TO_MPIA(SOL, TEMP2, numbr);
                      numbr++;
		   }else{
		      TEMP0 = SUB0I(FOURM, TEMP2);
                      ADD_TO_MPIA(SOL, TEMP0, numbr);
	              FREEMPI(TEMP0);
                      numbr++;
		   }
	           FREEMPI(TEMP2);
	           FREEMPI(TEMP3);
	           FREEMPI(TEMP1);
	   }
           
        }
	FREEMPI(MODULUS);
	FREEMPIA(SOLUTION);

	strcpy(buff, "repdefinite.out");
	outfile = fopen(buff, "w");

	TWOM = MULT_I(M, 2);
	for(i = 0; i < numbr; i++){
           N = COPYI(SOL->A[i]);
	   /*if(RSV(N, TWOM) > 0){
		   TEMP = N;
		   N = SUBI(FOURM, N);
		   FREEMPI(TEMP);
	   }*/
	   fprintf(outfile, "N = ");
	   FPRINTI(outfile, N);
	   fprintf(outfile, "\n");
	   if(EQUALI(N, TWOM)){
              FREEMPI(N);
              continue;
           }
	   /* now calculate L, where N^2-4ML=D */
	   TEMP1 = MULTI(N, N);
	   TEMP2 = SUBI(TEMP1, D);
	   L = INT0(TEMP2, FOURM);
	   FREEMPI(TEMP1);
	   FREEMPI(TEMP2);
           REDUCE_NEG0(M, N, L, &AA2, &BB2, &CC2, &R2, &S2, &T2, &U2, 0);
           r = EQUALI(AA1, AA2);
	   s = EQUALI(BB1, BB2);
           t = EQUALI(CC1, CC2);
	   FREEMPI(AA2);
	   FREEMPI(BB2);
	   FREEMPI(CC2);
	   if(r && s && t){
              TEMP1 = MULTI(R1, U2);
              TEMP2 = MULTI(S1, T2);
              ALPHA1 = SUBI(TEMP1, TEMP2);
	      FREEMPI(TEMP1);
	      FREEMPI(TEMP2);

              TEMP1 = MULTI(S1, R2);
              TEMP2 = MULTI(S2, R1);
              BETA1 = SUBI(TEMP1, TEMP2);
	      FREEMPI(TEMP1);
	      FREEMPI(TEMP2);

              TEMP1 = MULTI(U2, T1);
              TEMP2 = MULTI(U1, T2);
              GAMMA1 = SUBI(TEMP1, TEMP2);
	      FREEMPI(TEMP1);
	      FREEMPI(TEMP2);

              TEMP1 = MULTI(U1, R2);
              TEMP2 = MULTI(S2, T1);
              DELTA1 = SUBI(TEMP1, TEMP2);
	      FREEMPI(TEMP1);
	      FREEMPI(TEMP2);

	      FREEMPI(R2);
	      FREEMPI(S2);
	      FREEMPI(T2);
	      FREEMPI(U2);
	      if(print_flag){
                fprintf(outfile, "The unimodular transformation\n");
                fprintf(outfile, "x = ");
		FPRINTI(outfile, ALPHA1);
                fprintf(outfile, "X + ");
		FPRINTI(outfile, BETA1);
                fprintf(outfile, "Y\n");
                fprintf(outfile, "y = ");
		FPRINTI(outfile, GAMMA1);
                fprintf(outfile, "X + ");
		FPRINTI(outfile, DELTA1);
                fprintf(outfile, "Y\n");
                fprintf(outfile, "converts (");
		FPRINTI(outfile, A);
		fprintf(outfile, ",");
		FPRINTI(outfile, B);
		fprintf(outfile, ",");
		FPRINTI(outfile, C);
                fprintf(outfile, ") to ");
                fprintf(outfile, "(");
		FPRINTI(outfile, M);
		fprintf(outfile, ",");
		FPRINTI(outfile, N);
		fprintf(outfile, ",");
		FPRINTI(outfile, L);
		fprintf(outfile, ")\n");
	      }
	      FREEMPI(N);
	      FREEMPI(L);
	      fprintf(outfile, "solution: (");
	      FPRINTI(outfile, ALPHA1);
	      fprintf(outfile, ",");
	      FPRINTI(outfile, GAMMA1);
	      fprintf(outfile, ")\n");

	      X = MINUSI(ALPHA1);
	      Y = MINUSI(GAMMA1);
	      fprintf(outfile, "solution: (");
	      FPRINTI(outfile, X);
	      fprintf(outfile, ",");
	      FPRINTI(outfile, Y);
	      fprintf(outfile, ")\n");
	      solution_number = solution_number + 2;

	      if(D->D == 0 && D->S == -1 && D->V[0] == 4){ /* D = -4 */
                 T = ZEROI();
                 U = ONEI();
		 AUTOMORPH(A, B, D, T, U, X, Y, &XX, &YY);
	         fprintf(outfile, "solution: (");
	         FPRINTI(outfile, XX);
	         fprintf(outfile, ",");
	         FPRINTI(outfile, YY);
	         fprintf(outfile, ")\n");
		 FREEMPI(XX);
		 FREEMPI(YY);
		 AUTOMORPH(A, B, D, T, U, ALPHA1, GAMMA1, &XX, &YY);
	         fprintf(outfile, "solution: (");
	         FPRINTI(outfile, XX);
	         fprintf(outfile, ",");
	         FPRINTI(outfile, YY);
	         fprintf(outfile, ")\n");
		 FREEMPI(XX);
		 FREEMPI(YY);
		 FREEMPI(T);
		 FREEMPI(U);
	         solution_number = solution_number + 2;
              }
	      if(D->D == 0 && D->S == -1 && D->V[0] == 3){ /* D = -3 */
                 T = ONEI();
                 U = ONEI();
		 AUTOMORPH(A, B, D, T, U, X, Y, &XX, &YY);
	         fprintf(outfile, "solution: (");
	         FPRINTI(outfile, XX);
	         fprintf(outfile, ",");
	         FPRINTI(outfile, YY);
	         fprintf(outfile, ")\n");
		 FREEMPI(XX);
		 FREEMPI(YY);
		 AUTOMORPH(A, B, D, T, U, ALPHA1, GAMMA1, &XX, &YY);
		 FREEMPI(T);
	         fprintf(outfile, "solution: (");
	         FPRINTI(outfile, XX);
	         fprintf(outfile, ",");
	         FPRINTI(outfile, YY);
	         fprintf(outfile, ")\n");
		 FREEMPI(XX);
		 FREEMPI(YY);
                 T = MINUS_ONEI();
		 AUTOMORPH(A, B, D, T, U, X, Y, &XX, &YY);
	         fprintf(outfile, "solution: (");
	         FPRINTI(outfile, XX);
	         fprintf(outfile, ",");
	         FPRINTI(outfile, YY);
	         fprintf(outfile, ")\n");
		 FREEMPI(XX);
		 FREEMPI(YY);
		 AUTOMORPH(A, B, D, T, U, ALPHA1, GAMMA1, &XX, &YY);
	         fprintf(outfile, "solution: (");
	         FPRINTI(outfile, XX);
	         fprintf(outfile, ",");
	         FPRINTI(outfile, YY);
	         fprintf(outfile, ")\n");
                 FREEMPI(XX);
	         FREEMPI(YY);
		 FREEMPI(T);
		 FREEMPI(U);
		 solution_number = solution_number + 4;
	      }
              FREEMPI(X);
	      FREEMPI(Y);
	      fprintf(outfile, "\n");
	      FREEMPI(ALPHA1);
	      FREEMPI(BETA1);
	      FREEMPI(GAMMA1);
	      FREEMPI(DELTA1);
	   }else{
		continue;
	   }
        }
	FREEMPI(R1);
	FREEMPI(S1);
        FREEMPI(T1);
        FREEMPI(U1);
	FREEMPI(AA1);
	FREEMPI(BB1);
        FREEMPI(CC1);
	FREEMPI(TWOM);
	FREEMPI(FOURM);
	FREEMPI(D);
	FREEMPIA(SOL);
	fclose(outfile);
        return(solution_number);
}

void AUTOMORPH(MPI *A, MPI *B, MPI *D, MPI *T, MPI *U, MPI *X, MPI *Y, MPI **XX, MPI **YY){
	/* From Loo-Keng Hua, Introduction to Number Theory, 
	 * Theorem 4.2, pages 279-281 */
	MPI *TEMP0, *TEMP1, *TEMP2, *TEMP3, *TEMP4, *TEMP5, *TEMP6;
	MPI *T1, *T2, *T3, *T4, *T5, *T6, *TEMP;

	TEMP1 = MULT_I(A, 2);
	TEMP2 = MULTI(TEMP1, X);
	TEMP0 = MULTI(Y, B);
	TEMP3 = ADDI(TEMP2, TEMP0);
	TEMP4 = MULTI(TEMP3, U);
	TEMP5 = MULTI(Y, T);
	TEMP6 = ADDI(TEMP4, TEMP5);
	*YY = INT_(TEMP6, 2);

	T1 = MULTI(TEMP3, T);
	T2 = MULTI(D, Y);
	TEMP = T2;
	T2 = MULTI(T2, U);
	FREEMPI(TEMP);
	T3 = ADDI(T1, T2);
	TEMP = T3;
	T3 = INT_(T3, 2);
	FREEMPI(TEMP);
	T4 = MULTI(B, *YY);
	T5 = SUBI(T3, T4);
	T6 = MULT_I(A, 2);
	*XX = INT(T5, T6);

	FREEMPI(TEMP0);
	FREEMPI(TEMP1);
	FREEMPI(TEMP2);
	FREEMPI(TEMP3);
	FREEMPI(TEMP4);
	FREEMPI(TEMP5);
	FREEMPI(TEMP6);
	FREEMPI(T1);
	FREEMPI(T2);
	FREEMPI(T3);
	FREEMPI(T4);
	FREEMPI(T5);
	FREEMPI(T6);
}

MPI *REP_DEFINITEX(MPI *A, MPI *B, MPI *C, MPI *M, MPI *PRINT_FLAG){
	MPI *D, *TEMP1, *TEMP2, *TEMP3;
	USI i, print_flag;

	if(M->S <= 0){
		printf("M <= 0\n");
		return(NULL);
	}
	if(A->S <= 0){
		printf("A <= 0\n");
		return(NULL);
	}
	if(C->S <= 0){
		printf("C <= 0\n");
		return(NULL);
	}
	TEMP1 = MULTI(B, B);
	TEMP2 = MULTI(A, C);
	TEMP3 = MULT_I(TEMP2, 4);
	D = SUBI(TEMP1, TEMP3);
	FREEMPI(TEMP1);
	FREEMPI(TEMP2);
	FREEMPI(TEMP3);
	if(D->S >= 0){
		printf("B^2-4*A*C >= 0\n");
		FREEMPI(D);
		return(NULL);
	}
	print_flag = (USI)CONVERTI(PRINT_FLAG);
	if(print_flag != 0 && print_flag != 1){
		printf("print_flag != 0 or 1\n");
		FREEMPI(D);
		return(NULL);
	}
	i = REP_DEFINITE(A, B, C, M, print_flag);
	FREEMPI(D);
	return(CHANGE(i));
}