liptimer.c

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

/*
   touch lippar.h
   compile lip.c to produce lip.o
   compile and run timer.c (using lip.o), produces non-empty lippar.h
   recompile lip.c to produce new lip.o
*/


#include <stddef.h>
#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <signal.h>
#include <sys/time.h>
#include <sys/resource.h>
#include "lip.h"


/* do NUMBER timings of integers of size FIRST + i * INCR bits,
for 0 <= i < NUMBER */

long FIRST = 59;
long INCR = 120;
long NUMBER = 5;

#define ITER	10000		/* # of iterations for timings */
long SEED=11L;			/* for randomness */

/* do not worry about anything below this point */


#ifndef KAR_DEPTH
#define KAR_DEPTH 20
#endif
long K_M_C = 20;
long K_S_C = 20;

/*Function prototypes for static internal functions.*/


static void zhalt(
	char *c
	);

#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; \
}


/* for long division */
static double log10rad = -1.0;
static double epsilon;
static double fradix = (double)RADIX;
static double fudge = -1.0;
static double fudge1 = -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 convenience */
static long oner[] = {1, 1, 1};
static long glosho[] = {1, 1, 0};
static verylong one = &oner[1];
/* for karatsuba */
static verylong kar_mem[5*KAR_DEPTH];
static long kar_mem_initialized = 0;


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



/* SINGLE_MUL=0, PLAIN=0, KARAT=1 */
void
zmul3(
        verylong a,
        verylong b,
        verylong *c
        )
{       /* output not input */
        register long aneg;
        register long bneg;
        verylong olda;
        verylong oldb;

        if (ALLOCATE && (!a || !b))
        {
                zzero(c);
                return;
        }
        if (a == b)
        {
                zsq(a, c);
                return;
        }
        if (!kar_mem_initialized)
        {
                kar_mem_initialized = 1;
                for (aneg = (5 * KAR_DEPTH) - 1; aneg >= 0; aneg--)
                        kar_mem[aneg] = (verylong) 0;
        }
        olda = a;
        oldb = b;
        if (aneg = (*a < 0))
                a[0] = -a[0];
        if (bneg = (*b < 0))
                b[0] = -b[0];
        if (*a > *b)
                kar_mul3(a, b, c, (long) 0);
        else
                kar_mul3(b, a, c, (long) 0);
        if (aneg != bneg && ((*c)[1] || (*c)[0] != 1))
                (*c)[0] = -(*c)[0];
        if (aneg)
                olda[0] = -olda[0];
        if (bneg)
                oldb[0] = -oldb[0];
}


void
zsq3(
        verylong a,
        verylong *c
        )
{       /* output is not input */
        register long aneg;

        if (ALLOCATE && !a)
        {
                zzero(c);
                return;
        }
        if (!kar_mem_initialized)
        {
                kar_mem_initialized = 1;
                for (aneg = (5 * KAR_DEPTH) - 1; aneg >= 0; aneg--)
                        kar_mem[aneg] = (verylong) 0;
        }
        if (aneg = (*a < 0))
                a[0] = -a[0];
        kar_sq3(a, c, (long) 0);
        if (aneg)
                a[0] = -a[0];
}


#define zaddmulp3(  a,  b,  d,  t) { \
	register unsigned long aa= (*a) + (*t); \
        register unsigned long b1 = b & RADIXROOTM; \
        register unsigned long d1 = d & RADIXROOTM; \
        register unsigned long bd,b1d1,m; \
        register unsigned long ld = (d>>NBITSH); \
        register unsigned long lb = (b>>NBITSH); \
 \
        bd=lb*ld; \
        b1d1=b1*d1; \
        m=(lb+b1)*(ld+d1) - bd - b1d1; \
        aa += ( b1d1+ ((m&RADIXROOTM)<<NBITSH)); \
        (*t) = (aa >> NBITS) + bd + (m>>NBITSH); \
        (*a) = aa & RADIXM; \
}


