Annotation of OpenXM_contrib2/asir2000/engine/pari-mp.c, Revision 1.1
1.1 ! saito 1: /*
! 2: * $OpenXM: $
! 3: */
! 4: /* for f-itv.c */
! 5:
! 6: #include "genpari.h"
! 7: #include "itv-pari.h"
! 8:
! 9:
! 10: GEN
! 11: PariAddDown(GEN x, GEN y)
! 12: {
! 13: if(typ(x)==1) return (typ(y)==1) ? addii(x,y) : PariAddirDown(y,x);
! 14: return (typ(y)==1) ? PariAddirDown(y,x) : PariAddrrDown(x,y);
! 15: }
! 16:
! 17: GEN
! 18: PariAddirDown(GEN x, GEN y)
! 19: {
! 20: long l,e,ly,av,i,l1;
! 21: GEN z;
! 22:
! 23: if(!signe(x)) return rcopy(y);
! 24: if(!signe(y))
! 25: {
! 26: l=lgef(x)-(expo(y)>>TWOPOTBITS_IN_LONG);if((l<3)||(l>32767)) err(adder3);
! 27: z=cgetr(l);affir(x,z);return z;
! 28: }
! 29: else
! 30: {
! 31: e=expo(y)-expi(x);ly=lg(y);
! 32: if(e>0)
! 33: {
! 34: l=ly-(e>>TWOPOTBITS_IN_LONG);if(l<=2) return rcopy(y);
! 35: }
! 36: else
! 37: {
! 38: l=ly+((-e)>>TWOPOTBITS_IN_LONG)+1;if(l>32767) err(adder3);
! 39: }
! 40: av=avma;z=cgetr(l);affir(x,z);l1=av-avma;l=l1>>TWOPOTBYTES_IN_LONG;
! 41: z=PariAddrrDown(z,y);
! 42: for(i=lg(z)-1;i>=0;i--) z[i+l]=z[i];z+=l;avma+=l1;
! 43: }
! 44: return z;
! 45: }
! 46:
! 47: GEN
! 48: PariAddrrDown(GEN x, GEN y)
! 49: {
! 50: long sx=signe(x),sy=signe(y),lx=lg(x),ly=lg(y),lz,ex=expo(x),ey=expo(y),sz;
! 51: long av0=avma,e,l,i,d,m,flag,lp1,lp2,av,k,j,cex,f2;
! 52: GEN z,p1,p2;
! 53:
! 54: if(!sy)
! 55: {
! 56: if(!sx) {e=(ex>=ey)?ex:ey;z=cgetr(3);z[2]=0;z[1]=e+HIGHEXPOBIT;return z;}
! 57: e=ex-ey;
! 58: if(e<0) {z=cgetr(3);z[2]=0;z[1]=ey+HIGHEXPOBIT;return z;}
! 59: l=(e>>TWOPOTBITS_IN_LONG)+3;if(l>lx) l=lx;z=cgetr(l);
! 60: for(i=1;i<l;i++) z[i]=x[i];return z;
! 61: }
! 62: e=ey-ex;
! 63: if(!sx)
! 64: {
! 65: if(e<0) {z=cgetr(3);z[2]=0;z[1]=ex+HIGHEXPOBIT;return z;}
! 66: l=(e>>TWOPOTBITS_IN_LONG)+3;if(l>ly) l=ly;z=cgetr(l);
! 67: for(i=1;i<l;i++) z[i]=y[i];return z;
! 68: }
! 69: if(e)
! 70: {
! 71: if(e<0) {z=x;x=y;y=z;lz=lx;lx=ly;ly=lz;ey=ex;e= -e;sz=sx;sx=sy;sy=sz;}
! 72: d=(e>>TWOPOTBITS_IN_LONG);m=e&(BITS_IN_LONG-1);
! 73: if(d>=ly-2) return rcopy(y);
! 74: l=d+lx;
! 75: if(l>=ly)
! 76: {
! 77: flag=1;p1=cgetr(ly);lp1=ly;lp2=ly-d;
! 78: }
! 79: else
! 80: {
! 81: flag=0;p1=cgetr(l+1);lp2=lx+1;lp1=l+1;
! 82: }
! 83: av=avma;
! 84: if(m)
! 85: {
! 86: p2=cgetr(lp2);m=BITS_IN_LONG-m;
! 87: if(flag) {shiftl(x[lp2-1],m);k=hiremainder;}
! 88: else k=0;
! 89: for(i=lp2-1;i>=3;i--)
! 90: {
! 91: p2[i]=shiftl(x[i-1],m)+k;k=hiremainder;
! 92: }
! 93: p2[2]=k;
! 94: }
! 95: else p2=x;
! 96: }
! 97: else
! 98: {
! 99: l=(lx>ly)?ly:lx;p1=cgetr(l);av=avma;lp2=lp1=l;flag=2;p2=x;m=0;
! 100: }
! 101: if(sx==sy)
! 102: {
! 103: overflow=0;
! 104: if(m+flag) for(i=lp1-1,j=lp2-1;j>=2;i--,j--) p1[i]=addllx(p2[j],y[i]);
! 105: else
! 106: {
! 107: p1[lp1-1]=y[lp1-1];
! 108: for(i=lp1-2,j=lp2-2;j>=2;i--,j--) p1[i]=addllx(p2[j],y[i]);
! 109: }
! 110: if(overflow)
! 111: {
! 112: for(;(i>=2)&&(y[i]==(long)MAXULONG);i--) p1[i]=0;
! 113: if(i>=2) {cex=0;p1[i]=y[i]+1;while(i>=3) {i--;p1[i]=y[i];}}
! 114: else
! 115: {
! 116: cex=1;k=HIGHBIT;if(ey==(HIGHEXPOBIT-1)) err(adder4);
! 117: for(i=2;i<lp1;i++) {p1[i]=shiftlr(p1[i],1)+k;k=hiremainder;}
! 118:
! 119: if ( sx < 0 && hiremainder ) { /* |p1| + 1 */
! 120: overflow=1;
! 121: for(j=lg(p1);(overflow) && j>=2;j--) p1[j]=addllx(p1[j],0);
! 122: }
! 123:
! 124: }
! 125: }
! 126: else {cex=0;for(;i>=2;i--) p1[i]=y[i];}
! 127: p1[1]=evalsigne(sx)+ey+cex+HIGHEXPOBIT;
! 128: avma=av;return p1;
! 129: }
! 130: else
! 131: {
! 132: if(!e)
! 133: {
! 134: for(i=2;(i<l)&&(p2[i]==y[i]);i++);
! 135: if(i==l)
! 136: {
! 137: e=ex-((l-2)<<TWOPOTBITS_IN_LONG)+HIGHEXPOBIT;if(e<0) err(adder5);
! 138: if(e>EXPOBITS) err(adder4);
! 139: avma=av0;z=cgetr(3);z[2]=0;z[1]=e;return z;
! 140: }
! 141: else
! 142: {
! 143: f2=(((ulong)y[i])>((ulong)p2[i]))?1:0;
! 144: }
! 145: }
! 146: else f2=1;
! 147: if(f2)
! 148: {
! 149: overflow=0;
! 150: if(m+flag) for(i=lp1-1,j=lp2-1;j>=2;i--,j--) p1[i]=subllx(y[i],p2[j]);
! 151: else
! 152: {
! 153: p1[lp1-1]=y[lp1-1];
! 154: for(i=lp1-2,j=lp2-2;j>=2;i--,j--) p1[i]=subllx(y[i],p2[j]);
! 155: }
! 156: if(overflow)
! 157: {
! 158: for(;(i>=2)&&(!y[i]);i--) p1[i]=(long)MAXULONG;
! 159: p1[i]=y[i]-1;while(i>=3) {i--;p1[i]=y[i];}
! 160: }
! 161: else for(;i>=2;i--) p1[i]=y[i];
! 162: }
! 163: else
! 164: {
! 165: overflow=0;
! 166: if(m+flag) for(i=lp1-1;i>=2;i--) p1[i]=subllx(p2[i],y[i]);
! 167: else
! 168: {
! 169: p1[lp1-1]=subllx(0,y[lp1-1]);
! 170: for(i=lp1-2;i>=2;i--) p1[i]=subllx(p2[i],y[i]);
! 171: }
! 172: }
! 173: for(i=2;!p1[i];i++);j=i-2;avma=av+(j<<TWOPOTBYTES_IN_LONG);p1[j]=p1[0]-j;p1+=j;
! 174: m=bfffo(p1[2]);e=ey-(j<<TWOPOTBITS_IN_LONG)-m+HIGHEXPOBIT;
! 175: if(e<0) err(adder5);
! 176: p1[1]=f2 ? evalsigne(sy)+e : evalsigne(sx)+e;
! 177: if(m)
! 178: {
! 179: k=0;for(i=lp1-1-j;i>=2;i--) {p1[i]=shiftl(p1[i],m)+k;k=hiremainder;}
! 180: }
! 181: return p1;
! 182: }
! 183: }
! 184:
! 185: GEN
! 186: PariSubDown(GEN x, GEN y)
! 187: {
! 188: if(typ(x)==1) return (typ(y)==1) ? subii(x,y) : PariSubirDown(x,y);
! 189: return (typ(y)==1) ? PariSubriDown(x,y) : PariSubrrDown(x,y);
! 190: }
! 191:
! 192: GEN
! 193: PariSubrrDown(GEN x, GEN y)
! 194: {
! 195: long s=signe(y);
! 196: GEN z;
! 197:
! 198: if(x==y)
! 199: {
! 200: z=cgetr(3);z[2]=0;z[1]=HIGHEXPOBIT-(lg(x)<<TWOPOTBITS_IN_LONG);return z;
! 201: }
! 202: setsigne(y,-s);z=PariAddrrDown(x,y);setsigne(y,s);return z;
! 203: }
! 204:
! 205: GEN
! 206: PariSubirDown(GEN x, GEN y)
! 207: {
! 208: long s=signe(y);
! 209: GEN z;
! 210:
! 211: setsigne(y,-s);z=PariAddirDown(x,y);setsigne(y,s);return z;
! 212: }
! 213:
! 214: GEN
! 215: PariSubriDown(GEN x, GEN y)
! 216: {
! 217: long s=signe(y);
! 218: GEN z;
! 219:
! 220: setsigne(y,-s);z=PariAddirDown(y,x);setsigne(y,s);return z;
! 221: }
! 222:
! 223: GEN
! 224: PariMulDown(GEN x, GEN y)
! 225: {
! 226: if(typ(x)==1) return (typ(y)==1) ? mulii(x,y) : PariMulirDown(x,y);
! 227: return (typ(y)==1) ? PariMulirDown(y,x) : PariMulrrDown(x,y);
! 228: }
! 229:
! 230: GEN
! 231: PariMulrrDown(GEN x, GEN y)
! 232: {
! 233: long i,j,lx=lg(x),ly=lg(y),sx=signe(x),sy=signe(y),ex=expo(x),ey=expo(y);
! 234: long e,flag,garde,p1,p2,lz;
! 235: GEN z;
! 236:
! 237: e=ex+ey+HIGHEXPOBIT;if(e>=EXPOBITS) err(muler4);
! 238: if(e<0) err(muler5);
! 239: if((!sx)||(!sy)) {z=cgetr(3);z[2]=0;z[1]=e;return z;}
! 240: if(sy<0) sx= -sx;
! 241: if(lx>ly) {lz=ly;z=x;x=y;y=z;flag=1;} else {lz=lx;flag=(lx!=ly);}
! 242: z=cgetr(lz);z[1]=evalsigne(sx)+e;
! 243: if(flag) mulll(x[2],y[lz]);else hiremainder=0;
! 244: if(lz==3)
! 245: {
! 246: garde=flag ? addmul(x[2],y[2]) : mulll(x[2],y[2]);
! 247: if((long)hiremainder<0) {z[2]=hiremainder;z[1]++;}
! 248: else {z[2]=(garde<0)?(hiremainder<<1)+1:(hiremainder<<1);}
! 249: return z;
! 250: }
! 251: p1=x[lz-1];garde=hiremainder;
! 252: if(p1)
! 253: {
! 254: mulll(p1,y[3]);p2=addmul(p1,y[2]);
! 255: garde=addll(p2,garde);z[lz-1]=overflow+hiremainder;
! 256: }
! 257: else z[lz-1]=0;
! 258: for(j=lz-2;j>=3;j--)
! 259: {
! 260: p1=x[j];
! 261: if(p1)
! 262: {
! 263: mulll(p1,y[lz+2-j]);
! 264: p2=addmul(p1,y[lz+1-j]);
! 265: garde=addll(p2,garde);hiremainder+=overflow;
! 266: for(i=lz-j;i>=2;i--)
! 267: {
! 268: p2=addmul(p1,y[i]);
! 269: z[i+j-1]=addll(p2,z[i+j-1]);hiremainder+=overflow;
! 270: }
! 271: z[j]=hiremainder;
! 272: }
! 273: else z[j]=0;
! 274: }
! 275: p1=x[2];p2=mulll(p1,y[lz-1]);
! 276: garde=addll(p2,garde);hiremainder+=overflow;
! 277: for(i=lz-2;i>=2;i--)
! 278: {
! 279: p2=addmul(p1,y[i]);
! 280: z[i+1]=addll(p2,z[i+1]);hiremainder+=overflow;
! 281: }
! 282: z[2]=hiremainder;
! 283: if((long)hiremainder>0)
! 284: {
! 285: overflow=(garde<0)?1:0;
! 286: for(i=lz-1;i>=2;i--) {p1=z[i];z[i]=addllx(p1,p1);}
! 287: }
! 288: else z[1]++;
! 289: return z;
! 290: }
! 291:
! 292: GEN
! 293: PariMulirDown(GEN x, GEN y)
! 294: {
! 295: long sx=signe(x),sy,av,lz,ey,e,garde,p1,p2,i,j;
! 296: GEN z,temp;
! 297:
! 298: if(!sx) return gzero;
! 299: sy=signe(y);ey=expo(y);
! 300: if(!sy)
! 301: {
! 302: z=cgetr(3);z[2]=0;e=expi(x)+ey+HIGHEXPOBIT;if(e>EXPOBITS) err(muler6);
! 303: z[1]=e;return z;
! 304: }
! 305: lz=lg(y);if(sy<0) sx= -sx;
! 306: z=cgetr(lz);av=avma;
! 307: temp=cgetr(lz+1);affir(x,temp);x=y;y=temp;
! 308: e=expo(y)+ey+HIGHEXPOBIT;if(e>=EXPOBITS) err(muler4);
! 309: if(e<0) err(muler5);
! 310: z[1]=evalsigne(sx)+e;
! 311: mulll(x[2],y[lz]);
! 312: if(lz==3)
! 313: {
! 314: garde=addmul(x[2],y[2]);
! 315: if((long)hiremainder<0) {z[2]=hiremainder;z[1]++;}
! 316: else {z[2]=(garde<0)?(hiremainder<<1)+1:(hiremainder<<1);}
! 317: avma=av;return z;
! 318: }
! 319: garde=hiremainder;
! 320: p1=x[lz-1];mulll(p1,y[3]);p2=addmul(p1,y[2]);
! 321: garde=addll(p2,garde);z[lz-1]=overflow+hiremainder;
! 322: for(j=lz-2;j>=3;j--)
! 323: {
! 324: p1=x[j];mulll(p1,y[lz+2-j]);
! 325: p2=addmul(p1,y[lz+1-j]);
! 326: garde=addll(p2,garde);hiremainder+=overflow;
! 327: for(i=lz-j;i>=2;i--)
! 328: {
! 329: p2=addmul(p1,y[i]);
! 330: z[i+j-1]=addll(p2,z[i+j-1]);hiremainder+=overflow;
! 331: }
! 332: z[j]=hiremainder;
! 333: }
! 334: p1=x[2];p2=mulll(p1,y[lz-1]);
! 335: garde=addll(p2,garde);hiremainder+=overflow;
! 336: for(i=lz-2;i>=2;i--)
! 337: {
! 338: p2=addmul(p1,y[i]);
! 339: z[i+1]=addll(p2,z[i+1]);hiremainder+=overflow;
! 340: }
! 341: z[2]=hiremainder;
! 342: if((long)hiremainder>0)
! 343: {
! 344: overflow=(garde<0)?1:0;
! 345: for(i=lz-1;i>=2;i--) {p1=z[i];z[i]=addllx(p1,p1);}
! 346: }
! 347: else z[1]++;
! 348: avma=av;return z;
! 349: }
! 350:
! 351: #ifdef LONG_IS_32BIT
! 352: #define DIVCONVI 14
! 353: #endif
! 354:
! 355: #ifdef LONG_IS_64BIT
! 356: #define DIVCONVI 7
! 357: #endif
! 358:
! 359: GEN
! 360: PariDivDown(GEN x, GEN y)
! 361: {
! 362: if(typ(x)==1) return (typ(y)==1) ? divii(x,y) : PariDivirDown(x,y);
! 363: return (typ(y)==1) ? PariDivriDown(x,y) : PariDivrrDown(x,y);
! 364: }
! 365:
! 366: GEN
! 367: PariDivirDown(GEN x, GEN y)
! 368: {
! 369: GEN xr,z;
! 370: long av,ly;
! 371:
! 372: if(!signe(y)) err(diver5);
! 373: if(!signe(x)) return gzero;
! 374: ly=lg(y);z=cgetr(ly);av=avma;affir(x,xr=cgetr(ly+1));
! 375: xr=PariDivrrDown(xr,y);affrr(xr,z);avma=av;return z;
! 376: }
! 377:
! 378: GEN
! 379: PariDivriDown(GEN x, GEN y)
! 380: {
! 381: GEN yr,z;
! 382: long av,lx,ex,s=signe(y);
! 383:
! 384: if(!s) err(diver8);
! 385: if(!signe(x))
! 386: {
! 387: ex=expo(x)-expi(y)+HIGHEXPOBIT;
! 388: if(ex<0) err(diver12);
! 389: z=cgetr(3);z[1]=ex;z[2]=0;return z;
! 390: }
! 391: if((lg(y)==3)&&(y[2]>0)) return (s>0) ? divrs(x,y[2]) : divrs(x,-y[2]);
! 392: lx=lg(x);z=cgetr(lx);av=avma;affir(y,yr=cgetr(lx+1));
! 393: yr=PariDivrrDown(x,yr);affrr(yr,z);avma=av;return z;
! 394: }
! 395:
! 396: GEN
! 397: PariDivrrDown(GEN x, GEN y)
! 398: {
! 399: long sx=signe(x),sy=signe(y),lx,ly,lz,ex,ex1,i,z0;
! 400: ulong ldif,y0,y1,si,saux,qp,k,k3,k4,j;
! 401: GEN z;
! 402: if(!sy) err(diver9);
! 403: ex=expo(x)-expo(y)+HIGHEXPOBIT;
! 404: if(ex<=0) err(diver10);
! 405: if(ex>EXPOBITS) err(diver11);
! 406: if(!sx)
! 407: {
! 408: z=cgetr(3);z[1]=ex;z[2]=0;return z;
! 409: }
! 410: lx=lg(x);ly=lg(y);lz=(lx<=ly)?lx:ly;
! 411: z=cgetr(lz);if(sy<0) sx= -sx;
! 412: ex1=evalsigne(sx)+ex;
! 413: if(ly==3)
! 414: {
! 415: i=x[2];si=(lx>3)?x[3]:0;
! 416: if((ulong)i<(ulong)y[2])
! 417: {
! 418: hiremainder=i;z[2]=divll(si,y[2]);
! 419: z[1]=ex1-1;return z;
! 420: }
! 421: else
! 422: {
! 423: hiremainder=((ulong)i)>>1;
! 424: z[2]=(i&1)?divll((((ulong)si)>>1)|(HIGHBIT),y[2]):divll(((ulong)si)>>1,y[2]);
! 425: z[1]=ex1;return z;
! 426: }
! 427: }
! 428: z0= *z;*z=0;
! 429: for(i=2;i<=lz-1;i++) z[i-1]=x[i];
! 430: z[lz-1]=(lx>lz) ? x[lz] : 0;
! 431: ldif=ly-lz;if(!ldif) {y0=y[lz];y[lz]=0;}
! 432: if(ldif<=1) {y1=y[lz+1];y[lz+1]=0;}
! 433: si=y[2];saux=y[3];
! 434: for(i=0;i<lz-1;i++)
! 435: {
! 436: if(z[i]!=si)
! 437: {
! 438: if((ulong)z[i]>si)
! 439: {
! 440: overflow=0;
! 441: for(j=lz-i+1;j>=2;j--) z[i+j-2]=subllx(z[i+j-2],y[j]);
! 442: {z[i-1]++;for(j=i-1;j&&(!z[j]);j--) z[j-1]++;}
! 443: }
! 444: hiremainder=z[i];qp=divll(z[i+1],si);
! 445: overflow=0;k=hiremainder;
! 446: }
! 447: else
! 448: {
! 449: qp=(long)MAXULONG;k=addll(si,z[i+1]);
! 450: }
! 451: if(!overflow)
! 452: {
! 453: k3=subll(mulll(qp,saux),z[i+2]);k4=subllx(hiremainder,k);
! 454: while((!overflow)&&k4) {qp--;k3=subll(k3,saux);k4=subllx(k4,si);}
! 455: }
! 456: mulll(qp,y[lz+1-i]);
! 457: for(j=lz-i;j>=2;j--)
! 458: {
! 459: z[i+j-1]=subll(z[i+j-1],addmul(qp,y[j]));hiremainder+=overflow;
! 460: }
! 461: if((ulong)z[i]!=(ulong)hiremainder)
! 462: {
! 463: if((ulong)z[i]<(ulong)hiremainder)
! 464: {
! 465: overflow=0;qp--;
! 466: for(j=lz-i;j>=2;j--) z[i+j-1]=addllx(z[i+j-1],y[j]);
! 467: }
! 468: else
! 469: {
! 470: z[i]-=hiremainder;
! 471: while(z[i])
! 472: {
! 473: overflow=0;qp++;
! 474: if(!qp) {z[i-1]++;for(j=i-1;j&&(!z[j]);j--) z[j-1]++;}
! 475: for(j=lz-i;j>=2;j--) z[i+j-1]=subllx(z[i+j-1],y[j]);
! 476: z[i]-=overflow;
! 477: }
! 478: }
! 479: }
! 480: z[i]=qp;
! 481: }
! 482: if(!ldif) y[lz]=y0;if(ldif<=1) y[lz+1]=y1;
! 483: for(j=lz-1;j>=2;j--) z[j]=z[j-1];
! 484: if(*z)
! 485: {
! 486: k=HIGHBIT;
! 487: for(j=2;j<lz;j++) {z[j]=shiftlr(z[j],1)+k;k=hiremainder;}
! 488: }
! 489: else ex1--;
! 490: z[1]=ex1;*z=z0;return z;
! 491: }
! 492:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to pari-mp.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine