liptimer.c

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

	register long qq;
	long sa;
	long sb;
	long sq;
	verylong p;
	verylong pc;
	long sign;
	static verylong a = 0;
	static verylong b = 0;
	static verylong c = 0;
	static verylong d = 0;
	double btopinv;
	double aux;
	verylong q = *qqq;
	verylong r = *rrr;

#ifndef START
	if (fudge < 0)
		zstart0();
#endif

	if (ALLOCATE && !in_a)
	{
		zzero(qqq);
		zzero(rrr);
		return;
	}
	if ((!in_b) || (((sb = in_b[0]) == 1) && (!in_b[1])))
	{
		zhalt("division by zero in zdiv");
		return;
	}
	zcopy(in_a,&a);

	zcopy(in_b,&b);

	sign = (*a < 0 ? 2 : 0) | (*b < 0 ? 1 : 0);
	if (*a < 0)
		(*a) = (-(*a));
	sa = (*a);
	if (*b < 0)
		(*b) = (-(*b));
	sb = (*b);
	zsetlength(&c, (i = (sa > sb ? sa : sb) + 1), "in zdiv, locals\n");
	zsetlength(&d, i, "");
	zsetlength(&q, i, "in zdiv, third argument");
	*qqq = q;
	zsetlength(&r, sb + 1, "in zdiv, fourth argument");
	*rrr = r;
	p = &b[sb];
	if ((sb == 1) && (*p < RADIX))
		zintoz(zsdiv(a, *p, &q), &r);
	else
	{
		sq = sa - sb;	/* size of quotient */
		btopinv = (double) (*p) * fradix;
		if (sb > 1)
			btopinv += (*(--p));
		btopinv *= fradix;
		if (sb > 2)
			btopinv += (*(p - 1));
		btopinv = fudge / btopinv;
		p = &a[1];
		pc = &c[0];
		for (i = sa; i > 0; i--)
			*pc++ = *p++;
		p = pc;
		sa = 1 - sb;
		for (i = (-sq); i > 0; i--)
			*p++ = 0;
		*pc = 0;
		d[1] = 0;
		p = &d[sq + 1];
		for (i = sq; i >= 0; i--)
		{

			aux = fradix * (fradix * (*pc) + (*(pc - 1))) + 1.0;
			if (i > sa)
				aux += (*(pc - 2));
			qq = (long) (btopinv * aux + 0.5);
		/* dirty, but safe. On most machines (long)(btopinv*aux) */
		/* is correct, or one too big; on some however it becomes */
		/* too small. Could change zstart, but +0.5 and a while */
		/* instead of one single if is safer */
			if (qq > 0)
			{
				if (qq >= RADIX)
					qq = RADIX-1;
				zsubmul0(qq, &c[i], &b[0]);
				while (*pc < 0)
				{
					qq--;
					{
						zaddmulone(&c[i], &b[0]);
					}
				}
			}
			pc--;
			*p-- = qq;
		}
		sb--;
		while ((sb > 0) && (!(c[sb])))
			sb--;
		sb++;
		r[0] = sb;
		p = &r[1];
		pc = &c[0];
		for (i = sb; i > 0; i--)
			*p++ = *pc++;
		if (sq < 0)
		{
			q[0] = 1;
			q[1] = d[1];
		}
		else
		{
			sq++;
			while ((sq > 1) && (!(d[sq])))
				sq--;
			q[0] = sq;
			p = &q[1];
			pc = &d[1];
			for (i = sq; i > 0; i--)
				*p++ = *pc++;
		}
	}	/* non-trivial case */

	if (sign)
	{
		if (sign <= 2)
		{
			if (!(r[1]) && (r[0] == 1))
				q[0] = -(q[0]);
			else
			{
				zadd(q, one, &q);
				q[0] = -q[0];
				if (sign == 1)
					zsub(r, b, &r);
				else
					zsub(b, r, &r);
			}
		}
		else
			znegate(&r);
	}
	*qqq = q;
	*rrr = r;

}


#define zaddmulone1(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 & ((1L<<26)-1); \
                lams >>= 26; \
        } \
        *lama += lams; \
}


static
void zsubmul1(
        long ams,
        long *ama,
        long *amb
        )
{
        register long carry = (1L<<26);
        register long i = (*amb++);
        register double dams = (double) ((1L<<26)-ams);
        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+1)))*dams + 4503599627370496.0;
                lo = LO1(xx) & 0x3FFFFFF;
                hi = ((HI1(xx)<<6)|(LO1(xx)>>26)) & 0x3FFFFFF;
                lo = lo + (*ama) + carry;
                *ama = lo & 0x3FFFFFF;
                carry = hi + (lo >> 26);
                ama++;
                carry += ((1L<<26)-1) - (*(amb++));
                xx = yy;
        }
        lo = LO1(xx) & 0x3FFFFFF;
        hi = ((HI1(xx)<<6)|(LO1(xx)>>26)) & 0x3FFFFFF;
        lo = lo + (*ama) + carry;
        *ama = lo & 0x3FFFFFF;
        carry = hi + (lo >> 26);
        ama++;
        carry += ((1L<<26)-1) - (*amb);
        *ama += carry - (1L<<26);
}


