primes1.c

calc大数库

	}
}

MPI *SQROOTX(MPI *A, MPI *N, MPIA *Y, MPI **M, USI *l)
{
	MPI *tmp;

	if(N->S <= 0){
		printf("N <= 0\n");
		return(NULL);
	}
	else
	{
		tmp = SQROOT(A, N, Y, M, l);
		return(tmp);
	}
}

void CORNACCHIA(MPI *A, MPI *B, MPI *M)
/* Cornacchia's algorithm. See L'algorithme de Cornacchia, A. Nitaj,
 * Expositiones Mathematicae, 358-365.
 */
{
	MPI *A1, *A2, *Tmp1, *TT, *T1, *T2, *N;
	MPI *AA, *BB, *R, *Q1, *QQ, *MA, *MB; 
	MPIA ARRAY;
	USI ll, i, jj;
	int t1, t2;

	A1 = INVERSEM(A, M);	
	Tmp1 = MINUSI(B);
	A2 = MULTI(A1, Tmp1);
	FREEMPI(A1);
	FREEMPI(Tmp1);
	MA = INT0(M, A);
	MB = INT0(M, B);
	N = SQROOTX(A2, M, &ARRAY, &Tmp1, &ll);
	FREEMPI(N);
	FREEMPI(A2);
	FREEMPI(Tmp1);
	for( i = 0; i < ll; i++){
		TT = ZEROI();
		T1 = ONEI();
		AA = COPYI(M);
		BB = COPYI(ARRAY->A[i]);
		t2 = 0;
		jj = 0;
		R = COPYI(AA);
		T2 = ZEROI();
		while(t2 <= 0){
			Tmp1 = MULTI(R, R);
			t1 = RSV(Tmp1, MA);
			FREEMPI(Tmp1);
			if(t1 <= 0){
				printf("X=");
				PRINTI(R);
				printf(", Y=");
				if(T2->S == -1)
					T2->S = 1;
				PRINTI(T2);
				printf("\n");
				break;
			}
			jj++;
			if(jj == 1){
				FREEMPI(R);
				R = COPYI(BB);
				FREEMPI(T2);
				T2 = ONEI();
				t2 = RSV(T2, MB);
			}
			if(jj > 1){
				if(jj==2){
					FREEMPI(R);
					FREEMPI(T2);
				}
				R = MOD0(AA, BB);
				Q1 = INT0(AA, BB);
				FREEMPI(AA);
				AA = BB;
				BB = R;
				QQ = MINUSI(Q1);
				FREEMPI(Q1);
				T2 = MULTABC(TT, QQ, T1);
				FREEMPI(TT);
				FREEMPI(QQ);
				Tmp1= MULTI(T2, T2);
				t2 = RSV(Tmp1, MB);
				FREEMPI(Tmp1);
				TT = T1;
				T1 = T2;
			}
		}
		FREEMPI(TT);
		FREEMPI(T1);
		if(jj == 1){
			FREEMPI(T2);
			FREEMPI(R);
		}
		FREEMPI(AA);
		FREEMPI(BB);
	}
	FREEMPIA(ARRAY);
	FREEMPI(MA);
	FREEMPI(MB);
}


MPI *CORNACCHIAX(MPI *A, MPI *B, MPI *M)
{
	MPI *Tmp;
	int t;
	USI s;

	if (A->S <= 0)
	{
	  printf("A <= 0\n");
	  return NULL;
	}
	if (B->S <= 0)
	{
	  printf("B <= 0\n");
	  return NULL;
	}
	Tmp = ADD0I(A, B);
	t = RSV(Tmp, M);
	FREEMPI(Tmp);
	if (t > 0)
	{
	  printf("M < A + B\n");
	  return NULL;
	}
	Tmp = GCD(A, B);
	s = EQONEI(Tmp);
	FREEMPI(Tmp);
	if (s == 0)
	{
	  printf("gcd(A,B) > 1\n");
	  return NULL;
	}
	Tmp = GCD(A, M);
	s = EQONEI(Tmp);
	FREEMPI(Tmp);
	if (s == 0)
	{
	  printf("gcd(A,M) > 1\n");
	  return NULL;
	}
	
	CORNACCHIA(A, B, M);
	return(ONEI());
}

void GAUSS(MPI *A, MPI *B, MPI *C, MPI *N, MPI **alpha, MPI **gamma, MPI **M)
/* input: gcd(A,B,C)=1, |N|>1. 
 * output: (alpha,gamma), where A*alpha^2+B*alpha*gamma+C*gamma^2=M, gcd(M,N)=1.
 */
{
	MPI *ABSN, **QQ, *TMP1, *TMP2, *TMP3, *TMP4; 
	MPI **XX, **YY, *MM, *NN, *G;
	USI e, i;
	int s, t;

	ABSN = ABSI(N);
	e = FACTOR(ABSN);
	FREEMPI(ABSN);
	QQ = (MPI **)mmalloc(e * sizeof(MPI *));
	for(i = 0; i< scount; i++)
                QQ[i] = CHANGE(Q[i]);
	for(i = 0; i< s_count; i++)
                QQ[i + scount] = Q_[i];
	XX = (MPI **)mmalloc(e * sizeof(MPI *));
	YY = (MPI **)mmalloc(e * sizeof(MPI *));
	for(i = 0; i < e; i++){
		TMP1 = MOD(A, QQ[i]);
		TMP2 = MOD(C, QQ[i]);
		s = TMP1->S;
		t = TMP2->S;
		FREEMPI(TMP1);
		FREEMPI(TMP2);
		if(s){
			XX[i] = ONEI();
			YY[i] = ZEROI();
		}
		if(!s && t){
			XX[i] = ZEROI();
			YY[i] = ONEI();
		}
		if(!s && !t){
			XX[i] = ONEI();
			YY[i] = ONEI();
		}
	}
	TMP1 = CHINESE_ARRAY(XX, QQ, &MM, e);
	TMP2 = CHINESE_ARRAY(YY, QQ, &NN, e);
	FREEMPI(MM);
	FREEMPI(NN);
	G = GCD(TMP1, TMP2);
	*alpha = INT(TMP1, G);
	FREEMPI(TMP1);
	*gamma = INT(TMP2, G);
	FREEMPI(TMP2);
	FREEMPI(G);
	TMP1 = MULTI3(A, *alpha, *alpha);
	TMP2 = MULTI3(B, *alpha, *gamma);
	TMP3 = MULTI3(C, *gamma, *gamma);
	TMP4 = ADDI(TMP1, TMP2);
	*M = ADDI(TMP4, TMP3);
	FREEMPI(TMP1);
	FREEMPI(TMP2);
	FREEMPI(TMP3);
	FREEMPI(TMP4);
	for (i = 0; i < e; i++){
		FREEMPI(QQ[i]);
		FREEMPI(XX[i]);
		FREEMPI(YY[i]);
	}
	ffree((char *)QQ, e * sizeof(MPI *));
	ffree((char *)XX, e * sizeof(MPI *));
	ffree((char *)YY, e * sizeof(MPI *));
	return;
}

/*
#include <stdio.h>
#include "unistd.h"
*/

MPI *PRIME_GENERATOR(MPI *M, MPI *N)
/* lists the primes P satisfying M <= P <= N.  
 * Returns the number of such primes.
 */
{
	unsigned long j, k;
	MPI *C, *P, *Temp, *Temp1, *Temp2, *MM, *MMM, *NN, *ONE;
	unsigned int c=0;
	int t;
	char buff[20];
        FILE *outfile;

	ONE = ONEI();
	strcpy(buff, "primes.out");
	if(access(buff, R_OK) == 0)
		  unlink(buff);
	outfile = fopen(buff, "w");

	if(RSV(M, N)==1){
		Temp = COPYI(N);
		NN=COPYI(M);
		MM=Temp;
	}
	else{
		MM=COPYI(M);
		NN=COPYI(N);
	}
	if(MM->D==0 && MM->V[0]==2){
		c++;
		printf("2\n");
		fprintf(outfile, "2\n");
	}
	t = (MM->V[0])% 2;
	MMM = COPYI(MM);
	if(t == 0){
		Temp = MM;
		MM = ADD0I(MM, ONE);
		FREEMPI(Temp);
	}
	FREEMPI(ONE);
	P = COPYI(MM);
	while (RSV(P, NN) <= 0){
		j = 1;
		k = 1;
		while(1){
			Temp2 = CHANGE(PRIMES[k]);
			Temp1 = MULT_I(Temp2, (long)PRIMES[k]);
			FREEMPI(Temp2);
			t = RSV(P, Temp1);
			FREEMPI(Temp1);
			if(t < 0){
				break;
			}
			if(MOD0_(P, PRIMES[k])==0){
				j=0;
				break;
			}
			if(k >= 9591)
				break;
			k++;
           	}
		if(j == 1){
			c++;
			PRINTI(P); printf("\n");
			FPRINTI(outfile, P); fprintf(outfile, "\n");
		}
		Temp = P;
		P = ADD0_I(P, (USL)2);
		FREEMPI(Temp);
	}
	FREEMPI(P);
	printf("the number of primes in the range ");
	PRINTI(MMM); 
	printf(" to ");
	PRINTI(NN); 
	fprintf(outfile, "the number of primes in the range ");
	FPRINTI(outfile, MMM); 
	fprintf(outfile, " to ");
	FPRINTI(outfile, NN); 
	printf(" is %u\n",c);
	fprintf(outfile," is %u\n",c);
	fclose(outfile);
	C = CHANGE(c);
	FREEMPI(MM);
	FREEMPI(NN);
	FREEMPI(MMM);
	return (C);
}

MPI *PRIME_GENERATORX(MPI *M, MPI *N)
{
	MPI *BOUND, *ONE, *NUMBER;
	int s, t;

	BOUND = BUILDMPI(3);
	BOUND->V[0]  = 58368;
	BOUND->V[1]  = 21515;
	BOUND->V[2]  = 2;
	BOUND->S=1;
	if(RSV(N, BOUND) >= 0){
		printf("n>="); PRINTI(BOUND); printf("\n");
		FREEMPI(BOUND);
		return (NULL);
	}
	FREEMPI(BOUND);
	ONE = ONEI();
	s = RSV(M, ONE);
	t = RSV(N, ONE);
	FREEMPI(ONE);
	if(s <= 0){
		printf("m <= 1\n");
		return (NULL);
	}
	if(t <= 0){
		printf("n <= 1\n");
		return (NULL);
	}
	NUMBER = PRIME_GENERATOR(M, N);
	return (NUMBER);
}

MPI *SIGMAK(USI k, MPI *N)
/*
 *  Returns the sum of the kth powers of the divisors of N,
 *  returns NULL if factorization unsuccessful .
 *  Here 0 < k <= 10000 and N >= 1.
 */
{
	MPI *U, *V, *W, *Tmp;
	unsigned int i, s;

	if (EQONEI(N))
		return (ONEI());
	/* FREE_PRIMES = 0, so we can free the Q_[i] later */ 
	s = FACTOR(N);
	if (s == 0)
	{
		return ((MPI *)NULL);
	}
	U = ONEI();
	for (i = 0; i < scount; i++)
	{
		V = POWER_I((long)(Q[i]), k*(K[i] + 1));
		Tmp = V;
		V = SUB0_I(V, 1);
		FREEMPI(Tmp);
		W = POWER_I((long)(Q[i]), k);
		Tmp = W;
		W = SUB0_I(W, 1);
		FREEMPI(Tmp);
		Tmp = V;
		V = INT0(V, W);
		FREEMPI(Tmp);
		FREEMPI(W);
		Tmp = U;
		U = MULTI(U, V);
		FREEMPI(Tmp);
		FREEMPI(V);
	}
	for (i = 0; i < s_count; i++)
	{	V = POWERI(Q_[i], k*(K_[i] + 1));
		Tmp = V;
		V = SUB0_I(V, 1);
		FREEMPI(Tmp);
		W = POWERI(Q_[i], k);
		Tmp = W;
		W = SUB0_I(W, 1);
		FREEMPI(Tmp);
		Tmp = V;
		V = INT0(V, W);
		FREEMPI(W);
		FREEMPI(Tmp);
		Tmp = U;
		U = MULTI(U, V);
		FREEMPI(Tmp);
		FREEMPI(V);
	}
	for (i = 0; i < s_count; i++)
		FREEMPI(Q_[i]);
	return (U);
}

/*
 *  Returns the sum of the kth powers of the divisors of N,
 *  returns NULL if factorization unsuccessful .
 *  Here 0 < k <= 10000 and N >= 1.
 */
MPI *SIGMAKX(MPI *K, MPI *N){
	USI k;
	
	if (K->S <= 0){
	    printf("k <= 0\n");	
	    return(NULL);
	}
	if (N->S <= 0){
	    printf("n <= 0\n");	
	    return(NULL);
	}
	if (K->D > 0){
	    printf("k >= 2^16\n");	
	    return(NULL);
	}
	k = (USI)CONVERTI(K);
	if(k > 10000){
	    printf("k > 10000\n");	
	    return(NULL);
	}
	return(SIGMAK(k, N));
}


/* TAU(n) is Ramanujan's tau function. See T.M. Apostol, 'Modular functions
 * and Dirichlet series in number theory', 20-22, 140.
 * We assume n < 2^16
 */
MPI *TAU(USI n){
	USI i;
	MPI *S, *T, *TT, *TEMP1, *TEMP2, *TEMP3, *TEMP, *TEMPS;

	if(n == 1){
		return(ONEI());
	}
	S = ZEROI();
	for(i=1;i<n;i++){
		TT = CHANGE((USL)i);
		TEMP1 = SIGMAK(5, TT);
		FREEMPI(TT);
		TT = CHANGE((USL)(n-i));
		TEMP2 = SIGMAK(5, TT);
		FREEMPI(TT);
		TEMP = MULTI(TEMP1, TEMP2);
		FREEMPI(TEMP1);
		FREEMPI(TEMP2);
		TEMPS = S;
		S = ADD0I(S, TEMP);
		FREEMPI(TEMPS);
		FREEMPI(TEMP);
	}
	TT = CHANGE(n);
	TEMP1 = SIGMAK(11, TT);
	TEMP = TEMP1;
	TEMP1 = MULT_I(TEMP1, 65);
	FREEMPI(TEMP);

	TEMP2 = SIGMAK(5, TT);
	FREEMPI(TT);
	TEMP = TEMP2;
	TEMP2 = MULT_I(TEMP2, 691);
	FREEMPI(TEMP);

	TEMP3 = MULT_I(S, 252);
	TEMP = TEMP3;
	TEMP3 = MULT_I(TEMP3, 691);
	FREEMPI(TEMP);

	T = ADDI(TEMP1, TEMP2);
	FREEMPI(TEMP1);
	FREEMPI(TEMP2);
	TEMP = T;
	T = SUBI(T, TEMP3);
	FREEMPI(TEMP);
	FREEMPI(TEMP3);

	TEMP = T;
	T = INT_(T, 756);
	FREEMPI(TEMP);

	FREEMPI(S);
	return(T);
}

/* TAUX(N) is Ramanujan's tau function. See T.M. Apostol, 'Modular functions
 * and Dirichlet series in number theory', 20-22, 140.
 * We assume n < 2^16
 */
MPI *TAUX(MPI *N){
	USI n;

	if (N->S <= 0){
	    printf("n <= 0\n");	
	    return(NULL);
	}
	if (N->D > 0){
	    printf("n >= 2^16\n");	
	    return(NULL);
	}
	n = (USI)CONVERTI(N);
	return(TAU(n));
}

MPI *TAU_PRIMEPOWER(USI n, USI p){
	USI j, t, temp1, temp2;
	MPI *S, *T, *U, *TEMP, *TEMP1, *TEMP2;

	t = n/2;
	S = ZEROI();
	for(j = 0; j <= t;j++){
		temp1 = n - j;
		temp2 = temp1 - j;
		TEMP = BINOMIAL(temp1, temp2);
		TEMP1 = POWER_I(p, 11 * j);
		T = TAU(p);	
		TEMP2 = POWERI(T, temp2);
		FREEMPI(T);
		U = MULTI3(TEMP, TEMP1, TEMP2);
		FREEMPI(TEMP);
		FREEMPI(TEMP1);
		FREEMPI(TEMP2);
		if(j % 2){
			U->S = -(U->S);
		}
		TEMP = S;
		S = ADDI(S, U);
		FREEMPI(TEMP);
		FREEMPI(U);
	}	
	return(S);
}

MPI *TAU_PRIMEPOWERX(MPI *N, MPI *P){
	USI n, p, t;
	MPI *T;
	if (N->S <= 0){
	    printf("n <= 0\n");	
	    return(NULL);
	}
	if (N->D > 0){
	    printf("n >= 2^16\n");	
	    return(NULL);
	}
	n = (USI)CONVERTI(N);
	T = LUCAS(P);
	t = EQONEI(T);
	FREEMPI(T);
	if (t== 0){
	    printf("p is not prime\n");	
	    return(NULL);
	}

	p = (USI)CONVERTI(P);
	return(TAU_PRIMEPOWER(n, p));
}

MPI *TAU_COMPOSITE(USI n){
	USI i, d, *KGLOB, *QGLOB;
	MPI *U, *N, *TEMP, *T;

	U = ONEI();
	N = CHANGE(n);
	d = FACTOR(N);
	FREEMPI(N);
	KGLOB = (USI *)mmalloc(d * sizeof(USI));
	QGLOB = (USI *)mmalloc(d * sizeof(USI));
	for(i = 0; i < scount; i++){
		KGLOB[i] = K[i];
		QGLOB[i] = Q[i];
	}
	for(i = scount; i < d; i++){
		KGLOB[i] = K_[i];
		QGLOB[i] = CONVERTI(Q_[i]); /* n < 2^32 here */
	}
	for(i = 0; i < d; i++){
		TEMP = U;
		T = TAU_PRIMEPOWER(KGLOB[i], QGLOB[i]);
		U = MULTI(U, T);
		FREEMPI(TEMP);
		FREEMPI(T);
	}
	for (i = 0; i < s_count; i++)
		FREEMPI(Q_[i]);
	ffree(KGLOB, d * sizeof(USI));
	ffree(QGLOB, d * sizeof(USI));
	return(U);
}

MPI *TAU_COMPOSITEX(MPI *N){
	USI n;

	if (N->S <= 0){
	    printf("n <= 0\n");	
	    return(NULL);
	}
	if (N->D > 0){
	    printf("n >= 2^16\n");	
	    return(NULL);
	}
	n = (USI)CONVERTI(N);
	return(TAU_COMPOSITE(n));
}