#define zaddmul3(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--) \
        { \
                zaddmulp3(lama, *lamb, lams, &lamcarry); \
        /* Be careful, the last lama is unnormalized */ \
                lama++; \
                lamb++; \
        } \
        *lama += lamcarry; \
}


kar_mul3(
	verylong a,
	verylong b,
	verylong *c,
	long shi
	)
{
	register long al;
	register long hal;
	register long i;
	register long restoreb0 = b[0];
	register verylong pc;
	register long bbig = 1;
	verylong *a0;
	verylong *a1;
	verylong *a2;
	verylong *a3;
	verylong *a4;

	zsetlength(c, (hal = (al = a[0]) + (i = b[0])), "in kar_mul, third argument");
	if ((shi >= (5 * KAR_DEPTH)) || (al < K_M_C) || (i < K_M_C))
	{
		pc = &(*c)[1];
		for (i = hal; i > 0; i--)
			*pc++ = 0;
		pc = &(*c)[1];
		if (al <= *b)
			for (i = al; i; i--)
			{
				zaddmul3(*(++a), pc++, b);
			}
		else
			for (i = *b; i; i--)
			{
				zaddmul3(*(++b), pc++, a);
			}
		while ((hal > 1) && (!((*c)[hal])))
			hal--;
		(*c)[0] = hal;
		return;
	}
	a0 = &(kar_mem[shi]);
	a1 = &(kar_mem[shi + 1]);
	a2 = &(kar_mem[shi + 2]);
	a3 = &(kar_mem[shi + 3]);
	a4 = &(kar_mem[shi + 4]);
	hal = ((al + 1) >> 1);
	zsetlength(a0, al, "in kar_mul, locals\n");
	zsetlength(a1, al, "");
	zsetlength(a2, al, "");
	zsetlength(a3, al + hal, "");
	zsetlength(a4, al + 2, "");
	i = hal;
	while ((i > 1) && (!(a[i])))
		i--;
	a[0] = i;
	if (hal >= b[0])
		bbig = 0;
	else
	{
		i = hal;
		while ((i > 1) && (!(b[i])))
			i--;
		b[0] = i;
	}
	for (i = hal + 1; i <= al; i++)
		(*a1)[i - hal] = a[i];
	(*a1)[0] = al - hal;
	if (bbig)
	{
		for (i = hal + 1; i <= restoreb0; i++)
			(*a3)[i - hal] = b[i];
		(*a3)[0] = restoreb0 - hal;
	}
	kar_mul3(a, b, a4, shi + 5);
	zadd(a, (*a1), a0);
	a[0] = al;
	if (bbig)
	{
		kar_mul3((*a1), (*a3), c, shi + 5);
		zadd(b, (*a3), a2);
		b[0] = restoreb0;
		kar_mul3((*a0), (*a2), a3, shi + 5);
	}
	else
		kar_mul3((*a0), b, a3, shi + 5);
	zsubpos((*a3), (*a4), a3);
	if (bbig)
		zsubpos((*a3), *c, a3);
	zlshift((*a3), hal * NBITS, a3);
	hal <<= 1;
	if (bbig)
	{
		for (i = (*c)[0]; i; i--)
			(*c)[i + hal] = (*c)[i];
		for (i = hal; i > (*a4)[0]; i--)
			(*c)[i] = 0;
		for (; i; i--)
			(*c)[i] = (*a4)[i];
		(*c)[0] += hal;
	}
	else
	{
		for (i = (*a4)[0]; i >= 0; i--)
			(*c)[i] = (*a4)[i];
	}
	zadd(*c, (*a3), c);
}


#define zaddmulsq3(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--) \
        { \
                zaddmulp3(lsqa, *lsqb, lsqs, &lsqcarry); \
                lsqa++; \
                lsqb++; \
        } \
        *lsqa += lsqcarry; \
/* Be careful, the last lama is unnormalized */ \
}


#define zaddmulpsq3(_a, _b, _t) \
{ \
        register long lb = (_b); \
        register long b1 = (_b) & RADIXROOTM; \
        register long aa = *(_a) + b1 * b1; \
 \
        b1 = (b1 * (lb >>= NBITSH) << 1) + (aa >> NBITSH); \
        aa = (aa & RADIXROOTM) + ((b1 & RADIXROOTM) << NBITSH); \
        *(_t) = lb * lb + (b1 >> NBITSH) + (aa >> NBITS); \
        *(_a) = (aa & RADIXM); \
}



kar_sq3(
	verylong a,
	verylong *c,
	long shi
	)
{
	register long al;
	register long hal;
	register long i;
	register verylong pc;
	verylong *a0;
	verylong *a1;
	verylong *a2;

	zsetlength(c, (i = ((al = a[0]) << 1)), "in kar_sq, second argument");
	if ((shi >= (3 * KAR_DEPTH)) || (al < K_S_C))
	{
		register unsigned long uncar;
		long carry = 0;
		pc = &(*c)[1];
		for (; i > 0; i--)
			*pc++ = 0;
		for (hal = 1; hal <= al; hal++)
		{
			i += 2;
			{
				zaddmulsq3(al - hal, &((*c)[i]), &(a[hal]));
			}
			uncar = ((*c)[i - 1] << 1) + carry;
			(*c)[i - 1] = uncar & RADIXM;
			uncar = ((*c)[i] << 1) + (uncar >> NBITS);
			{
				zaddmulpsq3(&(*c)[i - 1], a[hal], &carry);
			}
			uncar += carry;
			carry = uncar >> NBITS;
			(*c)[i] = uncar & RADIXM;
		}
		while ((i > 1) && (!((*c)[i])))
			i--;
		(*c)[0] = i;
		return;
	}
	a0 = &(kar_mem[shi]);
	a1 = &(kar_mem[shi + 1]);
	a2 = &(kar_mem[shi + 2]);
	hal = ((al + 1) >> 1);
	zsetlength(a0, al + hal + 2, "in kar_sq, locals\n");
	zsetlength(a1, al + 2, "");
	zsetlength(a2, al, "");
	i = hal;
	while ((i > 1) && (!(a[i])))
		i--;
	a[0] = i;
	for (i = hal + 1; i <= al; i++)
		(*a0)[i - hal] = a[i];
	(*a0)[0] = al - hal;
	kar_sq3(a, a1, shi + 3);
	zadd(a, (*a0), a2);
	kar_sq3((*a0), c, shi + 3);
	a[0] = al;
	kar_sq3((*a2), a0, shi + 3);
	zsubpos((*a0), (*a1), a0);
	zsubpos((*a0), *c, a0);
	zlshift((*a0), hal * NBITS, a0);
	hal <<= 1;
	for (i = (*c)[0]; i; i--)
		(*c)[i + hal] = (*c)[i];
	for (i = hal; i > (*a1)[0]; i--)
		(*c)[i] = 0;
	for (; i; i--)
		(*c)[i] = (*a1)[i];
	(*c)[0] += hal;
	zadd(*c, (*a0), c);
}


/* SINGLE_MUL=0, PLAIN=0, KARAT=0 */

void
zmul0(
        verylong a,
        verylong b,
        verylong *c
        )
{       /* output not input */
        register long aneg;
        register long bneg;
        verylong olda;
        verylong oldb;

/*
printf("zmul0\n"); fflush(stdout);
*/
        if (ALLOCATE && (!a || !b))
        {
                zzero(c);
                return;
        }
        if (a == b)
        {
                zsq(a, c);
                return;
        }
        if (!kar_mem_initialized)
        {
                kar_mem_initialized = 1;
                for (aneg = (5 * KAR_DEPTH) - 1; aneg >= 0; aneg--)
                        kar_mem[aneg] = (verylong) 0;
        }
        olda = a;
        oldb = b;
        if (aneg = (*a < 0))
                a[0] = -a[0];
        if (bneg = (*b < 0))
                b[0] = -b[0];
        if (*a > *b)
                kar_mul0(a, b, c, (long) 0);
        else
                kar_mul0(b, a, c, (long) 0);
        if (aneg != bneg && ((*c)[1] || (*c)[0] != 1))
                (*c)[0] = -(*c)[0];
        if (aneg)
                olda[0] = -olda[0];
        if (bneg)
                oldb[0] = -oldb[0];
}


void
zsq0(
        verylong a,
        verylong *c
        )
{       /* output is not input */
        register long aneg;

        if (ALLOCATE && !a)
        {
                zzero(c);
                return;
        }
        if (!kar_mem_initialized)
        {
                kar_mem_initialized = 1;
                for (aneg = (5 * KAR_DEPTH) - 1; aneg >= 0; aneg--)
                        kar_mem[aneg] = (verylong) 0;
        }
        if (aneg = (*a < 0))
                a[0] = -a[0];
        kar_sq0(a, c, (long) 0);
        if (aneg)
                a[0] = -a[0];
}

static double fradix_inv0 = 1.0 / RADIX; /* Global constant */
#define zaddmulp0(_a, _b, _d, _t) \
{ \
        register at = *(_a) + *(_t); \
        register long aa = (at + (_b) * (_d)) & RADIXM; \
 \
        *(_t) = (long) (0.25 + fradix_inv0 * (((double) (at - aa)) \
                       + ((double) (_b)) * ((double) (_d)))); \
        *(_a) = aa; \
}


#define zaddmul0(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--) \
        { \
                zaddmulp0(lama, *lamb, lams, &lamcarry); \
        /* Be careful, the last lama is unnormalized */ \
                lama++; \
                lamb++; \
        } \
        *lama += lamcarry; \
}


kar_mul0(
	verylong a,
	verylong b,
	verylong *c,
	long shi
	)
{
	register long al;
	register long hal;
	register long i;
	register long restoreb0 = b[0];
	register verylong pc;
	register long bbig = 1;
	verylong *a0;
	verylong *a1;
	verylong *a2;
	verylong *a3;
	verylong *a4;

/*
printf("karmul0\n"); fflush(stdout);
*/
	zsetlength(c, (hal = (al = a[0]) + (i = b[0])), "in kar_mul, third argument");
	if ((shi >= (5 * KAR_DEPTH)) || (al < K_M_C) || (i < K_M_C))
	{
		pc = &(*c)[1];
		for (i = hal; i > 0; i--)
			*pc++ = 0;
		pc = &(*c)[1];
		if (al <= *b)
			for (i = al; i; i--)
			{
				zaddmul0(*(++a), pc++, b);
			}
		else
			for (i = *b; i; i--)
			{
				zaddmul0(*(++b), pc++, a);
			}
		while ((hal > 1) && (!((*c)[hal])))
			hal--;
		(*c)[0] = hal;
		return;
	}
	a0 = &(kar_mem[shi]);
	a1 = &(kar_mem[shi + 1]);
	a2 = &(kar_mem[shi + 2]);
	a3 = &(kar_mem[shi + 3]);
	a4 = &(kar_mem[shi + 4]);
	hal = ((al + 1) >> 1);
	zsetlength(a0, al, "in kar_mul, locals\n");
	zsetlength(a1, al, "");
	zsetlength(a2, al, "");
	zsetlength(a3, al + hal, "");
	zsetlength(a4, al + 2, "");
	i = hal;
	while ((i > 1) && (!(a[i])))
		i--;
	a[0] = i;
	if (hal >= b[0])
		bbig = 0;
	else
	{
		i = hal;
		while ((i > 1) && (!(b[i])))