liptimer.c

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

	{
		for (i = (*a4)[0]; i >= 0; i--)
			(*c)[i] = (*a4)[i];
	}
	zadd(*c, (*a3), c);
}


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


#define zaddmulpsq2(_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_sq2(
	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;
			{
				zaddmulsq2(al - hal, &((*c)[i]), &(a[hal]));
			}
			uncar = ((*c)[i - 1] << 1) + carry;
			(*c)[i - 1] = uncar & RADIXM;
			uncar = ((*c)[i] << 1) + (uncar >> NBITS);
			{
				zaddmulpsq2(&(*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_sq2(a, a1, shi + 3);
	zadd(a, (*a0), a2);
	kar_sq2((*a0), c, shi + 3);
	a[0] = al;
	kar_sq2((*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);
}


#define lower_radix1(in,out) { \
register long i=1,j=0,rb=0,nrb=0,bits; \
register verylong pin=in; \
long size; \
 \
zsetlength( (&out), (long) (30.0/26.0 * in[0] + 1),""); \
if ((size= *pin)==2) { \
        pin++; \
        out[1] = *pin & ((1L<<26)-1); \
        out[2] = *pin >> 26; \
        pin++; \
        out[2] |= (*pin & ((1L << (2*26-30))-1L)) << (30 - \
         26); \
        bits = *pin >> (2*26-30); \
        if (bits) { \
                out[3] = bits; \
                out[0]=3; \
                } \
        else \
                out[0]=2; \
        } \
else { \
        for (j=0; i <= size ; j++) { \
                if (nrb >= 26) { \
                        bits = rb & ((1L<<26)-1); \
                        rb >>= 26; \
                        nrb -= 26; \
                        } \
                else { \
                        bits = rb | ((((1L << (26 - nrb)) -1L) & \
                         pin[i]) << nrb); \
                        rb = pin[i] >> (26 - nrb); \
                        nrb += (30-26); \
                        i++; \
                        } \
                out[j+1]=bits; \
                } \
        if (rb) { \
                j++; \
                if (rb < (1L<<26)) \
                        out[j]= rb; \
                else { \
                        out[j]= rb & (( 1L << 26) -1); \
                        j++; \
                        out[j] = rb>>26; \
                        } \
                } \
        out[0]=j; \
        } \
}


void
zmul5(
        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_mul5(a, b, c, (long) 0);
        else
                kar_mul5(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
zsq5(
        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_sq5(a, c, (long) 0);
        if (aneg)
                a[0] = -a[0];
}

#define LO5(x) (((unsigned long *) &x)[0])
#define HI5(x) (((unsigned long *) &x)[1])

#define FetchHiLo5(hi,lo,x) \
{ \
   d_or_rep _l_xx; \
   _l_xx.d = x; \
   hi = _l_xx.rep[1]; \
   lo = _l_xx.rep[0]; \
}

void zaddmul5(
        long ams,
        long *ama,
        long *amb
        )
{
        register long carry = 0;
        register long i = (*amb++);
        register double dams = (double) ams;
        register double xx;
        register double yy;
        register unsigned long lo_wd, lo;
        register unsigned long hi_wd, hi;

        xx  =  ((double) (*amb))*dams + 4503599627370496.0;
        for (; i > 1; i--)
        {
                yy =  ((double) (*(++amb)))*dams +4503599627370496.0;
                FetchHiLo5(hi_wd, lo_wd, xx);
                lo = lo_wd & 0x3FFFFFF;
                hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
                lo = lo + (*ama) + carry;
                *ama = lo & 0x3FFFFFF;
                carry = hi + (lo >> 26);
                ama++;
                xx = yy;
        }

        FetchHiLo5(hi_wd, lo_wd, xx);
        lo = lo_wd & 0x3FFFFFF;
        hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
        lo = lo + (*ama) + carry;
        *ama = lo & 0x3FFFFFF;
        carry = hi + (lo >> 26);
        ama++;
        *ama += carry;
}


kar_mul5(
	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--)
			{
				zaddmul5(*(++a), pc++, b);
			}
		else
			for (i = *b; i; i--)
			{
				zaddmul5(*(++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_mul5(a, b, a4, shi + 5);
	zadd(a, (*a1), a0);
	a[0] = al;
	if (bbig)
	{
		kar_mul5((*a1), (*a3), c, shi + 5);
		zadd(b, (*a3), a2);
		b[0] = restoreb0;
		kar_mul5((*a0), (*a2), a3, shi + 5);
	}
	else
		kar_mul5((*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 ExtractHiLo5(Hi,Lo,x) \
{ \
double t=x+(4503599627370496.0); \
unsigned long *it = (unsigned long *)&t; \
Lo = it[0]; \
Hi = ((it[1]<<6)|(it[0]>>26)); \
Lo &= 0x3FFFFFF; \
Hi &= 0x3FFFFFF; \
}


#define zaddmulpsq5(a,b,t) \
{ \
   double __lx = ((double) (*a)) + ((double) b)*((double) b); \
   register long  __lhi = 0, __llo = 0;\
   ExtractHiLo5(__lhi,__llo,__lx);\
   (*a) = __llo;\
   (*t) = __lhi;\
}


void zaddmulsq5(
        long ams,
        long *ama,
        long *amb
        )
{
        register long carry = 0;
        register long i = ams;
        register double dams = (double) (*amb);
        double xx;
        register double yy;
        register unsigned long lo;
        register unsigned long hi;

        xx =  ((double) (*(++amb)))*dams + 4503599627370496.0;
        for (; i > 1; i--)
        {
                yy =  ((double) (*(++amb)))*dams + 4503599627370496.0;
                lo = LO5(xx) & 0x3FFFFFF;
                hi = ((HI5(xx)<<6)|(LO5(xx)>>26)) & 0x3FFFFFF;
                lo = lo + (*ama) + carry;
                *ama = lo & 0x3FFFFFF;
                carry = hi + (lo >> 26);
                ama++;
                xx = yy;
        }
        if (i==1)
        {
                lo = LO5(xx) & 0x3FFFFFF;
                hi = ((HI5(xx)<<6)|(LO5(xx)>>26)) & 0x3FFFFFF;
                lo = lo + (*ama) + carry;
                *ama = lo & 0x3FFFFFF;
                carry = hi + (lo >> 26);
                ama++;
        }
        *ama += carry;
}


kar_sq5(
	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;
			{
				zaddmulsq5(al - hal, &((*c)[i]), &(a[hal]));
			}
			uncar = ((*c)[i - 1] << 1) + carry;
			(*c)[i - 1] = uncar & RADIXM;
			uncar = ((*c)[i] << 1) + (uncar >> NBITS);
			{
				zaddmulpsq5(&(*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_sq5(a, a1, shi + 3);
	zadd(a, (*a0), a2);
	kar_sq5((*a0), c, shi + 3);
	a[0] = al;
	kar_sq5((*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);
}


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

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


void
zstart0()
{
        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
zdiv0(
	verylong in_a,
	verylong in_b,
	verylong *qqq,
	verylong *rrr
	)

{
	register long i;