reductio.c

calc大数库

		FREEMPI(UU);
		FREEMPI(VV);
		FREEMPI(F);
	        fclose(outfile);
		return(len);
	}else{
		printf("->");
		fprintf(outfile, "->");
		for(j=i;1;j++){
			TEMP1 = ADDI(F, UU);
			A = INT(TEMP1, VV);
			FREEMPI(TEMP1);

			TEMP1 = UU;
			TEMP2 = MULTI(A, VV);
			FREEMPI(A);
			UU = SUBI(TEMP2, UU);
			FREEMPI(TEMP1);
			FREEMPI(TEMP2);
			
			TEMP1 = VV;
			TEMP2 = MULTI(UU, UU);
			TEMP3 = SUBI(D, TEMP2);
			FREEMPI(TEMP2);
			VV = INT(TEMP3, VV);
			FREEMPI(TEMP3);

			TEMP2 = INT_(TEMP1, 2);
			X = MULT_I(TEMP2, (long)sign);
			FREEMPI(TEMP1);
			FREEMPI(TEMP2);

                        if(j>i){
			    printf("->");
			    fprintf(outfile, "->");
                        }
			printf("(");
			fprintf(outfile, "(");
			PRINTI(X);
			FPRINTI(outfile, X);
			printf(",");
			fprintf(outfile, ",");
			FREEMPI(X);
			PRINTI(UU);
			FPRINTI(outfile, UU);
			printf(",");
			fprintf(outfile, ",");
			TEMP1 = INT_(VV, 2);
			Y = MULT_I(TEMP1, (long)(-sign));
			FREEMPI(TEMP1);
			PRINTI(Y);
			FPRINTI(outfile, Y);
			FREEMPI(Y);
		        printf(")");
		        fprintf(outfile, ")");
			
			sign = -sign;
			if(COMPAREI(UU, U) == 0 && COMPAREI(VV, V) == 0){
				FREEMPI(UU);
				FREEMPI(VV);
				FREEMPI(F);
			        printf("\n");
			        fprintf(outfile, "\n");
	                        fclose(outfile);
				return(2*len);
			}
		}
	}

}

USI REDUCE_POS(MPI *A, MPI *B, MPI *C){
	MPI *D, *F, *ONE, *U, *V, *TEMP, *TEMP1, *TEMP2, *TEMP3, *TEMP4, *TEMP5;
	MPI *X, *Y, *AA, *BB, *CC;
        MPI *A_ORIG, *B_ORIG,*C_ORIG, *PP, *QQ, *RR, *SS, *HH;
        MPI *A_TMP, *B_TMP,*C_TMP, *T1, *T2, *T3, *T4;
	int r, s, sign;
	USI j, len, y;
	char buff[20];
	FILE *outfile;
	
        A_ORIG = COPYI(A);
        B_ORIG = COPYI(B);
        C_ORIG = COPYI(C);
        A_TMP = COPYI(A);
        B_TMP = COPYI(B);
        C_TMP = COPYI(C);
        PP = ONEI();
	QQ = ZEROI();
        RR = ZEROI();
        SS = ONEI();

        strcpy(buff, "reducepos.out");
	if(access(buff, R_OK) == 0)
	  unlink(buff);
	outfile = fopen(buff, "w");
	ONE = ONEI();
	AA = COPYI(A);
	BB = COPYI(B);
	CC = COPYI(C);
	TEMP1 = MULTI(BB, BB);
	TEMP2 = MULTI(AA, CC);
	TEMP3 = MULT_I(TEMP2, 4);
	FREEMPI(TEMP2);
	D = SUBI(TEMP1, TEMP3);
	FREEMPI(TEMP1);
	FREEMPI(TEMP3);
	F = BIG_MTHROOT(D, 2);
	printf("(");
	fprintf(outfile, "(");
	PRINTI(AA);
	FPRINTI(outfile, AA);
	printf(",");
	fprintf(outfile, ",");
	PRINTI(BB);
	FPRINTI(outfile, BB);
	printf(",");
	fprintf(outfile, ",");
	PRINTI(CC);
	FPRINTI(outfile, CC);
	printf(")->");
	fprintf(outfile, ")->");
	sign = 1;
	U = COPYI(BB);
	TEMP1 = MULT_I(CC, 2);
	V = MULT_I(TEMP1, (long)(-sign));
	FREEMPI(TEMP1);
	for(j=0; 1; j++){
		if(j){
			TEMP1 = MULTI(U, U);
			TEMP2 = SUBI(D, TEMP1);
			TEMP3 = MULT_I(V, 2);
			TEMP4 = INTI(TEMP2, TEMP3);
			X = MULT_I(TEMP4, (long)sign);
                        A_TMP = COPYI(TEMP4);
                        B_TMP = COPYI(U);
			FREEMPI(TEMP1);
			FREEMPI(TEMP2);
			FREEMPI(TEMP3);
			FREEMPI(TEMP4);
			
			TEMP1 = INT_(V, 2);
			TEMP2 = MINUSI(TEMP1);
			Y = MULT_I(TEMP2, (long)sign);
                        C_TMP = COPYI(Y);
			FREEMPI(TEMP1);
			FREEMPI(TEMP2);
                        if(j>1){
			   printf("->");
			   fprintf(outfile, "->");
                       }
			printf("(");
			fprintf(outfile, "(");
			PRINTI(X);
			FPRINTI(outfile, X);
			printf(",");
			fprintf(outfile, ",");
			PRINTI(U);
			FPRINTI(outfile, U);
			printf(",");
			fprintf(outfile, ",");
			PRINTI(Y);
			FPRINTI(outfile, Y);
			printf(")");
			fprintf(outfile, ")");
			FREEMPI(X);
			FREEMPI(Y);
		}
		sign = -sign;
		if(V->S >0 && U->S > 0 && COMPAREI(U, F) <= 0){
			TEMP1 = ADDI(U, V);
			r = COMPAREI(F, TEMP1);
			FREEMPI(TEMP1);

			TEMP2 = SUBI(V, U);
			s = COMPAREI(TEMP2, F);
			FREEMPI(TEMP2);

			if(r < 0 && s <= 0){
                                printf("\nUnimodular matrix (");
                                fprintf(outfile, "\nUnimodular matrix (");
                                PRINTI(PP);
                                FPRINTI(outfile, PP);
                                printf(",");
                                fprintf(outfile, ",");
                                PRINTI(QQ);
                                FPRINTI(outfile, QQ);
                                printf(",");
                                fprintf(outfile, ",");
                                PRINTI(RR);
                                FPRINTI(outfile, RR);
                                printf(",");
                                fprintf(outfile, ",");
                                PRINTI(SS);
                                FPRINTI(outfile, SS);
                                printf(") transforms (");
                                fprintf(outfile, ") transforms (");
                                PRINTI(A_ORIG);
                                FPRINTI(outfile, A_ORIG);
                                printf(",");
                                fprintf(outfile, ",");
                                PRINTI(B_ORIG);
                                FPRINTI(outfile, B_ORIG);
                                printf(",");
                                fprintf(outfile, ",");
                                PRINTI(C_ORIG);
                                FPRINTI(outfile, C_ORIG);
                                printf(") to the reduced form (");
                                fprintf(outfile, ") to the reduced form (");
                                if (j){
                                     PRINTI(A_TMP);
                                     FPRINTI(outfile, A_TMP);
                                     printf(",");
                                     fprintf(outfile, ",");
                                     PRINTI(B_TMP);
                                     FPRINTI(outfile, B_TMP);
                                     printf(",");
                                     fprintf(outfile, ",");
                                     PRINTI(C_TMP);
                                     FPRINTI(outfile, C_TMP);
                                     printf(")\n");
                                     fprintf(outfile, ")\n");
                                }else{
                                     PRINTI(A_ORIG);
                                     FPRINTI(outfile, A_ORIG);
                                     printf(",");
                                     fprintf(outfile, ",");
                                     PRINTI(B_ORIG);
                                     FPRINTI(outfile, B_ORIG);
                                     printf(",");
                                     fprintf(outfile, ",");
                                     PRINTI(C_ORIG);
                                     FPRINTI(outfile, C_ORIG);
                                     printf(")\n");
                                     fprintf(outfile, ")\n");
                                }
                                FREEMPI(A_ORIG);
                                FREEMPI(B_ORIG);
                                FREEMPI(C_ORIG);
                                FREEMPI(A_TMP);
                                FREEMPI(B_TMP);
                                FREEMPI(C_TMP);
                                FREEMPI(PP);
                                FREEMPI(QQ);
                                FREEMPI(RR);
                                FREEMPI(SS);
	                        fclose(outfile);
				len = CYCLE_PERIOD(D, U, V, j, sign);
	                        outfile = fopen(buff, "a");
				printf("cfrac has period length %u\n", llen);
				fprintf(outfile, "cfrac has period length %u\n", llen);
				printf("cycle of reduced forms has period length %u\n", len);
				fprintf(outfile, "cycle of reduced forms has period length %u\n", len);
				FREEMPI(ONE);
				FREEMPI(D);
				FREEMPI(F);
				FREEMPI(U);
				FREEMPI(V);
				FREEMPI(AA);
				FREEMPI(BB);
				FREEMPI(CC);
	                        fclose(outfile);
				return(len);
			}
		}
                FREEMPI(A_TMP);
                FREEMPI(B_TMP);
                FREEMPI(C_TMP);
		if(V->S > 0){
			TEMP1 = ADDI(F, U);
			X = INTI(TEMP1, V);
			FREEMPI(TEMP1);
		}else{
			TEMP1 = ADDI(F, U);
			TEMP2 = ADDI(TEMP1, ONE);
			FREEMPI(TEMP1);
			X = INTI(TEMP2, V);
			FREEMPI(TEMP2);
		}
                T1 = MINUSI(QQ);
                T3 = MINUSI(SS);
                y = j % 2;
                if(y == 0){
                   HH = COPYI(X);
                }else{
                   HH = MINUSI(X);
                }
                TEMP5 = MULTI(QQ, HH);
                T2 = ADDI(PP, TEMP5);
                FREEMPI(TEMP5);
                TEMP5 = MULTI(SS, HH);
                T4 = ADDI(RR, TEMP5);
                FREEMPI(TEMP5);
                FREEMPI(HH);
                TEMP = PP;
                PP = T1;
                FREEMPI(TEMP);
                TEMP = QQ;
                QQ = T2;
                FREEMPI(TEMP);
                TEMP = RR;
                RR = T3;
                FREEMPI(TEMP);
                TEMP = SS;
                SS = T4;
                FREEMPI(TEMP);
               
		TEMP1 = U;
		TEMP2 = MULTI(X, V);
		U = SUBI(TEMP2, U);
		FREEMPI(TEMP1);
		FREEMPI(TEMP2);
		FREEMPI(X);

		TEMP1 = V;
		TEMP2 = MULTI(U, U);
		TEMP3 = SUBI(D, TEMP2);
		V = INTI(TEMP3, V);
		FREEMPI(TEMP1);
		FREEMPI(TEMP2);
		FREEMPI(TEMP3);
	}
}

