primes1.c

calc大数库

{
	USL n;
	MPI *X;

	MPI *A = stackPop(S);
	MPI *P = stackPop(S);
	MPI *N = stackPop(S);
	n = CONVERTI(N);
	X = SQROOT1(A, P, n);
	printf("SQROOT of ");
	PRINTI(A);
        printf(" (mod ");
	PRINTI(P);
        printf("^%lu) is ", n);
	PRINTI(X);
	FREEMPI(X);
	FREEMPI(A);
	FREEMPI(P);
	FREEMPI(N);
	printf("\n");
	return;
}

USI sqroot2_number; /* global variable, used in SQROOT */

MPI *SQROOT2(MPI *A, USL n)
{
/* Solving x^2=A (mod 2^n), A odd.
 * returns -1 if not soluble, otherwise returns x.
 * In the case of solubility, changes global variable sqroot2_number to 1,2,4 
 * according as n= 1, 2  or > 2.
 */

	MPI *U, *V, *Tmp1, *Tmp2;
	USI i, t;
	USL s;

	if(n == 1){
		sqroot2_number = 1;
		return(ONEI());
	}
	if(n == 2){
		s = MOD_(A, 4);
		if (s != 1)
			return(MINUS_ONEI());
		else{
			sqroot2_number = 2;
			return(ONEI());
		}
	}
	s = MOD_(A, 8);
	if(s != 1)
		return(MINUS_ONEI());
	if(n == 3){
		sqroot2_number = 4;
		return(ONEI());
	}
	else{
		U = SQROOT2(A, n - 1);
		Tmp1 = MULTI(U, U);
                V = SUBI(Tmp1, A);
                FREEMPI(Tmp1);
		for(i = 0;i < n - 1;i++){
			Tmp1 = V;
			V = INT_(V, 2);
			FREEMPI(Tmp1);
		}
		Tmp2 = ONEI();
		Tmp1 = ADDI(V, Tmp2);
		FREEMPI(Tmp2);
		FREEMPI(V);
		t = (Tmp1->V[0]) % 2;
		FREEMPI(Tmp1);
		if (t == 0){
			Tmp1 = POWER_I(2, (USI)(n - 2));
			Tmp2 = U;
			U = ADDI(U, Tmp1);
			FREEMPI(Tmp1);
			FREEMPI(Tmp2);
		}
		return(U);
	}

}

void SQROOT2_W(Stack S)
{
	USL n;
	MPI *X;

	MPI *A = stackPop(S);
	MPI *N = stackPop(S);
	n = CONVERTI(N);
	X = SQROOT2(A, n);
	printf("SQROOT of ");
	PRINTI(A);
        printf(" (mod 2");
        printf("^%lu) is ", n);
	PRINTI(X);
	FREEMPI(X);
	FREEMPI(A);
	FREEMPI(N);
	printf("\n");
	return;
}

MPI *SQROOT3(MPI *A, MPI *P, USL n, MPI**EXPONENT)
/* Solving x^2=A (mod P^n), P a prime possibly dividing A.
 * Returns -1 if not soluble, otherwise returns x.
 * Also returns a "reduced moduli" explained as follows:
 * If P does not divide A, the story is well-known.
 * If P divides A, there are two cases:
 * (i) P^n divides A,
 * (ii) P^n does not divide A.
 * In case (i) there are two cases:
 * (a) n=2m+1, (b) n=2m.
 * In case (a), the solution is x=0 (mod P^(m+1)).
 * In case (b), the solution is x=0 (mod P^m).
 * (These are called reduced moduli.)
 * In case (i) gcd(A,P^n)=P^r, r<n. If r is odd, no solution.
 * If r=2m, write x=(P^m)*X, A=(p^2m)*A1, P not dividing A1. 
 * Then x^2=A (mod P^n) becomes X^2=A1 (mod P^(n-2m)).
 * If P is odd, this has 2 solutions X mod P^(n-2m) and we get two solutions
 * x (mod P^(n-m)). If P=2, we get
 *                            1 solution mod (2^(n-m)) if n-2m=1,
 *                            2 solutions mod (2^(n-m) if n-2m=2,
 *                            4 solutions mod (2^(n-m)) if n-2m>2.
 */
{
	MPI *D, *U, *V, *AA, *Tmp1, *T;
	USI r, m, s;
	int t;

/* The next two lines are important. eg. a call sqroot(431,10,&a[],&m,&l)
   sets global variable sqroot2_number=1, whereas a subsequent call to 
   sqroot(46,210,&a[],&m,&l) didn't call SQROOT2() and before the next fix
   was added, sqroot2_number remained equal to 1, causing problems as
   action is taken in SQROOT() if S[0] = 2 only if sqroot2_number > 0. */
	if ((P->D == 0) && P->V[0] == 2)
		sqroot2_number = 0;
	*EXPONENT = NULL;
	D = MOD(A, P);
	t = D->S;
	FREEMPI(D);
	if(t){
		*EXPONENT = POWERI(P, n);
		if ((P->D == 0) && P->V[0] == 2)
			return(SQROOT2(A, n));
		else
			return(SQROOT1(A, P, n));
	}
	else{
		U = POWERI(P, n);
		V = MOD(A, U);
		t = V->S;
		FREEMPI(U);
		FREEMPI(V);
		if(t == 0){
			if( (n % 2) == 0)
				*EXPONENT = POWERI(P, n/2);
			else
				*EXPONENT = POWERI(P, 1 + n/2);
			return(ZEROI());
		}
		r = 0;
		AA = COPYI(A);
		V = MOD(AA, P);
		while((V->S == 0) && (r <= n)){
			Tmp1 = AA;
			AA = INT(AA, P);
			FREEMPI(Tmp1);
			Tmp1 = V;
			V = MOD(AA, P);
			FREEMPI(Tmp1);
			r++;
		}
		FREEMPI(V);
		if(r % 2){
			FREEMPI(AA);
			return(MINUS_ONEI());
		}
		else{
			m = r / 2;
			s = n - r;
			*EXPONENT = POWERI(P, n - m);
			if((P->D == 0) && (P->V[0] == 2)){
				T = SQROOT2(AA, s);
				if(EQMINUSONEI(T)){
					FREEMPI(AA);
					return(T);
				}
				else{
					FREEMPI(T);
					U = POWERI(P, m);
					T = SQROOT2(AA, s);
					V = MULTI(U, T);
					FREEMPI(AA);
					FREEMPI(T);
					FREEMPI(U);
					return(V);
				}
			}
			else{
				T = SQROOT1(AA, P, s);
                                if(EQMINUSONEI(T)){
                                        FREEMPI(AA);
                                        return(T);
                                }
                                else{
                                        FREEMPI(T);
                                        U = POWERI(P, m);
                                        T = SQROOT1(AA, P, s);
                                        V = MULTI(U, T);
                                        FREEMPI(AA);
                                        FREEMPI(T);
                                        FREEMPI(U);
                                        return(V);
                                }

			}
		}
	}
}

void SQROOT3_W(Stack S)
{
	USL n;
	MPI *X, *EXPONENT;

	MPI *A = stackPop(S);
	MPI *P = stackPop(S);
	MPI *N = stackPop(S);
	n = CONVERTI(N);
	X = SQROOT3(A, P, n, &EXPONENT);
	printf("SQROOT of ");
	PRINTI(A);
        printf(" (mod ");
	PRINTI(P);
        printf("^%lu) is ", n);
	PRINTI(X);
	FREEMPI(X);
        printf(" (mod ");
	PRINTI(EXPONENT);
	FREEMPI(A);
	FREEMPI(P);
	FREEMPI(N);
	FREEMPI(EXPONENT);
	printf(")\n");
	return;
}

