lagrange.c

calc大数库

	int tt, ttt, u, n1 = 0;
	FILE *outfile;
        char buff[20];

	TWOA = MULT_I(A, 2);
        TWOC = MULT_I(C, 2);

/* first solve the congruence a\theta^2+b\theta+c=0 (mod |N|) */
	L = QUADRATIC(A,B,C,N,&SOL);
	if(L->S == 0){
		FREEMPI(L);
		FREEMPIA(SOL);
		execerror("There are no primitive solutions", "");
	}
	FSX = BUILDMPIA();
	FSY = BUILDMPIA();

	strcpy(buff, "binform.out");
	outfile = fopen(buff, "a");
 /* Let \Delta = B[1]^2-ac if B=2*B[1], else D = B^2-4AC.     */
	TMP1 = MULTI(B, B);
	FOUR = CHANGEI(4);
	TMP2 = MULTI3(FOUR, A, C);
	FREEMPI(FOUR);
	D = SUBI(TMP1, TMP2);
	FREEMPI(TMP1);
	FREEMPI(TMP2);
	t = (B->V[0]) % 2;
	if(t){
		printf("D=");
		fprintf(outfile, "D=");
		PRINTI(D);
		FPRINTI(outfile, D);
		printf("\n");
		fprintf(outfile, "\n");
		if(((D->D == 0) && (D->V[0] == 5)) && (((A->S) * (N ->S)) == -1))
				flag = 1;
	}
	else{
		TMP1 = D;
		D = INT_(D, 4);
		FREEMPI(TMP1);
		printf("Delta=");
		fprintf(outfile, "Delta=");
		PRINTI(D);
		FPRINTI(outfile, D);
		printf("\n");
		fprintf(outfile, "\n");
	}
	ABSN = ABSI(N);
	TWOABSN = MULT_I(ABSN, 2);
	tt = N->S;
	Q = MULTI(A, ABSN);
	TWOQ = MULT_I(Q, 2);
	MINUS2Q = MINUSI(TWOQ);
	MINUSQ = MINUSI(Q);
	ONE = ONEI();
	l = CONVERTI(L);
	FREEMPI(L);

	for (r = 0; r < l; r++){
		THETA = COPYI(SOL->A[r]);
		if(verbose){
			printf("THETA[%u]=", r);
			fprintf(outfile, "THETA[%u]=", r);
			PRINTI(THETA);
			FPRINTI(outfile, THETA);
			printf("\n");
			fprintf(outfile, "\n");
		}
		TMP1 = MULT_I(THETA, 2);
		TMP2 = MULTI(TMP1, A);
		FREEMPI(TMP1);
		n = ADDI(TMP2, B);
		FREEMPI(TMP2);
		P = INT_(n, 2);
		if(verbose){
			printf("n=");
			fprintf(outfile, "n=");
			PRINTI(n);
			FPRINTI(outfile, n);
			printf("\n");
			fprintf(outfile, "\n");
		}
		if (t == 0){
			P01 = MINUSI(P);
			P02 = COPYI(P);
			Q01 = COPYI(Q);
			Q02 = COPYI(MINUSQ);
		}
		else{
			P01 = MINUSI(n);
			P02 = COPYI(n);
			Q01 = COPYI(TWOQ);
			Q02 = COPYI(MINUS2Q);
		}
		period_length = SURD(D, ONE, P01, Q01, &AA1, &U1, &V1, &P1, &Q1, 1);
		period_length = SURD(D, ONE, P02, Q02, &AA2, &U2, &V2, &P2, &Q2, 1);
		if(period_length %2)
                      period_length = 2 * period_length;
		r1 = AA1->size;
		r2 = AA2->size;
		if(flag){
			s = r1 - period_length;
			if (s > 1){
				TMP1 = SUBI(P1->A[s - 1], P1->A[s - 2]);
				TMP2 = SUBI(Q1->A[s - 1], Q1->A[s - 2]);
				TMP3 = MULTI(ABSN, TMP1);
				TMP4 = MULTI(THETA, TMP2);
				X = ADDI(TMP3, TMP4);
				FREEMPI(TMP1);
				FREEMPI(TMP3);
				FREEMPI(TMP4);
			}
			else{
				TMP1 = ONEI();
				X = SUBI(P1->A[0], TMP1);
				FREEMPI(TMP1);
				TMP2 = COPYI(Q1->A[0]);
			}	
/* print from here */
		if(verbose){
			printf("Exceptional solution: s=%u, (X,Y) = (", s);
			fprintf(outfile, "Exceptional solution: s=%u, (X,Y) = (", s);
			PRINTI(X);
			FPRINTI(outfile, X);
			printf(", ");
			fprintf(outfile, ", ");
			PRINTI(TMP2);
			FPRINTI(outfile, TMP2);
			printf(")\n");
			fprintf(outfile, ")\n");
			GetReturn();
		}
/* to from here */
			ADD_TO_MPIA(FSX, X, n1);
			ADD_TO_MPIA(FSY, TMP2, n1);
			FLAGE = 1;
			FREEMPI(X);
			FREEMPI(TMP2);
		}

		for( k = 1; k < r1; k++){
			ttt = (k % 2) ? -1 : 1;
			d1 = (V1->A[k])->D;
			v1 = (V1->A[k])->V[0];
	/*printf("V1->A[%u]=", k); PRINTI(V1->A[k]);printf("\n"); */
	if((!t && ((d1 == 0) && (v1 == 1))) ||  (t && ((d1  == 0) && (v1 == 2)))){
		if(verbose){
			printf("processing omega:\n");
			fprintf(outfile, "processing omega:\n");
			printf("rcf:(");
			fprintf(outfile, "rcf:(");
			PRINTI(P01);
			FPRINTI(outfile, P01);
			printf("+sqrt(");
			fprintf(outfile, "+sqrt(");
			PRINTI(D);
			FPRINTI(outfile, D);
			printf("))/");
			fprintf(outfile, "))/");
			PRINTI(Q01);
			FPRINTI(outfile, Q01);
			printf("\n");
			fprintf(outfile, "\n");
			GetReturn();
		}
			if (!flag && ((V1->A[k])->S == ttt * tt)){/* calculate x=y\theta+|N|X */
				TMP1 = MULTI(ABSN, P1->A[k - 1]);
				TMP2 = MULTI(THETA, Q1->A[k - 1]);
				X = ADDI(TMP1, TMP2);
				FREEMPI(TMP1);
				FREEMPI(TMP2);
				Y = COPYI(Q1->A[k - 1]);
/* print from here */
			if(verbose){
				printf("(X1,Y1) = (");
				fprintf(outfile, "(X1,Y1) = (");
                                PRINTI(X);
                                FPRINTI(outfile, X);
				printf(", ");
				fprintf(outfile, ", ");
                                PRINTI(Y);
                                FPRINTI(outfile, Y);
                                printf(")\n");
                                fprintf(outfile, ")\n");
				GetReturn();
			}
/* to from here */
				ADD_TO_MPIA(FSX, X, n1);
				ADD_TO_MPIA(FSY, Y, n1);
				FREEMPI(X);
				FREEMPI(Y);
				FLAG1 = 1;
				break;
			}
		}
	}
		for( k = 1; k < r2; k++){
			ttt = (k % 2) ? 1 : -1;
			d2 = (V2->A[k])->D	;
			v2 = (V2->A[k])->V[0];
/*	printf("V2->A[%u]=", k); PRINTI(V2->A[k]);printf("\n");*/
	if((!t && ((d2 == 0) && (v2 == 1))) || (t && ((d2 == 0) && (v2  == 2)))){                       
			if(verbose){
				printf("processing omega_star\n");
				fprintf(outfile, "processing omega_star\n");
				printf("rcf:(");
				fprintf(outfile, "rcf:(");
				PRINTI(P02);
				FPRINTI(outfile, P02);
				printf("+sqrt(");
				fprintf(outfile, "+sqrt(");
				PRINTI(D);
				FPRINTI(outfile, D);
				printf("))/");
				fprintf(outfile, "))/");
				PRINTI(Q02);
				FPRINTI(outfile, Q02);
				printf("\n");
				fprintf(outfile, "\n");
				GetReturn();
			}
			if ((V2->A[k])->S == ttt * tt){/* calculate x=y\theta+|N|X */
				TMP1 = MULTI(ABSN, P2->A[k - 1]);
				TMP2 = MULTI(THETA, Q2->A[k - 1]);
				X = ADDI(TMP1, TMP2);
				FREEMPI(TMP1);
				FREEMPI(TMP2);
				Y = COPYI(Q2->A[k - 1]);
/* print from here */
			if(verbose){
				printf("(X2,Y2) = (");
				fprintf(outfile, "(X2,Y2) = (");
                                PRINTI(X);
                                FPRINTI(outfile, X);
				printf(", ");
				fprintf(outfile, ", ");
                                PRINTI(Y);
                                FPRINTI(outfile, Y);
                                printf(")\n");
                                fprintf(outfile, ")\n");
				GetReturn();
			}
/* to from here */
				if(!FLAG1 && !FLAGE){
					ADD_TO_MPIA(FSX, X, n1);
					ADD_TO_MPIA(FSY, Y, n1);
					FLAG2 = 1;
				}
				else{
					u = RSV(FSY->A[n1], Y);
					if(u > 0){
						TMP1 = FSX->A[n1];
						TMP2 = FSY->A[n1];
						FSX->A[n1] = COPYI(X);
						FSY->A[n1] = COPYI(Y);
						FREEMPI(TMP1);
						FREEMPI(TMP2);
					}
				}
				FREEMPI(X);
				FREEMPI(Y);
				break;
			}
		}
	}
	if(FLAGE || FLAG1 || FLAG2){
		n1++;
		FLAGE = 0;
		FLAG1 = 0;
		FLAG2 = 0;
	}
	FREEMPIA(AA1);
	FREEMPIA(AA2);
	FREEMPIA(U1);
	FREEMPIA(U2);
	FREEMPIA(V1);
	FREEMPIA(V2);
	FREEMPIA(P1);
	FREEMPIA(P2);
	FREEMPIA(Q1);
	FREEMPIA(Q2);
	FREEMPI(THETA);
	FREEMPI(P01);
	FREEMPI(P02);
	FREEMPI(Q01);
	FREEMPI(Q02);
	FREEMPI(P);
	FREEMPI(n);
}
      /*TMP1 = MULTI(A, N);
	t = EQONEI(TMP1);
	FREEMPI(TMP1);
	if(t && (n1 == 0)){
		ADD_TO_MPIA(FSX, ONE, n1);
		TMP1 = ZEROI();
		ADD_TO_MPIA(FSY, TMP1, n1);
		FREEMPI(TMP1);
		n1++;
	}*/
