📄 lip.c
字号:
long d,
verylong *bb
)
{
register long sa;
register long i;
register verylong pb;
register long bnegative = 0;
verylong b = *bb;
if (( d >= RADIX) || (-d >= RADIX)) {
STATIC verylong x = 0;
zintoz(d,&x);
zmulin(a,&x);
zcopy(x,bb);
FREESPACE(x);
return;
}
if (ALLOCATE && !a)
{
zzero(bb);
return;
}
if (ALLOCATE && !d)
{
zzero(bb);
return;
}
if ((sa = a[0]) < 0)
{
sa = (-sa);
if (d < 0)
d = (-d);
else
bnegative = 1;
}
else if ((bnegative = (d < 0)))
d = (-d);
zsetlength(&b, sa + 1, "in zsmul, third argument");
if (a == *bb) a = b;
*bb = b;
pb = &b[0];
for (i = sa; i >= 0; i--)
*pb++ = *a++;
b[0] = sa;
sa++;
*pb = 0;
zaddmul(d - 1, &b[1], &b[0]);
while ((sa > 1) && (!(b[sa])))
sa--;
b[0] = sa;
if (bnegative && (b[1] || b[0] != 1))
b[0] = (-b[0]);
}
static void
kar_mul(
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 < KAR_MUL_CROV) || (i < KAR_MUL_CROV))
{
pc = &(*c)[1];
for (i = hal; i > 0; i--)
*pc++ = 0;
pc = &(*c)[1];
if (al <= *b)
for (i = al; i; i--)
{
zaddmul(*(++a), pc++, b);
}
else
for (i = *b; i; i--)
{
zaddmul(*(++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_mul(a, b, a4, shi + 5);
zadd(a, (*a1), a0);
a[0] = al;
if (bbig)
{
kar_mul((*a1), (*a3), c, shi + 5);
zadd(b, (*a3), a2);
b[0] = restoreb0;
kar_mul((*a0), (*a2), a3, shi + 5);
}
else
kar_mul((*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);
}
void
zmul(
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_mul(a, b, c, (long) 0);
else
kar_mul(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];
}
static void
kar_sq(
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 < KAR_SQU_CROV))
{
register unsigned long uncar;
long carry = 0;
pc = &(*c)[1];
for (; i > 0; i--)
*pc++ = 0;
for (hal = 1; hal <= al; hal++)
{
i += 2;
{
zaddmulsq(al - hal, &((*c)[i]), &(a[hal]));
}
uncar = ((*c)[i - 1] << 1) + carry;
(*c)[i - 1] = uncar & RADIXM;
uncar = ((*c)[i] << 1) + (uncar >> NBITS);
{
zaddmulpsq(&(*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_sq(a, a1, shi + 3);
zadd(a, (*a0), a2);
kar_sq((*a0), c, shi + 3);
a[0] = al;
kar_sq((*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);
}
void
zsq(
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_sq(a, c, (long) 0);
if (aneg)
a[0] = -a[0];
}
void
zmul_plain(
verylong a,
verylong b,
verylong *cc
)
{
register long i;
register verylong pc;
register long sc;
verylong c = *cc;
long anegative;
long bnegative;
verylong olda;
verylong oldb;
if (ALLOCATE && (!a || !b))
{
zzero(cc);
return;
}
olda = a;
oldb = b;
if ((anegative = (*a < 0)))
a[0] = -a[0];
if (olda == oldb)
bnegative = anegative;
else if ((bnegative = (*b < 0)))
b[0] = -b[0];
zsetlength(&c, (sc = *a + *b), "in zmul_plain, third argument");
*cc = c;
pc = &c[1];
for (i = sc; i > 0; i--)
*pc++ = 0;
pc = &c[1];
if (*a <= *b)
for (i = *a; i; i--)
{
zaddmul(*(++a), pc++, b);
}
else
for (i = *b; i; i--)
{
zaddmul(*(++b), pc++, a);
}
while ((sc > 1) && (!(c[sc])))
sc--;
c[0] = sc;
if (anegative != bnegative && (c[1] || c[0] != 1))
c[0] = -c[0];
if (anegative)
olda[0] = -olda[0];
if (bnegative && oldb != olda)
oldb[0] = -oldb[0];
}
void
zmulin(
verylong b,
verylong *a
)
{
STATIC verylong mem = 0;
if (ALLOCATE && (!*a || !b))
{
zzero(a);
return;
}
zcopy(*a, &mem);
zmul(mem, b, a);
FREESPACE(mem);
}
void
zsq_plain(
verylong a,
verylong *cc
)
{
register long i;
register long sc;
register verylong pc;
register unsigned long uncar;
long carry = 0;
verylong c = *cc;
long anegative;
if (ALLOCATE && !a)
{
zzero(cc);
return;
}
if ((anegative = (*a)) < 0)
anegative = -anegative;
zsetlength(&c, (sc = 2 * anegative), "in zsq_plain, second argument");
*cc = c;
pc = &c[1];
for (i = sc; i > 0; i--)
*pc++ = 0;
for (sc = 1; sc <= anegative; sc++)
{
i += 2;
{
zaddmulsq(anegative - sc, &(c[i]), &(a[sc]));
}
uncar = (c[i - 1] << 1) + carry;
c[i - 1] = uncar & RADIXM;
uncar = (c[i] << 1) + (uncar >> NBITS);
{
zaddmulpsq(&c[i - 1], a[sc], &carry);
}
uncar += carry;
carry = uncar >> NBITS;
c[i] = uncar & RADIXM;
}
while ((i > 1) && (!(c[i])))
i--;
c[0] = i;
}
void
zsqin(
verylong *a
)
{
STATIC verylong mem = 0;
if (ALLOCATE && !*a)
{
zzero(a);
return;
}
zcopy(*a, &mem);
zsq(mem, a);
FREESPACE(mem);
}
#ifndef ALPHA50
long
zsdiv(
verylong a,
long d,
verylong *bb
)
{
register long sa;
verylong b = *bb;
#ifndef START
if (fudge < 0)
zstart();
#endif
if (ALLOCATE && !a)
{
zzero(bb);
return (0);
}
if (!d)
{
zhalt("division by zero in zsdiv");
return (0);
}
if ((sa = a[0]) < 0)
sa = (-sa);
zsetlength(&b, sa, "in zsdiv, third argument");
if (a == *bb) a = b;
*bb = b;
if ((d >= RADIX) || (d <= -RADIX))
{
STATIC verylong zd = 0;
STATIC verylong zb = 0;
zintoz(d, &zb);
zdiv(a, zb, &b, &zd);
*bb = b;
sa = ztoint(zd);
FREE2SPACE(zd,zb);
return (sa);
}
else
{
register long den = d;
register double deninv;
register long carry = 0;
register long i;
long flag = (*a < 0 ? 2 : 0) | (den < 0 ? 1 : 0);
if (den < 0)
den = -den;
deninv = (double)1/den;
if (a[sa] < den && sa > 1)
carry = a[sa--];
for (i = sa; i; i--)
zdiv21(carry, a[i], den, deninv, b[i]);
while ((sa > 1) && (!(b[sa])))
sa--;
b[0] = sa;
if (flag)
{
if (flag <= 2)
{
if (!carry)
znegate(&b);
else
{
zadd(b, one, &b);
b[0] = -b[0];
if (flag == 1)
carry = carry - den;
else
carry = den - carry;
*bb = b;
}
}
else
carry = -carry;
}
return (carry);
}
}
#endif
long
zsmod(
verylong a,
long b
)
{
STATIC verylong q = 0;
long y;
y= zsdiv(a, b, &q);
FREESPACE(q);
return y;
}
#define correct( q, x1, x2, x3, y1, y2, btopinv) { \
register long ix1=x1,ix2=x2,ix3=x3,iy1=y1,iy2=y2; \
long lprodlow=0,lprodhigh=0; \
\
zaddmulp(&lprodlow, q, iy2, &lprodhigh); \
if ((ix3 -= lprodlow) < 0) { \
ix3+= RADIX; \
ix2--; \
} \
lprodlow=0; \
zaddmulp(&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
zdiv(
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;
#ifdef ALPHA
double btopinv2;
#endif
double aux;
verylong q = *qqq;
verylong r = *rrr;
#ifndef START
if (fudge < 0)
zstart();
#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 + 1.0);
#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)) {
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)
{
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);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -