version 1.1, 2001/10/02 11:17:09 |
version 1.2, 2002/09/11 07:27:00 |
Line 30 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
Line 30 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
*/ |
*/ |
|
|
/* z2 := z1[imin..imax].f shifted left sh bits (feeding f from the right) */ |
/* z2 := z1[imin..imax].f shifted left sh bits (feeding f from the right) */ |
|
/* These macros work only for sh != 0 !!! */ |
#define shift_left2(z2,z1,imin,imax,f, sh,m) {\ |
#define shift_left2(z2,z1,imin,imax,f, sh,m) {\ |
register ulong _l,_k = ((ulong)f)>>m;\ |
register ulong _l,_k = ((ulong)f)>>m;\ |
GEN t1 = z1 + imax, t2 = z2 + imax, T = z1 + imin;\ |
GEN t1 = z1 + imax, t2 = z2 + imax, T = z1 + imin;\ |
Line 45 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
Line 46 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
shift_left2((z2),(z1),(imin),(imax),(f),(sh),(_m));\ |
shift_left2((z2),(z1),(imin),(imax),(f),(sh),(_m));\ |
} |
} |
|
|
/* z2 := f.z1[imin..imax-1] shifted right sh bits (feeding f from the left) */ |
#define shift_words_r(target,source,source_end,prepend, sh, sh_complement) {\ |
#define shift_right2(z2,z1,imin,imax,f, sh,m) {\ |
register ulong _k,_l = *source++;\ |
register GEN t1 = z1 + imin, t2 = z2 + imin, T = z1 + imax;\ |
*target++ = (_l>>(ulong)sh) | ((ulong)prepend<<(ulong)sh_complement);\ |
register ulong _k,_l = *t1++;\ |
while (source < source_end) {\ |
*t2++ = (_l>>(ulong)sh) | ((ulong)f<<(ulong)m);\ |
_k = _l<<(ulong)sh_complement; _l = *source++;\ |
while (t1 < T) {\ |
*target++ = (_l>>(ulong)sh) | _k;\ |
_k = _l<<(ulong)m; _l = *t1++;\ |
|
*t2++ = (_l>>(ulong)sh) | _k;\ |
|
}\ |
}\ |
} |
} |
|
#define shift_right2(z2,z1,imin,imax,f, sh,m) {\ |
|
register GEN s = (z1) + (imin), ta = (z2) + (imin), se = (z1) + (imax);\ |
|
shift_words_r(ta,s,se,(f),(sh),(m)); \ |
|
} |
|
/* z2 := f.z1[imin..imax-1] shifted right sh bits (feeding f from the left) */ |
#define shift_right(z2,z1,imin,imax,f, sh) {\ |
#define shift_right(z2,z1,imin,imax,f, sh) {\ |
register const ulong _m = BITS_IN_LONG - sh;\ |
register const ulong _m = BITS_IN_LONG - (sh);\ |
shift_right2((z2),(z1),(imin),(imax),(f),(sh),(_m));\ |
shift_right2((z2),(z1),(imin),(imax),(f),(sh),(_m));\ |
} |
} |
|
|
Line 95 addsispec(long s, GEN x, long nx) |
|
Line 99 addsispec(long s, GEN x, long nx) |
|
while (xd > x) *--zd = *--xd; |
while (xd > x) *--zd = *--xd; |
*--zd = evalsigne(1) | evallgefint(lz); |
*--zd = evalsigne(1) | evallgefint(lz); |
*--zd = evaltyp(t_INT) | evallg(lz); |
*--zd = evaltyp(t_INT) | evallg(lz); |
avma=(long)zd; return zd; |
avma=(gpmem_t)zd; return zd; |
} |
} |
|
|
#define swapspec(x,y, nx,ny) {long _a=nx;GEN _z=x; nx=ny; ny=_a; x=y; y=_z;} |
#define swapspec(x,y, nx,ny) {long _a=nx;GEN _z=x; nx=ny; ny=_a; x=y; y=_z;} |
Line 129 addiispec(GEN x, GEN y, long nx, long ny) |
|
Line 133 addiispec(GEN x, GEN y, long nx, long ny) |
|
while (xd > x) *--zd = *--xd; |
while (xd > x) *--zd = *--xd; |
*--zd = evalsigne(1) | evallgefint(lz); |
*--zd = evalsigne(1) | evallgefint(lz); |
*--zd = evaltyp(t_INT) | evallg(lz); |
*--zd = evaltyp(t_INT) | evallg(lz); |
avma=(long)zd; return zd; |
avma=(gpmem_t)zd; return zd; |
} |
} |
|
|
/* assume x >= y */ |
/* assume x >= y */ |
Line 158 subisspec(GEN x, long s, long nx) |
|
Line 162 subisspec(GEN x, long s, long nx) |
|
do *--zd = *--xd; while (xd > x); |
do *--zd = *--xd; while (xd > x); |
*--zd = evalsigne(1) | evallgefint(lz); |
*--zd = evalsigne(1) | evallgefint(lz); |
*--zd = evaltyp(t_INT) | evallg(lz); |
*--zd = evaltyp(t_INT) | evallg(lz); |
avma=(long)zd; return zd; |
avma=(gpmem_t)zd; return zd; |
} |
} |
|
|
/* assume x > y */ |
/* assume x > y */ |
Line 191 subiispec(GEN x, GEN y, long nx, long ny) |
|
Line 195 subiispec(GEN x, GEN y, long nx, long ny) |
|
do *--zd = *--xd; while (xd > x); |
do *--zd = *--xd; while (xd > x); |
*--zd = evalsigne(1) | evallgefint(lz); |
*--zd = evalsigne(1) | evallgefint(lz); |
*--zd = evaltyp(t_INT) | evallg(lz); |
*--zd = evaltyp(t_INT) | evallg(lz); |
avma=(long)zd; return zd; |
avma=(gpmem_t)zd; return zd; |
} |
} |
|
|
#ifndef __M68K__ |
#ifndef __M68K__ |
|
|
/* prototype of positive small ints */ |
/* prototype of positive small ints */ |
static long pos_s[] = { |
static long pos_s[] = { |
evaltyp(t_INT) | m_evallg(3), evalsigne(1) | evallgefint(3), 0 }; |
evaltyp(t_INT) | _evallg(3), evalsigne(1) | evallgefint(3), 0 }; |
|
|
/* prototype of negative small ints */ |
/* prototype of negative small ints */ |
static long neg_s[] = { |
static long neg_s[] = { |
evaltyp(t_INT) | m_evallg(3), evalsigne(-1) | evallgefint(3), 0 }; |
evaltyp(t_INT) | _evallg(3), evalsigne(-1) | evallgefint(3), 0 }; |
|
|
void |
void |
affir(GEN x, GEN y) |
affir(GEN x, GEN y) |
Line 213 affir(GEN x, GEN y) |
|
Line 217 affir(GEN x, GEN y) |
|
if (!s) |
if (!s) |
{ |
{ |
y[1] = evalexpo(-bit_accuracy(ly)); |
y[1] = evalexpo(-bit_accuracy(ly)); |
y[2] = 0; return; |
return; |
} |
} |
|
|
lx=lgefint(x); sh=bfffo(x[2]); |
lx = lgefint(x); sh = bfffo(x[2]); |
y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1); |
y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1); |
if (sh) |
if (sh) { |
{ |
|
if (lx>ly) |
if (lx>ly) |
{ shift_left(y,x,2,ly-1, x[ly],sh); } |
{ shift_left(y,x,2,ly-1, x[ly],sh); } |
else |
else |
Line 227 affir(GEN x, GEN y) |
|
Line 230 affir(GEN x, GEN y) |
|
for (i=lx; i<ly; i++) y[i]=0; |
for (i=lx; i<ly; i++) y[i]=0; |
shift_left(y,x,2,lx-1, 0,sh); |
shift_left(y,x,2,lx-1, 0,sh); |
} |
} |
return; |
|
} |
} |
|
else { |
if (lx>=ly) |
if (lx>=ly) |
for (i=2; i<ly; i++) y[i]=x[i]; |
for (i=2; i<ly; i++) y[i]=x[i]; |
else |
else |
{ |
{ |
for (i=2; i<lx; i++) y[i]=x[i]; |
for (i=2; i<lx; i++) y[i]=x[i]; |
for ( ; i<ly; i++) y[i]=0; |
for ( ; i<ly; i++) y[i]=0; |
|
} |
} |
} |
} |
} |
|
|
Line 244 affrr(GEN x, GEN y) |
|
Line 247 affrr(GEN x, GEN y) |
|
{ |
{ |
long lx,ly,i; |
long lx,ly,i; |
|
|
y[1] = x[1]; if (!signe(x)) { y[2]=0; return; } |
y[1] = x[1]; if (!signe(x)) return; |
|
|
lx=lg(x); ly=lg(y); |
lx=lg(x); ly=lg(y); |
if (lx>=ly) |
if (lx>=ly) |
Line 259 affrr(GEN x, GEN y) |
|
Line 262 affrr(GEN x, GEN y) |
|
GEN |
GEN |
shifti(GEN x, long n) |
shifti(GEN x, long n) |
{ |
{ |
|
return shifti3(x, n, 0); |
|
} |
|
|
|
GEN |
|
shifti3(GEN x, long n, long flag) |
|
{ |
const long s=signe(x); |
const long s=signe(x); |
long lx,ly,i,m; |
long lx,ly,i,m; |
GEN y; |
GEN y; |
Line 285 shifti(GEN x, long n) |
|
Line 294 shifti(GEN x, long n) |
|
} |
} |
else |
else |
{ |
{ |
|
long lyorig; |
|
|
|
if (s > 0) |
|
flag = 0; |
n = -n; |
n = -n; |
ly = lx - (n>>TWOPOTBITS_IN_LONG); |
ly = lyorig = lx - (n>>TWOPOTBITS_IN_LONG); |
if (ly<3) return gzero; |
if (ly<3) |
|
return flag ? stoi(-1) : gzero; |
y = new_chunk(ly); |
y = new_chunk(ly); |
m = n & (BITS_IN_LONG-1); |
m = n & (BITS_IN_LONG-1); |
if (!m) for (i=2; i<ly; i++) y[i]=x[i]; |
if (m) { |
else |
|
{ |
|
shift_right(y,x, 2,ly, 0,m); |
shift_right(y,x, 2,ly, 0,m); |
if (y[2] == 0) |
if (y[2] == 0) |
{ |
{ |
if (ly==3) { avma = (long)(y+3); return gzero; } |
if (ly==3) { avma = (gpmem_t)(y+3); return flag ? stoi(-1) : gzero; } |
ly--; avma = (long)(++y); |
ly--; avma = (gpmem_t)(++y); |
} |
} |
|
} else { |
|
for (i=2; i<ly; i++) y[i]=x[i]; |
} |
} |
|
/* With FLAG: round up instead of rounding down */ |
|
if (flag) { /* Check low bits of x */ |
|
i = lx; flag = 0; |
|
while (--i >= lyorig) |
|
if (x[i]) { flag = 1; break; } /* Need to increment y by 1 */ |
|
if (!flag && m) |
|
flag = x[lyorig - 1] & ((1<<m) - 1); |
|
} |
|
if (flag) { /* Increment y */ |
|
for (i = ly;;) |
|
{ |
|
if (--i < 2) |
|
{ /* Extend y on the left */ |
|
if (avma <= bot) err(errpile); |
|
avma = (gpmem_t)(--y); ly++; |
|
y[2] = 1; break; |
|
} |
|
if (++y[i]) break; |
|
/* Now propagate the bit into the next longword */ |
|
} |
|
} |
} |
} |
y[1] = evalsigne(s)|evallgefint(ly); |
y[1] = evalsigne(s)|evallgefint(ly); |
y[0] = evaltyp(t_INT)|evallg(ly); return y; |
y[0] = evaltyp(t_INT)|evallg(ly); return y; |
Line 524 addir(GEN x, GEN y) |
|
Line 559 addir(GEN x, GEN y) |
|
#else /* design issue: make 0.0 "absorbing" */ |
#else /* design issue: make 0.0 "absorbing" */ |
if (e>0) return rcopy(y); |
if (e>0) return rcopy(y); |
#endif |
#endif |
z = cgetr(3 + ((-e)>>TWOPOTBITS_IN_LONG)); |
return itor(x, 3 + ((-e)>>TWOPOTBITS_IN_LONG)); |
affir(x,z); return z; |
|
} |
} |
|
|
ly=lg(y); |
ly=lg(y); |
if (e>0) |
if (e > 0) |
{ |
{ |
l = ly - (e>>TWOPOTBITS_IN_LONG); |
l = ly - (e>>TWOPOTBITS_IN_LONG); |
if (l<3) return rcopy(y); |
if (l<3) return rcopy(y); |
} |
} |
else l = ly + ((-e)>>TWOPOTBITS_IN_LONG)+1; |
else l = ly + ((-e)>>TWOPOTBITS_IN_LONG)+1; |
z=cgetr(l); affir(x,z); y=addrr(z,y); |
y = addrr(itor(x,l), y); |
z = y+l; ly = lg(y); while (ly--) z[ly] = y[ly]; |
z = y+l; ly = lg(y); while (ly--) z[ly] = y[ly]; |
avma=(long)z; return z; |
avma=(gpmem_t)z; return z; |
} |
} |
|
|
GEN |
GEN |
Line 554 addrr(GEN x, GEN y) |
|
Line 588 addrr(GEN x, GEN y) |
|
if (!sx) |
if (!sx) |
{ |
{ |
if (e > 0) ex=ey; |
if (e > 0) ex=ey; |
z=cgetr(3); z[1]=evalexpo(ex); z[2]=0; return z; |
return realzero_bit(ex); |
} |
} |
if (e > 0) { z=cgetr(3); z[1]=evalexpo(ey); z[2]=0; return z; } |
if (e > 0) return realzero_bit(ey); |
lz = 3 + ((-e)>>TWOPOTBITS_IN_LONG); |
lz = 3 + ((-e)>>TWOPOTBITS_IN_LONG); |
lx = lg(x); if (lz>lx) lz=lx; |
lx = lg(x); if (lz>lx) lz=lx; |
z = cgetr(lz); while(--lz) z[lz]=x[lz]; |
z = cgetr(lz); while(--lz) z[lz]=x[lz]; |
Line 564 addrr(GEN x, GEN y) |
|
Line 598 addrr(GEN x, GEN y) |
|
} |
} |
if (!sx) |
if (!sx) |
{ |
{ |
if (e < 0) { z=cgetr(3); z[1]=evalexpo(ex); z[2]=0; return z; } |
if (e < 0) return realzero_bit(ex); |
lz = 3 + (e>>TWOPOTBITS_IN_LONG); |
lz = 3 + (e>>TWOPOTBITS_IN_LONG); |
ly = lg(y); if (lz>ly) lz=ly; |
ly = lg(y); if (lz>ly) lz=ly; |
z = cgetr(lz); while (--lz) z[lz]=y[lz]; |
z = cgetr(lz); while (--lz) z[lz]=y[lz]; |
Line 602 addrr(GEN x, GEN y) |
|
Line 636 addrr(GEN x, GEN y) |
|
|
|
if (sx==sy) |
if (sx==sy) |
{ /* addition */ |
{ /* addition */ |
i=lz-1; avma = (long)z; |
i=lz-1; avma = (gpmem_t)z; |
if (flag==0) { z[i] = y[i]; i--; } |
if (flag==0) { z[i] = y[i]; i--; } |
overflow=0; |
overflow=0; |
for (j=lx-1; j>=2; i--,j--) z[i] = addllx(x[j],y[i]); |
for (j=lx-1; j>=2; i--,j--) z[i] = addllx(x[j],y[i]); |
Line 613 addrr(GEN x, GEN y) |
|
Line 647 addrr(GEN x, GEN y) |
|
if (i==1) |
if (i==1) |
{ |
{ |
shift_right(z,z, 2,lz, 1,1); |
shift_right(z,z, 2,lz, 1,1); |
ey=evalexpo(ey+1); |
z[1] = evalsigne(sx) | evalexpo(ey+1); return z; |
z[1] = evalsigne(sx) | ey; return z; |
|
} |
} |
z[i] = y[i]+1; if (z[i--]) break; |
z[i] = y[i]+1; if (z[i--]) break; |
} |
} |
Line 629 addrr(GEN x, GEN y) |
|
Line 662 addrr(GEN x, GEN y) |
|
i=2; while (i<lx && x[i]==y[i]) i++; |
i=2; while (i<lx && x[i]==y[i]) i++; |
if (i==lx) |
if (i==lx) |
{ |
{ |
avma = (long)(z+lz); |
avma = (gpmem_t)(z+lz); |
e = evalexpo(ey - bit_accuracy(lx)); |
return realzero_bit(ey - bit_accuracy(lx)); |
z=cgetr(3); z[1]=e; z[2]=0; return z; |
|
} |
} |
f2 = ((ulong)y[i] > (ulong)x[i]); |
f2 = ((ulong)y[i] > (ulong)x[i]); |
} |
} |
Line 656 addrr(GEN x, GEN y) |
|
Line 688 addrr(GEN x, GEN y) |
|
|
|
x = z+2; i=0; while (!x[i]) i++; |
x = z+2; i=0; while (!x[i]) i++; |
lz -= i; z += i; m = bfffo(z[2]); |
lz -= i; z += i; m = bfffo(z[2]); |
e = evalexpo(ey - (m | (i<<TWOPOTBITS_IN_LONG))); |
|
if (m) shift_left(z,z,2,lz-1, 0,m); |
if (m) shift_left(z,z,2,lz-1, 0,m); |
z[1] = evalsigne(sx) | e; |
z[1] = evalsigne(sx) | evalexpo(ey - (m | (i<<TWOPOTBITS_IN_LONG))); |
z[0] = evaltyp(t_REAL) | evallg(lz); |
z[0] = evaltyp(t_REAL) | evallg(lz); |
avma = (long)z; return z; |
avma = (gpmem_t)z; return z; |
} |
} |
|
|
GEN |
GEN |
Line 719 mulsispec(long x, GEN y, long ny) |
|
Line 750 mulsispec(long x, GEN y, long ny) |
|
if (hiremainder) *--z = hiremainder; else lz--; |
if (hiremainder) *--z = hiremainder; else lz--; |
*--z = evalsigne(1) | evallgefint(lz); |
*--z = evalsigne(1) | evallgefint(lz); |
*--z = evaltyp(t_INT) | evallg(lz); |
*--z = evaltyp(t_INT) | evallg(lz); |
avma=(long)z; return z; |
avma=(gpmem_t)z; return z; |
} |
} |
|
|
GEN |
GEN |
Line 755 addsmulsi(long a, long b, GEN Y) |
|
Line 786 addsmulsi(long a, long b, GEN Y) |
|
if (hiremainder) *--z = hiremainder; else lz--; |
if (hiremainder) *--z = hiremainder; else lz--; |
*--z = evalsigne(1) | evallgefint(lz); |
*--z = evalsigne(1) | evallgefint(lz); |
*--z = evaltyp(t_INT) | evallg(lz); |
*--z = evaltyp(t_INT) | evallg(lz); |
avma=(long)z; return z; |
avma=(gpmem_t)z; return z; |
} |
} |
|
|
#ifndef __M68K__ |
#ifndef __M68K__ |
Line 784 mulsr(long x, GEN y) |
|
Line 815 mulsr(long x, GEN y) |
|
if (!s) |
if (!s) |
{ |
{ |
if (x<0) x = -x; |
if (x<0) x = -x; |
e = y[1] + (BITS_IN_LONG-1)-bfffo(x); |
e = expo(y) + (BITS_IN_LONG-1)-bfffo(x); |
if (e & ~EXPOBITS) err(muler2); |
return realzero_bit(e); |
z=cgetr(3); z[1]=e; z[2]=0; return z; |
|
} |
} |
if (x<0) { s = -s; x = -x; } |
if (x<0) { s = -s; x = -x; } |
if (x==1) { z=rcopy(y); setsigne(z,s); return z; } |
if (x==1) { z=rcopy(y); setsigne(z,s); return z; } |
Line 798 mulsr(long x, GEN y) |
|
Line 828 mulsr(long x, GEN y) |
|
|
|
sh = bfffo(hiremainder); m = BITS_IN_LONG-sh; |
sh = bfffo(hiremainder); m = BITS_IN_LONG-sh; |
if (sh) shift_left2(z,z, 2,lx-1, garde,sh,m); |
if (sh) shift_left2(z,z, 2,lx-1, garde,sh,m); |
e = evalexpo(m+e); |
z[1] = evalsigne(s) | evalexpo(m+e); return z; |
z[1] = evalsigne(s) | e; return z; |
|
} |
} |
|
|
GEN |
GEN |
Line 812 mulrr(GEN x, GEN y) |
|
Line 841 mulrr(GEN x, GEN y) |
|
LOCAL_HIREMAINDER; |
LOCAL_HIREMAINDER; |
LOCAL_OVERFLOW; |
LOCAL_OVERFLOW; |
|
|
e = evalexpo(expo(x)+expo(y)); |
e = expo(x)+expo(y); |
if (!sx || !sy) { z=cgetr(3); z[1]=e; z[2]=0; return z; } |
if (!sx || !sy) return realzero_bit(e); |
if (sy<0) sx = -sx; |
if (sy<0) sx = -sx; |
lz=lg(x); ly=lg(y); |
lz=lg(x); ly=lg(y); |
if (lz>ly) { lz=ly; z=x; x=y; y=z; flag=1; } else flag = (lz!=ly); |
if (lz>ly) { lz=ly; z=x; x=y; y=z; flag=1; } else flag = (lz!=ly); |
z=cgetr(lz); z[1] = evalsigne(sx) | e; |
z=cgetr(lz); |
if (lz==3) |
if (lz==3) |
{ |
{ |
if (flag) |
if (flag) |
Line 827 mulrr(GEN x, GEN y) |
|
Line 856 mulrr(GEN x, GEN y) |
|
} |
} |
else |
else |
garde = mulll(x[2],y[2]); |
garde = mulll(x[2],y[2]); |
if ((long)hiremainder<0) { z[2]=hiremainder; z[1]++; } |
if ((long)hiremainder<0) { z[2]=hiremainder; e++; } |
else z[2]=(hiremainder<<1) | (garde>>(BITS_IN_LONG-1)); |
else z[2]=(hiremainder<<1) | (garde>>(BITS_IN_LONG-1)); |
|
z[1] = evalsigne(sx)|evalexpo(e); |
return z; |
return z; |
} |
} |
|
|
Line 865 mulrr(GEN x, GEN y) |
|
Line 895 mulrr(GEN x, GEN y) |
|
z[i] = addll(addmul(p1,y1[i]), z[i]); |
z[i] = addll(addmul(p1,y1[i]), z[i]); |
} |
} |
z[2] = hiremainder+overflow; |
z[2] = hiremainder+overflow; |
if (z[2] < 0) z[1]++; |
if (z[2] < 0) e++; else shift_left(z,z,2,lzz,garde, 1); |
else |
z[1] = evalsigne(sx) | evalexpo(e); |
shift_left(z,z,2,lzz,garde, 1); |
|
return z; |
return z; |
} |
} |
|
|
Line 883 mulir(GEN x, GEN y) |
|
Line 912 mulir(GEN x, GEN y) |
|
if (!sx) return gzero; |
if (!sx) return gzero; |
if (!is_bigint(x)) return mulsr(itos(x),y); |
if (!is_bigint(x)) return mulsr(itos(x),y); |
sy=signe(y); ey=expo(y); |
sy=signe(y); ey=expo(y); |
if (!sy) |
if (!sy) return realzero_bit(expi(x)+ey); |
{ |
|
e = evalexpo(expi(x)+ey); |
|
z=cgetr(3); z[1]=e; z[2]=0; return z; |
|
} |
|
|
|
if (sy<0) sx = -sx; |
if (sy<0) sx = -sx; |
lz=lg(y); z=cgetr(lz); |
lz=lg(y); z=cgetr(lz); |
y1=cgetr(lz+1); |
y1 = itor(x, lz+1); x = y; y = y1; |
affir(x,y1); x=y; y=y1; |
e = expo(y)+ey; |
e = evalexpo(expo(y)+ey); |
|
z[1] = evalsigne(sx) | e; |
|
if (lz==3) |
if (lz==3) |
{ |
{ |
(void)mulll(x[2],y[3]); |
(void)mulll(x[2],y[3]); |
garde=addmul(x[2],y[2]); |
garde=addmul(x[2],y[2]); |
if ((long)hiremainder < 0) { z[2]=hiremainder; z[1]++; } |
if ((long)hiremainder < 0) { z[2]=hiremainder; e++; } |
else z[2]=(hiremainder<<1) | (garde>>(BITS_IN_LONG-1)); |
else z[2]=(hiremainder<<1) | (garde>>(BITS_IN_LONG-1)); |
avma=(long)z; return z; |
z[1] = evalsigne(sx) | evalexpo(e); |
|
avma=(gpmem_t)z; return z; |
} |
} |
|
|
(void)mulll(x[2],y[lz]); garde=hiremainder; |
(void)mulll(x[2],y[lz]); garde=hiremainder; |
Line 937 mulir(GEN x, GEN y) |
|
Line 961 mulir(GEN x, GEN y) |
|
z[i] = addll(addmul(p1,y1[i]), z[i]); |
z[i] = addll(addmul(p1,y1[i]), z[i]); |
} |
} |
z[2] = hiremainder+overflow; |
z[2] = hiremainder+overflow; |
if (z[2] < 0) z[1]++; |
if (z[2] < 0) e++; |
else |
else |
shift_left(z,z,2,lzz,garde, 1); |
shift_left(z,z,2,lzz,garde, 1); |
avma=(long)z; return z; |
z[1] = evalsigne(sx) | evalexpo(e); |
|
avma=(gpmem_t)z; return z; |
} |
} |
|
|
/* written by Bruno Haible following an idea of Robert Harley */ |
/* written by Bruno Haible following an idea of Robert Harley */ |
Line 974 modss(long x, long y) |
|
Line 999 modss(long x, long y) |
|
|
|
if (!y) err(moder1); |
if (!y) err(moder1); |
if (y<0) y=-y; |
if (y<0) y=-y; |
hiremainder=0; divll(labs(x),y); |
hiremainder=0; (void)divll(labs(x),y); |
if (!hiremainder) return gzero; |
if (!hiremainder) return gzero; |
return (x < 0) ? stoi(y-hiremainder) : stoi(hiremainder); |
return (x < 0) ? stoi(y-hiremainder) : stoi(hiremainder); |
} |
} |
Line 985 resss(long x, long y) |
|
Line 1010 resss(long x, long y) |
|
LOCAL_HIREMAINDER; |
LOCAL_HIREMAINDER; |
|
|
if (!y) err(reser1); |
if (!y) err(reser1); |
hiremainder=0; divll(labs(x),labs(y)); |
hiremainder=0; (void)divll(labs(x),labs(y)); |
return (x < 0) ? stoi(-((long)hiremainder)) : stoi(hiremainder); |
return (x < 0) ? stoi(-((long)hiremainder)) : stoi(hiremainder); |
} |
} |
|
|
Line 1044 umodiu(GEN y, ulong x) |
|
Line 1069 umodiu(GEN y, ulong x) |
|
GEN |
GEN |
modiu(GEN y, ulong x) { return utoi(umodiu(y,x)); } |
modiu(GEN y, ulong x) { return utoi(umodiu(y,x)); } |
|
|
|
/* return |y| \/ x */ |
|
GEN |
|
diviu(GEN y, ulong x) |
|
{ |
|
long sy=signe(y),ly,i; |
|
GEN z; |
|
LOCAL_HIREMAINDER; |
|
|
|
if (!x) err(diver4); |
|
if (!sy) { hiremainder=0; SAVE_HIREMAINDER; return gzero; } |
|
|
|
ly = lgefint(y); |
|
if (x <= (ulong)y[2]) hiremainder=0; |
|
else |
|
{ |
|
if (ly==3) { hiremainder=itou(y); SAVE_HIREMAINDER; return gzero; } |
|
hiremainder=y[2]; ly--; y++; |
|
} |
|
z = cgeti(ly); z[1] = evallgefint(ly) | evalsigne(1); |
|
for (i=2; i<ly; i++) z[i]=divll(y[i],x); |
|
SAVE_HIREMAINDER; return z; |
|
} |
|
|
#ifndef __M68K__ |
#ifndef __M68K__ |
|
|
GEN |
GEN |
Line 1096 divis(GEN y, long x) |
|
Line 1144 divis(GEN y, long x) |
|
GEN |
GEN |
divir(GEN x, GEN y) |
divir(GEN x, GEN y) |
{ |
{ |
GEN xr,z; |
GEN z; |
long av,ly; |
long ly; |
|
gpmem_t av; |
|
|
if (!signe(y)) err(diver5); |
if (!signe(y)) err(diver5); |
if (!signe(x)) return gzero; |
if (!signe(x)) return gzero; |
ly=lg(y); z=cgetr(ly); av=avma; xr=cgetr(ly+1); affir(x,xr); |
ly = lg(y); z = cgetr(ly); av = avma; |
affrr(divrr(xr,y),z); avma=av; return z; |
affrr(divrr(itor(x, ly+1), y), z); |
|
avma = av; return z; |
} |
} |
|
|
GEN |
GEN |
divri(GEN x, GEN y) |
divri(GEN x, GEN y) |
{ |
{ |
GEN yr,z; |
long lx, s = signe(y); |
long av,lx,s=signe(y); |
gpmem_t av; |
|
GEN z; |
|
|
if (!s) err(diver8); |
if (!s) err(diver8); |
if (!signe(x)) |
if (!signe(x)) return realzero_bit(expo(x) - expi(y)); |
{ |
|
const long ex = evalexpo(expo(x) - expi(y)); |
|
|
|
if (ex<0) err(diver12); |
|
z=cgetr(3); z[1]=ex; z[2]=0; return z; |
|
} |
|
if (!is_bigint(y)) return divrs(x, s>0? y[2]: -y[2]); |
if (!is_bigint(y)) return divrs(x, s>0? y[2]: -y[2]); |
|
|
lx=lg(x); z=cgetr(lx); |
lx = lg(x); z = cgetr(lx); av = avma; |
av=avma; yr=cgetr(lx+1); affir(y,yr); |
affrr(divrr(x, itor(y, lx+1)), z); |
affrr(divrr(x,yr),z); avma=av; return z; |
avma = av; return z; |
} |
} |
|
|
void |
void |
diviiz(GEN x, GEN y, GEN z) |
diviiz(GEN x, GEN y, GEN z) |
{ |
{ |
long av=avma,lz; |
gpmem_t av = avma; |
GEN xr,yr; |
if (typ(z) == t_INT) affii(divii(x,y), z); |
|
else { |
if (typ(z)==t_INT) { affii(divii(x,y),z); avma=av; return; } |
long lz = lg(z); |
lz=lg(z); xr=cgetr(lz); affir(x,xr); yr=cgetr(lz); affir(y,yr); |
affrr(divrr(itor(x,lz), itor(y,lz)), z); |
affrr(divrr(xr,yr),z); avma=av; |
} |
|
avma = av; |
} |
} |
|
|
void |
void |
mpdivz(GEN x, GEN y, GEN z) |
mpdivz(GEN x, GEN y, GEN z) |
{ |
{ |
long av=avma; |
gpmem_t av = avma; |
|
GEN r; |
|
|
if (typ(z)==t_INT) |
if (typ(z)==t_INT) |
{ |
{ |
if (typ(x)==t_REAL || typ(y)==t_REAL) err(divzer1); |
if (typ(x) == t_REAL || typ(y) == t_REAL) err(divzer1); |
affii(divii(x,y),z); avma=av; return; |
affii(divii(x,y), z); |
|
avma = av; return; |
} |
} |
if (typ(x)==t_INT) |
|
{ |
|
GEN xr,yr; |
|
long lz; |
|
|
|
if (typ(y)==t_REAL) { affrr(divir(x,y),z); avma=av; return; } |
if (typ(x) == t_INT) |
lz=lg(z); xr=cgetr(lz); affir(x,xr); yr=cgetr(lz); affir(y,yr); |
{ |
affrr(divrr(xr,yr),z); avma=av; return; |
if (typ(y) == t_REAL) |
|
r = divir(x,y); |
|
else |
|
{ |
|
long lz = lg(z); |
|
r = divrr(itor(x,lz), itor(y,lz)); |
|
} |
} |
} |
if (typ(y)==t_REAL) { affrr(divrr(x,y),z); avma=av; return; } |
else if (typ(y) == t_REAL) |
affrr(divri(x,y),z); avma=av; |
r = divrr(x,y); |
|
else |
|
r = divri(x,y); |
|
affrr(r, z); avma = av; |
} |
} |
|
|
GEN |
GEN |
divsr(long x, GEN y) |
divsr(long x, GEN y) |
{ |
{ |
long av,ly; |
gpmem_t av; |
GEN p1,z; |
long ly; |
|
GEN z; |
|
|
if (!signe(y)) err(diver3); |
if (!signe(y)) err(diver3); |
if (!x) return gzero; |
if (!x) return gzero; |
ly=lg(y); z=cgetr(ly); av=avma; |
ly = lg(y); z = cgetr(ly); av = avma; |
p1=cgetr(ly+1); affsr(x,p1); affrr(divrr(p1,y),z); |
affrr(divrr(stor(x,ly+1), y), z); |
avma=av; return z; |
avma = av; return z; |
} |
} |
|
|
GEN |
GEN |
Line 1182 modii(GEN x, GEN y) |
|
Line 1236 modii(GEN x, GEN y) |
|
case 1: return resii(x,y); |
case 1: return resii(x,y); |
default: |
default: |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
(void)new_chunk(lgefint(y)); |
(void)new_chunk(lgefint(y)); |
x = resii(x,y); avma=av; |
x = resii(x,y); avma=av; |
if (x==gzero) return x; |
if (x==gzero) return x; |
Line 1194 modii(GEN x, GEN y) |
|
Line 1248 modii(GEN x, GEN y) |
|
void |
void |
modiiz(GEN x, GEN y, GEN z) |
modiiz(GEN x, GEN y, GEN z) |
{ |
{ |
const long av = avma; |
const gpmem_t av = avma; |
affii(modii(x,y),z); avma=av; |
affii(modii(x,y),z); avma=av; |
} |
} |
|
|
GEN |
GEN |
divrs(GEN x, long y) |
divrs(GEN x, long y) |
{ |
{ |
long i,lx,ex,garde,sh,s=signe(x); |
long i,lx,garde,sh,s=signe(x); |
GEN z; |
GEN z; |
LOCAL_HIREMAINDER; |
LOCAL_HIREMAINDER; |
|
|
if (!y) err(diver6); |
if (!y) err(diver6); |
if (!s) |
if (!s) return realzero_bit(expo(x) - (BITS_IN_LONG-1)+bfffo(y)); |
{ |
|
z=cgetr(3); z[1] = x[1] - (BITS_IN_LONG-1) + bfffo(y); |
|
if (z[1]<0) err(diver7); |
|
z[2]=0; return z; |
|
} |
|
if (y<0) { s = -s; y = -y; } |
if (y<0) { s = -s; y = -y; } |
if (y==1) { z=rcopy(x); setsigne(z,s); return z; } |
if (y==1) { z=rcopy(x); setsigne(z,s); return z; } |
|
|
Line 1219 divrs(GEN x, long y) |
|
Line 1268 divrs(GEN x, long y) |
|
for (i=2; i<lx; i++) z[i]=divll(x[i],y); |
for (i=2; i<lx; i++) z[i]=divll(x[i],y); |
|
|
/* we may have hiremainder != 0 ==> garde */ |
/* we may have hiremainder != 0 ==> garde */ |
garde=divll(0,y); sh=bfffo(z[2]); ex=evalexpo(expo(x)-sh); |
garde=divll(0,y); sh=bfffo(z[2]); |
if (sh) shift_left(z,z, 2,lx-1, garde,sh); |
if (sh) shift_left(z,z, 2,lx-1, garde,sh); |
z[1] = evalsigne(s) | ex; |
z[1] = evalsigne(s) | evalexpo(expo(x)-sh); |
return z; |
return z; |
} |
} |
|
|
Line 1233 divrr(GEN x, GEN y) |
|
Line 1282 divrr(GEN x, GEN y) |
|
GEN z,x1; |
GEN z,x1; |
|
|
if (!sy) err(diver9); |
if (!sy) err(diver9); |
e = evalexpo(expo(x) - expo(y)); |
e = expo(x) - expo(y); |
if (!sx) { z=cgetr(3); z[1]=e; z[2]=0; return z; } |
if (!sx) return realzero_bit(e); |
if (sy<0) sx = -sx; |
if (sy<0) sx = -sx; |
e = evalsigne(sx) | e; |
|
lx=lg(x); ly=lg(y); |
lx=lg(x); ly=lg(y); |
if (ly==3) |
if (ly==3) |
{ |
{ |
Line 1248 divrr(GEN x, GEN y) |
|
Line 1296 divrr(GEN x, GEN y) |
|
l >>= 1; if (k&1) l |= HIGHBIT; |
l >>= 1; if (k&1) l |= HIGHBIT; |
k >>= 1; |
k >>= 1; |
} |
} |
z=cgetr(3); z[1]=e; |
z = cgetr(3); z[1] = evalsigne(sx) | evalexpo(e); |
hiremainder=k; z[2]=divll(l,y[2]); return z; |
hiremainder=k; z[2]=divll(l,y[2]); return z; |
} |
} |
|
|
lz = (lx<=ly)? lx: ly; |
lz = min(lx,ly); z = new_chunk(lz); |
x1 = (z=new_chunk(lz))-1; |
x1 = z-1; |
x1[1]=0; for (i=2; i<lz; i++) x1[i]=x[i]; |
x1[1]=0; for (i=2; i<lz; i++) x1[i]=x[i]; |
x1[lz] = (lx>lz)? x[lz]: 0; |
x1[lz] = (lx>lz)? x[lz]: 0; |
x=z; si=y[2]; saux=y[3]; |
x=z; si=y[2]; saux=y[3]; |
Line 1313 divrr(GEN x, GEN y) |
|
Line 1361 divrr(GEN x, GEN y) |
|
} |
} |
x1 = x-1; for (j=lz-1; j>=2; j--) x[j]=x1[j]; |
x1 = x-1; for (j=lz-1; j>=2; j--) x[j]=x1[j]; |
if (*x) { shift_right(x,x, 2,lz, 1,1); } else e--; |
if (*x) { shift_right(x,x, 2,lz, 1,1); } else e--; |
x[0]=evaltyp(t_REAL)|evallg(lz); |
x[0] = evaltyp(t_REAL)|evallg(lz); |
x[1]=e; return x; |
x[1] = evalsigne(sx) | evalexpo(e); |
|
return x; |
} |
} |
#endif /* !defined(__M68K__) */ |
#endif /* !defined(__M68K__) */ |
|
|
/* The following ones are not in mp.s (mulii is, with a different algorithm) */ |
/* The following ones are not in mp.s (mulii is, with a different algorithm) */ |
|
|
/* Integer division x / y: |
/* Integer division x / y: such that sign(r) = sign(x) |
* if z = ONLY_REM return remainder, otherwise return quotient |
* if z = ONLY_REM return remainder, otherwise return quotient |
* if z != NULL set *z to remainder |
* if z != NULL set *z to remainder |
* *z is the last object on stack (and thus can be disposed of with cgiv |
* *z is the last object on stack (and thus can be disposed of with cgiv |
|
|
dvmdii(GEN x, GEN y, GEN *z) |
dvmdii(GEN x, GEN y, GEN *z) |
{ |
{ |
long sx=signe(x),sy=signe(y); |
long sx=signe(x),sy=signe(y); |
long av,lx,ly,lz,i,j,sh,k,lq,lr; |
long lx, ly, lz, i, j, sh, k, lq, lr; |
|
gpmem_t av; |
ulong si,qp,saux, *xd,*rd,*qd; |
ulong si,qp,saux, *xd,*rd,*qd; |
GEN r,q,x1; |
GEN r,q,x1; |
|
|
Line 1447 DIVIDE: /* quotient is non-zero */ |
|
Line 1498 DIVIDE: /* quotient is non-zero */ |
|
while (lz--) *--qd = *--xd; |
while (lz--) *--qd = *--xd; |
*--qd = evalsigne(sy) | evallgefint(lq); |
*--qd = evalsigne(sy) | evallgefint(lq); |
*--qd = evaltyp(t_INT) | evallg(lq); |
*--qd = evaltyp(t_INT) | evallg(lq); |
avma = (long)qd; return (GEN)qd; |
avma = (gpmem_t)qd; return (GEN)qd; |
} |
} |
|
|
j=lq; while (j<lx && !x[j]) j++; |
j=lq; while (j<lx && !x[j]) j++; |
Line 1473 DIVIDE: /* quotient is non-zero */ |
|
Line 1524 DIVIDE: /* quotient is non-zero */ |
|
} |
} |
*--rd = evalsigne(sx) | evallgefint(lr); |
*--rd = evalsigne(sx) | evallgefint(lr); |
*--rd = evaltyp(t_INT) | evallg(lr); |
*--rd = evaltyp(t_INT) | evallg(lr); |
avma = (long)rd; return (GEN)rd; |
avma = (gpmem_t)rd; return (GEN)rd; |
} |
} |
|
|
lr = lz+2; |
lr = lz+2; |
Line 1517 DIVIDE: /* quotient is non-zero */ |
|
Line 1568 DIVIDE: /* quotient is non-zero */ |
|
while (lr--) *--qd = *--rd; |
while (lr--) *--qd = *--rd; |
*z = (GEN)qd; |
*z = (GEN)qd; |
} |
} |
avma = (long)qd; return q; |
avma = (gpmem_t)qd; return q; |
} |
} |
|
|
/* assume y > x > 0. return y mod x */ |
/* assume y > x > 0. return y mod x */ |
Line 1536 mppgcd_resiu(GEN y, ulong x) |
|
Line 1587 mppgcd_resiu(GEN y, ulong x) |
|
void |
void |
mppgcd_plus_minus(GEN x, GEN y, GEN res) |
mppgcd_plus_minus(GEN x, GEN y, GEN res) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
long lx = lgefint(x)-1; |
long lx = lgefint(x)-1; |
long ly = lgefint(y)-1, lt,m,i; |
long ly = lgefint(y)-1, lt,m,i; |
GEN t; |
GEN t; |
Line 1579 smulss(ulong x, ulong y, ulong *rem) |
|
Line 1630 smulss(ulong x, ulong y, ulong *rem) |
|
# define DIVCONVI 14 |
# define DIVCONVI 14 |
#endif |
#endif |
|
|
/* convert integer --> base 10^9 [assume x > 0] */ |
/* convert integer --> base 10^9 */ |
GEN |
GEN |
convi(GEN x) |
convi(GEN x) |
{ |
{ |
ulong av=avma, lim; |
gpmem_t av = avma, lim; |
long lz; |
long lz; |
GEN z,p1; |
GEN z,p1; |
|
|
lz = 3 + ((lgefint(x)-2)*15)/DIVCONVI; |
lz = 3 + ((lgefint(x)-2)*15)/DIVCONVI; |
z=new_chunk(lz); z[1] = -1; p1 = z+2; |
z = new_chunk(lz); z[1] = -1; p1 = z+2; |
lim = stack_lim(av,1); |
lim = stack_lim(av,1); |
for(;;) |
for(;;) |
{ |
{ |
x = divis(x,1000000000); *p1++ = hiremainder; |
x = diviu(x,1000000000); *p1++ = hiremainder; |
if (!signe(x)) { avma=av; return p1; } |
if (!signe(x)) { avma = av; return p1; } |
if (low_stack(lim, stack_lim(av,1))) x = gerepileuptoint((long)z,x); |
if (low_stack(lim, stack_lim(av,1))) x = gerepileuptoint((gpmem_t)z,x); |
} |
} |
} |
} |
#undef DIVCONVI |
#undef DIVCONVI |
|
|
GEN |
GEN |
truedvmdii(GEN x, GEN y, GEN *z) |
truedvmdii(GEN x, GEN y, GEN *z) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
GEN r, q = dvmdii(x,y,&r); /* assume that r is last on stack */ |
GEN r, q = dvmdii(x,y,&r); /* assume that r is last on stack */ |
GEN *gptr[2]; |
GEN *gptr[2]; |
|
|
Line 1653 truedvmdii(GEN x, GEN y, GEN *z) |
|
Line 1704 truedvmdii(GEN x, GEN y, GEN *z) |
|
} |
} |
|
|
if (z == ONLY_REM) |
if (z == ONLY_REM) |
{ |
{ /* r += sign(y) * y, that is |y| */ |
r = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2); |
r = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2); |
return gerepileuptoint(av, r); |
return gerepileuptoint(av, r); |
} |
} |
Line 1661 truedvmdii(GEN x, GEN y, GEN *z) |
|
Line 1712 truedvmdii(GEN x, GEN y, GEN *z) |
|
if (!z) return gerepileuptoint(av, q); |
if (!z) return gerepileuptoint(av, q); |
|
|
*z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2); |
*z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2); |
gptr[0]=&q; gptr[1]=z; gerepilemanysp(av,(long)r,gptr,2); |
gptr[0]=&q; gptr[1]=z; gerepilemanysp(av,(gpmem_t)r,gptr,2); |
return q; |
return q; |
} |
} |
|
|
|
/* Montgomery reduction. |
|
* N has k words, assume T >= 0 has less than 2k. |
|
* Return res := T / B^k mod N, where B = 2^BIL |
|
* such that 0 <= res < T/B^k + N and res has less than k words */ |
|
GEN |
|
red_montgomery(GEN T, GEN N, ulong inv) |
|
{ |
|
gpmem_t av; |
|
GEN Te,Td,Ne,Nd, scratch; |
|
ulong m,t,d,k = lgefint(N)-2; |
|
int carry; |
|
long i,j; |
|
LOCAL_HIREMAINDER; |
|
LOCAL_OVERFLOW; |
|
|
|
if (k == 0) return gzero; |
|
d = lgefint(T)-2; /* <= 2*k */ |
|
#ifdef DEBUG |
|
if (d > 2*k) err(bugparier,"red_montgomery"); |
|
#endif |
|
if (k == 1) |
|
{ /* as below, special cased for efficiency */ |
|
ulong n = (ulong)N[2]; |
|
t = (ulong)T[d+1]; |
|
m = t * inv; |
|
(void)addll(mulll(m, n), t); /* = 0 */ |
|
t = hiremainder + overflow; |
|
if (d == 2) |
|
{ |
|
t = addll((ulong)T[2], t); |
|
if (overflow) t -= n; /* t > n doesn't fit in 1 word */ |
|
} |
|
return utoi(t); |
|
} |
|
/* assume k >= 2 */ |
|
av = avma; scratch = new_chunk(k<<1); /* >= k + 2: result fits */ |
|
|
|
/* copy T to scratch space (pad with zeroes to 2k words) */ |
|
Td = (GEN)av; |
|
Te = T + (d+2); |
|
for (i=0; i < d ; i++) *--Td = *--Te; |
|
for ( ; i < (k<<1); i++) *--Td = 0; |
|
|
|
Te = (GEN)av; /* 1 beyond end of T mantissa */ |
|
Ne = N + k+2; /* 1 beyond end of N mantissa */ |
|
|
|
carry = 0; |
|
for (i=0; i<k; i++) /* set T := T/B nod N, k times */ |
|
{ |
|
Td = Te; /* one beyond end of (new) T mantissa */ |
|
Nd = Ne; |
|
m = *--Td * inv; /* solve T + m N = O(B) */ |
|
|
|
/* set T := (T + mN) / B */ |
|
Te = Td; |
|
(void)addll(mulll(m, *--Nd), *Td); /* = 0 */ |
|
for (j=1; j<k; j++) |
|
{ |
|
hiremainder += overflow; |
|
t = addll(addmul(m, *--Nd), *--Td); *Td = t; |
|
} |
|
overflow += hiremainder; |
|
t = addll(overflow, *--Td); *Td = t + carry; |
|
carry = (overflow || (carry && *Td == 0)); |
|
} |
|
if (carry) |
|
{ /* Td > N overflows (k+1 words), set Td := Td - N */ |
|
Td = Te; |
|
Nd = Ne; |
|
t = subll(*--Td, *--Nd); *Td = t; |
|
while (Td > scratch) { t = subllx(*--Td, *--Nd); *Td = t; } |
|
} |
|
|
|
/* copy result */ |
|
Td = (GEN)av; |
|
while (! *scratch) scratch++; /* strip leading zeroes */ |
|
while (Te > scratch) *--Td = *--Te; |
|
k = ((GEN)av - Td) + 2; |
|
*--Td = evalsigne(1) | evallgefint(k); |
|
*--Td = evaltyp(t_INT) | evallg(k); |
|
#ifdef DEBUG |
|
{ |
|
long l = lgefint(N)-2, s = BITS_IN_LONG*l; |
|
GEN R = shifti(gun, s); |
|
GEN res = resii(mulii(T, mpinvmod(R, N)), N); |
|
if (k > lgefint(N) |
|
|| !egalii(resii(Td,N),res) |
|
|| cmpii(Td, addii(shifti(T, -s), N)) >= 0) err(bugparier,"red_montgomery"); |
|
} |
|
#endif |
|
avma = (gpmem_t)Td; return Td; |
|
} |
|
|
/* EXACT INTEGER DIVISION */ |
/* EXACT INTEGER DIVISION */ |
|
|
/* Find c such that 1=c*b mod 2^BITS_IN_LONG, assuming b odd (unchecked) */ |
/* Find c such that 1=c*b mod 2^BITS_IN_LONG, assuming b odd (unchecked) */ |
static ulong |
ulong |
invrev(ulong b) |
invrev(ulong b) |
{ |
{ |
ulong x; |
ulong x; |
Line 1688 invrev(ulong b) |
|
Line 1832 invrev(ulong b) |
|
} |
} |
|
|
/* assume xy>0, y odd */ |
/* assume xy>0, y odd */ |
GEN |
GEN |
diviuexact(GEN x, ulong y) |
diviuexact(GEN x, ulong y) |
{ |
{ |
long i,lz,lx; |
long i,lz,lx; |
Line 1697 diviuexact(GEN x, ulong y) |
|
Line 1841 diviuexact(GEN x, ulong y) |
|
|
|
if (y == 1) return icopy(x); |
if (y == 1) return icopy(x); |
lx = lgefint(x); |
lx = lgefint(x); |
if (lx == 3) return stoi((ulong)x[2] / y); |
if (lx == 3) return utoi((ulong)x[2] / y); |
yinv = invrev(y); |
yinv = invrev(y); |
lz = (y <= (ulong)x[2]) ? lx : lx-1; |
lz = (y <= (ulong)x[2]) ? lx : lx-1; |
z = new_chunk(lz); |
z = new_chunk(lz); |
Line 1729 diviuexact(GEN x, ulong y) |
|
Line 1873 diviuexact(GEN x, ulong y) |
|
z += i-2; lz -= i-2; |
z += i-2; lz -= i-2; |
z[0] = evaltyp(t_INT)|evallg(lz); |
z[0] = evaltyp(t_INT)|evallg(lz); |
z[1] = evalsigne(1)|evallg(lz); |
z[1] = evalsigne(1)|evallg(lz); |
avma = (ulong)z; return z; |
avma = (gpmem_t)z; return z; |
} |
} |
|
|
/* Find z such that x=y*z, knowing that y | x (unchecked) |
/* Find z such that x=y*z, knowing that y | x (unchecked) |
Line 1738 diviuexact(GEN x, ulong y) |
|
Line 1882 diviuexact(GEN x, ulong y) |
|
GEN |
GEN |
diviiexact(GEN x, GEN y) |
diviiexact(GEN x, GEN y) |
{ |
{ |
long av,lx,ly,lz,vy,i,ii,sx = signe(x), sy = signe(y); |
long lx, ly, lz, vy, i, ii, sx = signe(x), sy = signe(y); |
|
gpmem_t av; |
ulong y0inv,q; |
ulong y0inv,q; |
GEN z; |
GEN z; |
|
|
Line 1758 diviiexact(GEN x, GEN y) |
|
Line 1903 diviiexact(GEN x, GEN y) |
|
avma = av; /* will erase our x,y when exiting */ |
avma = av; /* will erase our x,y when exiting */ |
/* now y is odd */ |
/* now y is odd */ |
ly = lgefint(y); |
ly = lgefint(y); |
if (ly == 3) |
if (ly == 3) |
{ |
{ |
x = diviuexact(x,(ulong)y[2]); |
x = diviuexact(x,(ulong)y[2]); |
if (signe(x)) setsigne(x,sx*sy); /* should have x != 0 at this point */ |
if (signe(x)) setsigne(x,sx*sy); /* should have x != 0 at this point */ |
Line 1783 diviiexact(GEN x, GEN y) |
|
Line 1928 diviiexact(GEN x, GEN y) |
|
/* x := x - q * y */ |
/* x := x - q * y */ |
(void)mulll(q,y[0]); limj = max(lx - lz, ii+3-ly); |
(void)mulll(q,y[0]); limj = max(lx - lz, ii+3-ly); |
{ /* update neither lowest word (could set it to 0) nor highest ones */ |
{ /* update neither lowest word (could set it to 0) nor highest ones */ |
register GEN x0 = x + (ii - 1), y0 = y - 1, xlim = x + limj; |
register GEN x0 = x + (ii - 1), y0 = y - 1, xlim = x + limj; |
for (; x0 >= xlim; x0--, y0--) |
for (; x0 >= xlim; x0--, y0--) |
{ |
{ |
*x0 = subll(*x0, addmul(q,*y0)); |
*x0 = subll(*x0, addmul(q,*y0)); |
Line 1810 diviiexact(GEN x, GEN y) |
|
Line 1955 diviiexact(GEN x, GEN y) |
|
z += i-2; lz -= (i-2); |
z += i-2; lz -= (i-2); |
z[0] = evaltyp(t_INT)|evallg(lz); |
z[0] = evaltyp(t_INT)|evallg(lz); |
z[1] = evalsigne(sx*sy)|evallg(lz); |
z[1] = evalsigne(sx*sy)|evallg(lz); |
avma = (ulong)z; return z; |
avma = (gpmem_t)z; return z; |
} |
} |
|
|
long |
long |
Line 1884 absr_cmp(GEN x, GEN y) |
|
Line 2029 absr_cmp(GEN x, GEN y) |
|
#define _sqri_l 47 |
#define _sqri_l 47 |
#define _muli_l 25 /* optimal on PII 350MHz + gcc 2.7.2.1 (gp-dyn) */ |
#define _muli_l 25 /* optimal on PII 350MHz + gcc 2.7.2.1 (gp-dyn) */ |
|
|
#if 0 /* for tunings */ |
#if 1 /* for tunings */ |
long KARATSUBA_SQRI_LIMIT = _sqri_l; |
long KARATSUBA_SQRI_LIMIT = _sqri_l; |
long KARATSUBA_MULI_LIMIT = _muli_l; |
long KARATSUBA_MULI_LIMIT = _muli_l; |
|
|
Line 1945 muliispec(GEN x, GEN y, long nx, long ny) |
|
Line 2090 muliispec(GEN x, GEN y, long nx, long ny) |
|
if (*zd == 0) { zd++; lz--; } /* normalize */ |
if (*zd == 0) { zd++; lz--; } /* normalize */ |
*--zd = evalsigne(1) | evallgefint(lz); |
*--zd = evalsigne(1) | evallgefint(lz); |
*--zd = evaltyp(t_INT) | evallg(lz); |
*--zd = evaltyp(t_INT) | evallg(lz); |
avma=(long)zd; return zd; |
avma=(gpmem_t)zd; return zd; |
} |
} |
|
|
#ifdef INLINE |
#ifdef INLINE |
Line 1967 sqrispec(GEN x, long nx) |
|
Line 2112 sqrispec(GEN x, long nx) |
|
*--zd = hiremainder; goto END; |
*--zd = hiremainder; goto END; |
} |
} |
xd = x + nx; |
xd = x + nx; |
|
|
/* compute double products --> zd */ |
/* compute double products --> zd */ |
p1 = *--xd; yd = xd; --zd; |
p1 = *--xd; yd = xd; --zd; |
*--zd = mulll(p1, *--yd); z2e = zd; |
*--zd = mulll(p1, *--yd); z2e = zd; |
|
|
if (*zd == 0) { zd++; lz--; } /* normalize */ |
if (*zd == 0) { zd++; lz--; } /* normalize */ |
*--zd = evalsigne(1) | evallgefint(lz); |
*--zd = evalsigne(1) | evallgefint(lz); |
*--zd = evaltyp(t_INT) | evallg(lz); |
*--zd = evaltyp(t_INT) | evallg(lz); |
avma=(long)zd; return zd; |
avma=(gpmem_t)zd; return zd; |
} |
} |
|
|
/* return (x shifted left d words) + y. Assume d > 0, x > 0 and y >= 0 */ |
/* return (x shifted left d words) + y. Assume d > 0, x > 0 and y >= 0 */ |
|
|
quickmulii(GEN a, GEN b, long na, long nb) |
quickmulii(GEN a, GEN b, long na, long nb) |
{ |
{ |
GEN a0,c,c0; |
GEN a0,c,c0; |
long av,n0,n0a,i; |
long n0, n0a, i; |
|
gpmem_t av; |
|
|
if (na < nb) swapspec(a,b, na,nb); |
if (na < nb) swapspec(a,b, na,nb); |
if (nb == 1) return mulsispec(*b, a, na); |
if (nb == 1) return mulsispec(*b, a, na); |
|
|
resiimul(GEN x, GEN sy) |
resiimul(GEN x, GEN sy) |
{ |
{ |
GEN r, q, y = (GEN)sy[1], invy; |
GEN r, q, y = (GEN)sy[1], invy; |
long av = avma, k; |
long k; |
|
gpmem_t av = avma; |
|
|
k = cmpii(x, y); |
k = cmpii(x, y); |
if (k <= 0) return k? icopy(x): gzero; |
if (k <= 0) return k? icopy(x): gzero; |
Line 2140 resmod2n(GEN x, long n) |
|
Line 2287 resmod2n(GEN x, long n) |
|
{ |
{ |
long hi,l,k,lx,ly; |
long hi,l,k,lx,ly; |
GEN z, xd, zd; |
GEN z, xd, zd; |
|
|
if (!signe(x) || !n) return gzero; |
if (!signe(x) || !n) return gzero; |
|
|
l = n & (BITS_IN_LONG-1); /* n % BITS_IN_LONG */ |
l = n & (BITS_IN_LONG-1); /* n % BITS_IN_LONG */ |
k = n >> TWOPOTBITS_IN_LONG; /* n / BITS_IN_LONG */ |
k = n >> TWOPOTBITS_IN_LONG; /* n / BITS_IN_LONG */ |
lx = lgefint(x); |
lx = lgefint(x); |
if (lx < k+3) return icopy(x); |
if (lx < k+3) return icopy(x); |
|
|
quicksqri(GEN a, long na) |
quicksqri(GEN a, long na) |
{ |
{ |
GEN a0,c; |
GEN a0,c; |
long av,n0,n0a,i; |
long n0, n0a, i; |
|
gpmem_t av; |
|
|
if (na < KARATSUBA_SQRI_LIMIT) return sqrispec(a,na); |
if (na < KARATSUBA_SQRI_LIMIT) return sqrispec(a,na); |
i=(na>>1); n0=na-i; na=i; |
i=(na>>1); n0=na-i; na=i; |
Line 2216 karamulrr1(GEN x, GEN y) |
|
Line 2364 karamulrr1(GEN x, GEN y) |
|
/* ensure that lg(y) >= lg(x) */ |
/* ensure that lg(y) >= lg(x) */ |
if (lx>ly) { lx=ly; z=x; x=y; y=z; flag=1; } else flag = (lx!=ly); |
if (lx>ly) { lx=ly; z=x; x=y; y=z; flag=1; } else flag = (lx!=ly); |
if (lx < MULRR_LIMIT) return mulrr(x,y); |
if (lx < MULRR_LIMIT) return mulrr(x,y); |
sx=signe(x); sy=signe(y); |
sx=signe(x); sy=signe(y); e = expo(x)+expo(y); |
e = evalexpo(expo(x)+expo(y)); |
if (!sx || !sy) return realzero_bit(e); |
if (!sx || !sy) { z=cgetr(3); z[2]=0; z[1]=e; return z; } |
|
if (sy<0) sx = -sx; |
if (sy<0) sx = -sx; |
ly=lx+flag; z=cgetr(lx); |
ly=lx+flag; z=cgetr(lx); |
lz2 = (lx>>1); lz3 = lx-lz2; |
lz2 = (lx>>1); lz3 = lx-lz2; |
Line 2246 karamulrr1(GEN x, GEN y) |
|
Line 2393 karamulrr1(GEN x, GEN y) |
|
garde = (hi[lx+2] << 1); |
garde = (hi[lx+2] << 1); |
shift_left(z,hi,2,lx+1, garde, 1); |
shift_left(z,hi,2,lx+1, garde, 1); |
} |
} |
z[1]=evalsigne(sx) | e; |
z[1]=evalsigne(sx) | evalexpo(e); |
if (garde < 0) |
if (garde < 0) |
{ /* round to nearest */ |
{ /* round to nearest */ |
i=lx+2; do z[--i]++; while (z[i]==0); |
i=lx+2; do z[--i]++; while (z[i]==0); |
if (i==1) z[2]=HIGHBIT; |
if (i==1) z[2]=HIGHBIT; |
} |
} |
avma=(long)z; return z; |
avma=(gpmem_t)z; return z; |
} |
} |
|
|
GEN |
GEN |
Line 2264 karamulrr2(GEN x, GEN y) |
|
Line 2411 karamulrr2(GEN x, GEN y) |
|
if (lx>ly) { lx=ly; z=x; x=y; y=z; flag=1; } else flag = (lx!=ly); |
if (lx>ly) { lx=ly; z=x; x=y; y=z; flag=1; } else flag = (lx!=ly); |
if (lx < MULRR2_LIMIT) return mulrr(x,y); |
if (lx < MULRR2_LIMIT) return mulrr(x,y); |
ly=lx+flag; sx=signe(x); sy=signe(y); |
ly=lx+flag; sx=signe(x); sy=signe(y); |
e = evalexpo(expo(x)+expo(y)); |
e = expo(x)+expo(y); |
if (!sx || !sy) { z=cgetr(3); z[2]=0; z[1]=e; return z; } |
if (!sx || !sy) return realzero_bit(e); |
if (sy<0) sx = -sx; |
if (sy<0) sx = -sx; |
z=cgetr(lx); |
z=cgetr(lx); |
hi=quickmulii(y+2,x+2,ly-2,lx-2); |
hi=quickmulii(y+2,x+2,ly-2,lx-2); |
Line 2279 karamulrr2(GEN x, GEN y) |
|
Line 2426 karamulrr2(GEN x, GEN y) |
|
garde = (hi[lx] << 1); |
garde = (hi[lx] << 1); |
shift_left(z,hi,2,lx-1, garde, 1); |
shift_left(z,hi,2,lx-1, garde, 1); |
} |
} |
z[1]=evalsigne(sx) | e; |
z[1]=evalsigne(sx) | evalexpo(e); |
if (garde < 0) |
if (garde < 0) |
{ /* round to nearest */ |
{ /* round to nearest */ |
i=lx; do z[--i]++; while (z[i]==0); |
i=lx; do z[--i]++; while (z[i]==0); |
if (i==1) z[2]=HIGHBIT; |
if (i==1) z[2]=HIGHBIT; |
} |
} |
avma=(long)z; return z; |
avma=(gpmem_t)z; return z; |
} |
} |
|
|
GEN |
GEN |
|
|
karamulir(GEN x, GEN y, long flag) |
karamulir(GEN x, GEN y, long flag) |
{ |
{ |
long sx=signe(x),lz,i; |
long sx=signe(x),lz,i; |
GEN z,temp,z1; |
GEN z, z1; |
|
|
if (!sx) return gzero; |
if (!sx) return gzero; |
lz=lg(y); z=cgetr(lz); |
lz = lg(y); z = cgetr(lz); |
temp=cgetr(lz+1); affir(x,temp); |
z1 = karamulrr(itor(x, lz+1), y, flag); |
z1=karamulrr(temp,y,flag); |
|
for (i=1; i<lz; i++) z[i]=z1[i]; |
for (i=1; i<lz; i++) z[i]=z1[i]; |
avma=(long)z; return z; |
avma=(gpmem_t)z; return z; |
} |
} |
#endif |
#endif |
|
|
#ifdef LONG_IS_64BIT |
#ifdef LONG_IS_64BIT |
|
long |
|
expodb(double x) |
|
{ |
|
union { double f; ulong i; } fi; |
|
const int mant_len = 52; /* mantissa bits (excl. hidden bit) */ |
|
const int exp_mid = 0x3ff;/* exponent bias */ |
|
|
#if PARI_BYTE_ORDER == LITTLE_ENDIAN_64 || PARI_BYTE_ORDER == BIG_ENDIAN_64 |
if (x==0.) return -1023; |
#else |
fi.f = x; |
error... unknown machine |
return ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid; |
#endif |
} |
|
|
GEN |
GEN |
dbltor(double x) |
dbltor(double x) |
{ |
{ |
GEN z = cgetr(3); |
GEN z; |
long e; |
long e; |
union { double f; ulong i; } fi; |
union { double f; ulong i; } fi; |
const int mant_len = 52; /* mantissa bits (excl. hidden bit) */ |
const int mant_len = 52; /* mantissa bits (excl. hidden bit) */ |
Line 2333 dbltor(double x) |
|
Line 2485 dbltor(double x) |
|
const int expo_len = 11; /* number of bits of exponent */ |
const int expo_len = 11; /* number of bits of exponent */ |
LOCAL_HIREMAINDER; |
LOCAL_HIREMAINDER; |
|
|
if (x==0) { z[1]=evalexpo(-308); z[2]=0; return z; } |
if (x==0.) return realzero_bit(-1023); |
fi.f = x; |
fi.f = x; z = cgetr(DEFAULTPREC); |
e = ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid; |
e = ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid; |
z[1] = evalexpo(e) | evalsigne(x<0? -1: 1); |
z[1] = evalexpo(e) | evalsigne(x<0? -1: 1); |
z[2] = (fi.i << expo_len) | HIGHBIT; |
z[2] = (fi.i << expo_len) | HIGHBIT; |
|
|
return fi.f; |
return fi.f; |
} |
} |
|
|
#else |
#else /* LONG_IS_64BIT */ |
|
|
#if PARI_BYTE_ORDER == LITTLE_ENDIAN |
#if PARI_DOUBLE_FORMAT == 1 |
# define INDEX0 1 |
# define INDEX0 1 |
# define INDEX1 0 |
# define INDEX1 0 |
#elif PARI_BYTE_ORDER == BIG_ENDIAN |
#elif PARI_DOUBLE_FORMAT == 0 |
# define INDEX0 0 |
# define INDEX0 0 |
# define INDEX1 1 |
# define INDEX1 1 |
#else |
|
error... unknown machine |
|
#endif |
#endif |
|
|
|
long |
|
expodb(double x) |
|
{ |
|
union { double f; ulong i[2]; } fi; |
|
const int mant_len = 52; /* mantissa bits (excl. hidden bit) */ |
|
const int exp_mid = 0x3ff;/* exponent bias */ |
|
const int shift = mant_len-32; |
|
|
|
if (x==0.) return -1023; |
|
fi.f = x; |
|
{ |
|
const ulong a = fi.i[INDEX0]; |
|
return ((a & (HIGHBIT-1)) >> shift) - exp_mid; |
|
} |
|
} |
|
|
GEN |
GEN |
dbltor(double x) |
dbltor(double x) |
{ |
{ |
Line 2385 dbltor(double x) |
|
Line 2551 dbltor(double x) |
|
union { double f; ulong i[2]; } fi; |
union { double f; ulong i[2]; } fi; |
const int mant_len = 52; /* mantissa bits (excl. hidden bit) */ |
const int mant_len = 52; /* mantissa bits (excl. hidden bit) */ |
const int exp_mid = 0x3ff;/* exponent bias */ |
const int exp_mid = 0x3ff;/* exponent bias */ |
const int shift = mant_len-32; |
|
const int expo_len = 11; /* number of bits of exponent */ |
const int expo_len = 11; /* number of bits of exponent */ |
|
const int shift = mant_len-32; |
|
|
if (x==0) { z=cgetr(3); z[1]=evalexpo(-308); z[2]=0; return z; } |
if (x==0.) return realzero_bit(-1023); |
fi.f = x; z=cgetr(4); |
fi.f = x; z=cgetr(DEFAULTPREC); |
{ |
{ |
const ulong a = fi.i[INDEX0]; |
const ulong a = fi.i[INDEX0]; |
const ulong b = fi.i[INDEX1]; |
const ulong b = fi.i[INDEX1]; |
|
|
union { double f; ulong i[2]; } fi; |
union { double f; ulong i[2]; } fi; |
const int mant_len = 52; /* mantissa bits (excl. hidden bit) */ |
const int mant_len = 52; /* mantissa bits (excl. hidden bit) */ |
const int exp_mid = 0x3ff;/* exponent bias */ |
const int exp_mid = 0x3ff;/* exponent bias */ |
const int shift = mant_len-32; |
|
const int expo_len = 11; /* number of bits of exponent */ |
const int expo_len = 11; /* number of bits of exponent */ |
|
const int shift = mant_len-32; |
|
|
if (typ(x)==t_INT && !s) return 0.0; |
if (typ(x)==t_INT && !s) return 0.0; |
if (typ(x)!=t_REAL) err(typeer,"rtodbl"); |
if (typ(x)!=t_REAL) err(typeer,"rtodbl"); |
|
|
fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len); |
fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len); |
return fi.f; |
return fi.f; |
} |
} |
#endif |
#endif /* LONG_IS_64BIT */ |
|
|
/* Old cgiv without reference count (which was not used anyway) |
/* Old cgiv without reference count (which was not used anyway) |
* Should be a macro. |
* Should be a macro. |
|
|
cgiv(GEN x) |
cgiv(GEN x) |
{ |
{ |
if (x == (GEN) avma) |
if (x == (GEN) avma) |
avma = (long) (x+lg(x)); |
avma = (gpmem_t) (x+lg(x)); |
} |
} |
|
|
/********************************************************************/ |
/********************************************************************/ |
|
|
* |
* |
* Note that these routines do not know and do not need to know about the |
* Note that these routines do not know and do not need to know about the |
* PARI stack. |
* PARI stack. |
* |
* |
* Added 2001Jan15: |
* Added 2001Jan15: |
* rgcduu() is a variant of xxgcduu() which does not have f (the effect is |
* rgcduu() is a variant of xxgcduu() which does not have f (the effect is |
* that of f=0), but instead has a ulong vmax parameter, for use in rational |
* that of f=0), but instead has a ulong vmax parameter, for use in rational |
Line 2747 rgcduu(ulong d, ulong d1, ulong vmax, |
|
Line 2913 rgcduu(ulong d, ulong d1, ulong vmax, |
|
* as many partial quotients of d/d1 as possible, and assigns to [u,u1;v,v1] |
* as many partial quotients of d/d1 as possible, and assigns to [u,u1;v,v1] |
* the product of all the [0, 1; 1, q_j] thus obtained, where the leftmost |
* the product of all the [0, 1; 1, q_j] thus obtained, where the leftmost |
* factor arises from the quotient of the first division step. |
* factor arises from the quotient of the first division step. |
* |
* |
* For use in rational reconstruction, input param vmax can be given a |
* For use in rational reconstruction, input param vmax can be given a |
* nonzero value. In this case, we will return early as soon as v1 > vmax |
* nonzero value. In this case, we will return early as soon as v1 > vmax |
* (i.e. it is guaranteed that v <= vmax). --2001Jan15 |
* (i.e. it is guaranteed that v <= vmax). --2001Jan15 |
Line 3282 lgcdii(ulong* d, ulong* d1, |
|
Line 3448 lgcdii(ulong* d, ulong* d1, |
|
*u = xu1; *u1 = xu; *v = xv1; *v1 = xv; |
*u = xu1; *u1 = xu; *v = xv1; *v1 = xv; |
break; |
break; |
} |
} |
|
|
res++; /* commit dd, xu1, xv1 */ |
res++; /* commit dd, xu1, xv1 */ |
dd = tmpd; xu1 = tmpu; xv1 = tmpv; |
dd = tmpd; xu1 = tmpu; xv1 = tmpv; |
#ifdef DEBUG_LEHMER |
#ifdef DEBUG_LEHMER |
Line 3303 lgcdii(ulong* d, ulong* d1, |
|
Line 3469 lgcdii(ulong* d, ulong* d1, |
|
*================================== |
*================================== |
* If a is invertible, return 1, and set res = a^{ -1 } |
* If a is invertible, return 1, and set res = a^{ -1 } |
* Otherwise, return 0, and set res = gcd(a,b) |
* Otherwise, return 0, and set res = gcd(a,b) |
* |
* |
* This is sufficiently different from bezout() to be implemented separately |
* This is sufficiently different from bezout() to be implemented separately |
* instead of having a bunch of extra conditionals in a single function body |
* instead of having a bunch of extra conditionals in a single function body |
* to meet both purposes. |
* to meet both purposes. |
|
|
invmod(GEN a, GEN b, GEN *res) |
invmod(GEN a, GEN b, GEN *res) |
{ |
{ |
GEN v,v1,d,d1,q,r; |
GEN v,v1,d,d1,q,r; |
long av,av1,lim; |
gpmem_t av, av1, lim; |
long s; |
long s; |
ulong g; |
ulong g; |
ulong xu,xu1,xv,xv1; /* Lehmer stage recurrence matrix */ |
ulong xu,xu1,xv,xv1; /* Lehmer stage recurrence matrix */ |
Line 3481 bezout(GEN a, GEN b, GEN *pu, GEN *pv) |
|
Line 3647 bezout(GEN a, GEN b, GEN *pu, GEN *pv) |
|
{ |
{ |
GEN t,u,u1,v,v1,d,d1,q,r; |
GEN t,u,u1,v,v1,d,d1,q,r; |
GEN *pt; |
GEN *pt; |
long av,av1,lim; |
gpmem_t av, av1, lim; |
long s; |
long s; |
int sa, sb; |
int sa, sb; |
ulong g; |
ulong g; |
Line 3802 cbezout(long a,long b,long *uu,long *vv) |
|
Line 3968 cbezout(long a,long b,long *uu,long *vv) |
|
AUTHOR = {Collins, George E. and Encarnaci{\'o}n, Mark J.}, |
AUTHOR = {Collins, George E. and Encarnaci{\'o}n, Mark J.}, |
TITLE = {Efficient rational number reconstruction}, |
TITLE = {Efficient rational number reconstruction}, |
JOURNAL = {J. Symbolic Comput.}, |
JOURNAL = {J. Symbolic Comput.}, |
FJOURNAL = {Journal of Symbolic Computation}, |
|
VOLUME = {20}, |
VOLUME = {20}, |
YEAR = {1995}, |
YEAR = {1995}, |
NUMBER = {3}, |
NUMBER = {3}, |
PAGES = {287--297}, |
PAGES = {287--297}, |
ISSN = {0747-7171}, |
|
MRCLASS = {11Y16 (68Q40)}, |
|
MRNUMBER = {97c:11116}, |
|
MRREVIEWER = {Maurice Mignotte}, |
|
} |
} |
* Preprint available from: |
* Preprint available from: |
* ftp://ftp.risc.uni-linz.ac.at/pub/techreports/1994/94-64.ps.gz |
* ftp://ftp.risc.uni-linz.ac.at/pub/techreports/1994/94-64.ps.gz |
*/ |
*/ |
|
|
/* #define DEBUG_RATLIFT */ |
/* #define DEBUG_RATLIFT */ |
|
|
ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bmax) |
ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bmax) |
{ |
{ |
GEN d,d1,v,v1,q,r; |
GEN d,d1,v,v1,q,r; |
ulong av = avma, av1, lim; |
gpmem_t av = avma, av1, lim; |
long lb,lr,lbb,lbr,s,s0; |
long lb,lr,lbb,lbr,s,s0; |
ulong vmax; |
ulong vmax; |
ulong xu,xu1,xv,xv1; /* Lehmer stage recurrence matrix */ |
ulong xu,xu1,xv,xv1; /* Lehmer stage recurrence matrix */ |
Line 3904 ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bm |
|
Line 4065 ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bm |
|
* Furthermore, at the start of the loop body we have in fact |
* Furthermore, at the start of the loop body we have in fact |
* (c') 0 <= v < v1 <= bmax, d > d1 > amax >= 0, |
* (c') 0 <= v < v1 <= bmax, d > d1 > amax >= 0, |
* (and are never done already). |
* (and are never done already). |
* |
* |
* Main loop is similar to those of invmod() and bezout(), except for |
* Main loop is similar to those of invmod() and bezout(), except for |
* having to determine appropriate vmax bounds, and checking termination |
* having to determine appropriate vmax bounds, and checking termination |
* conditions. The signe(d1) condition is only for paranoia |
* conditions. The signe(d1) condition is only for paranoia |
|
|
long l = lgefint(a); |
long l = lgefint(a); |
ulong x,y,z,k,m; |
ulong x,y,z,k,m; |
int L, s; |
int L, s; |
|
|
if (l <= 3) return l==2? 0: usqrtsafe((ulong)a[2]); |
if (l <= 3) return l==2? 0: usqrtsafe((ulong)a[2]); |
s = bfffo(a[2]); |
s = bfffo(a[2]); |
if (s > 1) |
if (s > 1) |
|
|
} |
} |
while (x < y); |
while (x < y); |
return y; |
return y; |
|
} |
|
|
|
/* target should point to a buffer of source_end - source + 1 ulongs. |
|
|
|
fills this buffer by bits of ulongs in source..source_end-1 shifted |
|
right sh units; the "most significant" sh bits of the result are |
|
set to be the least significant sh bits of prepend. |
|
|
|
The ordering of bits in this bitmap is the same as for t_INT. |
|
|
|
sh should not exceed BITS_IN_LONG. |
|
*/ |
|
void |
|
shift_r(ulong *target, ulong *source, ulong *source_end, ulong prepend, ulong sh) |
|
{ |
|
if (sh) { /* shift_words_r() works */ |
|
register ulong sh_complement = BITS_IN_LONG - sh; |
|
|
|
shift_words_r(target, source, source_end, prepend, sh, sh_complement); |
|
} else { |
|
int i; |
|
|
|
for (i=0; i < source_end - source; i++) |
|
target[i] = source[i]; |
|
} |
} |
} |