void
zstart1()
{
        double local_one = (double) 1;
        double local_half = 1 / (double) 2;
        epsilon = local_one;
        fudge1 = 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 != fudge1)
        {
                epsilon = epsilon * local_half;
                fudge1 = local_one + epsilon;
        }
        epsilon += epsilon;
        if ((epsilon * (1L<<26)) > local_one)
        {
                fprintf(stderr,
                        "zstart failure: recompile with smaller NBITS\n");
                exit(0);
        }
        epsilon *= 3;
        fudge1 = 67108864.0 + epsilon * 67108864.0;
}


void
zdiv1(
	verylong in_a,
	verylong in_b,
	verylong *qqq,
	verylong *rrr
	)

{
	register long i;
	register long qq;
	long sa;
	long sb;
	long sq;
	verylong p;
	verylong pc;
	long sign;
	static verylong a = 0;
	static verylong b = 0;
	static verylong c = 0;
	static verylong d = 0;
	double btopinv;
	double aux;
	verylong q = *qqq;
	verylong r = *rrr;

#ifndef START
	if (fudge1 < 0)
		zstart1();
#endif

	if (ALLOCATE && !in_a)
	{
		zzero(qqq);
		zzero(rrr);
		return;
	}
	if ((!in_b) || (((sb = in_b[0]) == 1) && (!in_b[1])))
	{
		zhalt("division by zero in zdiv");
		return;
	}
	zcopy(in_a,&a);

	zcopy(in_b,&b);

	sign = (*a < 0 ? 2 : 0) | (*b < 0 ? 1 : 0);
	if (*a < 0)
		(*a) = (-(*a));
	sa = (*a);
	if (*b < 0)
		(*b) = (-(*b));
	sb = (*b);
	zsetlength(&c, (i = (sa > sb ? sa : sb) + 1), "in zdiv, locals\n");
	zsetlength(&d, i, "");
	zsetlength(&q, i, "in zdiv, third argument");
	*qqq = q;
	zsetlength(&r, sb + 1, "in zdiv, fourth argument");
	*rrr = r;
	p = &b[sb];
	if ((sb == 1) && (*p < (1L<<26)))
		zintoz(zsdiv(a, *p, &q), &r);
	else
	{
		sq = sa - sb;	/* size of quotient */
		btopinv = (double) (*p) * 67108864.0;
		if (sb > 1)
			btopinv += (*(--p));
		btopinv *= 67108864.0;
		if (sb > 2)
			btopinv += (*(p - 1));
		btopinv = fudge1 / btopinv;
		p = &a[1];
		pc = &c[0];
		for (i = sa; i > 0; i--)
			*pc++ = *p++;
		p = pc;
		sa = 1 - sb;
		for (i = (-sq); i > 0; i--)
			*p++ = 0;
		*pc = 0;
		d[1] = 0;
		p = &d[sq + 1];
		for (i = sq; i >= 0; i--)
		{

			aux = 67108864.0 * (67108864.0 * (*pc) + (*(pc - 1))) + 1.0;
			if (i > sa)
				aux += (*(pc - 2));
			qq = (long) (btopinv * aux + 0.5);
		/* dirty, but safe. On most machines (long)(btopinv*aux) */
		/* is correct, or one too big; on some however it becomes */
		/* too small. Could change zstart, but +0.5 and a while */
		/* instead of one single if is safer */
			if (qq > 0)
			{
				if (qq >= (1L<<26))
					qq = (1L<<26)-1;
				zsubmul1(qq, &c[i], &b[0]);
				while (*pc < 0)
				{
					qq--;
					{
						zaddmulone1(&c[i], &b[0]);
					}
				}
			}
			pc--;
			*p-- = qq;
		}
		sb--;
		while ((sb > 0) && (!(c[sb])))
			sb--;
		sb++;
		r[0] = sb;
		p = &r[1];
		pc = &c[0];
		for (i = sb; i > 0; i--)
			*p++ = *pc++;
		if (sq < 0)
		{
			q[0] = 1;
			q[1] = d[1];
		}
		else
		{
			sq++;
			while ((sq > 1) && (!(d[sq])))
				sq--;
			q[0] = sq;
			p = &q[1];
			pc = &d[1];
			for (i = sq; i > 0; i--)
				*p++ = *pc++;
		}
	}	/* non-trivial case */

	if (sign)
	{
		if (sign <= 2)
		{
			if (!(r[1]) && (r[0] == 1))
				q[0] = -(q[0]);
			else
			{
				zadd1(q, one, &q);
				q[0] = -q[0];
				if (sign == 1)
					zsub1(r, b, &r);
				else
					zsub1(b, r, &r);
			}
		}
		else
			znegate(&r);
	}
	*qqq = q;
	*rrr = r;

}


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

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

#define correct2( q, x1, x2, x3, y1, y2, btopinv) { \
        register long ix1=x1,ix2=x2,ix3=x3,iy1=y1,iy2=y2; \
        long lprodlow=0,lprodhigh=0; \
\
        zaddmulp2(&lprodlow, q, iy2, &lprodhigh); \
        if ((ix3 -= lprodlow) < 0) { \
                ix3+= RADIX; \
                ix2--; \
        } \
        lprodlow=0; \
        zaddmulp2(&lprodlow, q, iy1, &lprodhigh); \
        if ((ix2 -= lprodlow) < 0) { \
                ix2 += RADIX; \
                --ix1; \
        } \
        if (ix1 < lprodhigh) q--; \
        else if ((ix1 -= lprodhigh)) { \
                q += (long) ((fradix * (fradix * ix1 + ix2))*btopinv ); \
        } \
        else { \
                q += (long) ((fradix * ix2 + ix3)*btopinv ); \
        } \
}

#ifndef ALPHA50
void
zdiv2(
	verylong in_a,
	verylong in_b,
	verylong *qqq,
	verylong *rrr
	)

{
	register long i;
	register long qq;
	long sa;
	long sb;
	long sq;
	verylong p;
	verylong pc;
	long sign;
	static verylong a = 0;
	static verylong b = 0;
	static verylong c = 0;
	static verylong d = 0;
	double btopinv;
	double aux;
	verylong q = *qqq;
	verylong r = *rrr;
#ifdef ALPHA
        double btopinv2;
#endif

#ifndef START
	if (fudge < 0)
		zstart0();
#endif

	if (ALLOCATE && !in_a)
	{
		zzero(qqq);
		zzero(rrr);
		return;
	}
	if ((!in_b) || (((sb = in_b[0]) == 1) && (!in_b[1])))
	{
		zhalt("division by zero in zdiv");
		return;
	}
	zcopy(in_a,&a);

	zcopy(in_b,&b);

	sign = (*a < 0 ? 2 : 0) | (*b < 0 ? 1 : 0);
	if (*a < 0)
		(*a) = (-(*a));
	sa = (*a);
	if (*b < 0)
		(*b) = (-(*b));
	sb = (*b);
	zsetlength(&c, (i = (sa > sb ? sa : sb) + 1), "in zdiv, locals\n");
	zsetlength(&d, i, "");
	zsetlength(&q, i, "in zdiv, third argument");
	*qqq = q;
	zsetlength(&r, sb + 1, "in zdiv, fourth argument");
	*rrr = r;
	p = &b[sb];
	if ((sb == 1) && (*p < RADIX))
		zintoz(zsdiv(a, *p, &q), &r);
	else
	{
		sq = sa - sb;	/* size of quotient */
		btopinv = (double) (*p) * fradix;
		if (sb > 1)
			btopinv += (*(--p));
		btopinv *= fradix;
#ifdef ALPHA
                btopinv2=btopinv;
                btopinv2 = fudge2/btopinv2;
#else
                if (sb > 2)
                        btopinv += (*(p - 1));
#endif
		btopinv = fudge / btopinv;
		p = &a[1];
		pc = &c[0];
		for (i = sa; i > 0; i--)
			*pc++ = *p++;
		p = pc;
		sa = 1 - sb;
		for (i = (-sq); i > 0; i--)
			*p++ = 0;
		*pc = 0;
		d[1] = 0;
		p = &d[sq + 1];
		for (i = sq; i >= 0; i--)
		{

			aux = fradix * (fradix * (*pc) + (*(pc - 1))) + 1.0;
#ifndef ALPHA
                        if (i > sa)
                                aux += (*(pc - 2));
#endif
#ifdef ALPHA
                        qq = (long) (btopinv2 * aux + 0.5);
#else
                        qq = (long) (btopinv * aux + 0.5);
                /* dirty, but safe. On most machines (long)(btopinv*aux) */
                /* is correct, or one too big; on some however it becomes */
                /* too small. Could change zstart, but +0.5 and a while */
                /* instead of one single if is safer */
#endif
			if (qq > 0)
			{
#ifdef ALPHA
                                if (qq > (1L<<48)) {
                                        correct2(qq,*pc,*(pc-1),(i>sa) ? *(pc-2):0,b[sb],b[sb-1],btopinv);
                                        if (qq >= RADIX)
                                                qq = RADIX-1;
                                }
#else
                                if (qq >= RADIX)
                                        qq = RADIX-1;
#endif
				zsubmul2(qq, &c[i], &b[0]);
				while (*pc < 0)
				{
					qq--;
					{
						zaddmulone(&c[i], &b[0]);
					}
				}
			}
			pc--;
			*p-- = qq;
		}
		sb--;
		while ((sb > 0) && (!(c[sb])))
			sb--;
		sb++;
		r[0] = sb;
		p = &r[1];
		pc = &c[0];
		for (i = sb; i > 0; i--)
			*p++ = *pc++;
		if (sq < 0)
		{
			q[0] = 1;
			q[1] = d[1];
		}
		else
		{
			sq++;
			while ((sq > 1) && (!(d[sq])))
				sq--;
			q[0] = sq;
			p = &q[1];
			pc = &d[1];
			for (i = sq; i > 0; i--)
				*p++ = *pc++;
		}
	}	/* non-trivial case */

	if (sign)
	{
		if (sign