reductio.c

calc大数库

/* reduction.c */
#ifdef _WIN32
#include "unistd_DOS.h"
#else
#include <unistd.h>
#endif
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "integer.h"
#include "fun.h"

extern unsigned long PRIME[];
unsigned int llen;

unsigned int SQUAREFREE_TEST1000(USL n){
	USI i;
	USL temp;
	for(i=0; i <= 167; i++){
		temp = PRIME[i]*PRIME[i];
		if(n % temp == 0)
			return(0);
	}
	return(1);
}

USI global_count;
MPIA GLOBALARRAY_U;
MPIA GLOBALARRAY_V;

unsigned int PERIOD(MPI *D, MPI *U, MPI *V){
	USI j, k, cycle_length;
	MPI *F, *R, *S, *A, *tmp, *tmp1, *tmp2;

	F = BIG_MTHROOT(D, 2);
	k=global_count; /* initially set to 0 in POS() below */
	R = COPYI(U); 
	S = COPYI(V);
	for(j = k; 1; j++){
		tmp = ADD0I(F, R); 
		A = INT0(tmp, S); 
		FREEMPI(tmp);
		tmp = MULTI(A, S); 
		FREEMPI(A);
		tmp1 = R; 
		R = SUB0I(tmp, R);
                FREEMPI(tmp); 
		FREEMPI(tmp1);
                tmp = S; 
		tmp1 = MULTI(R, R);
                tmp2 = SUB0I(D, tmp1); 
		S = INT0(tmp2, S);
                ADD_TO_MPIA(GLOBALARRAY_U, R, j);
                ADD_TO_MPIA(GLOBALARRAY_V, S, j);
                FREEMPI(tmp); 
		FREEMPI(tmp1); 
                FREEMPI(tmp2);
 	
		if (EQUALI(U, R) && EQUALI(V, S)){
                        FREEMPI(R); 
			FREEMPI(S); 
			FREEMPI(F);
			break;
                }
                if(j == R0){
                        FREEMPI(R); 
			FREEMPI(S); 
			FREEMPI(F);
                        execerror("j = R0", "");
                }
	}
	cycle_length = j+1-k;
	global_count=global_count+cycle_length;
        return(cycle_length);
}
	
unsigned int CHECK_UVARRAYS(MPI *U, MPI *V){
	USI i;
	for(i = 0; i < global_count; i++){
		if(EQUALI(U, (GLOBALARRAY_U)->A[i]) && EQUALI(V, (GLOBALARRAY_V)->A[i])){
			return(0);
                }
        }
        return(1);
}

USI POS(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, *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);
        printf("Creating a complete list of reduced forms of discriminant %lu\n", 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 */
				        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;
					    }
					    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-");PRINTI(D);printf("*y^2=-4 has a solution\n");
        }else{
              printf("x^2-");PRINTI(D);printf("*y^2=-4 has no solution\n");
        }
        printf("h(%lu)=%u\n", d, class_number);
	FREEMPI(ONE);
	FREEMPI(F);
	FREEMPI(G);
	FREEMPI(DD);
	FREEMPIA(GLOBALARRAY_U);
        FREEMPIA(GLOBALARRAY_V);
        return(class_number);
}

MPI *POSX(MPI *D){
	USI class_number, w;
	USL d;
	int t;	
	MPI *ONE, *N, *TEMP;

	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);
	}
	w = SQUAREFREE_TEST1000(d);
	if(w == 0){
                printf("D is not squarefree\n");
                return(NULL);
        }else{
		class_number = POS(D);
		N = CHANGE(class_number);
		return(N);
	}
}

