i5m.c

calc大数库

/* i5m.c */
/* modular arithmetic */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "integer.h"
#include "fun.h"

extern unsigned int reciprocal[][Z0];
extern long int nettbytes;
void DUMPMPI(MPI *Mptr, char *name)
{
	unsigned int i;
	if (Mptr == NULL) {
	  printf("NULL\n");
	  return;
	}
	  
	printf("%s:", name);
	printf("Mptr->S = %d, Mptr->D = %u, ", Mptr->S, Mptr->D);
	for (i = 0; i <= Mptr->D; i++)
		printf("Mptr->V[%u] = %lu, ", i, Mptr->V[i]);
	printf("\n");
	return;
}

MPI *ADDM(MPI *Aptr, MPI *Bptr, MPI *Mptr)
/*
 * *Cptr = (*Aptr + *Bptr) mod (*Mptr). 
 * Here 0 <= *Aptr, *Bptr, *Cptr < *Mptr. 
 */
{
	MPI *Cptr, *TempI;

	Cptr = ADD0I(Aptr, Bptr);
	if (RSV(Cptr, Mptr) >= 0)
	{
		TempI = Cptr;
		Cptr = SUB0I(Cptr, Mptr);
		FREEMPI(TempI);
	}
	return (Cptr);
}

unsigned long ADDm(USL a, USL b, USL m)
/*
 * returns a + b mod m, if m > 0, where 0 <= a,b < m < 2^32.
 * a + b in GF(4) if m = 0.
 */
{
	unsigned long c;

	if (a == 0)
		return b;
	if (b == 0)
		return a;
	if (m)
	{
		c = (a < m - b) ? a + b : a - (m - b);
		return c;
	}
	else
	{
		if (a == 1 && b == 1)
			return (USL)0;
		else if (a == 1 && b == 2)
			return (USL)3;
		else if (a == 1 && b == 3)
			return (USL)2;
		else if (a == 2 && b == 1)
			return (USL)3;
		else if (a == 2 && b == 2)
			return (USL)0;
		else if (a == 2 && b == 3)
			return (USL)1;
		else if (a == 3 && b == 1)
			return (USL)3;
		else if (a == 3 && b == 2)
			return (USL)1;
		else  /*(a == 3 && b == 3)*/
			return (USL)0;
	}
}

MPI *MINUSM(MPI *Aptr, MPI *Mptr)
/*
 * *Bptr = -(*Aptr) mod (*Mptr). 
 * Here 0 <= *Aptr, *Bptr < *Mptr. 
 */
{

	if (Aptr->S == 0)
		return ZEROI();
	else
		return SUB0I(Mptr, Aptr);
}

unsigned long MINUSm(USL a, USL m)
/*
 * returns -a mod m if m > 0; here 0 <= a < m < 2^32.
 * -a  = a in GF(4) if m = 0.
 */
{
	if (m)
	{
		if (a == 0)
			return (USL)0;
		else
			return (m - a);
	}
	else
		return a;
}

MPI *SUBM(MPI *Aptr, MPI *Bptr, MPI *Mptr)
/*
 * *Cptr = (*Aptr - *Bptr) mod (*Mptr). 
 * Here 0 <= *Aptr, *Bptr, *Cptr < *Mptr.
 */
{
	MPI *Cptr, *TempI;

	Cptr = MINUSM(Bptr, Mptr);
	TempI = Cptr;
	Cptr = ADDM(Aptr, Cptr, Mptr);
	FREEMPI(TempI);
	return (Cptr);
}

unsigned long SUBm(USL a, USL b, USL m)
/*
 * returns a - b mod m if m > 0; here 0 <= a, b < m < 2^32.
 * a - b = a + b in GF(4) if m = 0.
 */
{
	if (m)
	{
		if (a >= b)
			return (a - b);
		else
		{
			return (MINUSm(b - a, m));
		}
	}
	else
		return ADDm(a, b, (USL)0);
}

MPI *MULTM(MPI *Aptr, MPI *Bptr, MPI *Mptr)
/*
 * *Cptr = (*Aptr * *Bptr) mod (*Mptr).
 * Here 0 <= *Aptr, *Bptr, *Cptr < *Mptr.
 */
{
	MPI *TempI, *Cptr;

	Cptr = MULTI(Aptr, Bptr);
	TempI = Cptr;
	Cptr = MOD0(Cptr, Mptr);
	FREEMPI(TempI);
	return (Cptr);
}

unsigned long MULTm(USL a, USL b, USL m)
/*
 * returns a * b mod m  if m > 0; here 0 <= a, b < m < 2^32.
 * a * b in GF(4) if m = 0.
 */
{
	if (a == 0 || b == 0)
		return (USL)0;
	if (m)
		return (RUSSm(0, a, b, m));
	else
	{
		if (a == 1)
			return b;
		else if (b == 1)
			return a;
		else if (a == 2 && b == 2)
			return (USL)3;
		else if (a == 2 && b == 3)
			return (USL)1;
		else if (a == 3 && b == 2)
			return (USL)1;
		else  /*(a == 3 && b == 3)*/
			return (USL)2;
	}
}

MPI *INVERSEM(MPI *Nptr, MPI *Mptr)
/*
 * here gcd(*Nptr, *Mptr) = 1, 1 <= *Nptr < *Mptr. 
 * *Kptr satisfies 1 <= *Kptr < *Mptr with (*Kptr) * (*Nptr) = 1  mod *Mptr.
 * See Knuth, p. 325. 
 */
{
	MPI *A, *B, *C, *H, *L, *Q, *Kptr;
	MPI *TempI;
	unsigned long i;

	if (EQONEI(Nptr))
		return ONEI();
	A = COPYI(Mptr);
	B = COPYI(Nptr);
	C = MOD(A, B);
	Kptr = ONEI();
	L = ZEROI();
	i = 1;
	while (C->S > 0)
	{
		Q = INT0(A, B);
		FREEMPI(A);
		A = B;
		B = C;
		C = MOD0(A, B);
		TempI = Q;
		Q = MULTI(Q, Kptr);
		FREEMPI(TempI);
		H = ADD0I(L, Q);
		FREEMPI(Q);
		FREEMPI(L);
		L = Kptr;
		Kptr = H;
		i = 1 - i;
	}
	FREEMPI(L);
	FREEMPI(A);
	FREEMPI(B);
	FREEMPI(C);
	if (i == 0)
	{
		TempI = Kptr;
		Kptr = SUB0I(Mptr, Kptr);
		FREEMPI(TempI);
	}
	return (Kptr);
}

unsigned long GCDm(USL m, USL n)
/*
 * returns the gcd of m and n.
 */
{
	unsigned long c;

	if (n == 0)
		return m;
        c = m % n;
        while (c)
	{
		m = n;
                n = c;
                c = m % n;
        }
	return n;
}   

unsigned long INVERSEm(USL n, USL m)
/*
 *  If 1 <= n < m < 2^32, gcd(n, m) = 1; returns k, 1 <= k < m with
 * n * k congruent to 1 mod m. 
 * returns n^{-1} in GF(4) if m = 0.
 */
{
	unsigned long a, b, c, h, i, k, l, q;
	unsigned int t;

	if (n == 1)
		return ((USL)1);
	if (m)
	{
		if (n == m - 1)
			return (m - 1);
		if (m <= Z0 + 1)
		{
			t =reciprocal[m - 2][n - 1];
			return (USL)t;
		}
		a = m;
		b = n;
		c = a % b;	
		l = 0;
		k = 1;
		i = 1;
		while (c)
		{
			q = a / b;
			a = b;
			b = c;
			c = a % b;
			h = l + q * k; /* no overflow - see Knuth p. 595 */
			l = k;
			k = h;
			i = 1 - i;
		}
		if (i)
			return (k);
		else
			return (m - k);
		return (k);
	}
	else
	{
		if (n == 2)
			return (USL)3;
		else /* (n == 3)*/
			return (USL)2;
	}
}

MPI *DIVM(MPI *Aptr, MPI *Bptr, MPI *Mptr)
/*
 * *Cptr = (*Aptr / *Bptr) mod (*Mptr). 
 * Here 0 <= *Aptr, *Bptr, *Cptr < *Mptr, gcd(*Bptr, *Mptr) = 1.
 */
{
	MPI *K, *Cptr;

	K = INVERSEM(Bptr, Mptr);
	Cptr = MULTM(Aptr, K, Mptr);
	FREEMPI(K);
	return (Cptr);
}

unsigned long DIVm(USL a, USL b, USL m)
/*
 * returns a / b mod m if m > 0, a / b in GF(4) if m = 0.
 * here 0 <= a, b < m < R0, gcd(b, m) = 1 if m > 0, b != 0 if m = 0.
 */
{
	unsigned long k;

	k = INVERSEm(b, m);
	return (MULTm(a, k, m));
}

unsigned long POWERm(USL a, MPI *Bptr, USL m)
/*
 * returns a^ *Bptr mod m.
 * here 0 <= a < m < R0 and 0 <= *Bptr.
 */
{
	unsigned long x, z;
	MPI *Y, *TempI;
	
	x = a;
	Y = COPYI(Bptr);
	z = 1;
	while (Y->S)
        {
		while (!(Y->V[0] & 1))
		{
			TempI = Y;
			Y = INT0_(Y, (USL)2);
			FREEMPI(TempI);
			x = MULTm(x, x, m);
		} 
		TempI = Y;
		Y = SUB0_I(Y, (USL)1);
		FREEMPI(TempI);
		z = MULTm(z, x, m);
	}
	FREEMPI(Y);
	return (z);
}

unsigned long POWER_m(USL a, USL y, USL m)
/*
 * returns a^ y mod m. 
 * Here 0 <= a < m < R0 and 0 <= y is an unsigned int.
 */ 
{
	unsigned long x, z;
	
	x = a;
	z = 1;
	while (y)
        {
		while (y % 2 == 0)
	        {
			y = y / 2;
			x = MULTm(x, x, m);
		} 
		y = y - 1;
		z = MULTm(z, x, m);
	}
	return (z);
}

MPI *MPOWER(MPI *Aptr, MPI *Bptr, MPI *Cptr)
/*
 * *Aptr, *Bptr, *Cptr, *Eptr  are MPI'S, *Cptr > 0, *Bptr >= 0,
 * Eptr = (*Aptr)^ *Bptr (mod *Cptr).
 */ 
{
	MPI *X, *Y, *Z, *TempI;
	
	X = MOD(Aptr, Cptr);
	Y = COPYI(Bptr);
	Z = ONEI();
	while (Y->S)
        {
		while (!(Y->V[0] & (USL)1))
	        {
			TempI = Y;
			Y = INT0_(Y, (USL)2);
			FREEMPI(TempI);
			TempI = X;
			X = MULTI(X, X);
			FREEMPI(TempI);
			TempI = X;
			X = MOD0(X, Cptr);
			FREEMPI(TempI);
		} 
		TempI = Y;
		Y = SUB0_I(Y, (USL)1);
		FREEMPI(TempI);
		TempI = Z;
		Z = MULTM(Z, X, Cptr);
		FREEMPI(TempI);
	}
	FREEMPI(X);
	FREEMPI(Y);
	return (Z);
}

MPI *MPOWER_M(MPI *Aptr, USL b, MPI *Cptr)
/*
 * *Aptr, *Cptr, *Eptr  are MPI'S, *Cptr > 0, b an unsigned int.
 * *Eptr = (*Aptr)^ b (mod *Cptr). 
 */ 
{
	MPI *X, *Z, *TempI;
	unsigned long y;
	
	X = MOD(Aptr, Cptr);
	y = b;
	Z = ONEI();
	while (y)
        {
		while (y % 2 == 0)
	        {
			y = y / 2;
			TempI = X;
			X = MULTI(X, X);
			FREEMPI(TempI);
			TempI = X;
			X = MOD0(X, Cptr);
			FREEMPI(TempI);
		} 
		y = y - 1;
		TempI = Z;
		Z = MULTM(Z, X, Cptr);
		FREEMPI(TempI);
	}
	FREEMPI(X);
	return (Z);