MPI *REDUCE_POSX(MPI *A, MPI *B, MPI *C){
	MPI *D, *F, *G, *TEMP, *TEMP1, *TEMP2, *TEMP3;
	USI len;
	int t;
	
	TEMP1 = MULTI(B, B);
	TEMP2 = MULTI(A, C);
	TEMP3 = MULT_I(TEMP2, 4);
	FREEMPI(TEMP2);
	D = SUBI(TEMP1, TEMP3);
	FREEMPI(TEMP1);
	FREEMPI(TEMP3);

	F = BIG_MTHROOT(D, 2);
	G = MULTI(F, F);
	t = COMPAREI(G, D);
	FREEMPI(F);
	FREEMPI(G);
	if(t == 0){
		printf("B^2-4AC is a perfect square\n");
		FREEMPI(D);
		return(NULL);
	}
	t = D->S;
	if(t <= 0){
		printf("B^2-4*A*C <= 0\n");
		FREEMPI(D);
		return(NULL);
	}
	TEMP = POWER_I (10, 6);
	t = RSV(D, TEMP);
	FREEMPI(TEMP);
	if(t >= 0){
		printf("D >= 10^6\n");
		FREEMPI(D);
		return(NULL);
	}
	FREEMPI(D);
	len = REDUCE_POS(A, B, C);
	TEMP1 = CHANGE((USL)len);
	return(TEMP1);
}

USI POS0(MPI *D){
	USI a, b, e, t, z, parity_flag, class_number, x;
	USL d, f; /* d will be restricted to d < 10^6. */
	MPI *A, *B, *C, *TEMP, *F, *G, *H, *I, *J, *K, *L;
	MPI *ONE, *R;
	int sign, sign1;

	GLOBALARRAY_U = BUILDMPIA();
	GLOBALARRAY_V = BUILDMPIA();
	class_number = 0;
	ONE = ONEI();
	global_count = 0;
	parity_flag = 0;
	d = CONVERTI(D);

	if((d-1) % 4 != 0){
		e=2;
	}else{
		e=1;
	}
        printf("Creating a complete list of reduced forms of discriminant %lu\n", d);
	F = BIG_MTHROOT(D, 2);
	f = CONVERTI(F);
	G = INT0_(F, 2);
        for(a=1;a<=f;a++){
			A = CHANGE(a);
			I = MULT_I(A, 2);
			J = MULT_I(A, 4);
                for(b=e;b<=f;b=b+2){
			B = CHANGE(b);
			TEMP = MULTI(B, B);
			H = SUBI(TEMP, D);
			FREEMPI(TEMP);
			TEMP = MOD(H, J);
			sign = TEMP->S;
			FREEMPI(TEMP);
			if(sign == 0){
			   TEMP = SUBI(F, I);
			   sign1 = COMPAREI(TEMP, B);
			   FREEMPI(TEMP);
	   if((RSV(A, G) <= 0 && sign1 < 0) || (RSV(A, G) > 0 && sign1 >= 0)){
				    R = INT(H, J);
		                    TEMP = ABSI(R);
				    K = GCD(A, B);
				    L = GCD(K, TEMP);
				    FREEMPI(K);
				    FREEMPI(TEMP);
				    t = EQONEI(L);
				    FREEMPI(L);
				    if(t){
					TEMP = MINUSI(H);
					C = INT0(TEMP, J); /* C > 0 */
				        FREEMPI(TEMP);
					TEMP = MULT_I(C, 2);
					x = CHECK_UVARRAYS(B, TEMP);
					if(x){
					    class_number=class_number+1;
					    z=PERIOD(D,B,TEMP);
					    if(z % 2){
					       parity_flag = 1;
					    }
					    printf("[%u]: (",class_number);
				 	    PRINTI(A);
					    printf(", ");
				 	    PRINTI(B);
					    printf(", ");
				 	    PRINTI(R);
					    printf(")\n");
					}
				        FREEMPI(TEMP);
				        FREEMPI(C);
				    }
				    FREEMPI(R);
				}
			}
			FREEMPI(H); 
			FREEMPI(B);
		}
		FREEMPI(A);
		FREEMPI(I);
		FREEMPI(J);
	}

	if(parity_flag){
              printf("x^2-%lu*y^2=-4 has a solution\n", d);
        }else{
              printf("x^2-%lu*y^2=-4 has no solution\n", d);
	      class_number = 2*class_number;
              printf("There are reduced forms (-a,b,-c) as well\n");
        }
        printf("h(%lu)=%u\n", d, class_number);
	FREEMPI(ONE);
	FREEMPI(F);
	FREEMPI(G);
	FREEMPIA(GLOBALARRAY_U);
        FREEMPIA(GLOBALARRAY_V);
        return(class_number);
}

MPI *POS0X(MPI *D){
	USI class_number;
	USL d;
	int t;	
	MPI *ONE, *N, *TEMP, *F, *G;

	d = CONVERTI(D);
	ONE = ONEI();
	t = COMPAREI(D, ONE);
	FREEMPI(ONE);
	if(t <= 0){
                printf(" D <= 1\n");
		return(NULL);
	}
	TEMP = POWER_I (10, 6);
	t = RSV(D, TEMP);
	FREEMPI(TEMP);
	if(t >= 0){
		printf("D >= 10^6\n");
		return(NULL);
	}
	F = BIG_MTHROOT(D, 2);
	G = MULTI(F, F);
	t = COMPAREI(G, D);
	FREEMPI(F);
	FREEMPI(G);
	if(t == 0){
		printf("B^2-4AC is a perfect square\n");
		return(NULL);
        }else{
		class_number = POS0(D);
		N = CHANGE(class_number);
		return(N);
	}
}

USL TABLENEG(MPI *M, MPI *N){
/* The following gives a table of h(-d) for all squarefree d
 * in the range M<=d<=N<10^6. 
 */
	MPI *DD, *ONE, *ZERO;
	USL count, m, n, d, r;
	USI h, t;
	long dd;
	char buff[20];
	FILE *outfile;

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

        strcpy(buff, "tableneg.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{
                r = d % 4;
		dd = (long)d;
                if(r != 3){
		     DD = CHANGEL((long)(-4) * d);
                }else{

		     DD = CHANGEL(-dd);
                }
		ONE = ONEI();
		ZERO = ZEROI();
                t = NEG(DD, ONE, ZERO);
		FREEMPI(DD);
		FREEMPI(ONE);
		FREEMPI(ZERO);
		count++;
                printf("h(%ld) = %u\n",-d, t);
                fprintf(outfile, "h(%ld) = %u\n",-d, t);
            }
	}
	if(count==0){
		printf("no squarefree d in the range [%lu,%lu]\n",m,n);
	}
	fclose(outfile);
	return (count);
}

MPI *TABLENEGX(MPI *M, MPI *N){
	MPI *LIMIT, *COUNT;
	USL count;

	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);
	if(M->S <= 0 || N->S <= 0){
		printf("M < 1 or N < 1\n");
		return (NULL);
	}
	count = TABLENEG(M, N);
	COUNT = CHANGE(count);
	return(COUNT);
}

USI POS1(MPI *D, int *norm){
/* D is squarefree. This function performs Lagrange's method on all
 * reduced quadratic irrationals (b+\sqrt(Disc))/2|c|, where 4*c divides 
 * Disc-b^2, Disc being the discriminant. The class-number h(D) of Q(sqrt(D)
 * is calculated, as well as the norm of the fundamental unit.
 * For use in TABLEPOS(M,N) below.
 */

	USI a, b, e, t, z, parity_flag, class_number, x;
	USL d, f; /* d will be restricted to d < 10^6. */
	MPI *A, *B, *C, *TEMP, *F, *G, *H, *I, *J, *K, *L;
	MPI *ONE, *R, *DD;
	int sign, sign1;

	GLOBALARRAY_U = BUILDMPIA();
	GLOBALARRAY_V = BUILDMPIA();
	class_number = 0;
	ONE = ONEI();
	global_count = 0;
	parity_flag = 0;
	d = CONVERTI(D);
        /* creates a fundamental discriminant */
	if((d-1) % 4 != 0){
		d=4*d;
		e=2;
	}else{
		e=1;
	}
	DD = CHANGE(d);
	F = BIG_MTHROOT(DD, 2);
	f = CONVERTI(F);
	G = INT0_(F, 2);
        for(a=1;a<=f;a++){
			A = CHANGE(a);
			I = MULT_I(A, 2);
			J = MULT_I(A, 4);
                for(b=e;b<=f;b=b+2){
			B = CHANGE(b);
			TEMP = MULTI(B, B);
			H = SUBI(TEMP, DD);
			FREEMPI(TEMP);
			TEMP = MOD(H, J);
			sign = TEMP->S;
			FREEMPI(TEMP);
			if(sign == 0){
			   TEMP = SUBI(F, I);
			   sign1 = COMPAREI(TEMP, B);
			   FREEMPI(TEMP);
	   if((RSV(A, G) <= 0 && sign1 < 0) || (RSV(A, G) > 0 && sign1 >= 0)){
				    R = INT(H, J);
		                    TEMP = ABSI(R);
				    K = GCD(A, B);
				    L = GCD(K, TEMP);
				    FREEMPI(K);
				    FREEMPI(TEMP);
				    t = EQONEI(L);
				    FREEMPI(L);
				    if(t){
					TEMP = MINUSI(H);
					C = INT0(TEMP, J); /* C > 0 */