lip.c

应用密码学这本书的源代码

 \
	lmmd = (zninv * lmmi) & RADIXM; \
	for (lmmi = *lmmpb++; lmmi > 0; lmmi--) \
	{ \
		zaddmulp(lmmpa, *lmmpb, lmmd, &zmmtemp); \
		lmmpa++; \
		*lmmpb++; \
	} \
	if (((*lmmpa += zmmtemp) & RADIX) > 0) \
	{ \
		(*lmmpa++) &= RADIXM; \
		(*lmmpa)++; \
	} \
}

#define zaddmul(ams, ama, amb) \
{ \
	register long lami; \
	register long lams = (ams); \
	register verylong lama = (ama); \
	register verylong lamb = (amb); \
	long lamcarry = 0; \
 \
	for (lami = (*lamb++); lami > 0; lami--) \
	{ \
		zaddmulp(lama, *lamb, lams, &lamcarry); \
	/* Be careful, the last lama is unnormalized */ \
		lama++; \
		lamb++; \
	} \
        *lama += lamcarry; \
}


#define zaddmulsq(sql, sqa, sqb) \
{ \
	register long lsqi = (sql); \
	register long lsqs = *(sqb); \
	register verylong lsqa = (sqa); \
	register verylong lsqb = (sqb); \
	long lsqcarry = 0; \
 \
	lsqb++; \
	for (; lsqi > 0; lsqi--) \
	{ \
		zaddmulp(lsqa, *lsqb, lsqs, &lsqcarry); \
		lsqa++; \
		lsqb++; \
	} \
	*lsqa += lsqcarry; \
/* Be careful, the last lama is unnormalized */ \
}

static void
zsubmul(
	long r,
	verylong a,
	verylong b
	)
{
	register long rd = RADIX - r;
	register long i;
	long carry = RADIX;

	for (i = (*b++); i > 0; i--)
	{
		zaddmulp(a, *b, rd, &carry);
		a++;
		carry += RADIXM - (*b++);
	}
	*a += carry - RADIX;	/* unnormalized */
}
#endif


#ifdef ALPHA50
#define zaldivaddmulp(  a,  b,  d,  t) { \
        register long b1 = b & ALPHA50RADIXROOTM; \
        register long d1 = d & ALPHA50RADIXROOTM; \
        register long bd,b1d1,m,aa= (*a) + (*t); \
        register long ld = (d>>ALPHA50NBITSH); \
        register long lb = (b>>ALPHA50NBITSH); \
 \
        bd=lb*ld; \
        b1d1=b1*d1; \
        m=(lb+b1)*(ld+d1) - bd - b1d1; \
	aa += ( b1d1+ ((m&ALPHA50RADIXROOTM)<<ALPHA50NBITSH)); \
        (*t) = (aa >> ALPHA50NBITS) + bd + (m>>ALPHA50NBITSH); \
        (*a) = aa & ALPHA50RADIXM; \
}
#endif


#define zaddmulone(ama, amb) \
{ \
	register long lami; \
	register long lams = 0; \
	register verylong lama = (ama); \
	register verylong lamb = (amb); \
 \
	lams = 0; \
	for (lami = (*lamb++); lami > 0; lami--) \
	{ \
		lams += (*lama + *lamb++); \
		*lama++ = lams & RADIXM; \
		lams >>= NBITS; \
	} \
	*lama += lams; \
}

#define zalphaaddmulone(ama, amb) \
{ \
        register long lami; \
        register long lams = 0; \
        register verylong lama = (ama); \
        register verylong lamb = (amb); \
 \
        lams = 0; \
        for (lami = (*lamb++); lami > 0; lami--) \
        { \
                lams += (*lama + *lamb++); \
                *lama++ = lams & ALPHA50RADIXM; \
                lams >>= ALPHA50NBITS; \
        } \
        *lama += lams; \
}


/*
	zdiv21 returns quot, numhigh so

	quot = (numhigh*RADIX + numlow)/denom;
	numhigh  = (numhigh*RADIX + numlow)%denom;

Assumes 0 <= numhigh < denom < RADIX and 0 <= numlow < RADIX.
*/


#ifdef ALPHA

#define zdiv21(numhigh, numlow, denom, deninv, quot) \
{ \
	register long lq21 = (long) (((fradix * (double) (numhigh)) \
			 + (double) (numlow)) * (deninv)); \
	register long lr21; \
	register long dd; \
	long lprodhigh = 0; \
	long lprodlow = 0; \
 \
	zaddmulp(&lprodlow, lq21, denom, &lprodhigh); \
	dd = numhigh-lprodhigh; \
	lprodhigh=0; /* this line here because of compiler error */ \
	if (dd > 1) { \
		register long correction = (((double)dd)*(double)RADIX + numlow - lprodlow)*deninv; \
		lq21 += correction; \
		lprodlow = 0; \
		zaddmulp(&lprodlow, lq21, denom, &lprodhigh); \
		dd = numhigh-lprodhigh; \
	} \
	else if (dd < -1) { \
		register long correction = (((double)dd)*(double)RADIX + numlow - lprodlow)*deninv; \
		lq21 += correction; \
		lprodlow = 0; \
		zaddmulp(&lprodlow, lq21, denom, &lprodhigh); \
		dd = numhigh-lprodhigh; \
	} \
	lr21 = RADIX * dd + (numlow - lprodlow); \
	if (dd=lr21*(deninv)) { \
		lq21 += dd; \
		lr21 -= denom*dd; \
	} \
	if (lr21 < 0) \
	{ \
		while (lr21 < 0) { \
			lq21--; \
			lr21 += denom; \
		}; \
	} \
	else \
	{ \
		while (lr21 >= denom) \
		{ \
			lr21 -= denom; \
			lq21++; \
		}; \
	} \
	quot = lq21; \
	numhigh = lr21; \
}
#else

#ifdef ALPHA50

#define zdiv21(numhigh, numlow, denom, deninv, quot) \
{ \
        register long lr21; \
        register long lq21 = (long) (((alpha50fradix * (double) (numhigh)) \
                         + (double) (numlow)) * (deninv)); \
        long lprodhigh = 0; \
        long lprodlow = 0; \
 \
        zaldivaddmulp(&lprodlow, lq21, denom, &lprodhigh); \
        lr21 = ALPHA50RADIX *(numhigh - lprodhigh) + (numlow - lprodlow); \
        if (lr21 < 0) \
        { \
                do \
                { \
                        lq21--; \
                } while ((lr21 += denom) < 0); \
        } \
        else \
        { \
                while (lr21 >= denom) \
                { \
                        lr21 -= denom; \
                        lq21++; \
                }; \
        } \
        quot = lq21; \
        numhigh = lr21; \
}

#else

#ifndef WIN32
#define zdiv21(numhigh, numlow, denom, deninv, quot) \
{ \
        register long lr21; \
        register long lq21 = (long) (((fradix * (double) (numhigh)) \
                         + (double) (numlow)) * (deninv)); \
        long lprodhigh = 0; \
        long lprodlow = 0; \
 \
        zaddmulp(&lprodlow, lq21, denom, &lprodhigh); \
        lr21 = RADIX * (numhigh - lprodhigh) + (numlow - lprodlow); \
        if (lr21 < 0) \
        { \
                do \
                { \
                        lq21--; \
                } while ((lr21 += denom) < 0); \
        } \
        else \
        { \
                while (lr21 >= denom) \
                { \
                        lr21 -= denom; \
                        lq21++; \
                }; \
        } \
        quot = lq21; \
        numhigh = lr21; \
}
#else

#define zdiv21(numhigh, numlow, denom, deninv, quot) \
{ \
        register long lr21; \
        register long lq21 = (long) (((fradix * (double) (numhigh)) \
                         + (double) (numlow)) * (deninv)); \
 /* Following works in many two's complement architectures. */ \
        lr21 = (numhigh << NBITS) + numlow - lq21 * denom; \
        if (lr21 < 0) \
        { \
                do \
                { \
                        lq21--; \
                } while ((lr21 += denom) < 0); \
        } \
        else \
        { \
                while (lr21 >= denom) \
                { \
                        lr21 -= denom; \
                        lq21++; \
                }; \
        } \
        quot = lq21; \
        numhigh = lr21; \
}

#endif


#endif
#endif


long
zmulmods(
	long a,
	long b,
	long n
	)
{
	register long nn=n;
	register long na= (a >= nn||a<0) ? (a % nn) : a;
	register long nb= (b >= nn||b<0 ) ? (b % nn) : b;
	register long lqmul = (long) (((double)na) * ((double)nb) / ((double)nn));
#ifdef PLAIN_OR_KARAT
	register long lr;
	long lprodhigh1 = 0;
	long lprodhigh2 = 0;
	long lprodlow1 = 0;
	long lprodlow2 = 0;

	zaddmulp(&lprodlow1, na, nb, &lprodhigh1);
	zaddmulp(&lprodlow2, lqmul, nn, &lprodhigh2);
# ifdef ALPHA
	lr = lprodhigh1 - lprodhigh2;
	if (lr > 2) {
		double correction= (((double)lr)*(double)RADIX + lprodlow1 - lprodlow2)/((double) nn);
		lqmul -= (long) correction;
		lprodhigh1 = 0; lprodhigh2 = 0; lprodlow1 = 0; lprodlow2 = 0;
		zaddmulp(&lprodlow1, na, nb, &lprodhigh1);
		zaddmulp(&lprodlow2, lqmul, nn, &lprodhigh2);
	}
	else if (lr < -2) {
		double correction= (((double)lr)*(double)RADIX + lprodlow1 - lprodlow2)/((double) nn);

		lqmul = lqmul+(long) correction;
		lprodhigh1 = 0; lprodhigh2 = 0; lprodlow1 = 0; lprodlow2 = 0;
		zaddmulp(&lprodlow1, na, nb, &lprodhigh1);
		zaddmulp(&lprodlow2, lqmul, nn, &lprodhigh2);
	}
# endif
	lr = RADIX * (lprodhigh1 - lprodhigh2) + (lprodlow1 - lprodlow2);
#else
/* Many machines compute the following modulo 2^32, which is OK */
	register long lr = na * nb - lqmul * nn;
#endif
	while (lr >= nn)
		lr -= nn;
	while (lr < 0)
		lr +=  nn;
	return lr;
}


#define MulLo(rres,a,b) \
{ \
   register double _x = ((double) a) * ((double) b); \
   double _t = _x+(4503599627370496.0); \
   register long *_lo = (long *)&_t; \
   rres =  _lo[LO_WD];\
}

long
zmulmod26(
	long a,
	long b,
	long n,
	double bninv
	)
{
        register long lqmul = (long)(((double)a)*bninv);
        register long p1;
	register long p2;
        MulLo(p1,a,b);
        MulLo(p2,lqmul,n);
        p1 -= p2;
        if (p1 < 0)
           p1 += n; 
        else if (p1 >= n)
	   p1 -= n;
	return(p1);
} 

extern void *calloc();

/* global variables */

/* for long division */
static double log10rad = -1.0;
static double log16rad = -1.0;
static double epsilon = 0.0;
static double fradix = (double)RADIX;
static double fudge = -1.0;
#ifdef ALPHA
static double fudge2 = -1.0;
#endif
#ifdef ALPHA50
static double alpha50fudge = -1.0;
static double alpha50fradix = (double) ALPHA50RADIX;
#endif
/* for random generator */
static verylong zseed = 0;
static verylong zranp = 0;
static verylong zprroot = 0;
/* for small prime genaration */
static short *lowsieve = 0;
static short *movesieve = 0;
static long pindex = 0;
static long pshift = -1;
static long lastp = 0;
/* for convenience */
static long oner[] = {1, 1, 1};
static long glosho[] = {1, 1, 0};
static verylong one = &oner[1];
/* for m_ary exponentiation */
static verylong **exp_odd_powers = 0;
/* for karatsuba */
static verylong kar_mem[5*KAR_DEPTH];
static long kar_mem_initialized = 0;

#ifndef WIN32
#include <sys/times.h>
#include <limits.h>
#ifndef CLK_TCK
#define CLK_TCK 60
#endif
#endif

double 
getutime()
{
#ifdef OKWIN32

      FILETIME create_ft,exit_ft,kernel_ft,user_ft;
      GetProcessTimes(GetCurrentProcess(),
                      &create_ft,&exit_ft,&kernel_ft,&user_ft);

      return ((double)user_ft.dwLowDateTime)/1e7;

#elif WIN32
	return((double)(clock()/CLK_TCK));
#else
	struct tms tbuf; times(&tbuf); return(1.0*tbuf.tms_utime/CLK_TCK);
#endif
}

double 
getstime()
{
#ifdef OKWIN32
      FILETIME create_ft,exit_ft,kernel_ft,user_ft;
      GetProcessTimes(GetCurrentProcess(),
                      &create_ft,&exit_ft,&kernel_ft,&user_ft);

      return ((double)kernel_ft.dwLowDateTime)/1e7;

#elif WIN32
	return((double)(clock()/CLK_TCK));
#else
	struct tms tbuf; times(&tbuf); return(1.0*tbuf.tms_stime/CLK_TCK);
#endif
}

double 
gettime()
{
#ifdef OKWIN32
      FILETIME create_ft,exit_ft,kernel_ft,user_ft;
      GetProcessTimes(GetCurrentProcess(),
                      &create_ft,&exit_ft,&kernel_ft,&user_ft);

      return ((double)(kernel_ft.dwLowDateTime+user_ft.dwLowDateTime))/1e7;
#elif WIN32
	return((double)(clock()/CLK_TCK));
#else
  struct tms tbuf;  times(&tbuf);
  return(1.0*(tbuf.tms_utime+tbuf.tms_stime)/CLK_TCK);
#endif
}

static void
hidetimer(
	FILE *f,
	long what
	)
{
	static double keep_time = 0.0;
	if (what)
	{
		fprintf(f,"%8.5f sec.\n",gettime()-keep_time);
		fflush(f);
	}
	else
		keep_time = gettime();
}

void
starttime()
{
	hidetimer(stderr,0);
}

void
printtime(
	FILE *f
	)
{
	hidetimer(f,1);
}

static void
zhalt(
        char *c
	)
{
#ifdef NOHALT
	fprintf(stderr,"error:\n   %s\ncontinue...\n",c);
	fflush(stderr);
#else
	fprintf(stderr,"fatal error:\n   %s\nexit...\n",c);
	fflush(stderr);
	(void)exit((int)0);
#endif
}

void
zsetlength(
	verylong *v,
	long len,
	char *str
	)
{
	verylong x = *v;

	if (x)
	{
		if (len <= x[-1])
			return;
#if (PRT_REALLOC>0)
			fprintf(stderr,"%s reallocating to %ld\n", str, len);
			fflush(stderr);
#endif
		x[-1] = len;
		if (!(x = (verylong)realloc((void*)(&(x[-1])), (size_t)(len + 2) * SIZEOFLONG)))
		{
			fprintf(stderr,"%d bytes realloc failed\n", ((int)len + 2) * SIZEOFLONG);
			zhalt("reallocation failed in zsetlength");
		}
	}
	else if (len >= 0)
	{
		if (len < SIZE)
			len = SIZE;
#if (PRT_REALLOC>0)
			fprintf(stderr,"%s allocating to %ld\n", str, len);
			fflush(stderr);
#endif
		if (!(x = (verylong)calloc((size_t)(len + 2), (size_t)SIZEOFLONG)))
		{
			fprintf(stderr,"%d bytes calloc failed\n", ((int)len + 2) * SIZEOFLONG);
			zhalt("allocation failed in zsetlength");
		}
		x[0] = len;
		x[1] = 1;
		x[2] = 0;
	}
	else
		zhalt("negative size allocation in zsetlength");
	x++;
	*v = x;
}

void
zfree(
	verylong *x
	)
{
	if (!(*x))
		return;
	{
		verylong y = (*x - 1);
		free((void*)y);
		*x = 0;
		return;
	}
}

double 
zdoub(
	verylong n
	)
{
	double res;
	register long i;

	if (ALLOCATE && !n)
		return ((double) 0);
	if ((i = n[0]) < 0)
		i = -i;
	res = (double) (n[i--]);
	for (; i; i--)
		res = res * fradix + (double) (n[i]);
	if (n[0] > 0)
		return (res);
	return (-res);
}

void
zstart()
{
	double local_one = (double) 1;
	double local_half = 1 / (double) 2;
	epsilon = local_one;
	fudge = local_one + epsilon;

 /* the following loop is sick, but we have to: */
 /* comparing one to one + epsilon does not have */
 /* the desired effect on some machines, because */
 /* one + epsilon is Not cast to double (as it should) */
	while (local_one != fudge)
	{
		epsilon = epsilon * local_half;
		fudge = local_one + epsilon;
	}
	epsilon += epsilon;
#ifndef ALPHA_OR_ALPHA50
	if ((epsilon * RADIX) > local_one)
	{
		fprintf(stderr,
			"zstart failure: recompile with smaller NBITS\n");
		exit(0);
	}
#endif
	epsilon *= 3;
	fudge = fradix + epsilon * fradix;
#ifdef ALPHA
	fudge2 = fradix - epsilon * fradix;
#endif
#ifdef ALPHA50
	alpha50fudge = alpha50fradix + epsilon * alpha50fradix;
#endif
}


void
zzero(
	verylong *aa
	)
{
	if (!(*aa)) zsetlength(aa, (long)SIZE, "in zzero, first argument");
	(*aa)[0] = (long) 1;
	(*aa)[1] = (long) 0;
}

void
zone(
	verylong *aa
	)
{
	if (!(*aa)) zsetlength(aa, (long)SIZE, "in zone, first argument");
	(*aa)[0] = (long) 1;
	(*aa)[1] = (long) 1;
}

void
zcopy(
	verylong a,
	verylong *bb
	)
{
	register long i;
	verylong b = *bb;

	if (ALLOCATE && !a)
	{
		zzero(bb);
		return;
	}
	if (a != b)
	{
		if ((i = *a) < 0)
			i = (-i);
		zsetlength(&b, i, "in zcopy, second argument");
		*bb = b;
		for (; i >= 0; i--)
			*b++ = *a++;
	}
}

void
zintoz(
	long d,
	verylong *aa
	)
{
	register long i = 0;
	register long anegative = 0;
	verylong a = *aa;

	if (!(*aa)) zsetlength(&a, (long)SIZE, "in zintoz, second argument");
	*aa = a;
	if (d < 0)
	{
		anegative = 1;
		d = -d;
	}
	a[1] = 0;
	while (d > 0)
	{
		a[++i] = d & RADIXM;
		d >>= NBITS;
	}
	if (i > 0)
		a[0] = i;
	else
		a[0] = 1;
	if (anegative)
		a[0] = (-a[0]);
}

void
zuintoz(
	unsigned long d,
	verylong *aa
	)
{
	register long i = 0;
	verylong a = *aa;

	if (!(*aa)) zsetlength(&a, (long)SIZE, "in zuintoz, second argument");
	*aa = a;
	a[1] = 0;
	while (d)
	{
		a[++i] = (long)(d & RADIXM);
		d >>= NBITS;
	}
	if (i > 0)
		a[0] = i;
	else
		a[0] = 1;
}

long 
ztoint(
	verylong a
	)
{
	register long d;
	register long sa;

	if (ALLOCATE && !a)
		return (0);
	if ((sa = *a) < 0)
		sa = -sa;
	d = *(a += sa);
	while (--sa)
	{
		d <<= NBITS;
		d += *(--a);
	}
	if ((*(--a)) < 0)
		return (-d);
	return (d);
}

unsigned long 
ztouint(
	verylong a
	)
{
	register unsigned long d;
	register long sa;

	if (ALLOCATE && !a)
		return (0);
	if ((sa = *a) < 0)