if(FLAG){
	if(n1 == 0)
		printf("There are no primitive solutions (x,y)!\n");
	else{
		for (i = 0; i < n1; i++){
			printf("solution (");
			fprintf(outfile, "solution (");
			PRINTI(FSX->A[i]);
			FPRINTI(outfile, FSX->A[i]);
			printf(",");
			fprintf(outfile, ",");
			PRINTI(FSY->A[i]);
			FPRINTI(outfile, FSY->A[i]);
			printf(")");
			fprintf(outfile, ")");
			x = FSX->A[i];
			y = FSY->A[i];
			Z = EUCLIDI(x, y, &DELTA, &BETA);
			FREEMPI(Z);
			BETA->S = -(BETA->S);
			/* x(DELTA)-y(BETA)=1 */
			TMP1 = MULTI3(TWOA, x, BETA);
			TMP2 = MULTI3(TWOC, y, DELTA);
			TMP5 = MULTI3(B, x, DELTA);
			TMP6 = MULTI3(B, BETA, y);
			/* now add TMP1, TMP2, TPM5, TMP6 */
			SUM = ADDI(TMP1, TMP2);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			FREEMPI(BETA);
			FREEMPI(DELTA);
			TMP1 = SUM;
			SUM = ADDI(SUM, TMP5);
			FREEMPI(TMP1);
			FREEMPI(TMP5);
			TMP1 = SUM;
			SUM = ADDI(SUM, TMP6);
			FREEMPI(TMP1);
			FREEMPI(TMP6);
			TMP1 = SUM;
			SUM = MOD(SUM, TWOABSN);
			FREEMPI(TMP1);
			if(RSV(SUM, ABSN) > 0){
				TMP1 = SUM;
				SUM = SUBI(SUM, TWOABSN);
				FREEMPI(TMP1);
			}
                        MODULUS = MULT_I(ABSN, 2);
			printf(": n = ");
			fprintf(outfile, ": n = ");
			PRINTI(SUM);
			FPRINTI(outfile, SUM);
			printf("(mod ");
			PRINTI(MODULUS);
			printf(")\n");
			fprintf(outfile, "(mod ");
			FPRINTI(outfile, MODULUS);
			fprintf(outfile, ")\n");
			FREEMPI(SUM);
			FREEMPI(MODULUS);
			GetReturn();
		}
	}
}
	FREEMPIA(SOL);
	FREEMPI(D);
	FREEMPI(Q);
	FREEMPI(ABSN);
	FREEMPI(TWOABSN);
	FREEMPI(TWOQ);
	FREEMPI(TWOA);
	FREEMPI(TWOC);
	FREEMPI(MINUSQ);
	FREEMPI(MINUS2Q);
	FREEMPI(ONE);
	*N1 = n1;
	*FSX1 = FSX;
	*FSY1 = FSY;
	fclose(outfile);
	return;
}

void BINARYFORM1(MPI *AA, MPI *BB, MPI *CC, MPI *NN, MPI *VERB)
/* We reduce to the case GCD(AA,NN)=1 and then use BINARYFORM.
 * The n of p.3 of paper http://www.maths.uq.edu.au/~krm/binary.pdf
 * is also printed.
 * Verbose output if VERB = 1; otherwise 0;
 * Completed 21/12/00. 
 */ 
{
	MPI *a, *b, *c, *x, *y, *TWOA, *TWOC, *MODULUS;
	MPI *G, *alpha, *beta, *M, *gamma, *delta, *Z;
	MPI *ABSN, *TWOABSN, *BETA, *DELTA, *X, *D, *FOUR;
	MPI *TMP1, *TMP2, *TMP3, *TMP4, *TMP5, *TMP6, *SUM;
	MPI *A, *B, *C, *N;
	USI t, length, i, FLAG = 0, verb;
	MPIA FSX1, FSY1;
	int s;
	FILE *outfile;
        char buff[20];
		strcpy(buff, "binform.out");
		outfile = fopen(buff, "w");

/* Check that NN is non-zero */
	s = NN->S;
	if (s == 0)
		{
		  printf("NN = 0\n");
		  return;
		}
/* Check to see d=gcd(AA,BB,CC) divides NN and if so
 * to replace AA,BB,CC,NN by AA/d,BB/d,CC/d,NN/d.
 */
	TMP1 = GCD(AA,BB);
	TMP2 = GCD(TMP1,CC);
	TMP3 = MOD(NN, TMP2);
	s = TMP3->S;
	FREEMPI(TMP1);
	FREEMPI(TMP3);
	if (s){
	        printf("gcd(A,B,C) does not divide N\n");
	 	FREEMPI(TMP2);
	        return;
	}
	else{
		A = INT(AA, TMP2);
		B = INT(BB, TMP2);
		C = INT(CC, TMP2);
		N = INT(NN, TMP2);
	}
	FREEMPI(TMP2);

/* Check that B^2-4AC > 0 */
	TMP1 = MULTI(B, B);
        FOUR = CHANGEI(4);
        TMP2 = MULTI3(FOUR, A, C);
        FREEMPI(FOUR);
        D = SUBI(TMP1, TMP2);
	s = D->S;
        FREEMPI(TMP1);
        FREEMPI(TMP2);
	if (s <= 0)
	{
	  printf("D <= 0\n");
          FREEMPI(D);
          FREEMPI(A);
          FREEMPI(B);
          FREEMPI(C);
          FREEMPI(N);
	  return;
	}
/* Check that B^2-4AC is not a perfect square */
	X = BIG_MTHROOT(D, 2);
	G = MULTI(X, X);
	t = EQUALI(D, G);
	FREEMPI(G);
	FREEMPI(X);
        FREEMPI(D);
	if (t)
	{
	  printf("D is a perfect square\n");
          FREEMPI(A);
          FREEMPI(B);
          FREEMPI(C);
          FREEMPI(N);
	  return;
	}

	TWOA = MULT_I(A, 2);
	TWOC = MULT_I(C, 2);
	ABSN = ABSI(N);
	TWOABSN = MULT_I(ABSN, 2);
	G = GCD(A,N);
	t = EQONEI(G);
	FREEMPI(G);
	if(t){
		a = COPYI(A);
		b = COPYI(B);
		c = COPYI(C);
		FLAG = 1;
	}
	else{
		GAUSS(A, B, C, N, &alpha, &gamma, &M);
		Z = EUCLIDI(alpha, gamma, &delta, &beta);
		FREEMPI(Z);
		beta->S = -(beta->S);
		a = M;
		strcpy(buff, "binform.out");
		outfile = fopen(buff, "w");

		/* alpha * delta - beta * gamma = 1 */
		if(VERB->S){
			printf("alpha=");PRINTI(alpha);printf("\n");
			fprintf(outfile, "alpha=");FPRINTI(outfile, alpha);fprintf(outfile, "\n");
			printf("beta=");PRINTI(beta);printf("\n");
			fprintf(outfile, "beta=");FPRINTI(outfile, beta);fprintf(outfile, "\n");
			printf("gamma=");PRINTI(gamma);printf("\n");
			fprintf(outfile, "gamma=");FPRINTI(outfile, gamma);fprintf(outfile, "\n");
			printf("delta=");PRINTI(delta);printf("\n");
			fprintf(outfile, "delta=");FPRINTI(outfile, delta);fprintf(outfile, "\n");
			GetReturn();
		}

		/* now calculate a, b, c, when A*x^2+B*x*y+C*y^2 is
		   transformed into a*X^2+b*X*Y+c*Y^2 under the
		   transformation x=alpha*X+beta*Y, y=gamma*X+delta*Y */

		TMP1 = MULTI3(A, beta, beta);
		TMP2 = MULTI3(B, beta, delta);
		TMP3 = MULTI3(C, delta, delta);
		TMP4 = ADDI(TMP1, TMP2);
		c = ADDI(TMP4, TMP3);
		FREEMPI(TMP1);
		FREEMPI(TMP2);
		FREEMPI(TMP3);
		FREEMPI(TMP4);
		TMP1 = MULTI3(TWOA, alpha, beta);
		TMP2 = MULTI3(TWOC, gamma, delta);
		TMP5 = MULTI3(B, alpha, delta);
		TMP6 = MULTI3(B, beta, gamma);
		/* now add TMP1, TMP2, TPM5, TMP6 */
		SUM = ADDI(TMP1, TMP2);
		FREEMPI(TMP1);
		FREEMPI(TMP2);
		TMP1 = SUM;
		SUM = ADDI(SUM, TMP5);
		FREEMPI(TMP1);
		FREEMPI(TMP5);
		TMP1 = SUM;
		SUM = ADDI(SUM, TMP6);
		FREEMPI(TMP1);
		FREEMPI(TMP6);
		b = COPYI(SUM);
		FREEMPI(SUM);
		/* Now we have gcd(a,N)=1 */
	}
	verb = CONVERTI(VERB);
	fclose(outfile);
        BINARYFORM(a, b, c, N, &FSX1, &FSY1, &length, FLAG, verb);
	
	outfile = fopen(buff, "a");
if(!t){
	for(i = 0; i < length; i++){
		printf("i=%u: ",i);
		fprintf(outfile, "i=%u: ",i);
		TMP1 = MULTI(alpha, FSX1->A[i]);
		TMP2 = MULTI(beta, FSY1->A[i]);
		TMP3 = MULTI(gamma, FSX1->A[i]);
		TMP4 = MULTI(delta, FSY1->A[i]);
		x = ADDI(TMP1, TMP2);
		FREEMPI(TMP1);
		FREEMPI(TMP2);
		y = ADDI(TMP3, TMP4);
		FREEMPI(TMP3);
                FREEMPI(TMP4);
		printf("Solution (");
		fprintf(outfile, "Solution (");
		PRINTI(x);
		FPRINTI(outfile, x);
		printf(",");
		fprintf(outfile, ",");
		PRINTI(y);
		FPRINTI(outfile, y);
		printf("): ");
		fprintf(outfile, "): ");
		Z = EUCLIDI(x, y, &DELTA, &BETA);
		FREEMPI(Z);
		BETA->S = -(BETA->S);
		TMP1 = MULTI3(TWOA, x, BETA);
		TMP2 = MULTI3(TWOC, y, DELTA);
		TMP5 = MULTI3(B, x, DELTA);
		TMP6 = MULTI3(B, BETA, y);
		/* now add TMP1, TMP2, TPM5, TMP6 */
		SUM = ADDI(TMP1, TMP2);
		FREEMPI(TMP1);
		FREEMPI(TMP2);
		FREEMPI(BETA);
		FREEMPI(DELTA);
		TMP1 = SUM;
		SUM = ADDI(SUM, TMP5);
		FREEMPI(TMP1);
		FREEMPI(TMP5);
		TMP1 = SUM;
		SUM = ADDI(SUM, TMP6);
		FREEMPI(TMP1);
		FREEMPI(TMP6);
		TMP1 = SUM;
		SUM = MOD(SUM, TWOABSN);
		FREEMPI(TMP1);
		if(RSV(SUM, ABSN) > 0){
			TMP1 = SUM;
			SUM = SUBI(SUM, TWOABSN);
			FREEMPI(TMP1);
		}
                MODULUS = MULT_I(ABSN, 2);
		printf(": n = ");
		fprintf(outfile, ": n = ");
		PRINTI(SUM);
		FPRINTI(outfile, SUM);
		printf("(mod ");
		PRINTI(MODULUS);
		printf(")\n");
		fprintf(outfile, "(mod ");
		FPRINTI(outfile, MODULUS);
		fprintf(outfile, ")\n");
		FREEMPI(SUM);
		FREEMPI(MODULUS);
		FREEMPI(x);
		FREEMPI(y);
	}
	FREEMPI(alpha);
	FREEMPI(beta);
	FREEMPI(gamma);
	FREEMPI(delta);
}
	FREEMPIA(FSX1);
	FREEMPIA(FSY1);
	FREEMPI(a);
	FREEMPI(b);
	FREEMPI(c);
	FREEMPI(TWOA);
	FREEMPI(TWOC);
	FREEMPI(ABSN);
	FREEMPI(TWOABSN);
        FREEMPI(A);
        FREEMPI(B);
        FREEMPI(C);
        FREEMPI(N);
	fclose(outfile);
	
	return;
}