func.c

calc大数库

	}
	M = MOD0(B, N);
	t = M->S;
	FREEMPI(M);
	if (t == 0)
	{
	  printf("N divides base\n");
	  return;
	}
	b = CONVERTI(B);	
	r = MILLER(N, b);	
	PRINTI(N);
	if (r)
		printf(" passes Miller's test to base %lu\n", b);
	else
		printf(" fails Miller's test to base %lu\n", b);
	return;
}

MPI *DIVISORX(MPI *N)
/* Returns the number of divisors of N. 
 * Returns NULL on failure.  Handled by wrapper */
{
	MPI *M;

	if (N->S <= 0)
	{
	  printf("N <= 0\n");
	  return NULL;
	}
	M = DIVISOR(N);
	if (M == NULL)
	{
	  printf("Factorization unsuccessful\n");
	  return NULL;
	}
	return (M);
}

MPI *MOBIUSX(MPI *N)
/* Returns  Mu(N). */
{
	MPI *M;

	if (N->S <= 0)
	{
	  printf("N <= 0\n");
	  return NULL;
	}
	M = MOBIUS(N);
	if (M == NULL)
	{
	  printf("Factorization unsuccessful\n");
	  return NULL;
	}
	return (M);
}

MPI *EULERX(MPI *N)
/* Returns Euler's function phi(N). */
{
	MPI *M;

	if (N->S <= 0)
	{
	  printf("N <= 0\n");
	  return NULL;
	}
	M = EULER(N);
	if (M == NULL)
	{
	  printf("Factorization unsuccessful\n");
	  return NULL;
	}
	return (M);
}

MPI *SIGMAX(MPI *N)
/* Returns sigma(N), the sum of the divisors of N. */
{
	MPI *M;

	if (N->S <= 0)
	{
	  printf("N <= 0\n");
	  return NULL;
	}
	M = SIGMA(N);
	if (M == NULL)
	{
	  printf("Factorization unsuccessful\n");
	  return NULL;
	}
	return (M);
}

MPI *LPRIMROOTX(MPI *P)
/* Returns the least primitive root mod P. */
{
	MPI *M;

	if (P->S <= 0 || (P->D == 0 && P->V[0] == 1))
	{
	  printf("P <= 1\n");
	  return NULL;
	}
	M = LPRIMROOT(P);
	if (M->S  == 0)
	{
	  printf("factorization of P-1 unsuccessful\n");
	  return NULL;
	}
	return (M);
}

MPI *ORDERPX(MPI *A, MPI *P)
/*
 * Returns the order of A mod P, a prime.
 */
{
	if (P->S <= 0 || (P->D == 0 && P->V[0] == 1))
	{
	  printf("P <= 1\n");
	  return NULL;
	}
	return (ORDERP(A, P));
}

MPI *ORDERQX(MPI *A, MPI *P, MPI *N)
/*
 * Returns the order of A mod P^N, P a prime.
 */
{
	unsigned int n;

	if (N->S <= 0 || N->D || (N->D == 0 && (N->V[0] > 1000)))
	{
	  printf("N <= 0 or N > 1000\n");
	  return NULL;
	}
	if (P->S <= 0 || (P->D == 0 && P->V[0] == 1))
	{
	  printf("P <= 1\n");
	  return NULL;
	}
	n = (USI)(CONVERTI(N));
	return (ORDERQ(A, P, n));
}

MPI *ORDERMX(MPI *A, MPI *M)
/*
 * Returns the order of A mod M.
 */
{
	MPI *G;
	unsigned int t;

	if (M->S <= 0 || (M->D == 0 && M->V[0] == 1))
	{
	  printf("M <= 1\n");
	  return NULL;
	}
	G = GCD(A, M);
	t = EQONEI(G);
	FREEMPI(G);
	if (!t)
	{
	  printf("GCD(A, M) > 1\n");
	  return NULL;
	}
	return (ORDERM(A, M));
}

MPI *LUCASUX(MPI *N, MPI *Q, MPI *M)
/*
 * returns U_n (mod m), where U_{k+1}=U_k-qU_{k-1}, U_0=0,U_1=1.
 * D=1-4q != 0.
 * We use the recurrence relations:
 * U_{k+1}=U_k-qU_{k-1},
 * U_{2k}=-2qU_{k-1}U_k+U_kU_k,
 * U_{2k-1}=U_kU_k-qU_{k-1}U_{k-1}
 * See Niven, Zuckermann and Montgomery, p.202.
 * Also D. Bressoud, Factorization and Primality Testing, p.179-191 and 
 * C. Pomerance, Kent State Lecture Notes, 19-22.
 */
{
	if (M->S <= 0 || (M->D == 0 && M->V[0] == 1))
	{
		FREEMPI(N); FREEMPI(Q); FREEMPI(M);
		execerror("M <= 1", "");
	}
	return (LUCASU(N, Q, M));
}

MPI *LUCASBX(MPI *N)
/*
 * Let N > 1 be odd. Then LUCASB(N) determines an integer D(!=-3) of the form
 * 4m+1,  such that the Jacobi symbol j(D,N)= -1. Then with Q=(1-d)/4, the Lucas
 * pseudoprime test examines L=U_{(N+1)/2}mod N. If L=0, N is a Lucas probable 
 * prime, otherwise N is composite. See "The Pseudoprimes to 25.10^9", 
 * Mathematics of computation, 35 (1980) 1003-1026. At the end of this paper 
 * it is conjectured that if N is a strong base 2 pseudoprime and a Lucas 
 * probable prime, then N is in fact a prime. A $30 prize is offered for a  
 * counterexample.
 * if LUCASB(N) = 1, then N is a Lucas probable prime, else N is composite.
 * Returns NULL if N is a perfect square.
 */
{
	if (N->S <= 0 || (N->D == 0 && N->V[0] == 1))
	{
		FREEMPI(N);
		execerror("N <= 1", "");
	}
	if ((N->V[0]) % 2 == 0)
	{
		FREEMPI(N);
		execerror("N is even", "");
	}
	return (LUCASB(N));
}

MPI *LUCASX(MPI *N)
/* Here N is odd and > 1.
 * If LUCAS(N) returns 1, then N is a strong base 2 pseudoprime and a Lucas
 * probable prime; if LUCAS(N) returns 0, then N is composite.
 * See "The Pseudoprimes to 25.10^9", Mathematics of computation, 35 (1980)
 * 1003-1026. At the end of this paper it is conjectured that if N is a strong
 * base 2 pseudoprime and a Lucas probable prime, then N is in fact a prime. 
 * A $30 prize is offered for a counterexample.
 */
{
	if (N->S <= 0 || (N->D == 0 && N->V[0] == 1))
	{
	  printf("N <= 1\n");
	  return NULL;
	}
	if ((N->V[0]) % 2 == 0)
	{
	  printf("N is even\n");
	  return NULL;
	}
	return (LUCAS(N));
}

MPI *LENGTHX(MPI *N)
/* Returns the number of decimal digits of N. */
{
	unsigned long z;

	z = LENGTHI(N);
	if (N->S < 0)
		z--;
	return (CHANGE(z));
}

MPI *MULT32X(MPI *X, MPI *Y)
{
	USL x, y, z;

	x = CONVERTI(X);
	y = CONVERTI(Y);
	z = MULT32(x, y);
	return (CHANGE(z));
}

