func.c

calc大数库

/* func.c */
#include <stdio.h>
#include "integer.h"
#include "fun.h"

unsigned int VERBOSE;
unsigned int FREE_PRIMES;
extern unsigned long PRIME[];
unsigned int ECMAX;

MPI *FACTORX(MPI *N)
{
	unsigned int z;

	if (N->S <= 0 || (N->D == 0 && N->V[0] == 1))
	{
	  printf("N <= 1\n");
	  return NULL;
	}
	printf("enter the number of elliptic curves to be used:");
	scanf("%u", &ECMAX);
	getchar();
	VERBOSE = 1;
	FREE_PRIMES = 1;
	z = FACTOR(N);
	VERBOSE = 0;
	FREE_PRIMES = 0;
	if (z == 0)
	{
	  printf("Factorization unsuccessful\n");
	  return NULL;
	}
	return (CHANGE(z));
}

/* Returns NULL on failure.  Wrapper handles this. */
MPI *Nextprime(MPI *N)
{
	MPI *tmp;
	if (N->S < 0 )
	{
	  printf("N <= 0\n");
	  return NULL;
	}
	tmp = NEXT_PRIMEX(N);
	return(tmp);
}

MPI *JACOB(MPI *M, MPI *N)
/*
 * Evaluates the Jacobi symbol (M/N); N odd, N > 0.
 */
{
	int t;
	MPI *MM;

	if (MOD0_(N, (USL)2) == 0)
	{
	  printf("N is even\n");
	  return NULL;
	}
	if (N->S < 0)
	{
	  printf("N is negative\n");
	  return NULL;
	}
	MM = MOD(M, N);
	if (N->D == 0)
		t = JACOBI(CONVERTI(MM), CONVERTI(N));
	else	
		t = JACOBIB(MM, N);
	FREEMPI(MM);
	if (t == 0)
		return(ZEROI());
	else if (t == 1)
		return( ONEI());
	else
		return( MINUS_ONEI());
}


MPI *GCD_ARRAYVX(MPIA M, MPIA *Y)
/*
 * Returns d=gcd(M[0],...,M[N-1]) and an array Y[] of MPI's such that
 * d = M[0]*Y[0]+...+M[N-1]*Y[N-1]. Here N > 1.
 */
{
	MPI *Z;
	MPIA V;

	Z = GCD_ARRAYV(M, &V);
	*Y = V;	
	return (Z);
}
MPI *SQRTMX(MPI *x, MPI *p)
/*
 * Calculates sqrt(x) (mod p) using "A simple and fast probabilistic algorithm
 * for computing square roots modulo a prime number", I.E.E.E. Trans. Inform.
 * Theory, IT-32, 1986, 846-847, R. Peralta.
 * Here x is a quadratic residue mod p. x can be negative.
 * Returns NULL on failure.  This is taken care of by wrapper. 
 */
{
	unsigned long X, P, Z;
	MPI *M, *N=NULL, *T;
	int t;

	if (p->S <= 0 || (p->D == 0 && p->V[0] <= 2))
	{
	  printf("p <= 2\n");
	}
	T = LUCAS(p);
	t = T->S;
	FREEMPI(T);
	if (!t)
	{
	  printf("2nd argument is not a prime\n");
	  return N;
	}
	if (JACOBIB(x, p) != 1)
	{
	  printf("x is not a quadratic residue mod p\n");
	  return N;
	}
	if (p->D == 0)
	{
		M = MOD(x, p);
		X = CONVERTI(M);
		P = CONVERTI(p);
		Z = SQRTm(X, P);
		FREEMPI(M);
		return(CHANGE(Z));
	}
	else
	{
		N = SQRTM(x, p);
		return (N);
	}
}

MPI *CONGRX(MPI *A, MPI *B, MPI *M, MPI **N)
/*
 * Returns the least solution (mod N) of the congruence AX=B(mod M), where
 * N = M / gcd(A, M).
 * returns NULL if function fails. Handled by wrapper CONGRX_W.
 */
{
	MPI *Z=NULL;

	if (M->S <= 0 || (M->D == 0 && M->V[0] == 1))
	{
	  printf("M <= 1\n");
	  *N = ZEROI();
	  return Z;
	}
	Z = CONGR(A, B, M, N);
	if (Z == NULL)
	{
	  printf("the congruence has no solution\n");
	  *N = ZEROI();
	}
	return (Z);
}

MPI *CHINESEX(MPI *A, MPI *B, MPI *M, MPI *N, MPI **Mptr)
/*
 * Returns the solution mod *Mptr=lcm[M,N] of the congruences X = A (mod M)
 * and X = B (mod N), if soluble.
 * Returns NULL if it failed.  This case is handled by wrapper CHINESEX_W.
 */
{
	MPI *Z=NULL;
		
	if (M->S <= 0 || (M->D == 0 && M->V[0] == 1))
	{
	  printf("M <= 1\n");
	  *Mptr = ZEROI();
	  return(Z);
	}
	if (N->S <= 0 || (N->D == 0 && N->V[0] == 1))
	{
	  printf("N <= 1\n");
	  *Mptr = ZEROI();
	  return(Z);
	}
	Z = CHINESE(A, B, M, N, Mptr);
	if (Z == NULL)
	{
	  printf("the system has no solution\n");
	  *Mptr = ZEROI();
	}
	  return (Z);
}

MPI *CHINESE_ARRAYX(MPIA A, MPIA M, MPI **Mptr)
/*
 * Returns the solution mod *Mptr=lcm[M[0],...,M[n-1] of the congruences
 * X = A[i] (mod M[i]),0 <= i < N, if soluble; otherwise returns NULL.
 */
{
	unsigned int i, n;
	MPI *Z;

	n = ((A->size >= M->size) ? M->size: A->size);
	for (i = 0; i < n; i++)
	{
		if (M->A[i]->S <= 0 || (M->A[i]->D == 0 && M->A[i]->V[0] == 1))
		{
		  *Mptr = ZEROI();
		  printf("M[i] <= 1\n");
		  return NULL;
		}
	}
	Z = CHINESE_ARRAY(A->A, M->A, Mptr, n);
	if (Z == NULL)
	{
	  *Mptr = ZEROI();
	  printf("the system has no solution\n");
	  return NULL;

	}
	return (Z);
}

void MTHROOTX(MPI *Aptr, MPI *Bptr, MPI *M, MPI *R)
/*
 * *Aptr and *Bptr are positive MPI'S.
 * The mthroot of *Aptr/(*Bptr) is computed to R decimal places, R >= 0 ;
 * M, R are integers, 0 < M * R < R0 * R0.
 */
{
	unsigned int m, r;

	if (Aptr->S < 0 || Bptr->S <= 0 )
	{
	  printf("Numerator or denominator <= 0\n");
	  return;
	}
	if (M->S <= 0 || M->D >= 1 || M->V[0] == 1)
	{
	  printf("m <= 1 or m >= R0\n");
	  return;
	}
	if (R->S < 0 || R->D >= 1)
	{
	  printf("r < 0 or r >= R0\n");
	  return;
	}
	m = CONVERTI(M);
	r = CONVERTI(R);
	MTHROOT(Aptr, Bptr, m, r);
	return ;
}

MPI *BIG_MTHROOTX(MPI *Aptr, MPI *M)
/*
 * *Eptr, the integer part of the Mth root of the positive MPI *Aptr,
 * 1 < M < R0, is obtained by Newton's method, using the integer part function.
 * Returns NULL on failure.  Handled by wrapper function. 
 */
{
	unsigned int m;
	if (Aptr->S <= 0)
	{
	  printf("argument <= 0\n");
	  return NULL;
	}
	if (M->S <= 0 || M->D >= 1 || M->V[0] == 1)
	{
	  printf("m <= 1 or m >= R\n");
	  return NULL;
	}
	m = CONVERTI(M);
	return (BIG_MTHROOT(Aptr, m));
}

MPI  *FUND_UNITX(MPI *D, MPI **Xptr, MPI **Yptr)
/*
 * This is a program for finding the fundamental unit of Q(sqrt(D)).
 * The algorithm is based on K. Rosen, Elementary number theory
 * and its applications, p382, B.A. Venkov, Elementary Number theory, p.62 
 * and D. Knuth, Art of computer programming, Vol.2, p359, with Pohst's trick 
 * of using half the period.
 * w=(1+sqrt(D))/2 if D=1 (mod 4), w=sqrt(D) otherwise.
 * The norm of the fundamental unit is returned.
 *
 * Returns NULL on failure.  Wrapper function handles this.
 */
{
	MPI *G, *X, *tmp;
	unsigned int t;

	if (D->S <= 0)
	{
	  printf("D <= 0\n");
	  *Xptr = ZEROI();
	  *Yptr = ZEROI();
	  return NULL;
	}
	X = BIG_MTHROOT(D, 2);
	G = MULTI(X, X);
	t = EQUALI(D, G);
	FREEMPI(G);
	FREEMPI(X);
	if (t)
	{
		printf("D is a perfect square\n");
		*Xptr = ZEROI();
		*Yptr = ZEROI();
		return NULL;
	}
	tmp = FUND_UNIT(D, Xptr, Yptr);
	return (tmp);
}

MPI  *PELX(MPI *D, MPI *E, MPI **Xptr, MPI **Yptr)
/*
 * This is a program for finding the least solution of Pell's equation
 * x*x - D*y*y = +-1.
 * The algorithm is based on K. Rosen, Elementary number theory
 * and its applications, p382, B.A. Venkov, Elementary Number theory, p.62 
 * and D. Knuth, Art of computer programming, Vol.2, p359, with Pohst's trick 
 * of using half the period.
 * The norm of the least solution is returned.
 * The partial quotients are printed iff E->S != 0.
 *
 * Returns NULL on failure.  Wrapper handles this.
 */
{
	MPI *G, *X, *tmp;
	unsigned int t;

	if (D->S <= 0)
	{
	  printf("D <= 0\n");
	  *Xptr = ZEROI();
	  *Yptr = ZEROI();
	  return NULL;
	}
	X = BIG_MTHROOT(D, 2);
	G = MULTI(X, X);
	t = EQUALI(D, G);
	FREEMPI(G);
	FREEMPI(X);
	if (t)
	{
	  printf("D is a perfect square\n");
	  *Xptr = ZEROI();
	  *Yptr = ZEROI();
	  return NULL;
	}
	tmp = PEL(D, E, Xptr, Yptr);
	return (tmp);
}

MPI *SURDX(MPI *D, MPI *T, MPI *U, MPI *V, MPIA *A_ARRAY, MPIA *U_ARRAY, MPIA *V_ARRAY, MPIA *P_ARRAY, MPIA *Q_ARRAY)
/*
 * This function uses the continued fraction algorithm expansion in K. Rosen,
 * Elementary Number theory and its applications,p.379-381 and Knuth's
 * The art of computer programming, Vol.2, p. 359. It locates the first complete
 * quotient that is reduced and then uses the function REDUCED(D,U,V,i) above
 * to locate and return the period of the continued fraction expansion.
 *
 * Returns NULL on failure.  Wrapper handles this. 
 */
{
	MPI *G, *X;
	unsigned int t;
	unsigned long f;

	if (D->S <= 0)
	{
	  printf("D <= 0\n");
	  return NULL;
	}
	X = BIG_MTHROOT(D, 2);
	G = MULTI(X, X);
	f = EQUALI(D, G);
	FREEMPI(G);
	FREEMPI(X);
	if (f)
	{
	  printf("D is a perfect square\n");
	  return NULL;
	}
	if (V->S == 0)
	{
	  printf("V = 0\n");
	  return NULL;
	}
	t = SURD(D, T, U, V, A_ARRAY, U_ARRAY, V_ARRAY, P_ARRAY, Q_ARRAY, 0);
	return (CHANGE((USL)t));
}

MPI *MPOWERX(MPI *Aptr, MPI *Bptr, MPI *Cptr)
/*
 * *Aptr, *Bptr, *Cptr, *Eptr  are MPI'S, *Cptr > 0, *Bptr >= 0,
 * (*Aptr)^ *Bptr (mod *Cptr) is returned.
 * 
 * Returns NULL on failure.  Handled by wrapper.
 */ 
{
	if (Cptr->S <= 0)
	{
	  printf("Modulus <= 0\n");
	  return NULL;
	}
	if (Cptr->D == 0 && Cptr->V[0] == 1)
	{
	  printf("Modulus = 1\n");
	  return NULL;
	}
	return (MPOWER(Aptr, Bptr, Cptr));
}

MPI *INVERSEMX(MPI *A, MPI *M)
/*
 * Returns the inverse of A mod M.
 */
{
	MPI *temp, *Z;
	unsigned int t;

	if (M->S <= 0 || (M->D == 0 && M->V[0] == 1))
	{
	  printf("M <= 1\n");
	  return NULL;
	}
	temp = GCD(A, M);
	t = EQONEI(temp);
	FREEMPI(temp);
	if (!t)
	{
	  printf("A is not relatively prime to M\n");
	  return NULL;
	}
	temp = MOD(A, M);
	Z = INVERSEM(temp, M);
	FREEMPI(temp);
	return (Z);
}

void MILLERX(MPI *N, MPI *B)
/*
 * *N is odd, > 1, and does not divide b, 0 < b < R0.
 *  If *N passes Miller's test for base *B, *N is likely to be prime.
 * If *N fails Miller's test for base *B, then *N is composite.
 */
{
	unsigned long b, r;
	int t;
	MPI *M;

	if (N->S <= 0 || (N->D == 0 && N->V[0] == 1))
	{
	  printf("N <= 1\n");
	  return;
	}
	if (B->S <= 0 || B->D >= 1 || B->V[0] == 1)
	{
	  printf("BASE <= 1 or BASE >= R0\n");
	  return;
	}
	if (N->V[0] %2 == 0)
	{
	  printf("N is even\n");
	  return;