USI NEG(MPI *D, MPI *FLAG, MPI *TABLE_FLAG){
/*
 * This is Henri Cohen's Algorithm 5.3.5, p. 228.
 * for finding the class number h(D) of binary quadratic forms
 * of discriminant D, when D<0.
 * Here D=0(mod 4) or 1(mod 4).
 * If flag=1, we print only the primitive forms.
 * h(D) is returned in each case.
 * Davenport's Higher Arithmetic has a table of forms, which
 * lists the imprimitive ones with an asterisk.
 * If d is the discriminant of an imaginary quadratic field K,
 * then the primitive forms class-number h(D) is also the class number of K.
 * We print forms only when TABLE_FLAG is nonzero.
 */
	MPI *A, *B, *BB, *GG, *GGG, *TEMP, *TEMP1, *ABSD, *ONE;
	MPI *Q, *QQ, *T;
	int t;
	USI g, h;
	USL temp;

	ONE = ONEI();
	h = 0;
	g = 1;
	if(TABLE_FLAG->S){
		if(FLAG->S == 1){
		    printf("determining primitive forms\n");
		}else{
		    printf("determining primitive and imprimitive forms\n");
		}
	}
	temp = MOD_(D, 4);
	if(temp == 0){
                B = ZEROI();
	}else{
                B = ONEI();
	}
	ABSD = MINUSI(D);
	TEMP = INT0_(ABSD, 3);
	BB = BIG_MTHROOT(TEMP, 2);
	FREEMPI(TEMP);
	Q = NULL;
	A = NULL;
	while(RSV(B, BB)<=0){
		TEMP = MULTI(B, B);
		TEMP1 = ADD0I(TEMP, ABSD);
		FREEMPI(TEMP);
		Q = INT0_(TEMP1, 4);
		FREEMPI(TEMP1);
		A = COPYI(B);
		if(RSV(A, ONE) <= 0){
			TEMP = A;
			A = ONEI();
			FREEMPI(TEMP);
		}
		QQ = BIG_MTHROOT(Q, 2);
		while(RSV(A, QQ) <= 0){
			TEMP = MOD0(Q, A);
			t = TEMP->S;
			FREEMPI(TEMP);
			if(t == 0){
				T = INT0(Q, A);
				if(FLAG->S == 1){
					GG = GCD(A, B);
					GGG = GCD(GG, T);
					FREEMPI(GG);
					if(RSV(GGG, ONE) > 0){
						g = 0;
					}
					FREEMPI(GGG);
				}
				if(g == 1){
					if(RSV(A, B)==0 || RSV(A, T)==0 || B->S == 0){
						if(TABLE_FLAG->S){
							printf("(");
							PRINTI(A);
							printf(",");
							PRINTI(B);
							printf(",");
							PRINTI(T);
							printf(")\n");
						}
							h = h + 1;
						}else{
						if(TABLE_FLAG->S){
							printf("(");
							PRINTI(A);
							printf(",");
							PRINTI(B);
							printf(",");
							PRINTI (T);
							printf(")\n");

							printf("(");
							PRINTI(A);
							printf(",");
							TEMP = MINUSI(B);
							PRINTI(TEMP);
							FREEMPI(TEMP);
							printf(",");
							PRINTI(T);
							printf(")\n");
						}
							h = h + 2;
						}
					}else{
						g=1;
					}
				FREEMPI(T);
			}
			TEMP = A;
			A = ADD0_I(A, 1);
			FREEMPI (TEMP);
		}
		FREEMPI(A);
		FREEMPI(QQ);
		FREEMPI(Q);
		TEMP = B;
		B = ADD0_I(B, 2);
		FREEMPI(TEMP);
	}
	FREEMPI(B);
	FREEMPI(BB);
	FREEMPI(ABSD);
	FREEMPI(ONE);
	return(h);
}

MPI *NEGX(MPI *D, MPI *FLAG){
	MPI *ZERO, *H, *TEMP, *ONE;
	USL temp;
	USI h;
	int s, t;

	ZERO = ZEROI();
	ONE = ONEI();
	t = COMPAREI(D, ZERO);
	if(t >= 0){
		printf("D >= 0\n");
		FREEMPI(ZERO);
		FREEMPI(ONE);
		return(NULL);
	}
	TEMP = POWER_I (10, 6);
	s = RSV(D, TEMP);
	FREEMPI(TEMP);
	if(s >= 0){
		printf("|D| >= 10^6\n");
		FREEMPI(ZERO);
		FREEMPI(ONE);
		return(NULL);
	}
	t = COMPAREI(FLAG, ONE);
	if(t > 0){
		printf("flag > 1\n");
		FREEMPI(ONE);
		FREEMPI(ZERO);
		return(NULL);
	}
	t = COMPAREI(FLAG, ZERO);
	FREEMPI(ZERO);
	if(t < 0){
		printf("flag < 0\n");
		FREEMPI(ONE);
		return(NULL);
	}
	temp = MOD_(D, 4);
	if(temp == 2 || temp ==3){
		printf("D is not congruent to 0 or 1 (mod 4)\n");
		FREEMPI(ONE);
		return(NULL);
	}
	h = NEG(D, FLAG, ONE);
	FREEMPI(ONE);
	H = CHANGE((USL)h);
	return(H);
}

USI REDUCE_NEG(MPI *A, MPI *B, MPI *C){
/* This is Gauss' algorithm for reducing a positive definite binary
 * quadratic form. 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 (A,B,C)
 * satisfies -A<B<=A, C>=A, with B>=0 if C=A.
 * 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){/* -A < B <= A <= C */
		r = COMPAREI(A, C);
		if(t < 0 || (r == 0 && B->S >= 0)){/* A < C or A=C and B>=0 */
			printf("(");
			PRINTI(A);
			printf(",");
			PRINTI(B);
			printf(",");
			PRINTI(C);
			printf(") is reduced\n");
                	printf("Transforming matrix:1,0,0,1\n");
			FREEMPI(D);
                        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++;
		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);
	}
	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);
}

MPI *REDUCE_NEGX(MPI *A, MPI *B, MPI *C){
	MPI *D, *TEMP1, *TEMP2, *TEMP3;
	USI i;

	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);
	}
	i = REDUCE_NEG(A, B, C);
	FREEMPI(D);
	return(CHANGE(i));
}

USI CYCLE_PERIOD(MPI *D, MPI *U, MPI *V, USI i, int sign){
	MPI *A, *F, *UU, *VV, *X, *Y, *TEMP1, *TEMP2, *TEMP3;
	USI len, j;
	char buff[20];
	FILE *outfile;

        strcpy(buff, "reducepos.out");
	/*if(access(buff, R_OK) == 0)
	  unlink(buff);*/
	outfile = fopen(buff, "a");
	F = BIG_MTHROOT(D, 2);
	UU = COPYI(U);
	VV = COPYI(V);
	printf("\ncycle:\n->");
	fprintf(outfile, "\ncycle:\n->");
	/* i is created by reduce(a,b,c) below and indexes the ith
   convergent of the (U+sqrt(D))/V created there) */

	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);
		VV = INT(TEMP3, VV);
		FREEMPI(TEMP2);
		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){
			break;
		}
	}
	len = j + 1 - i;
        llen = len;

	if(len % 2 == 0){