MPI *NEXTPRIMEAPX(MPI *A, MPI *B, MPI *M)
/*
 * Finds the first Lucas probable prime P, A | P - B, M <= P.
 * Here A is even, B is odd, A > B >= 1, gcd(A,B) = 1, M > 1.
 */
{
	MPI *tmp1;
	unsigned int s;

	if (A->V[0] % 2)
	{
	  printf("A is odd\n");
	  return NULL;
	}
	if (B->V[0] % 2 == 0)
	{
	  printf("B is even\n");
	  return NULL;
	}
	if (A->S <= 0 || (A->D == 0 && A->V[0] == 1))
	{
	  printf("A <= 1\n");
	  return NULL;
	}
	if (B->S <= 0)
	{
	  printf("B <= 0\n");
	  return NULL;
	}
	if (M->S <= 0 || (M->D == 0 && M->V[0] == 1))
	{
	  printf("M <= 1\n");
	  return NULL;
	}
	if (RSV(A, B) <= 0)
	{
	  printf("A <= B\n");
	  return NULL;
	}
	tmp1 = GCD(A, B);
	s = EQONEI(tmp1);
	FREEMPI(tmp1);
	if (!s)
	{
	  printf("GCD(A,B) > 1\n");
	  return NULL;
	}
	return (NEXTPRIMEAP(A, B, M));
}
void SHALLIT()
{
	MPI *K, *N, *tmp, *L, *M;
	unsigned int i;

	K = ONEI();
	L = ONEI();
	VERBOSE=1;
	for (i = 2; i <= 20; i++)
	{
		M = CHANGE(4 * i - 2);
		tmp = MULTI(M, L);
		FREEMPI(M);
		N = ADD0I(tmp, K);
		FREEMPI(tmp);
		FACTOR(N);
		printf("completed K[%u]\n", i);
		printf("is the factorization of K[%u] = ", i);
		PRINTI(N);
		printf("\n");
		FREEMPI(K);
		K = L;
		L = N;
	}
	VERBOSE=0;
	FREEMPI(K);
	FREEMPI(L);
}

MPI *EFACTORX(MPI *N, MPI *M, MPI *P)
/* Elliptic curve factoring. 
 * Warning: 2329 is the number of array elements in PRIMEPOWERS[] 
 * in primepow.h
 */
{
	MPI *Z;
	unsigned long m, p;

	if (N->S <= 0 || (N->D == 0 && N->V[0] == 1))
	{
	  printf("N <= 1\n");
	  return NULL;
	}
	if (M->S <= 0 || (M->D == 0 && M->V[0] <= 10))
	{
	  printf("M <= 0 or M <= 10\n");
	  return NULL;
	}
	if ((M->D == 0 && (M->V[0] > 2328)) || M->D)
	{
	  printf("M > 2328\n");
	  return NULL;
	}
	if (P->S <= 0 || P->D >= 2)
	{
	  printf("P <= 0 or P >= 2^32\n");
	  return NULL;
	}
	m = CONVERTI(M);
	p = CONVERTI(P);
	Z = EFACTOR(N, m, p);
	if (Z == NULL)
	{
	  printf("Factorization unsuccessful\n");
	  return NULL;
	}
	return (Z);
}
void ABS_NEAREST_INTRX()
{
	unsigned int u;
	MPI *tempI, *G, *S;
	MPR *R;

	R = BUILDMPR();
	printf("enter the numerator:");
	R->N = INPUTI(&u);
	Flush();
	printf("enter the denominator (> 0):");
	R->D = INPUTI(&u);
	Flush();
	G = GCD(R->N, R->D);
	tempI = R->N;
	R->N = INT(R->N, G);
	FREEMPI(tempI);
	tempI = R->D;
	R->D = INT(R->D, G);
	FREEMPI(tempI);
	FREEMPI(G);
	S = ABS_NEAREST_INTR(R);
	printf("abs-nearest-int of "); PRINTR(R); printf("=");  PRINTI(S);
        printf("\n");
	FREEMPI(S);
	FREEMPR(R);
	return;
}


MPI *FIBONACCI(USI n)
/*
 * Returns the nth Fibonacci number.
 */
{
	MPI *A, *B, *C;
	USI i;

	C = NULL;
	if (n == 0)
		return(ZEROI());
	if (n == 1)
		return(ONEI());
	A = ZEROI();
	B = ONEI();
	for (i = 2; i <= n; i++)
	{
		C = ADD0I(A, B);
		FREEMPI(A);
		A = B;
		B = C;
	}
	FREEMPI(A);
	return(C);
}
MPI *LUCASS(USI n)
/*
 * Returns the nth Lucas number.
 */
{
	MPI *A, *B, *C, *THREE;
	USI i;

	C = NULL;
	THREE = BUILDMPI(1);
	THREE->S = 1;
	THREE->V[0] = 3;
	if (n == 0)
		return(ONEI());
	if (n == 1)
		return(THREE);
	A = ZEROI();
	B = ONEI();
	for (i = 2; i <= n; i++)
	{
		C = ADD0I(A, B);
		FREEMPI(A);
		A = B;
		B = C;
	}
	FREEMPI(A);
	return(C);
}