MPI *SQROOT(MPI *A, MPI *N, MPIA *SOL, MPI **MODULUS, USI *lptr)
/* Returns an array SOL consisting of the solutions of x^2=A (mod N), 
 * with 0<=x<=N/2 where N = prod QQ[i]^KK[i].
 * Returns the number of solutions (mod N) and their MODULUS if soluble,
 * otherwise returns -1.
 * The algorithm uses the Chinese remainder theorem after solving mod
 * QQ[i]^KK[i], i=1,...,e=omega(N). Note N is even if and only if QQ[0]=2.
 * Finally we join all solutions together using the CRT.
 * Some care is needed to exhibit all solutions. We operate on the D[]
 * array below, as follows:
 * l= number of primes QQ[i] dividing N, for which P^a=QQ[i]^KK[i] does not 
 * divide N, while jj is the remaining number.
 * D[] is the array formed by a single solution each of x^2=A mod P^a
 * SS is the product of the reduced moduli for these l primes.
 * OM is the product of the reduced moduli for the jj primes.
 * S[] is the array formed by the moduli P^a.
 * If l>1, we produce all 2^l l-tuples (e[0]*D[0],...,e[l-1]*D[l-1]),
 * where e[i]=1 or -1, if s is odd or D[0]=2 and "sqroot2_number">1, but only 
 * 2^(l-1) * (l-1)-tuples (e[0]*D[0],...,e[l-2]*D[l-2]), if SS is even and 
 * "sqroot2_number"=1.
 * It is convenient to use the CRT with D[l]=om, S[l]=0, even when OM=1,
 * so as to produce all solutions.
 * If SS is even, we place D[0] at the end of the D[] array, making it D[l-1].
 * This is done so as to conveniently operate on the initial part of the D[]
 * array, where the initial S[i] are only odd prime powers.
 * The program aborts if the number of prime factors of N which do not divide A
 * is > 31. (Not an easy program to get right. Finished 12th April 2000.)
 */
{
	USI l, e, jj, i, r, j, k, t, *KK, rr, ii; 
	MPI *C, *OM, **D, **D1, **DD, **DD1, **S, *E, *X, *AA, *B, *Y, *SS;
	MPI *Tmp1, *Tmp2, *Tmp3, *SO, *Z, **QQ, *F, *X1;
	int tt;

	if(EQONEI(N)){
			*SOL = BUILDMPIA(); 
			Y = ZEROI();
			ADD_TO_MPIA(*SOL, Y, 0);
			FREEMPI(Y);
			*MODULUS = ONEI();
			*lptr = 1;
			return(ONEI());
	}
	e = FACTOR(N);

	QQ = (MPI **)mmalloc(e * sizeof(MPI *));
	KK = (USI *)mmalloc(e * sizeof(USI));
	for(i = 0; i< scount; i++){
		QQ[i] = CHANGE(Q[i]);
		KK[i] = K[i];
	}
	for(i = 0; i< s_count; i++){
		QQ[i + scount] = Q_[i];
		KK[i + scount] = K_[i];
	}
	l = 0;
	jj = 0;
	OM = ONEI();
	SS = ONEI();
	D = (MPI **)mmalloc(e * sizeof(MPI *));
	S = (MPI **)mmalloc((e+1) * sizeof(MPI *));
	for(i = 0;i < e;i++){
		C = SQROOT3(A, QQ[i], KK[i], &E);
		t = EQMINUSONEI(C);
		tt = C->S;
		if(t){
			FREEMPI(OM);
			FREEMPI(E);
			FREEMPI(SS);
			FREEMPI(C);
			for (ii = 0; ii < e; ii++)
				FREEMPI(QQ[ii]);
			ffree((char *)D, e * sizeof(MPI *));
			ffree((char *)S, (e+1) * sizeof(MPI *));
			ffree((char *)QQ, e * sizeof(MPI *));
			ffree((char *)KK, e * sizeof(USI));
			*SOL = BUILDMPIA(); 
			ADD_TO_MPIA(*SOL, NULL, 0);
			*MODULUS = (MPI *)NULL;
			*lptr = 0;
			return(MINUS_ONEI());
		}
		else if (tt == 0){
			Tmp1 = OM;
			OM = MULTI(OM, E);
			FREEMPI(Tmp1);
			FREEMPI(E); /* added 10/12.00 */
			FREEMPI(C);
		}
		else{
			D[l] = C;
			S[l] = E;
			Tmp1 = SS;
			SS = MULTI(SS, S[l]);
			FREEMPI(Tmp1);
			l++;
		}
	}
	for (i = 0; i < e; i++)
		FREEMPI(QQ[i]);
	ffree((char *) QQ, e * sizeof(MPI *));
	ffree((char *) KK, e * sizeof(USI));
	if(l)
		SO = MULTI(S[0], OM);
	else
		SO = NULL;
	if(l == 0){
		X = INT(N, OM);
		FREEMPI(SS);
		ffree((char *) D, e * sizeof(MPI *));
		ffree((char *) S, (e+1) * sizeof(MPI *));
		*SOL = BUILDMPIA(); 
		Y = ZEROI();
		ADD_TO_MPIA(*SOL, Y, 0);
		FREEMPI(Y);
		*MODULUS = OM;
		*lptr = 1;
		return(X);
	}
	if(l == 1){
		if(!EQONEI(OM)){
			Tmp1 = ZEROI();
			X = CHINESE(D[0], Tmp1, S[0], OM, &Tmp2);
			FREEMPI(Tmp2);
			FREEMPI(Tmp1);
		}
		else
			X = COPYI(D[0]);
		if((S[0]->V[0]) % 2){
			Tmp1 = MULT_I(N, 2);
			Y = INT(Tmp1, SO);
			FREEMPI(Tmp1);
			FREEMPI(SS);
			FREEMPI(OM);
			FREEMPI(D[0]);
			ffree((char *) D, e * sizeof(MPI *));
			for(i=0;i<l;i++) /* changed 4/1/01 */
				FREEMPI(S[i]);
			ffree((char *) S, (e+1) * sizeof(MPI *));
			*SOL = BUILDMPIA(); 
			ADD_TO_MPIA(*SOL, X, 0);
			FREEMPI(X);
			*MODULUS = SO;
			*lptr = 1;
			return(Y);
		}
		else{/* S[0] is even */
			t = sqroot2_number;
			if(t == 4){
				Tmp2 = INT_(S[0], 2);
				Tmp3 = ADDI(D[0], Tmp2);
				FREEMPI(Tmp2);
				Y = MOD(Tmp3, S[0]);
				FREEMPI(Tmp3);
				if(!EQONEI(OM)){
					Tmp1 = ZEROI();
					X1 = CHINESE(Y, Tmp1, S[0], OM, &Tmp2);
					FREEMPI(Tmp1);
					FREEMPI(Tmp2);
				}
				else
					X1 = COPYI(Y);
				FREEMPI(Y);
				FREEMPI(D[0]);
				FREEMPI(SS);
				FREEMPI(OM);
				for(i=0;i<l;i++)/* 20/12/00 - had e not l */
					FREEMPI(S[i]);
				ffree((char *) D, e * sizeof(MPI *));
				ffree((char *) S, (e+1) * sizeof(MPI *));
				Y = MULT_I(N, 4);
				Z = INT0(Y, SO);
				FREEMPI(Y);
				*SOL = BUILDMPIA(); 
				ADD_TO_MPIA(*SOL, X, 0);
				ADD_TO_MPIA(*SOL, X1, 1);
				*MODULUS = SO;
				FREEMPI(X);
				FREEMPI(X1);
				*lptr = 2; /* bugfix: used to have =1, 4/Jul/00 */
				return(Z);
			}
			else if (t == 2){
				Y = MULT_I(N, 2);
				Z = INT0(Y, SO);
				FREEMPI(D[0]);
				FREEMPI(SS);
				FREEMPI(OM);
				FREEMPI(Y);
				ffree((char *) D, e * sizeof(MPI *));
				for(i=0;i<l;i++)/* 20/12/00 - had e not l */
					FREEMPI(S[i]);
				ffree((char *) S, (e+1) * sizeof(MPI *));
				*SOL = BUILDMPIA(); 
				ADD_TO_MPIA(*SOL, X, 0);
				FREEMPI(X);
				*MODULUS = SO;
				*lptr = 1;
				return(Z);
			}
			else{
				Z = INT0(N, SO);
				FREEMPI(D[0]);
				FREEMPI(SS);
				FREEMPI(OM);
				ffree((char *) D, e * sizeof(MPI *));
				for(i=0;i<e;i++)
					FREEMPI(S[i]);
				ffree((char *) S, (e+1) * sizeof(MPI *));
				*SOL = BUILDMPIA(); 
				ADD_TO_MPIA(*SOL, X, 0);
				FREEMPI(X);
				*MODULUS = SO;
				*lptr = 1;
				return(Z);
			}
		}
	}
	if(l > 31)
		execerror("l in SQROOT >31", "");
	F = POWER_I(2, l);
	k = (USI)CONVERTI(F);
	D1 = (MPI **)mmalloc((l + 1)* sizeof(MPI *));
	S[l] = OM;
	D1[l] = ZEROI();
	FREEMPI(SO);
	SO = MULTI(SS, OM);
	FREEMPI(SS);
	rr = 0;
	*SOL = BUILDMPIA(); 
	if((N->V[0]) % 2){
		for(i = 0;i < k;i++){
			j = i;
			for(r = 0;r < l;r++){
				t = j % 2;
				if (t)
					D1[r] = MINUSI(D[r]);
				else
					D1[r] = COPYI(D[r]);
				j = j / 2;
			}
			X = CHINESE_ARRAY(D1, S, &Tmp2, l + 1);
			for(r = 0;r < l;r++)
				FREEMPI(D1[r]);
			Tmp1 = MULT_I(X, 2);
			if(RSV(Tmp1, Tmp2) <= 0){
				ADD_TO_MPIA(*SOL, X, rr);
				rr++;
			}
			FREEMPI(X);
			FREEMPI(Tmp1);
			FREEMPI(Tmp2);
		}
		for(i = 0;i < l;i++){
			FREEMPI(D[i]);
			FREEMPI(S[i]);
		}
		FREEMPI(D1[l]);
		FREEMPI(S[l]);
		ffree((char *) D, e * sizeof(MPI *));
		ffree((char *) S, (e+1) * sizeof(MPI *));
		ffree((char *) D1, (l + 1) * sizeof(MPI *));
		Tmp1 = MULTI(F, N);
		FREEMPI(F);
		Z = INT0(Tmp1, SO);
		FREEMPI(Tmp1);
		*MODULUS = SO;
		*lptr = rr;
		return(Z);
	}
	else{/* N is even */
		DD = (MPI **)mmalloc((l + 1) * sizeof(MPI *));
		DD1 = (MPI **)mmalloc((l + 1) * sizeof(MPI *));
		if(sqroot2_number == 4)
			DD1[l] = ZEROI();
		AA = D[0];
		B = S[0];
		for(r = 1;r < l;r++){
			D[r - 1] = D[r];
			S[r - 1] = S[r];
		}
		D[l - 1] = AA;
		S[l - 1] = B;
		if(sqroot2_number == 4){
			for(r = 0;r < l - 1;r++)
				DD[r] = COPYI(D[r]);	
			Tmp2 = INT_(S[l - 1], 2);
			Tmp3 = ADDI(D[l - 1], Tmp2);
			FREEMPI(Tmp2);
			DD[l - 1] = MOD(Tmp3, S[l - 1]);
			FREEMPI(Tmp3);
		}
		if(sqroot2_number == 1)
			k = k / 2;
			
		for(i = 0;i < k;i++){
			j = i;
			for(r = 0;r < l;r++){
				t = j % 2;
				if (t){
					D1[r] = MINUSI(D[r]);
					if(sqroot2_number == 4)
						DD1[r] = MINUSI(DD[r]);
				}
				else{
					D1[r] = COPYI(D[r]);
					if(sqroot2_number == 4)
						DD1[r] = COPYI(DD[r]);
				}
				j = j / 2;
			}
			X = CHINESE_ARRAY(D1, S, &Tmp2, l + 1);
			Tmp1 = MULT_I(X, 2);
			if(RSV(Tmp1, Tmp2) <= 0){
				ADD_TO_MPIA(*SOL, X, rr);
				rr++;
			}
			FREEMPI(Tmp1);
			FREEMPI(X);
			FREEMPI(Tmp2);
			if(sqroot2_number == 4){
				X = CHINESE_ARRAY(DD1, S, &Tmp2, l + 1);
				Tmp1 = MULT_I(X, 2);
				if(RSV(Tmp1, Tmp2) <= 0){
					ADD_TO_MPIA(*SOL, X, rr);
					rr++;
				}
				FREEMPI(X);
				FREEMPI(Tmp1);
				FREEMPI(Tmp2);
			}
			for(r = 0;r < l;r++){
				FREEMPI(D1[r]);
				if(sqroot2_number == 4)
					FREEMPI(DD1[r]);
			}
		}
		for(i = 0;i < l;i++){
			FREEMPI(D[i]);
			FREEMPI(S[i]);
		}
		FREEMPI(D1[l]);
		if(sqroot2_number == 4){
			FREEMPI(DD1[l]);
			for(i = 0;i < l;i++)
				FREEMPI(DD[i]);
		}
		FREEMPI(S[l]);
		ffree((char *) D, e * sizeof(MPI *));
		ffree((char *) S, (e+1) * sizeof(MPI *));
		ffree((char *) D1, (l + 1) * sizeof(MPI *));
		ffree((char *) DD1, (l + 1) * sizeof(MPI *));
		ffree((char *) DD, (l + 1) * sizeof(MPI *));
		Tmp1 = MULTI(F, N);
		if(sqroot2_number == 1){
			Tmp2 = Tmp1;
			Tmp1 = INT0_(Tmp1, 2);
			FREEMPI(Tmp2);
		}
		if(sqroot2_number == 4){
			Tmp2 = Tmp1;
			Tmp1 = MULT_I(Tmp1, 2);
			FREEMPI(Tmp2);
		}
		FREEMPI(F);
		Z = INT0(Tmp1, SO);
		FREEMPI(Tmp1);
		*MODULUS = SO;
		*lptr = rr;
		return(Z);