📄 lip.c
字号:
q[0] = -q[0];
if (sign == 1)
zsub(r, b, &r);
else
zsub(b, r, &r);
}
}
else
znegate(&r);
}
*qqq = q;
*rrr = r;
FREE2SPACE(a,b); FREE2SPACE(c,d);
}
void
zmod(
verylong in_a,
verylong in_b,
verylong *rr
)
{
register long i;
register long qq;
verylong r = *rr;
STATIC verylong a = 0;
STATIC verylong b = 0;
STATIC verylong c = 0;
long sa;
long sb;
long sq;
verylong p;
verylong pc;
long sign;
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);
}
long
zmontred(
verylong a,
verylong *cc
)
{
register long i;
register long j;
register long a0;
long neg = 0;
long subcnt = 0;
verylong c = *cc;
STATIC verylong x = 0;
verylong px;
if (!zn)
{
zhalt("undefined Montgomery modulus in zmontred");
return(-1);
}
if (ALLOCATE && (!a))
{
zzero(cc);
return(-1);
}
i = (zntop << 1) + 1;
a0 = a[0];
if (a0 < 0) {
a0 = -a0;
neg = 1;
}
zsetlength(&x, i, "in zmontred, local");
zsetlength(&c, zntop+1, "in zmontred, second argument");
*cc = c;
for (j=1; j <= a0; j++) x[j] = a[j];
for (; j <= i; j++) x[j] = (long)0;
x[0] = a0;
/*
printf("copy a: ");
for (j=0;j<=i;j++) printf("%ld ",x[j]);
*/
px = &x[1];
for (j = zntop; j > 0; j--) zmmulp(px++);
/*
printf("redu a: ");
for (j=0;j<=i;j++) printf("%ld ",x[j]);
*/
for (j=1; j <= zntop+1; j++) c[j] = *(px++);
/*
printf(" : ");
for (j=0;j<=zntop+1;j++) printf("%ld ",c[j]);
*/
i = zntop+1;
while ((i > 1) && (!(c[i]))) i--;
c[0] = i;
while (zcompare(c, zn) >= 0) {
zsubpos(c, zn, &c);
subcnt ++;
}
if (neg && (!ziszero(c))) zsub(zn,c,&c);
*cc = c;
FREESPACE(x);
return(subcnt);
}
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--)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -