lip.c

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

	double btopinv;
#ifdef ALPHA
        double btopinv2;
#endif
	double aux;

/*printf("in zmod: "); zwrite(in_a); printf(" mod "); zwriteln(in_b); fflush(stdout);
*/
#ifndef START
	if (fudge < 0)
		zstart();
#endif
	if (ALLOCATE && !in_a)
	{
		zzero(rr);
		return;
	}
	if ((!in_b) || (((sb = in_b[0]) == 1) && (!in_b[1])))
	{
		zhalt("division by zero in zmod");
		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, (sa > sb ? sa : sb) + 1, "in zmod, local");
	zsetlength(&r, sb + 1, "in zmod, third argument");
	*rr = r;
	p = &b[sb];
	if ((sb == 1) && (*p < RADIX))
		zintoz(zsdiv(a, *p, &c), &r);
	else
	{
		sq = sa - sb;
		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;
		for (i = sq; i >= 0; i--)
		{
			aux = fradix * (fradix * (*pc) + (*(pc - 1))) + 1.0;
#ifdef ALPHA
			qq = (long) (btopinv2 * aux + 1.0);
#else
			if (i > sa)
				aux += (*(pc - 2));
			qq = (long) (btopinv * aux + 0.5);

		/* see comment in zdiv */
#endif
			if (qq > 0)
			{
#ifdef ALPHA
				if (qq > (1L<<48)) {
					correct(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
				zsubmul(qq, &c[i], &b[0]);
				while (*pc < 0)
				{
					{
						zaddmulone(&c[i], &b[0]);
					}
				}
			}
			pc--;
		}	/* loop on i */
		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++;
	}	/* non-trivial case */

	if (sign)
	{
		if ((sign <= 2) && (!((r[0] == 1) && !(r[1]))))
		{
			if (sign == 1)
				zsub(r, b, &r);
			else
				zsub(b, r, &r);
		}
		else
			znegate(&r);	/* negating r==0 doesn't hurt */
	}
	*rr = r;
	FREE2SPACE(a,b); FREESPACE(c);

}
#endif


void
zaddmod(
	verylong a,
	verylong b,
	verylong n,
	verylong *c
	)
{
	if (ALLOCATE && !n)
	{
		zhalt("modulus zero in zaddmod");
		return;
	}
	if (ALLOCATE && !a)
	{
		if (b)
			zcopy(b, c);
		else
			zzero(c);
		return;
	}
	if (ALLOCATE && !b)
	{
		zcopy(a, c);
		return;
	}
	zadd(a, b, c);
	if (zcompare(*c, n) >= 0)
		zsubpos(*c, n, c);

}

void
zsubmod(
	verylong a,
	verylong b,
	verylong n,
	verylong *c
	)
{
	if (ALLOCATE && !n)
	{
		zhalt("modulus zero in zsubmod");
		return;
	}
	if (ALLOCATE && !b)
	{
		if (a)
			zcopy(a, c);
		else
			zzero(c);
		return;
	}
	if (ALLOCATE && !a)
	{
		if (b[1] || (b[0] != 1))
			zsubpos(n, b, c);
		else
			zzero(c);
		return;
	}
	zsub(a, b, c);
	if ((*c)[0] < 0)
		zadd(*c, n, c);
}

void
zsmulmod(
	verylong a,
	long d,
	verylong n,
	verylong *c
	)
{
	STATIC verylong mem = 0;

	if (ALLOCATE && !n)
	{
		zhalt("modulus zero in zsmulmod");
		return;
	}
	if (ALLOCATE && (!a || !d))
	{
		zzero(c);
		return;
	}
	zsmul(a, d, &mem);
	zmod(mem, n, c);
	FREESPACE(mem);
}

void
zmulmod(
	verylong a,
	verylong b,
	verylong n,
	verylong *c
	)
{
	STATIC verylong mem = 0;

	if (ALLOCATE && !n)
	{
		zhalt("modulus zero in zmulmod");
		return;
	}
	if (ALLOCATE && (!a || !b))
	{
		zzero(c);
		return;
	}
	zsetlength(&mem, a[0] + b[0] + 2, "in zmulmod, local");
	zmul(a, b, &mem);
	zmod(mem, n, c);
	FREESPACE(mem);
}

void
zmulinmod(
	verylong a,
	verylong *b,
	verylong n
	)
{
	if (ALLOCATE && !n)
	{
		zhalt("modulus zero in zmulinmod");
		return;
	}
	if (ALLOCATE && (!a || !(*b)))
	{
		zzero(b);
		return;
	}
	zmulin(a, b);
	zmod(*b, n, b);
}

void
zsqmod(
	verylong a,
	verylong n,
	verylong *c
	)
{
	STATIC verylong mem = 0;

	if (ALLOCATE && !n)
	{
		zhalt("modulus zero in zsqmod");
		return;
	}
	if (ALLOCATE && !a)
	{
		zzero(c);
		return;
	}
	zsetlength(&mem, 2 * a[0] + 2, "in zsqmod, local");
	zsq(a, &mem);
	zmod(mem, n, c);
	FREESPACE(mem);
}

void
zsqinmod(
	verylong *a,
	verylong n
	)
{
	if (ALLOCATE && !n)
	{
		zhalt("modulus zero in zsqinmod");
		return;
	}
	if (ALLOCATE && !(*a))
	{
		zzero(a);
		return;
	}
	zsqin(a);
	zmod(*a, n, a);
}

void
zdivmod(
	verylong a,
	verylong b,
	verylong n,
	verylong *c
	)
{
	STATIC verylong g = 0;
	STATIC verylong na = 0;
	STATIC verylong nb = 0;
	STATIC verylong r = 0;

	if (ALLOCATE && !n)
	{
		zhalt("modulus zero in zdivmod");
		return;
	}
	if (ALLOCATE && !a)
	{
		zzero(c);
		return;
	}
	if (ALLOCATE && !b)
	{
		zhalt("division by zero in zdivmod");
		return;
	}
	zgcd(a, b, &g);
	if ((g[1] == 1) && (g[0] == 1))
	{
		if (zinv(b, n, &g))
		{
			zhalt("undefined quotient in zdivmod");
			zcopy(g,c);
		}
		else
			zmulmod(a, g, n, c);
	}
	else
	{
		zdiv(a, g, &na, &r);
		zdiv(b, g, &nb, &r);
		if (zinv(nb, n, &g))
		{
			zhalt("undefined quotient in zdivmod");
                        zcopy(g,c);
                }
		else
			zmulmod(na, g, n, c);
	}
	FREE2SPACE(g,na); FREE2SPACE(nb,r);
}

void
zinvmod(
	verylong a,
	verylong n,
	verylong *c
	)
{
	if (ALLOCATE && !n)
	{
		zhalt("modulus zero in zinvmod");
		return;
	}
	if (ALLOCATE && !a)
	{
		zhalt("division by zero in zinvmod");
		return;
	}
	if (zinv(a, n, c))
		zhalt("undefined inverse in zinvmod");
}

static void
zmstartint(
	verylong n
	)
{
	STATIC verylong x = 0;
	long i;

	if (!n || !(n[1] & 1) || ((zntop = n[0]) < 1)
			|| ((zntop == 1) && (n[1] <= 1)))
	{
		zhalt("even or negative modulus in zmstart");
		return;
	}
	if (n[zntop] >= (RADIX >> 1))
		zntop++;
	zsetlength(&x, (long) (zntop + 1), "in zmstart, local");
	zsetlength(&zn, (long) (zntop + 1), "in zmstart, globals\n");
	zsetlength(&zr, (long) (zntop + 1), "");
	zsetlength(&zrr, (long) (zntop + 1), "");
	zsetlength(&zrrr, (long) (zntop + 1), "");
	zsetlength(&znm, (long) (zntop + 1), "");
	if (zcompare(n, zn))
	{
		zcopy(n, &zn);
#ifdef PLAIN_OR_KARAT
		zninv1 = zinvs(RADIX-zn[1],RADIX);
		zninv2 = zninv1 >> NBITSH;
		zninv1 &= RADIXROOTM;
#else
		zninv = zinvs(RADIX-zn[1],RADIX);
#endif
		for (i = 1; i <= zntop; i++)
			x[i] = 0;
		x[zntop + 1] = 1;
		x[0] = zntop + 1;
		zdiv(x, zn, &zrr, &zr);
		zsqmod(zr, zn, &zrr);
		zmulmod(zr, zrr, zn, &zrrr);
		zsubmod(zn, zr, zn, &znm);
	}
	FREESPACE(x);
}

static void
zmkeep(
	verylong n
	)
{
	if (znotinternal)
		zcopy(zn, &zoldzn);
	zmstartint(n);
}

static void
zmback(
	)
{
	if (znotinternal)
		zmstartint(zoldzn);
}

void
zmstart(
	verylong n
	)
{
	zmstartint(n);
	znotinternal = 1;
}
void
zmfree(
	)
{
	znotinternal = 0;
}

void
ztom(
	verylong a,
	verylong *bb
	)
{
	if (!zn)
	{
		zhalt("undefined Montgomery modulus in ztom");
		return;
	}
	if (ALLOCATE && !a)
	{
		zzero(bb);
		return;
	}
	if ((zscompare(a,0)<0) || (zcompare(zn, a) < 0))
	{
		zmod(a, zn, bb);
		zmontmul(*bb, zrr, bb);
	}
	else
		zmontmul(a, zrr, bb);
}

void
zmtoz(
	verylong a,
	verylong *bb
	)
{
	if (!zn)
	{
		zhalt("undefined Montgomery modulus in zmtoz");
		return;
	}
	if (ALLOCATE && !a)
	{
		zzero(bb);
		return;
	}
	zmontmul(a, one, bb);
}

void
zmontadd(
	verylong a,
	verylong b,
	verylong *c
	)
{
	if (!zn)
	{
		zhalt("undefined Montgomery modulus in zmontadd");
		return;
	}
	zaddmod(a, b, zn, c);
}

void
zmontsub(
	verylong a,
	verylong b,
	verylong *c
	)
{
	if (!zn)
	{
		zhalt("undefined Montgomery modulus in zmontsub");
		return;
	}
	zsubmod(a, b, zn, c);
}

void
zsmontmul(
	verylong a,
	long d,
	verylong *c
	)
{
	if (!zn)
	{
		zhalt("undefined Montgomery modulus in zsmontmul");
		return;
	}
	zsmulmod(a, d, zn, c);
}

void
zmontmul(
	verylong a,
	verylong b,
	verylong *cc
	)
{
	register long i;
	register long j;
	verylong c = *cc;
	STATIC verylong x = 0;
	verylong px;
	verylong pc;

	if (!zn)
	{
		zhalt("undefined Montgomery modulus in zmontmul");
		return;
	}
	if (ALLOCATE && (!a || !b))
	{
		zzero(cc);
		return;
	}
	zsetlength(&x, (i = (zntop << 1) + 1), "in zmontmul, local");
	zsetlength(&c, zntop, "in zmontmul, third argument");
	if (a == *cc) a = c;
	if (b == *cc) b = c;
	*cc = c;
        if (!kar_mem_initialized)
        {
                kar_mem_initialized = 1;
                for (j = (5 * KAR_DEPTH) - 1; j >= 0; j--)
                        kar_mem[j] = (verylong) 0;
        }
	if (*a <= *b)
		kar_mul(b, a, &x, (long) 0);
	else
		kar_mul(a, b, &x, (long) 0);
	for (; i > x[0]; i--)
		x[i] = (long) 0;
	px = &x[1];
	for (i = zntop; i > 0; i--)
		zmmulp(px++);
	pc = &c[1];
	for (i = zntop; i > 0; i--)
		*pc++ = *px++;
	i = zntop;
	while ((i > 1) && (!(*--pc)))
		i--;
	c[0] = i;
	if (zcompare(c, zn) >= 0)
		zsubpos(c, zn, &c);
	*cc = c;
	FREESPACE(x);
}

void
zmontsq(
	verylong a,
	verylong *cc
	)
{
	register long i;
	register long j;
	verylong c = *cc;
	STATIC verylong x = 0;
	verylong px;
	verylong pc;

	if (!zn)
	{
		zhalt("undefined Montgomery modulus in zmontsq");
		return;
	}
	if (ALLOCATE && !a)
	{
		zzero(cc);
		return;
	}
	zsetlength(&x, (i = (zntop << 1) + 1), "in zmontsq, local");
	zsetlength(&c, zntop, "in zmontsq, third argument");
	if (a == *cc) a = c;
	*cc = c;
        if (!kar_mem_initialized)
        {
                kar_mem_initialized = 1;
                for (j = (5 * KAR_DEPTH) - 1; j >= 0; j--)
                        kar_mem[j] = (verylong) 0;
        }
	kar_sq(a, &x, (long) 0);
	for (; i > x[0]; i--)
		x[i] = (long) 0;
	px = &x[1];
	for (i = zntop; i > 0; i--)
		zmmulp(px++);
	pc = &c[1];
	for (i = zntop; i > 0; i--)
		*pc++ = *px++;
	i = zntop;
	while ((i > 1) && (!(*--pc)))
		i--;
	c[0] = i;
	if (zcompare(c, zn) >= 0)
		zsubpos(c, zn, &c);
	*cc = c;
	FREESPACE(x);
}

void
zmontdiv(
	verylong a,
	verylong b,
	verylong *c
	)
{
	STATIC verylong g = 0;
	STATIC verylong na = 0;
	STATIC verylong nb = 0;
	STATIC verylong r = 0;

	if (!zn)
	{
		zhalt("undefined Montgomery modulus in zmontdiv");
		return;
	}
	if (ALLOCATE && !a)
	{
		zzero(c);
		return;
	}
	if (ALLOCATE && !b)
	{
		zhalt("division by zero in zmontdiv");
		return;
	}
	zgcd(a, b, &g);
	if ((g[1] == 1) && (g[0] == 1))
	{
		if (zinv(b, zn, &g))
		{
			zhalt("undefined quotient in zmontdiv");
			zcopy(g,c);
			FREESPACE(g);
			return;
		}
		else
			zmontmul(g, a, &g);
	}
	else
	{
		zdiv(a, g, &na, &r);
		zdiv(b, g, &nb, &r);
		if (zinv(nb, zn, &g))
		{
			zhalt("undefined quotient in zmontdiv");
			zcopy(g,c);
                        return;
		}
		else
			zmontmul(g, na, &g);
	}
	zmontmul(g, zrrr, c);
	FREE2SPACE(g,na); FREE2SPACE(nb,r);
}

void
zmontinv(
	verylong a,
	verylong *c
	)
{
	if (!zn)
	{
		zhalt("undefined Montgomery modulus in zmontinv");
		return;
	}
        if (ALLOCATE && !a)
	{
                zhalt("division by zero in zmontinv");
		return;
	}
	if (zinv(a, zn, c))
	{
		zhalt("undefined inverse in zmontinv");
		return;
	}
	zmontmul(*c, zrrr, c);
}


void
zrstarts(
	long s
	)
{
	register long i;
	register long ones;
	register long res;
	STATIC verylong three = 0;

	zintoz(s, &zseed);
	if (!three)
	{
		zintoz((long) 3, &three);
		ones = 107 / NBITS;
		res = 107 % NBITS;
		zsetlength(&zranp, ones + 1, "in zrstarts, zranp");
		for (i = ones; i; i--)
			zranp[i] = RADIXM;
		if (res)
			zranp[++ones] = (RADIXM >> (NBITS - res));
		zra