This is mpi.c in view mode; [Download] [Up]
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~ ~*/
/*~ OPERATIONS DE BASE (NOYAU) ~*/
/*~ ~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
/* This file was modified by W. Schelter to be suitable for optimization
and inlining of assembler for maximum speed
*/
#include "config.h"
#include "genpari.h"
#include "arith.h"
GEN mulsi(x,y)
plong x;
GEN y;
{ TEMPVARS
plong s=signe(y),ly=lgef(y),i;
GEN z,zp,yp; ulong hiremainder;
if((!x)||(!s)) return gzero;
if(x<0) {s= -s;
x= -x;
if (x < 0) /* -2^31 */
{return mulii(stoi(1<<31),y);}
}
z=cgeti(ly+1);
hiremainder=0;
MP_START_LOW(yp,y,ly); MP_START_LOW(zp,z,ly+1);
i = MP_COUNT_LG(ly);
WHILE_COUNT(--i) { MP_NEXT_UP(zp) = addmul(x,MP_NEXT_UP(yp));}
if(hiremainder) {MP_NEXT_UP(zp)=hiremainder;
setlgef(z,ly+1);}
else {avma+=4;z[1]=z[0]-1;z++;setlgef(z,ly);}
setsigne(z,s);return z;
}
int expi(x)
GEN x;
{
plong lx=x[1]&0xffff;
return lx==2 ? -8388608 : ((lx-2)<<5)-bfffo(x[2])-1;
}
GEN addsi(x,y)
plong x;
GEN y;
{
plong sx,sy,ly,p,i;
ulong overflow;
GEN z; TEMPVARS
if(!x) return icopy(y);
sy=signe(y);if(!sy) return stoi(x);
if(x<0) {sx= -1;
x= -x;
if (x < 0) /* x=-2^31 */
return addii(MOST_NEGS,y);
} else sx=1;
ly=lgef(y);
if(sx==sy)
{
p=addll(x,y[ly-1]);
if(overflow)
{
z=cgeti(ly+1);z[ly]=p;
for(i=ly-1;(i>2)&&(y[i-1]==0xffffffff);i--) z[i]=0;
if(i>2)
{
z[i]=y[i-1]+1;i--;while(i>=3) {z[i]=y[i-1];i--;}
z[2]=z[1]=z[0]-1;z++;avma+=4;
}
else {z[2]=1;z[1]=z[0];}
}
else
{
z=cgeti(ly);z[ly-1]=p;for(i=1;i<ly-1;i++) z[i]=y[i];
}
setsigne(z,sx);
}
else
{
if(ly==3)
{
if((ulong)y[2]>(ulong)x)
{
z=cgeti(3);z[1]=(sy<<24)+3;z[2]=y[2]-x;return z;
}
if(y[2]==x) return gzero;
z=cgeti(3);z[1]=((-sy)<<24)+3;z[2]=x-y[2];return z;
}
p=subll(y[ly-1],x);
if(overflow)
{
z=cgeti(ly);z[ly-1]=p;
for(i=ly-2;!(y[i]);i--) z[i]=0xffffffff;
z[i]=y[i]-1;
if((i>2)||z[i]) {i--;for(;i>=1;i--) z[i]=y[i];}
else
{
z[2]=z[1]=z[0]-1;z++;avma+=4;setsigne(z,sy);
}
}
else
{
z=cgeti(ly);z[ly-1]=p;for(i=1;i<ly-1;i++) z[i]=y[i];
}
}
return z;
}
GEN addii(x,y)
GEN x,y;
{
plong sx,sy,sz,lx,ly,i,j,p;
GEN z,xp,yp,zp,xpp,xhigh; TEMPVARS
ulong overflow;
lx = lgef(x);
ly = lgef(y);
if (lx< ly) {z=x;x=y;y=z; sx=lx; lx=ly; ly=sx;}
/* lx>=ly so sx==0 ==> sy==0 */
if (0 == (sy=signe(y))) return icopy(x);
sx = signe(x);
if(sx==sy)
{
z=cgeti(lx+1);overflow=0;
MP_START_LOW(zp,z,lx+1);
MP_START_LOW(xp,x,lx);
MP_START_LOW(yp,y,ly);
#ifdef QUICK_LOOP
i = ly - 2;
QUICK_LOOP(i,ADDXCC);
#else
i = MP_COUNT_LG(ly);
WHILE_COUNT(--i)
{ADDLLX(MP_NEXT_UP(xp),MP_NEXT_UP(yp),MP_NEXT_UP(zp));}
#endif
if(overflow)
{ GEN xhigh = &MP_HIGH(x,lx);
again:
{ GEN xpp = &MP_NEXT_UP(xp);
if (xpp >= xhigh)
{ if (*xpp == 0xffffffff)
{ MP_NEXT_UP(zp)=0;
goto again;}
else
{ MP_NEXT_UP(zp) = *xpp + 1;
while ((xpp = &MP_NEXT_UP(xp)) >= xhigh)
{ MP_NEXT_UP(zp) = *xpp ;}
z[1]=z[0]-1;z[2]=x[1];z++;avma+=4;}}
else {z[2]=1;z[1]=x[1]+1;}
}}
else
{ j = COUNT(lx - ly);
WHILE_COUNT( --j)
{ MP_NEXT_UP(zp) = MP_NEXT_UP(xp);}
z[1]=z[0]-1;z[2]=x[1];z++;avma+=4;
}
}
else
{
if(lx==ly)
/* we have to compare x and y */
{ j = MP_COUNT_LG(lx);
MP_START_HIGH(xp,x,lx);
MP_START_HIGH(yp,y,lx);
WHILE_COUNT(--j)
{ ulong tx = MP_NEXT_DOWN(xp);
ulong ty = MP_NEXT_DOWN(yp);
if ( ty > tx)
{z=x;x=y;y=z;sz=sx;sx=sy;sy=sz;
goto DIFFER;}
else
if ( tx > ty)
{goto DIFFER;}}
SAME: return gzero;
DIFFER:;
}
z=cgeti(lx);overflow=0;
MP_START_LOW(xp,x,lx);MP_START_LOW(yp,y,ly);MP_START_LOW(zp,z,lx);
i = MP_COUNT_LG(ly);
#ifdef QUICK_LOOP
i = ly - 2;
QUICK_LOOP(i,SUBXCC);
#else
i = MP_COUNT_LG(ly);
WHILE_COUNT(--i)
{SUBLLX(MP_NEXT_UP(xp),MP_NEXT_UP(yp),MP_NEXT_UP(zp));}
#endif
if(overflow)
{ ulong tx ;
while((tx=MP_NEXT_UP(xp)) == 0)
MP_NEXT_UP(zp) = 0xffffffff;
if (xp >= (xhigh = &MP_HIGH(x,lx)))
{ MP_NEXT_UP(zp) = tx -1;
while ((xpp = &MP_NEXT_UP(xp)) >= xhigh)
{ MP_NEXT_UP(zp) = *xpp;}}
}
else
{ i = COUNT(lx - ly);
WHILE_COUNT(--i)
MP_NEXT_UP(zp) = MP_NEXT_UP(xp);
}
if(z[2]) z[1]=x[1];
else
{ zp = &z[3];
while (*zp ==0){zp++;} /* x was != y by above */
zp -= 2;
i = zp - z;
zp[1] = (zp[0] = z[0]-i);
z = zp;
setsigne(z,sx);
avma+=(i<<2);
}
}
return z;
}
GEN mulss(x,y)
plong x,y;
{
plong s,p1;
GEN z; ulong hiremainder;
if((!x)||(!y)) return gzero;
s=1;
if(x<0)
{s= -1;
x= -x;
if (x<0)
return mulsi(y,stoi(x));
}
if(y<0) {s= -s;
y= -y;
if(y<0)
return mulsi((s > 0 ? x : -x),ABS_MOST_NEGS);
}
p1=mulll(x,y);
if(hiremainder) {z=cgeti(4);z[2]=hiremainder;z[3]=p1;}
else {z=cgeti(3);z[2]=p1;}
z[1]=z[0];setsigne(z,s);return z;
}
GEN mulii(x,y)
GEN x,y;
{
plong i,j,lx=lgef(x),ly=lgef(y),sx,sy,lz,p1,p2;
GEN z; TEMPVARS
GEN zz,yy,zp,xx;
GEN ylow;
ulong hiremainder;
ulong overflow;
sx=signe(x);if(!sx) return gzero;
sy=signe(y);if(!sy) return gzero;
if(sy<0) sx= -sx;
if(lx>ly) {z=x;x=y;y=z;lz=lx;lx=ly;ly=lz;}
lz=lx+ly-2;if(lz>=0x10000) err(muler1);
z=cgeti(lz);z[1]=z[0];setsigne(z,sx);
MP_START_LOW(xx,x,lx);
p1 = MP_NEXT_UP(xx);
hiremainder=0;
i = COUNT(ly-2);
MP_START_LOW(yy,y,ly);
MP_START_LOW(zz,z,lz);
WHILE_COUNT (--i)
{ MP_NEXT_UP(zz) = addmul(p1,MP_NEXT_UP(yy));}
MP_NEXT_UP(zz) = hiremainder;
/* restart zz one above bottom */
MP_START_LOW(zz,z,lz);
MP_START_LOW(ylow,y,ly);
ly = COUNT(ly - MP_CODE_WORDS);
lx -= MP_CODE_WORDS;
while (--lx > 0) /* one less iteration first term of x, already used */
{ plong tem;
register plong p11;
p11 = MP_NEXT_UP(xx);
i = ly;
yy = ylow;
zp = &MP_NEXT_UP(zz); /* *zp = second from low word of z first time through */
tem = 0;
/* ZerO is just a 68k kludge to getit to keep 0 in a reg during this loop*/
#undef ZERO
#define ZERO ZerO
{ int ZerO = 0;
WHILE_COUNT(--i)
{ p2 = MP_NEXT_UP(yy);
p2 = mulul(p2,p11,hiremainder);
MP_NEXT_UP(zp);
p2 = add_carry(p2,*zp,hiremainder);
p2 = add_carry(p2,tem,hiremainder);
*zp = p2;
tem = hiremainder;
}
}
MP_NEXT_UP(zp) = hiremainder;
#undef ZERO
#define ZERO 0
}
if(!MP_HIGH(z,lz))
{ /* shift header one along decreasing lg and lgef */
z[2]=z[1]-1;z[1]=z[0]-1;z++;avma+=4;
}
return z;
}
GEN confrac(x)
GEN x;
{
plong lx=lg(x),ex= -expo(x)-1,ex1,av=avma,ly,ey;
plong lr,nbdec,k,i,j; ulong hiremainder;
GEN y,res; TEMPVARS
ey=((lx-2)<<5)+ex;ly=(ey+63)>>5;y=cgeti(ly);ex1=ex>>5; /* 95 dans mp.s faux? */
for(i=0;i<ex1;i++) y[i]=0;
ex&=31;
if(!ex) for(j=2;j<lx;j++) y[i++]=x[j];
else
{
k=0;
for(j=2;j<lx;j++) {y[i++]=shiftlr(x[j],ex)+k;k=hiremainder;}
y[ly-2]=k;
}
y[ly-1]=0;
nbdec=ey*0.30103+1;lr=(nbdec+17)/9;res=cgeti(lr);
*res=nbdec;
for(j=1;j<lr;j++)
{
hiremainder=0;
for(i=ly-1;i>=0;i--) y[i]=addmul(y[i],1000000000);
res[j]=hiremainder;
}
avma=av;return res;
}
/* x/y : uses hiremainder for return */
GEN divss(x,y)
plong x,y;
{
plong p1;
if(!y) err(diver1);
if (x == (1<<31))
return divis(stoi(x),y);
hiremainder=0;p1=divll((ulong)abs(x),(ulong)abs(y));
if(y<0) {hiremainder= -((plong)hiremainder);p1= -p1;}
if(x<0) p1= -p1;
return stoi(p1);
}
GEN modss(x,y)
plong x,y;
{
plong y1; ulong hiremainder;
if(!y) err(moder1);
if (x == (1<<31))
return modis(stoi(x),y);
hiremainder=0;divll(abs(x),y1=abs(y));
if(!hiremainder) return gzero;
return (((plong)hiremainder)<0) ? stoi(y1-hiremainder) : stoi(hiremainder);
}
GEN resss(x,y)
plong x,y;
{ ulong hiremainder;
if(!y) err(reser1);
hiremainder=0;divll(abs(x),abs(y));
return (y<0) ? stoi(-((plong)hiremainder)) : stoi(hiremainder);
}
/* uses hiremainder for return */
GEN divsi(x,y)
plong x;
GEN y;
{
plong s=signe(y),ly=lgef(y),p1;
if(!s) err(diver2);
if((!x)||(ly>3)||(y[2]<0)) {hiremainder=x;return gzero;}
if (x== 1<<31)
return divii(stoi(x),y);
hiremainder=0;p1=divll(abs(x),y[2]);
if(signe(y)<0) {hiremainder= -((plong)hiremainder);p1= -p1;}
if(x<0) p1= -p1;
return stoi(p1);
}
/* this uses the GLOBAL hiremainder to return its remainder
We cannot make it a local.
*/
GEN divis(y,x)
plong x;
GEN y;
{ ulong hi;
plong s=signe(y),ly=lgef(y),i,d;
GEN z;
if(!x) err(diver4);
if(!s) {hiremainder=0;return gzero;}
if(x<0) {s= -s;x= -x;
if (x < 0)
return divii(y,stoi(x));
}
if((ulong)x>(ulong)y[2])
{
if(ly==3) {hiremainder=itos(y);return gzero;}
else {z=cgeti(ly-1);d=1;hi=y[2];}
}
else {z=cgeti(ly);d=0;hi=0;}
for(i=d+2;i<ly;i++) z[i-d]=divul(y[i],x,hi);
z[1]=z[0];setsigne(z,s);
hiremainder= (s < 0 ? -((plong)hi) : hi) ;
return z;
}
/*
# entree : x pointe sur i2 de type I (dividende)
# y pointe sur i1 de type I (diviseur)
# z contient un pointeur sur le reste si l'on
# veut a la fois q et r, 0 si l'on ne veut que le
# quotient, -1 si l'on ne veut que le reste
# sortie : resultat est q si celui-ci est attendu, et sinon r
# r. z pointe sur r si q et r sont attendus
# (toutes les zones sont creees)
# remarque : il s'agit de la 'fausse division' ; le reste est
# du signe du dividende
*/
GEN dvmdii(x,y,z)
GEN x,y,*z;
{
plong av,lx,ly,lz,i,j,dec,sh,k,k1,sx=signe(x),sy=signe(y);
plong saux,k3,k4,av1,flk4;
ulong si,qp;
GEN p1,p2,p3,p4,xp,p1p,p2p,pp; TEMPVARS
ulong hiremainder,overflow;
if(!sy) err(dvmer1);
if(!sx)
{
if(((plong)z==0xffffffff)||((plong)z==0)) return gzero;
*z=gzero;return gzero;
}
lx=lgef(x);ly=lgef(y);lz=lx-ly;
if(lz<0)
{
if((plong)z==0xffffffff) return icopy(x);
if(z==0) return gzero;
*z=icopy(x);return gzero;
}
av=avma;if(sx<0) sy= -sy;
if(ly==3)
{ int lgp1;
si=y[2];
MP_START_HIGH(xp,x,lx);
if (si > (ulong)MP_HIGH(x,lx))
{ lgp1=lx-1; hiremainder= MP_NEXT_DOWN(xp);}
else { lgp1=lx; hiremainder=0;}
p1 = cgeti(lgp1); i = MP_COUNT_LG(lgp1);
MP_START_HIGH(p1p,p1,lgp1);
WHILE_COUNT(--i) { MP_NEXT_DOWN(p1p) = divll(MP_NEXT_DOWN(xp),si);}
if((plong)z==0xffffffff)
{
avma=av;if(!hiremainder) return gzero;
p2=cgeti(3);p2[1]=(sx<<24)+3;p2[2]=hiremainder;return p2;
}
if(lgp1!= 2) {p1[1]=p1[0];setsigne(p1,sy);} else {avma=av;p1=gzero;}
if(z==0) return p1;
if(!hiremainder) *z=gzero;
else {p2=cgeti(3);p2[1]=(sx<<24)+3;p2[2]=hiremainder;*z=p2;}
return p1;
}
else
{
p1=cgeti(lx);
sh=bfffo(y[2]);
if(sh)
{ GEN p2p,yp;
MP_START_HIGH(yp,y,ly);
p2=cgeti(ly);
k=shiftl(MP_NEXT_DOWN(yp),sh);
MP_START_HIGH(p2p,p2,ly);
i = MP_COUNT_LG(ly-1);
WHILE_COUNT(--i)
{
k1=shiftl(MP_NEXT_DOWN(yp),sh);
MP_NEXT_DOWN(p2p) = k + hiremainder;
k = k1;
}
MP_NEXT_DOWN(p2p) = k ; k=0;
MP_START_HIGH(xp,x,lx);
MP_START_HIGH(p1p,p1,lx);
MP_NEXT_UP(p1p) ; /* yes go out of range !! */
i = MP_COUNT_LG(lx);
WHILE_COUNT (--i)
{ k1 = shiftl(MP_NEXT_DOWN(xp),sh);
MP_NEXT_DOWN(p1p) = k + hiremainder; k = k1;
}
MP_NEXT_DOWN(p1p) = k;
}
else {
MP_START_HIGH(xp,x,lx);
MP_START_HIGH(p1p,p1,lx);
MP_NEXT_UP(p1p) ; /* yes go out of range !! */
MP_NEXT_DOWN(p1p) = 0;
j = MP_COUNT_LG(lx);
WHILE_COUNT (-- j)
{ MP_NEXT_DOWN(p1p) = MP_NEXT_DOWN(xp);}
p2 = y;}
si=p2[2];saux=p2[3];
MP_START_HIGH(p1p,p1,lx); MP_NEXT_UP(p1p) ; /* out of bound */
i = COUNT(lz+1);
WHILE_COUNT(--i)
{ GEN pp;
if(MP_NEXT_DOWN(p1p)==si)
{ /* Using fact that next_down does post increment */
qp=0xffffffff;k=addll(si,*p1p);
}
else
{
hiremainder=p1p[-1];qp=divll(*p1p,si);
overflow=0;k=hiremainder;
}
if(!overflow)
{
/* k1=mulll(qp,saux);k3=subll(k1,p1p[1]);k+=overflow;
flk4=((ulong)hiremainder>(ulong)k);k4=subll(hiremainder,k);
while(flk4) {qp--;k3=subll(k3,saux);
k4-=overflow;flk4=((ulong)k4>(ulong)si);
k4=subll(k4,si);}
*/
k1=mulll(qp,saux);k3=subll(k1,p1p[1]);
k4=subllx(hiremainder,k);
while((!overflow)&&k4) {qp--;k3=subll(k3,saux);k4=subllx(k4,si);}
}
hiremainder=0;
j = MP_COUNT_LG(ly);
MP_START_LOW(pp,p1p,ly-2);
MP_START_LOW(p2p,p2,ly);
WHILE_COUNT(--j)
{ GEN ppp;
k1=addmul(qp,MP_NEXT_UP(p2p));
ppp = &MP_NEXT_UP(pp);
*ppp =subll(*ppp,k1);hiremainder+=overflow;
}
if((ulong)p1p[-1]<(ulong)hiremainder)
{
overflow=0;qp--;
j = MP_COUNT_LG(ly);
MP_START_LOW(pp,p1p,ly-2);
MP_START_LOW(p2p,p2,ly);
WHILE_COUNT(--j){ GEN ppp = &MP_NEXT_UP(pp);
ADDLLX(*ppp,MP_NEXT_UP(p2p),*ppp);}
}
p1p[-1] = qp;
}
av1=avma;
if((plong)z!=0xffffffff)
{ulong lgp3 = lz + 2;
MP_START_LOW(p1p,p1,lgp3);
if (p1[1]) {lgp3++;}
else if (lz==0) sy=0;
p3 = cgeti(lgp3);
MP_START_LOW(pp,p3,lgp3);
j = MP_COUNT_LG(lgp3);
WHILE_COUNT(--j)
{MP_NEXT_UP(pp) = MP_NEXT_UP(p1p) ;}
if(lgp3<3) {p3[1]=2;} else {p3[1]=p3[0];setsigne(p3,sy);}
}
if(z!=0)
{
for(j=lz+2;(j<lx)&&(!p1[j]);j++);
if(j==lx) p4=icopy(gzero);
else
{
p4=cgeti(lx-j+2);p4[1]=p4[0];
if(!sh) for(i=2;j<lx;j++,i++) p4[i]=p1[j];
else
{
hiremainder=0;k1=shiftlr(p1[j++],sh);k=hiremainder;
if(k1) {p4[2]=k1;dec=1;} else
{p4[1]=p4[0]-1;p4++;avma+=4;p4[1]=p4[0];dec=0;}
for(i=2+dec;j<lx;j++,i++)
{
p4[i]=shiftlr(p1[j],sh)+k;k=hiremainder;
}
}
setsigne(p4,sx);
}
}
if((plong)z==0xffffffff) return gerepile(av,av1,p4);
if((plong)z==0) return gerepile(av,av1,p3);
dec=lpile(av,av1,0)>>2;*z=p4+dec;return p3+dec;
}
}
/* machines which provide an inline version of mulul need
to provide a function for calls where that inlining can't take place
*/
#ifdef NEED_MULUL3
ulong
mulul3(a,b,c)
ulong a,b,*c;
{ return mulul(a,b,*c);}
#endif
#ifdef NEED_DIVUL3
ulong
divul3(a,b,c)
ulong a,b,*c;
{ return divul(a,b,*c);}
#endif
/*
;;- Local variables:
;;- version-control:t
;;- End:
*/
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.