Annotation of OpenXM_contrib/pari/src/kernel/none/mp.c, Revision 1.1
1.1 ! maekawa 1: /***********************************************************************/
! 2: /** **/
! 3: /** MULTIPRECISION KERNEL **/
! 4: /** **/
! 5: /***********************************************************************/
! 6: /* $Id: mp.c,v 1.1.1.1 1999/09/16 13:47:57 karim Exp $ */
! 7: /* most of the routines in this file are commented out in 68k */
! 8: /* version (#ifdef __M68K__) since they are defined in mp.s */
! 9: #include "pari.h"
! 10:
! 11: /* z2 := z1[imin..imax].f shifted left sh bits (feeding f from the right) */
! 12: #define shift_left2(z2,z1,imin,imax,f, sh,m) {\
! 13: register ulong _i,_l,_k = ((ulong)f)>>m;\
! 14: for (_i=imax; _i>imin; _i--) {\
! 15: _l = z1[_i];\
! 16: z2[_i] = (_l<<(ulong)sh) | _k;\
! 17: _k = _l>>(ulong)m;\
! 18: }\
! 19: z2[imin]=(z1[imin]<<(ulong)sh) | _k;\
! 20: }
! 21: #define shift_left(z2,z1,imin,imax,f, sh) {\
! 22: register const ulong _m = BITS_IN_LONG - sh;\
! 23: shift_left2((z2),(z1),(imin),(imax),(f),(sh),(_m));\
! 24: }
! 25:
! 26: /* z2 := f.z1[imin..imax-1] shifted right sh bits (feeding f from the left) */
! 27: #define shift_right2(z2,z1,imin,imax,f, sh,m) {\
! 28: register ulong _i,_k,_l = z1[2];\
! 29: z2[imin] = (_l>>(ulong)sh) | ((ulong)f<<(ulong)m);\
! 30: for (_i=imin+1; _i<imax; _i++) {\
! 31: _k = _l<<(ulong)m; _l = z1[_i];\
! 32: z2[_i] = (_l>>(ulong)sh) | _k;\
! 33: }\
! 34: }
! 35: #define shift_right(z2,z1,imin,imax,f, sh) {\
! 36: register const ulong _m = BITS_IN_LONG - sh;\
! 37: shift_right2((z2),(z1),(imin),(imax),(f),(sh),(_m));\
! 38: }
! 39:
! 40: #ifdef INLINE
! 41: INLINE
! 42: #endif
! 43: GEN
! 44: addsispec(long s, GEN x, long nx)
! 45: {
! 46: GEN xd, zd = (GEN)avma;
! 47: long lz;
! 48: LOCAL_OVERFLOW;
! 49:
! 50: lz = nx+3; (void)new_chunk(lz);
! 51: xd = x + nx;
! 52: *--zd = addll(*--xd, s);
! 53: if (overflow)
! 54: for(;;)
! 55: {
! 56: if (xd == x) { *--zd = 1; break; } /* enlarge z */
! 57: *--zd = ((ulong)*--xd) + 1;
! 58: if (*zd) { lz--; break; }
! 59: }
! 60: else lz--;
! 61: while (xd > x) *--zd = *--xd;
! 62: *--zd = evalsigne(1) | evallgefint(lz);
! 63: *--zd = evaltyp(t_INT) | evallg(lz);
! 64: avma=(long)zd; return zd;
! 65: }
! 66:
! 67: #define swapspec(x,y, nx,ny) {long _a=nx;GEN _z=x; nx=ny; ny=_a; x=y; y=_z;}
! 68:
! 69: #ifdef INLINE
! 70: INLINE
! 71: #endif
! 72: GEN
! 73: addiispec(GEN x, GEN y, long nx, long ny)
! 74: {
! 75: GEN xd,yd,zd;
! 76: long lz;
! 77: LOCAL_OVERFLOW;
! 78:
! 79: if (nx < ny) swapspec(x,y, nx,ny);
! 80: if (ny == 1) return addsispec(*y,x,nx);
! 81: zd = (GEN)avma;
! 82: lz = nx+3; (void)new_chunk(lz);
! 83: xd = x + nx;
! 84: yd = y + ny;
! 85: *--zd = addll(*--xd, *--yd);
! 86: while (yd > y) *--zd = addllx(*--xd, *--yd);
! 87: if (overflow)
! 88: for(;;)
! 89: {
! 90: if (xd == x) { *--zd = 1; break; } /* enlarge z */
! 91: *--zd = ((ulong)*--xd) + 1;
! 92: if (*zd) { lz--; break; }
! 93: }
! 94: else lz--;
! 95: while (xd > x) *--zd = *--xd;
! 96: *--zd = evalsigne(1) | evallgefint(lz);
! 97: *--zd = evaltyp(t_INT) | evallg(lz);
! 98: avma=(long)zd; return zd;
! 99: }
! 100:
! 101: /* assume x >= y */
! 102: #ifdef INLINE
! 103: INLINE
! 104: #endif
! 105: GEN
! 106: subisspec(GEN x, long s, long nx)
! 107: {
! 108: GEN xd, zd = (GEN)avma;
! 109: long lz;
! 110: LOCAL_OVERFLOW;
! 111:
! 112: lz = nx+2; (void)new_chunk(lz);
! 113: xd = x + nx;
! 114: *--zd = subll(*--xd, s);
! 115: if (overflow)
! 116: for(;;)
! 117: {
! 118: *--zd = ((ulong)*--xd) - 1;
! 119: if (*xd) break;
! 120: }
! 121: if (xd == x)
! 122: while (*zd == 0) { zd++; lz--; } /* shorten z */
! 123: else
! 124: do *--zd = *--xd; while (xd > x);
! 125: *--zd = evalsigne(1) | evallgefint(lz);
! 126: *--zd = evaltyp(t_INT) | evallg(lz);
! 127: avma=(long)zd; return zd;
! 128: }
! 129:
! 130: /* assume x >= y */
! 131: #ifdef INLINE
! 132: INLINE
! 133: #endif
! 134: GEN
! 135: subiispec(GEN x, GEN y, long nx, long ny)
! 136: {
! 137: GEN xd,yd,zd;
! 138: long lz;
! 139: LOCAL_OVERFLOW;
! 140:
! 141: if (ny==1) return subisspec(x,*y,nx);
! 142: zd = (GEN)avma;
! 143: lz = nx+2; (void)new_chunk(lz);
! 144: xd = x + nx;
! 145: yd = y + ny;
! 146: *--zd = subll(*--xd, *--yd);
! 147: while (yd > y) *--zd = subllx(*--xd, *--yd);
! 148: if (overflow)
! 149: for(;;)
! 150: {
! 151: *--zd = ((ulong)*--xd) - 1;
! 152: if (*xd) break;
! 153: }
! 154: if (xd == x)
! 155: while (*zd == 0) { zd++; lz--; } /* shorten z */
! 156: else
! 157: do *--zd = *--xd; while (xd > x);
! 158: *--zd = evalsigne(1) | evallgefint(lz);
! 159: *--zd = evaltyp(t_INT) | evallg(lz);
! 160: avma=(long)zd; return zd;
! 161: }
! 162:
! 163: #ifndef __M68K__
! 164:
! 165: /* prototype of positive small ints */
! 166: static long pos_s[] = {
! 167: evaltyp(t_INT) | m_evallg(3), evalsigne(1) | evallgefint(3), 0 };
! 168:
! 169: /* prototype of negative small ints */
! 170: static long neg_s[] = {
! 171: evaltyp(t_INT) | m_evallg(3), evalsigne(-1) | evallgefint(3), 0 };
! 172:
! 173: void
! 174: affir(GEN x, GEN y)
! 175: {
! 176: const long s=signe(x),ly=lg(y);
! 177: long lx,sh,i;
! 178:
! 179: if (!s)
! 180: {
! 181: y[1] = evalexpo(-bit_accuracy(ly));
! 182: y[2] = 0; return;
! 183: }
! 184:
! 185: lx=lgefint(x); sh=bfffo(x[2]);
! 186: y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
! 187: if (sh)
! 188: {
! 189: if (lx>ly)
! 190: { shift_left(y,x,2,ly-1, x[ly],sh); }
! 191: else
! 192: {
! 193: for (i=lx; i<ly; i++) y[i]=0;
! 194: shift_left(y,x,2,lx-1, 0,sh);
! 195: }
! 196: return;
! 197: }
! 198:
! 199: if (lx>=ly)
! 200: for (i=2; i<ly; i++) y[i]=x[i];
! 201: else
! 202: {
! 203: for (i=2; i<lx; i++) y[i]=x[i];
! 204: for ( ; i<ly; i++) y[i]=0;
! 205: }
! 206: }
! 207:
! 208: void
! 209: affrr(GEN x, GEN y)
! 210: {
! 211: long lx,ly,i;
! 212:
! 213: y[1] = x[1]; if (!signe(x)) { y[2]=0; return; }
! 214:
! 215: lx=lg(x); ly=lg(y);
! 216: if (lx>=ly)
! 217: for (i=2; i<ly; i++) y[i]=x[i];
! 218: else
! 219: {
! 220: for (i=2; i<lx; i++) y[i]=x[i];
! 221: for ( ; i<ly; i++) y[i]=0;
! 222: }
! 223: }
! 224:
! 225: GEN
! 226: shifti(GEN x, long n)
! 227: {
! 228: const long s=signe(x);
! 229: long lx,ly,i,m;
! 230: GEN y;
! 231:
! 232: if (!s) return gzero;
! 233: if (!n) return icopy(x);
! 234: lx = lgefint(x);
! 235: if (n > 0)
! 236: {
! 237: GEN z = (GEN)avma;
! 238: long d = n>>TWOPOTBITS_IN_LONG;
! 239:
! 240: ly = lx+d; y = new_chunk(ly);
! 241: for ( ; d; d--) *--z = 0;
! 242: m = n & (BITS_IN_LONG-1);
! 243: if (!m) for (i=2; i<lx; i++) y[i]=x[i];
! 244: else
! 245: {
! 246: register const ulong sh = BITS_IN_LONG - m;
! 247: shift_left2(y,x, 2,lx-1, 0,m,sh);
! 248: i = ((ulong)x[2]) >> sh;
! 249: if (i) { ly++; y = new_chunk(1); y[2] = i; }
! 250: }
! 251: }
! 252: else
! 253: {
! 254: n = -n;
! 255: ly = lx - (n>>TWOPOTBITS_IN_LONG);
! 256: if (ly<3) return gzero;
! 257: y = new_chunk(ly);
! 258: m = n & (BITS_IN_LONG-1);
! 259: if (!m) for (i=2; i<ly; i++) y[i]=x[i];
! 260: else
! 261: {
! 262: shift_right(y,x, 2,ly, 0,m);
! 263: if (y[2] == 0)
! 264: {
! 265: if (ly==3) { avma = (long)(y+3); return gzero; }
! 266: ly--; avma = (long)(++y);
! 267: }
! 268: }
! 269: }
! 270: y[1] = evalsigne(s)|evallgefint(ly);
! 271: y[0] = evaltyp(t_INT)|evallg(ly); return y;
! 272: }
! 273:
! 274: GEN
! 275: mptrunc(GEN x)
! 276: {
! 277: long d,e,i,s,m;
! 278: GEN y;
! 279:
! 280: if (typ(x)==t_INT) return icopy(x);
! 281: if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gzero;
! 282: d = (e>>TWOPOTBITS_IN_LONG) + 3;
! 283: m = e & (BITS_IN_LONG-1);
! 284: if (d > lg(x)) err(truer2);
! 285:
! 286: y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
! 287: if (++m == BITS_IN_LONG)
! 288: for (i=2; i<d; i++) y[i]=x[i];
! 289: else
! 290: {
! 291: const register ulong sh = BITS_IN_LONG - m;
! 292: shift_right2(y,x, 2,d,0, sh,m);
! 293: }
! 294: return y;
! 295: }
! 296:
! 297: /* integral part */
! 298: GEN
! 299: mpent(GEN x)
! 300: {
! 301: long d,e,i,lx,m;
! 302: GEN y;
! 303:
! 304: if (typ(x)==t_INT) return icopy(x);
! 305: if (signe(x) >= 0) return mptrunc(x);
! 306: if ((e=expo(x)) < 0) return stoi(-1);
! 307: d = (e>>TWOPOTBITS_IN_LONG) + 3;
! 308: m = e & (BITS_IN_LONG-1);
! 309: lx=lg(x); if (d>lx) err(truer2);
! 310: y = new_chunk(d);
! 311: if (++m == BITS_IN_LONG)
! 312: {
! 313: for (i=2; i<d; i++) y[i]=x[i];
! 314: i=d; while (i<lx && !x[i]) i++;
! 315: if (i==lx) goto END;
! 316: }
! 317: else
! 318: {
! 319: const register ulong sh = BITS_IN_LONG - m;
! 320: shift_right2(y,x, 2,d,0, sh,m);
! 321: if (x[d-1]<<m == 0)
! 322: {
! 323: i=d; while (i<lx && !x[i]) i++;
! 324: if (i==lx) goto END;
! 325: }
! 326: }
! 327: /* set y:=y+1 */
! 328: for (i=d-1; i>=2; i--) { y[i]++; if (y[i]) goto END; }
! 329: y=new_chunk(1); y[2]=1; d++;
! 330: END:
! 331: y[1] = evalsigne(-1) | evallgefint(d);
! 332: y[0] = evaltyp(t_INT) | evallg(d); return y;
! 333: }
! 334:
! 335: int
! 336: cmpsi(long x, GEN y)
! 337: {
! 338: ulong p;
! 339:
! 340: if (!x) return -signe(y);
! 341:
! 342: if (x>0)
! 343: {
! 344: if (signe(y)<=0) return 1;
! 345: if (lgefint(y)>3) return -1;
! 346: p=y[2]; if (p == (ulong)x) return 0;
! 347: return p < (ulong)x ? 1 : -1;
! 348: }
! 349:
! 350: if (signe(y)>=0) return -1;
! 351: if (lgefint(y)>3) return 1;
! 352: p=y[2]; if (p == (ulong)-x) return 0;
! 353: return p < (ulong)(-x) ? -1 : 1;
! 354: }
! 355:
! 356: int
! 357: cmpii(GEN x, GEN y)
! 358: {
! 359: const long sx = signe(x), sy = signe(y);
! 360: long lx,ly,i;
! 361:
! 362: if (sx<sy) return -1;
! 363: if (sx>sy) return 1;
! 364: if (!sx) return 0;
! 365:
! 366: lx=lgefint(x); ly=lgefint(y);
! 367: if (lx>ly) return sx;
! 368: if (lx<ly) return -sx;
! 369: i=2; while (i<lx && x[i]==y[i]) i++;
! 370: if (i==lx) return 0;
! 371: return ((ulong)x[i] > (ulong)y[i]) ? sx : -sx;
! 372: }
! 373:
! 374: int
! 375: cmprr(GEN x, GEN y)
! 376: {
! 377: const long sx = signe(x), sy = signe(y);
! 378: long ex,ey,lx,ly,lz,i;
! 379:
! 380: if (sx<sy) return -1;
! 381: if (sx>sy) return 1;
! 382: if (!sx) return 0;
! 383:
! 384: ex=expo(x); ey=expo(y);
! 385: if (ex>ey) return sx;
! 386: if (ex<ey) return -sx;
! 387:
! 388: lx=lg(x); ly=lg(y); lz = (lx<ly)?lx:ly;
! 389: i=2; while (i<lz && x[i]==y[i]) i++;
! 390: if (i<lz) return ((ulong)x[i] > (ulong)y[i]) ? sx : -sx;
! 391: if (lx>=ly)
! 392: {
! 393: while (i<lx && !x[i]) i++;
! 394: return (i==lx) ? 0 : sx;
! 395: }
! 396: while (i<ly && !y[i]) i++;
! 397: return (i==ly) ? 0 : -sx;
! 398: }
! 399:
! 400: GEN
! 401: addss(long x, long y)
! 402: {
! 403: if (!x) return stoi(y);
! 404: if (x>0) { pos_s[2] = x; return addsi(y,pos_s); }
! 405: neg_s[2] = -x; return addsi(y,neg_s);
! 406: }
! 407:
! 408: GEN
! 409: addsi(long x, GEN y)
! 410: {
! 411: long sx,sy,ly;
! 412: GEN z;
! 413:
! 414: if (!x) return icopy(y);
! 415: sy=signe(y); if (!sy) return stoi(x);
! 416: if (x<0) { sx=-1; x=-x; } else sx=1;
! 417: if (sx==sy)
! 418: {
! 419: z = addsispec(x,y+2, lgefint(y)-2);
! 420: setsigne(z,sy); return z;
! 421: }
! 422: ly=lgefint(y);
! 423: if (ly==3)
! 424: {
! 425: const long d = y[2] - x;
! 426: if (!d) return gzero;
! 427: z=cgeti(3);
! 428: if (y[2] < 0 || d > 0) {
! 429: z[1] = evalsigne(sy) | evallgefint(3);
! 430: z[2] = d;
! 431: }
! 432: else {
! 433: z[1] = evalsigne(-sy) | evallgefint(3);
! 434: z[2] =-d;
! 435: }
! 436: return z;
! 437: }
! 438: z = subisspec(y+2,x, ly-2);
! 439: setsigne(z,sy); return z;
! 440: }
! 441:
! 442: GEN
! 443: addii(GEN x, GEN y)
! 444: {
! 445: long sx = signe(x);
! 446: long sy = signe(y);
! 447: long lx,ly;
! 448: GEN z;
! 449:
! 450: if (!sx) return sy? icopy(y): gzero;
! 451: if (!sy) return icopy(x);
! 452: lx=lgefint(x); ly=lgefint(y);
! 453:
! 454: if (sx==sy)
! 455: z = addiispec(x+2,y+2,lx-2,ly-2);
! 456: else
! 457: { /* sx != sy */
! 458: long i = lx - ly;
! 459: if (i==0)
! 460: {
! 461: i = absi_cmp(x,y);
! 462: if (!i) return gzero;
! 463: }
! 464: if (i<0) { sx=sy; swapspec(x,y, lx,ly); } /* ensure |x| >= |y| */
! 465: z = subiispec(x+2,y+2,lx-2,ly-2);
! 466: }
! 467: setsigne(z,sx); return z;
! 468: }
! 469:
! 470: GEN
! 471: addsr(long x, GEN y)
! 472: {
! 473: if (!x) return rcopy(y);
! 474: if (x>0) { pos_s[2]=x; return addir(pos_s,y); }
! 475: neg_s[2] = -x; return addir(neg_s,y);
! 476: }
! 477:
! 478: GEN
! 479: addir(GEN x, GEN y)
! 480: {
! 481: long e,l,ly;
! 482: GEN z;
! 483:
! 484: if (!signe(x)) return rcopy(y);
! 485: if (!signe(y))
! 486: {
! 487: l=lgefint(x)-(expo(y)>>TWOPOTBITS_IN_LONG);
! 488: if (l<3) err(adder3);
! 489: z=cgetr(l); affir(x,z); return z;
! 490: }
! 491:
! 492: e = expo(y)-expi(x); ly=lg(y);
! 493: if (e>0)
! 494: {
! 495: l = ly - (e>>TWOPOTBITS_IN_LONG);
! 496: if (l<3) return rcopy(y);
! 497: }
! 498: else l = ly + ((-e)>>TWOPOTBITS_IN_LONG)+1;
! 499: z=cgetr(l); affir(x,z); y=addrr(z,y);
! 500: z = y+l; ly = lg(y); while (ly--) z[ly] = y[ly];
! 501: avma=(long)z; return z;
! 502: }
! 503:
! 504: GEN
! 505: addrr(GEN x, GEN y)
! 506: {
! 507: long sx=signe(x),sy=signe(y),ex=expo(x),ey=expo(y);
! 508: long e,m,flag,i,j,f2,lx,ly,lz;
! 509: GEN z;
! 510: LOCAL_OVERFLOW;
! 511:
! 512: e = ey-ex;
! 513: if (!sy)
! 514: {
! 515: if (!sx)
! 516: {
! 517: if (e > 0) ex=ey;
! 518: z=cgetr(3); z[1]=evalexpo(ex); z[2]=0; return z;
! 519: }
! 520: if (e > 0) { z=cgetr(3); z[1]=evalexpo(ey); z[2]=0; return z; }
! 521: lz = 3 + ((-e)>>TWOPOTBITS_IN_LONG);
! 522: lx = lg(x); if (lz>lx) lz=lx;
! 523: z = cgetr(lz); while(--lz) z[lz]=x[lz];
! 524: return z;
! 525: }
! 526: if (!sx)
! 527: {
! 528: if (e < 0) { z=cgetr(3); z[1]=evalexpo(ex); z[2]=0; return z; }
! 529: lz = 3 + (e>>TWOPOTBITS_IN_LONG);
! 530: ly = lg(y); if (lz>ly) lz=ly;
! 531: z = cgetr(lz); while (--lz) z[lz]=y[lz];
! 532: return z;
! 533: }
! 534:
! 535: if (e < 0) { z=x; x=y; y=z; ey=ex; i=sx; sx=sy; sy=i; e=-e; }
! 536: /* now ey >= ex */
! 537: lx=lg(x); ly=lg(y);
! 538: if (e)
! 539: {
! 540: long d = e >> TWOPOTBITS_IN_LONG, l = ly-d;
! 541: if (l > lx) { flag=0; lz=lx+d+1; }
! 542: else if (l > 2) { flag=1; lz=ly; lx=l; }
! 543: else return rcopy(y);
! 544: m = e & (BITS_IN_LONG-1);
! 545: }
! 546: else
! 547: {
! 548: if (lx > ly) lx = ly;
! 549: flag=2; lz=lx; m=0;
! 550: }
! 551: z = cgetr(lz);
! 552: if (m)
! 553: { /* shift right x m bits */
! 554: const ulong sh = BITS_IN_LONG-m;
! 555: GEN p1 = x; x = new_chunk(lx+1);
! 556: shift_right2(x,p1,2,lx, 0,m,sh);
! 557: if (flag==0)
! 558: {
! 559: x[lx] = p1[lx-1] << sh;
! 560: if (x[lx]) { flag = 3; lx++; }
! 561: }
! 562: }
! 563:
! 564: if (sx==sy)
! 565: { /* addition */
! 566: i=lz-1; avma = (long)z;
! 567: if (flag==0) { z[i] = y[i]; i--; }
! 568: overflow=0;
! 569: for (j=lx-1; j>=2; i--,j--) z[i] = addllx(x[j],y[i]);
! 570:
! 571: if (overflow)
! 572: for (;;)
! 573: {
! 574: if (i==1)
! 575: {
! 576: shift_right(z,z, 2,lz, 1,1);
! 577: ey=evalexpo(ey+1); if (ey & ~EXPOBITS) err(adder4);
! 578: z[1] = evalsigne(sx) | ey; return z;
! 579: }
! 580: z[i] = y[i]+1; if (z[i--]) break;
! 581: }
! 582: for (; i>=2; i--) z[i]=y[i];
! 583: z[1] = evalsigne(sx) | evalexpo(ey); return z;
! 584: }
! 585:
! 586: /* subtraction */
! 587: if (e) f2 = 1;
! 588: else
! 589: {
! 590: i=2; while (i<lx && x[i]==y[i]) i++;
! 591: if (i==lx)
! 592: {
! 593: avma = (long)(z+lz);
! 594: e = evalexpo(ey - bit_accuracy(lx));
! 595: if (e & ~EXPOBITS) err(adder4);
! 596: z=cgetr(3); z[1]=e; z[2]=0; return z;
! 597: }
! 598: f2 = ((ulong)y[i] > (ulong)x[i]);
! 599: }
! 600:
! 601: /* result is non-zero. f2 = (y > x) */
! 602: i=lz-1;
! 603: if (f2)
! 604: {
! 605: if (flag==0) { z[i] = y[i]; i--; }
! 606: j=lx-1; z[i] = subll(y[i],x[j]); i--; j--;
! 607: for (; j>=2; i--,j--) z[i] = subllx(y[i],x[j]);
! 608: if (overflow)
! 609: for (;;) { z[i] = y[i]-1; if (y[i--]) break; }
! 610: for (; i>=2; i--) z[i]=y[i];
! 611: sx = sy;
! 612: }
! 613: else
! 614: {
! 615: if (flag==0) { z[i] = -y[i]; i--; overflow=1; } else overflow=0;
! 616: for (; i>=2; i--) z[i]=subllx(x[i],y[i]);
! 617: }
! 618:
! 619: x = z+2; i=0; while (!x[i]) i++;
! 620: lz -= i; z += i; m = bfffo(z[2]);
! 621: e = evalexpo(ey - (m | (i<<TWOPOTBITS_IN_LONG)));
! 622: if (e & ~EXPOBITS) err(adder5);
! 623: if (m) shift_left(z,z,2,lz-1, 0,m);
! 624: z[1] = evalsigne(sx) | e;
! 625: z[0] = evaltyp(t_REAL) | evallg(lz);
! 626: avma = (long)z; return z;
! 627: }
! 628:
! 629: GEN
! 630: mulss(long x, long y)
! 631: {
! 632: long s,p1;
! 633: GEN z;
! 634: LOCAL_HIREMAINDER;
! 635:
! 636: if (!x || !y) return gzero;
! 637: if (x<0) { s = -1; x = -x; } else s=1;
! 638: if (y<0) { s = -s; y = -y; }
! 639: p1 = mulll(x,y);
! 640: if (hiremainder)
! 641: {
! 642: z=cgeti(4); z[1] = evalsigne(s) | evallgefint(4);
! 643: z[2]=hiremainder; z[3]=p1; return z;
! 644: }
! 645: z=cgeti(3); z[1] = evalsigne(s) | evallgefint(3);
! 646: z[2]=p1; return z;
! 647: }
! 648: #endif
! 649:
! 650: GEN
! 651: muluu(ulong x, ulong y)
! 652: {
! 653: long p1;
! 654: GEN z;
! 655: LOCAL_HIREMAINDER;
! 656:
! 657: if (!x || !y) return gzero;
! 658: p1 = mulll(x,y);
! 659: if (hiremainder)
! 660: {
! 661: z=cgeti(4); z[1] = evalsigne(1) | evallgefint(4);
! 662: z[2]=hiremainder; z[3]=p1; return z;
! 663: }
! 664: z=cgeti(3); z[1] = evalsigne(1) | evallgefint(3);
! 665: z[2]=p1; return z;
! 666: }
! 667:
! 668: /* assume ny > 0 */
! 669: #ifdef INLINE
! 670: INLINE
! 671: #endif
! 672: GEN
! 673: mulsispec(long x, GEN y, long ny)
! 674: {
! 675: GEN yd, z = (GEN)avma;
! 676: long lz = ny+3;
! 677: LOCAL_HIREMAINDER;
! 678:
! 679: (void)new_chunk(lz);
! 680: yd = y + ny; *--z = mulll(x, *--yd);
! 681: while (yd > y) *--z = addmul(x,*--yd);
! 682: if (hiremainder) *--z = hiremainder; else lz--;
! 683: *--z = evalsigne(1) | evallgefint(lz);
! 684: *--z = evaltyp(t_INT) | evallg(lz);
! 685: avma=(long)z; return z;
! 686: }
! 687:
! 688: GEN
! 689: mului(ulong x, GEN y)
! 690: {
! 691: long s = signe(y);
! 692: GEN z;
! 693:
! 694: if (!s || !x) return gzero;
! 695: z = mulsispec(x, y+2, lgefint(y)-2);
! 696: setsigne(z,s); return z;
! 697: }
! 698:
! 699: #ifndef __M68K__
! 700:
! 701: GEN
! 702: mulsi(long x, GEN y)
! 703: {
! 704: long s = signe(y);
! 705: GEN z;
! 706:
! 707: if (!s || !x) return gzero;
! 708: if (x<0) { s = -s; x = -x; }
! 709: z = mulsispec(x, y+2, lgefint(y)-2);
! 710: setsigne(z,s); return z;
! 711: }
! 712:
! 713: GEN
! 714: mulsr(long x, GEN y)
! 715: {
! 716: long lx,i,s,garde,e,sh,m;
! 717: GEN z;
! 718: LOCAL_HIREMAINDER;
! 719:
! 720: if (!x) return gzero;
! 721: s = signe(y);
! 722: if (!s)
! 723: {
! 724: if (x<0) x = -x;
! 725: e = y[1] + (BITS_IN_LONG-1)-bfffo(x);
! 726: if (e & ~EXPOBITS) err(muler2);
! 727: z=cgetr(3); z[1]=e; z[2]=0; return z;
! 728: }
! 729: if (x<0) { s = -s; x = -x; }
! 730: if (x==1) { z=rcopy(y); setsigne(z,s); return z; }
! 731:
! 732: lx = lg(y); e = expo(y);
! 733: z=cgetr(lx); y--; garde=mulll(x,y[lx]);
! 734: for (i=lx-1; i>=3; i--) z[i]=addmul(x,y[i]);
! 735: z[2]=hiremainder;
! 736:
! 737: sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;
! 738: if (sh) shift_left2(z,z, 2,lx-1, garde,sh,m);
! 739: e = evalexpo(m+e); if (e & ~EXPOBITS) err(muler2);
! 740: z[1] = evalsigne(s) | e; return z;
! 741: }
! 742:
! 743: GEN
! 744: mulrr(GEN x, GEN y)
! 745: {
! 746: long sx = signe(x), sy = signe(y);
! 747: long i,j,ly,lz,lzz,e,flag,p1;
! 748: ulong garde;
! 749: GEN z, y1;
! 750: LOCAL_HIREMAINDER;
! 751: LOCAL_OVERFLOW;
! 752:
! 753: e = evalexpo(expo(x)+expo(y)); if (e & ~EXPOBITS) err(muler4);
! 754: if (!sx || !sy) { z=cgetr(3); z[1]=e; z[2]=0; return z; }
! 755: if (sy<0) sx = -sx;
! 756: lz=lg(x); ly=lg(y);
! 757: if (lz>ly) { lz=ly; z=x; x=y; y=z; flag=1; } else flag = (lz!=ly);
! 758: z=cgetr(lz); z[1] = evalsigne(sx) | e;
! 759: if (lz==3)
! 760: {
! 761: if (flag)
! 762: {
! 763: (void)mulll(x[2],y[3]);
! 764: garde = addmul(x[2],y[2]);
! 765: }
! 766: else
! 767: garde = mulll(x[2],y[2]);
! 768: if ((long)hiremainder<0) { z[2]=hiremainder; z[1]++; }
! 769: else z[2]=(hiremainder<<1) | (garde>>(BITS_IN_LONG-1));
! 770: return z;
! 771: }
! 772:
! 773: if (flag) { (void)mulll(x[2],y[lz]); garde = hiremainder; } else garde=0;
! 774: lzz=lz-1; p1=x[lzz];
! 775: if (p1)
! 776: {
! 777: (void)mulll(p1,y[3]);
! 778: garde = addll(addmul(p1,y[2]), garde);
! 779: z[lzz] = overflow+hiremainder;
! 780: }
! 781: else z[lzz]=0;
! 782: for (j=lz-2, y1=y-j; j>=3; j--)
! 783: {
! 784: p1 = x[j]; y1++;
! 785: if (p1)
! 786: {
! 787: (void)mulll(p1,y1[lz+1]);
! 788: garde = addll(addmul(p1,y1[lz]), garde);
! 789: for (i=lzz; i>j; i--)
! 790: {
! 791: hiremainder += overflow;
! 792: z[i] = addll(addmul(p1,y1[i]), z[i]);
! 793: }
! 794: z[j] = hiremainder+overflow;
! 795: }
! 796: else z[j]=0;
! 797: }
! 798: p1=x[2]; y1++;
! 799: garde = addll(mulll(p1,y1[lz]), garde);
! 800: for (i=lzz; i>2; i--)
! 801: {
! 802: hiremainder += overflow;
! 803: z[i] = addll(addmul(p1,y1[i]), z[i]);
! 804: }
! 805: z[2] = hiremainder+overflow;
! 806: if (z[2] < 0) z[1]++;
! 807: else
! 808: shift_left(z,z,2,lzz,garde, 1);
! 809: return z;
! 810: }
! 811:
! 812: GEN
! 813: mulir(GEN x, GEN y)
! 814: {
! 815: long sx=signe(x),sy,lz,lzz,ey,e,p1,i,j;
! 816: ulong garde;
! 817: GEN z,y1;
! 818: LOCAL_HIREMAINDER;
! 819: LOCAL_OVERFLOW;
! 820:
! 821: if (!sx) return gzero;
! 822: if (!is_bigint(x)) return mulsr(itos(x),y);
! 823: sy=signe(y); ey=expo(y);
! 824: if (!sy)
! 825: {
! 826: e = evalexpo(expi(x)+ey);
! 827: if (e & ~EXPOBITS) err(muler6);
! 828: z=cgetr(3); z[1]=e; z[2]=0; return z;
! 829: }
! 830:
! 831: if (sy<0) sx = -sx;
! 832: lz=lg(y); z=cgetr(lz);
! 833: y1=cgetr(lz+1);
! 834: affir(x,y1); x=y; y=y1;
! 835: e = evalexpo(expo(y)+ey);
! 836: if (e & ~EXPOBITS) err(muler4);
! 837: z[1] = evalsigne(sx) | e;
! 838: if (lz==3)
! 839: {
! 840: (void)mulll(x[2],y[3]);
! 841: garde=addmul(x[2],y[2]);
! 842: if ((long)hiremainder < 0) { z[2]=hiremainder; z[1]++; }
! 843: else z[2]=(hiremainder<<1) | (garde>>(BITS_IN_LONG-1));
! 844: avma=(long)z; return z;
! 845: }
! 846:
! 847: (void)mulll(x[2],y[lz]); garde=hiremainder;
! 848: lzz=lz-1; p1=x[lzz];
! 849: if (p1)
! 850: {
! 851: (void)mulll(p1,y[3]);
! 852: garde=addll(addmul(p1,y[2]),garde);
! 853: z[lzz] = overflow+hiremainder;
! 854: }
! 855: else z[lzz]=0;
! 856: for (j=lz-2, y1=y-j; j>=3; j--)
! 857: {
! 858: p1=x[j]; y1++;
! 859: if (p1)
! 860: {
! 861: (void)mulll(p1,y1[lz+1]);
! 862: garde = addll(addmul(p1,y1[lz]), garde);
! 863: for (i=lzz; i>j; i--)
! 864: {
! 865: hiremainder += overflow;
! 866: z[i] = addll(addmul(p1,y1[i]), z[i]);
! 867: }
! 868: z[j] = hiremainder+overflow;
! 869: }
! 870: else z[j]=0;
! 871: }
! 872: p1=x[2]; y1++;
! 873: garde = addll(mulll(p1,y1[lz]), garde);
! 874: for (i=lzz; i>2; i--)
! 875: {
! 876: hiremainder += overflow;
! 877: z[i] = addll(addmul(p1,y1[i]), z[i]);
! 878: }
! 879: z[2] = hiremainder+overflow;
! 880: if (z[2] < 0) z[1]++;
! 881: else
! 882: shift_left(z,z,2,lzz,garde, 1);
! 883: avma=(long)z; return z;
! 884: }
! 885:
! 886: /* written by Bruno Haible following an idea of Robert Harley */
! 887: long
! 888: vals(ulong z)
! 889: {
! 890: static char tab[64]={-1,0,1,12,2,6,-1,13,3,-1,7,-1,-1,-1,-1,14,10,4,-1,-1,8,-1,-1,25,-1,-1,-1,-1,-1,21,27,15,31,11,5,-1,-1,-1,-1,-1,9,-1,-1,24,-1,-1,20,26,30,-1,-1,-1,-1,23,-1,19,29,-1,22,18,28,17,16,-1};
! 891: #ifdef LONG_IS_64BIT
! 892: long s;
! 893: #endif
! 894:
! 895: if (!z) return -1;
! 896: #ifdef LONG_IS_64BIT
! 897: if (! (z&0xffffffff)) { s = 32; z >>=32; } else s = 0;
! 898: #endif
! 899: z |= ~z + 1;
! 900: z += z << 4;
! 901: z += z << 6;
! 902: z ^= z << 16; /* or z -= z<<16 */
! 903: #ifdef LONG_IS_64BIT
! 904: return s + tab[(z&0xffffffff)>>26];
! 905: #else
! 906: return tab[z>>26];
! 907: #endif
! 908: }
! 909:
! 910: GEN
! 911: modss(long x, long y)
! 912: {
! 913: LOCAL_HIREMAINDER;
! 914:
! 915: if (!y) err(moder1);
! 916: if (y<0) y=-y;
! 917: hiremainder=0; divll(labs(x),y);
! 918: if (!hiremainder) return gzero;
! 919: return (x < 0) ? stoi(y-hiremainder) : stoi(hiremainder);
! 920: }
! 921:
! 922: GEN
! 923: resss(long x, long y)
! 924: {
! 925: LOCAL_HIREMAINDER;
! 926:
! 927: if (!y) err(reser1);
! 928: hiremainder=0; divll(labs(x),labs(y));
! 929: return (x < 0) ? stoi(-((long)hiremainder)) : stoi(hiremainder);
! 930: }
! 931:
! 932: GEN
! 933: divsi(long x, GEN y)
! 934: {
! 935: long p1, s = signe(y);
! 936: LOCAL_HIREMAINDER;
! 937:
! 938: if (!s) err(diver2);
! 939: if (!x || lgefint(y)>3 || ((long)y[2])<0)
! 940: {
! 941: hiremainder=x; SAVE_HIREMAINDER; return gzero;
! 942: }
! 943: hiremainder=0; p1=divll(labs(x),y[2]);
! 944: if (x<0) { hiremainder = -((long)hiremainder); p1 = -p1; }
! 945: if (s<0) p1 = -p1;
! 946: SAVE_HIREMAINDER; return stoi(p1);
! 947: }
! 948: #endif
! 949:
! 950: GEN
! 951: modui(ulong x, GEN y)
! 952: {
! 953: LOCAL_HIREMAINDER;
! 954:
! 955: if (!signe(y)) err(diver2);
! 956: if (!x || lgefint(y)>3) hiremainder=x;
! 957: else
! 958: {
! 959: hiremainder=0; (void)divll(x,y[2]);
! 960: }
! 961: return utoi(hiremainder);
! 962: }
! 963:
! 964: GEN
! 965: modiu(GEN y, ulong x)
! 966: {
! 967: long sy=signe(y),ly,i;
! 968: LOCAL_HIREMAINDER;
! 969:
! 970: if (!x) err(diver4);
! 971: if (!sy) return gzero;
! 972: ly = lgefint(y);
! 973: if (x <= (ulong)y[2]) hiremainder=0;
! 974: else
! 975: {
! 976: if (ly==3) return utoi((sy > 0)? (ulong)y[2]: x - (ulong)y[2]);
! 977: hiremainder=y[2]; ly--; y++;
! 978: }
! 979: for (i=2; i<ly; i++) (void)divll(y[i],x);
! 980: return utoi((sy > 0)? hiremainder: x - hiremainder);
! 981: }
! 982:
! 983: #ifndef __M68K__
! 984:
! 985: GEN
! 986: modsi(long x, GEN y)
! 987: {
! 988: long s = signe(y);
! 989: GEN p1;
! 990: LOCAL_HIREMAINDER;
! 991:
! 992: if (!s) err(diver2);
! 993: if (!x || lgefint(y)>3 || ((long)y[2])<0) hiremainder=x;
! 994: else
! 995: {
! 996: hiremainder=0; (void)divll(labs(x),y[2]);
! 997: if (x<0) hiremainder = -((long)hiremainder);
! 998: }
! 999: if (!hiremainder) return gzero;
! 1000: if (x>0) return stoi(hiremainder);
! 1001: if (s<0)
! 1002: { setsigne(y,1); p1=addsi(hiremainder,y); setsigne(y,-1); }
! 1003: else
! 1004: p1=addsi(hiremainder,y);
! 1005: return p1;
! 1006: }
! 1007:
! 1008: GEN
! 1009: divis(GEN y, long x)
! 1010: {
! 1011: long sy=signe(y),ly,s,i;
! 1012: GEN z;
! 1013: LOCAL_HIREMAINDER;
! 1014:
! 1015: if (!x) err(diver4);
! 1016: if (!sy) { hiremainder=0; SAVE_HIREMAINDER; return gzero; }
! 1017: if (x<0) { s = -sy; x = -x; } else s=sy;
! 1018:
! 1019: ly = lgefint(y);
! 1020: if ((ulong)x <= (ulong)y[2]) hiremainder=0;
! 1021: else
! 1022: {
! 1023: if (ly==3) { hiremainder=itos(y); SAVE_HIREMAINDER; return gzero; }
! 1024: hiremainder=y[2]; ly--; y++;
! 1025: }
! 1026: z = cgeti(ly); z[1] = evallgefint(ly) | evalsigne(s);
! 1027: for (i=2; i<ly; i++) z[i]=divll(y[i],x);
! 1028: if (sy<0) hiremainder = - ((long)hiremainder);
! 1029: SAVE_HIREMAINDER; return z;
! 1030: }
! 1031:
! 1032: GEN
! 1033: divir(GEN x, GEN y)
! 1034: {
! 1035: GEN xr,z;
! 1036: long av,ly;
! 1037:
! 1038: if (!signe(y)) err(diver5);
! 1039: if (!signe(x)) return gzero;
! 1040: ly=lg(y); z=cgetr(ly); av=avma; xr=cgetr(ly+1); affir(x,xr);
! 1041: affrr(divrr(xr,y),z); avma=av; return z;
! 1042: }
! 1043:
! 1044: GEN
! 1045: divri(GEN x, GEN y)
! 1046: {
! 1047: GEN yr,z;
! 1048: long av,lx,s=signe(y);
! 1049:
! 1050: if (!s) err(diver8);
! 1051: if (!signe(x))
! 1052: {
! 1053: const long ex = evalexpo(expo(x) - expi(y));
! 1054:
! 1055: if (ex<0) err(diver12);
! 1056: z=cgetr(3); z[1]=ex; z[2]=0; return z;
! 1057: }
! 1058: if (!is_bigint(y)) return divrs(x, s>0? y[2]: -y[2]);
! 1059:
! 1060: lx=lg(x); z=cgetr(lx);
! 1061: av=avma; yr=cgetr(lx+1); affir(y,yr);
! 1062: affrr(divrr(x,yr),z); avma=av; return z;
! 1063: }
! 1064:
! 1065: void
! 1066: diviiz(GEN x, GEN y, GEN z)
! 1067: {
! 1068: long av=avma,lz;
! 1069: GEN xr,yr;
! 1070:
! 1071: if (typ(z)==t_INT) { affii(divii(x,y),z); avma=av; return; }
! 1072: lz=lg(z); xr=cgetr(lz); affir(x,xr); yr=cgetr(lz); affir(y,yr);
! 1073: affrr(divrr(xr,yr),z); avma=av;
! 1074: }
! 1075:
! 1076: void
! 1077: mpdivz(GEN x, GEN y, GEN z)
! 1078: {
! 1079: long av=avma;
! 1080:
! 1081: if (typ(z)==t_INT)
! 1082: {
! 1083: if (typ(x)==t_REAL || typ(y)==t_REAL) err(divzer1);
! 1084: affii(divii(x,y),z); avma=av; return;
! 1085: }
! 1086: if (typ(x)==t_INT)
! 1087: {
! 1088: GEN xr,yr;
! 1089: long lz;
! 1090:
! 1091: if (typ(y)==t_REAL) { affrr(divir(x,y),z); avma=av; return; }
! 1092: lz=lg(z); xr=cgetr(lz); affir(x,xr); yr=cgetr(lz); affir(y,yr);
! 1093: affrr(divrr(xr,yr),z); avma=av; return;
! 1094: }
! 1095: if (typ(y)==t_REAL) { affrr(divrr(x,y),z); avma=av; return; }
! 1096: affrr(divri(x,y),z); avma=av;
! 1097: }
! 1098:
! 1099: GEN
! 1100: divsr(long x, GEN y)
! 1101: {
! 1102: long av,ly;
! 1103: GEN p1,z;
! 1104:
! 1105: if (!signe(y)) err(diver3);
! 1106: if (!x) return gzero;
! 1107: ly=lg(y); z=cgetr(ly); av=avma;
! 1108: p1=cgetr(ly+1); affsr(x,p1); affrr(divrr(p1,y),z);
! 1109: avma=av; return z;
! 1110: }
! 1111:
! 1112: GEN
! 1113: modii(GEN x, GEN y)
! 1114: {
! 1115: switch(signe(x))
! 1116: {
! 1117: case 0: return gzero;
! 1118: case 1: return resii(x,y);
! 1119: default:
! 1120: {
! 1121: long av = avma;
! 1122: (void)new_chunk(lgefint(y));
! 1123: x = resii(x,y); avma=av;
! 1124: if (x==gzero) return x;
! 1125: return subiispec(y+2,x+2,lgefint(y)-2,lgefint(x)-2);
! 1126: }
! 1127: }
! 1128: }
! 1129:
! 1130: void
! 1131: modiiz(GEN x, GEN y, GEN z)
! 1132: {
! 1133: const long av = avma;
! 1134: affii(modii(x,y),z); avma=av;
! 1135: }
! 1136:
! 1137: GEN
! 1138: divrs(GEN x, long y)
! 1139: {
! 1140: long i,lx,ex,garde,sh,s=signe(x);
! 1141: GEN z;
! 1142: LOCAL_HIREMAINDER;
! 1143:
! 1144: if (!y) err(diver6);
! 1145: if (!s)
! 1146: {
! 1147: z=cgetr(3); z[1] = x[1] - (BITS_IN_LONG-1) + bfffo(y);
! 1148: if (z[1]<0) err(diver7);
! 1149: z[2]=0; return z;
! 1150: }
! 1151: if (y<0) { s = -s; y = -y; }
! 1152: if (y==1) { z=rcopy(x); setsigne(z,s); return z; }
! 1153:
! 1154: z=cgetr(lx=lg(x)); hiremainder=0;
! 1155: for (i=2; i<lx; i++) z[i]=divll(x[i],y);
! 1156:
! 1157: /* we may have hiremainder != 0 ==> garde */
! 1158: garde=divll(0,y); sh=bfffo(z[2]); ex=evalexpo(expo(x)-sh);
! 1159: if (ex & ~EXPOBITS) err(diver7);
! 1160: if (sh) shift_left(z,z, 2,lx-1, garde,sh);
! 1161: z[1] = evalsigne(s) | ex;
! 1162: return z;
! 1163: }
! 1164:
! 1165: GEN
! 1166: divrr(GEN x, GEN y)
! 1167: {
! 1168: long sx=signe(x), sy=signe(y), lx,ly,lz,e,i,j;
! 1169: ulong si,saux;
! 1170: GEN z,x1;
! 1171:
! 1172: if (!sy) err(diver9);
! 1173: e = evalexpo(expo(x) - expo(y));
! 1174: if (e & ~EXPOBITS) err(diver11);
! 1175: if (!sx) { z=cgetr(3); z[1]=e; z[2]=0; return z; }
! 1176: if (sy<0) sx = -sx;
! 1177: e = evalsigne(sx) | e;
! 1178: lx=lg(x); ly=lg(y);
! 1179: if (ly==3)
! 1180: {
! 1181: ulong k = x[2], l = (lx>3)? x[3]: 0;
! 1182: LOCAL_HIREMAINDER;
! 1183: if (k < (ulong)y[2]) e--;
! 1184: else
! 1185: {
! 1186: l >>= 1; if (k&1) l |= HIGHBIT;
! 1187: k >>= 1;
! 1188: }
! 1189: z=cgetr(3); z[1]=e;
! 1190: hiremainder=k; z[2]=divll(l,y[2]); return z;
! 1191: }
! 1192:
! 1193: lz = (lx<=ly)? lx: ly;
! 1194: x1 = (z=new_chunk(lz))-1;
! 1195: x1[1]=0; for (i=2; i<lz; i++) x1[i]=x[i];
! 1196: x1[lz] = (lx>lz)? x[lz]: 0;
! 1197: x=z; si=y[2]; saux=y[3];
! 1198: for (i=0; i<lz-1; i++)
! 1199: { /* x1 = x + (i-1) */
! 1200: ulong k,qp;
! 1201: LOCAL_HIREMAINDER;
! 1202: LOCAL_OVERFLOW;
! 1203: if ((ulong)x1[1] == si)
! 1204: {
! 1205: qp = MAXULONG; k=addll(si,x1[2]);
! 1206: }
! 1207: else
! 1208: {
! 1209: if ((ulong)x1[1] > si) /* can't happen if i=0 */
! 1210: {
! 1211: GEN y1 = y+1;
! 1212: overflow=0;
! 1213: for (j=lz-i; j>0; j--) x1[j] = subllx(x1[j],y1[j]);
! 1214: j=i; do x[--j]++; while (j && !x[j]);
! 1215: }
! 1216: hiremainder=x1[1]; overflow=0;
! 1217: qp=divll(x1[2],si); k=hiremainder;
! 1218: }
! 1219: if (!overflow)
! 1220: {
! 1221: long k3 = subll(mulll(qp,saux), x1[3]);
! 1222: long k4 = subllx(hiremainder,k);
! 1223: while (!overflow && k4) { qp--; k3=subll(k3,saux); k4=subllx(k4,si); }
! 1224: }
! 1225: j = lz-i+1;
! 1226: if (j<ly) (void)mulll(qp,y[j]); else { hiremainder=0 ; j=ly; }
! 1227: for (j--; j>1; j--)
! 1228: {
! 1229: x1[j] = subll(x1[j], addmul(qp,y[j]));
! 1230: hiremainder += overflow;
! 1231: }
! 1232: if (x1[1] != hiremainder)
! 1233: {
! 1234: if ((ulong)x1[1] < hiremainder)
! 1235: {
! 1236: qp--;
! 1237: overflow=0; for (j=lz-i; j>1; j--) x1[j]=addllx(x1[j], y[j]);
! 1238: }
! 1239: else
! 1240: {
! 1241: x1[1] -= hiremainder;
! 1242: while (x1[1])
! 1243: {
! 1244: qp++; if (!qp) { j=i; do x[--j]++; while (j && !x[j]); }
! 1245: overflow=0; for (j=lz-i; j>1; j--) x1[j]=subllx(x1[j],y[j]);
! 1246: x1[1] -= overflow;
! 1247: }
! 1248: }
! 1249: }
! 1250: x1[1]=qp; x1++;
! 1251: }
! 1252: x1 = x-1; for (j=lz-1; j>=2; j--) x[j]=x1[j];
! 1253: if (*x) { shift_right(x,x, 2,lz, 1,1); } else e--;
! 1254: x[0]=evaltyp(t_REAL)|evallg(lz);
! 1255: x[1]=e; return x;
! 1256: }
! 1257: #endif /* !defined(__M68K__) */
! 1258: /* The following ones are not in mp.s (mulii is, with a different algorithm) */
! 1259:
! 1260: /* Integer division x / y:
! 1261: * if z = ONLY_REM return remainder, otherwise return quotient
! 1262: * if z != NULL set *z to remainder
! 1263: * *z is the last object on stack (and thus can be disposed of with cgiv
! 1264: * instead of gerepile)
! 1265: * If *z is zero, we put gzero here and no copy.
! 1266: * space needed: lx + ly
! 1267: */
! 1268: GEN
! 1269: dvmdii(GEN x, GEN y, GEN *z)
! 1270: {
! 1271: long sx=signe(x),sy=signe(y);
! 1272: long av,lx,ly,lz,i,j,sh,k,lq,lr;
! 1273: ulong si,qp,saux, *xd,*rd,*qd;
! 1274: GEN r,q,x1;
! 1275:
! 1276: if (!sy) err(dvmer1);
! 1277: if (!sx)
! 1278: {
! 1279: if (!z || z == ONLY_REM) return gzero;
! 1280: *z=gzero; return gzero;
! 1281: }
! 1282: lx=lgefint(x);
! 1283: ly=lgefint(y); lz=lx-ly;
! 1284: if (lz <= 0)
! 1285: {
! 1286: if (lz == 0)
! 1287: {
! 1288: for (i=2; i<lx; i++)
! 1289: if (x[i] != y[i])
! 1290: {
! 1291: if ((ulong)x[i] > (ulong)y[i]) goto DIVIDE;
! 1292: goto TRIVIAL;
! 1293: }
! 1294: if (z == ONLY_REM) return gzero;
! 1295: if (z) *z = gzero;
! 1296: if (sx < 0) sy = -sy;
! 1297: return stoi(sy);
! 1298: }
! 1299: TRIVIAL:
! 1300: if (z == ONLY_REM) return icopy(x);
! 1301: if (z) *z = icopy(x);
! 1302: return gzero;
! 1303: }
! 1304: DIVIDE: /* quotient is non-zero */
! 1305: av=avma; if (sx<0) sy = -sy;
! 1306: if (ly==3)
! 1307: {
! 1308: LOCAL_HIREMAINDER;
! 1309: si = y[2];
! 1310: if (si <= (ulong)x[2]) hiremainder=0;
! 1311: else
! 1312: {
! 1313: hiremainder = x[2]; lx--; x++;
! 1314: }
! 1315: q = new_chunk(lx); for (i=2; i<lx; i++) q[i]=divll(x[i],si);
! 1316: if (z == ONLY_REM)
! 1317: {
! 1318: avma=av; if (!hiremainder) return gzero;
! 1319: r=cgeti(3);
! 1320: r[1] = evalsigne(sx) | evallgefint(3);
! 1321: r[2]=hiremainder; return r;
! 1322: }
! 1323: q[1] = evalsigne(sy) | evallgefint(lx);
! 1324: q[0] = evaltyp(t_INT) | evallg(lx);
! 1325: if (!z) return q;
! 1326: if (!hiremainder) { *z=gzero; return q; }
! 1327: r=cgeti(3);
! 1328: r[1] = evalsigne(sx) | evallgefint(3);
! 1329: r[2] = hiremainder; *z=r; return q;
! 1330: }
! 1331:
! 1332: x1 = new_chunk(lx); sh = bfffo(y[2]);
! 1333: if (sh)
! 1334: { /* normalize so that highbit(y) = 1 (shift left x and y by sh bits)*/
! 1335: register const ulong m = BITS_IN_LONG - sh;
! 1336: r = new_chunk(ly);
! 1337: shift_left2(r, y,2,ly-1, 0,sh,m); y = r;
! 1338: shift_left2(x1,x,2,lx-1, 0,sh,m);
! 1339: x1[1] = ((ulong)x[2]) >> m;
! 1340: }
! 1341: else
! 1342: {
! 1343: x1[1]=0; for (j=2; j<lx; j++) x1[j]=x[j];
! 1344: }
! 1345: x=x1; si=y[2]; saux=y[3];
! 1346: for (i=0; i<=lz; i++)
! 1347: { /* x1 = x + i */
! 1348: LOCAL_HIREMAINDER;
! 1349: LOCAL_OVERFLOW;
! 1350: if ((ulong)x1[1] == si)
! 1351: {
! 1352: qp = MAXULONG; k=addll(si,x1[2]);
! 1353: }
! 1354: else
! 1355: {
! 1356: hiremainder=x1[1]; overflow=0;
! 1357: qp=divll(x1[2],si); k=hiremainder;
! 1358: }
! 1359: if (!overflow)
! 1360: {
! 1361: long k3 = subll(mulll(qp,saux), x1[3]);
! 1362: long k4 = subllx(hiremainder,k);
! 1363: while (!overflow && k4) { qp--; k3=subll(k3,saux); k4=subllx(k4,si); }
! 1364: }
! 1365: hiremainder=0;
! 1366: for (j=ly-1; j>1; j--)
! 1367: {
! 1368: x1[j] = subll(x1[j], addmul(qp,y[j]));
! 1369: hiremainder+=overflow;
! 1370: }
! 1371: if ((ulong)x1[1] < hiremainder)
! 1372: {
! 1373: overflow=0; qp--;
! 1374: for (j=ly-1; j>1; j--) x1[j] = addllx(x1[j],y[j]);
! 1375: }
! 1376: x1[1]=qp; x1++;
! 1377: }
! 1378:
! 1379: lq = lz+2;
! 1380: if (!z)
! 1381: {
! 1382: qd = (ulong*)av;
! 1383: xd = (ulong*)(x + lq);
! 1384: if (x[1]) { lz++; lq++; }
! 1385: while (lz--) *--qd = *--xd;
! 1386: *--qd = evalsigne(sy) | evallgefint(lq);
! 1387: *--qd = evaltyp(t_INT) | evallg(lq);
! 1388: avma = (long)qd; return (GEN)qd;
! 1389: }
! 1390:
! 1391: j=lq; while (j<lx && !x[j]) j++;
! 1392: if (z == ONLY_REM)
! 1393: {
! 1394: lz = lx-j; if (lz==0) { avma = av; return gzero; }
! 1395: rd = (ulong*)av; lr = lz+2;
! 1396: xd = (ulong*)(x + lx);
! 1397: if (!sh) while (lz--) *--rd = *--xd;
! 1398: else
! 1399: { /* shift remainder right by sh bits */
! 1400: const ulong shl = BITS_IN_LONG - sh;
! 1401: ulong l;
! 1402: xd--;
! 1403: while (--lz) /* fill r[3..] */
! 1404: {
! 1405: l = *xd >> sh;
! 1406: *--rd = l | (*--xd << shl);
! 1407: }
! 1408: l = *xd >> sh;
! 1409: if (l) *--rd = l; else lr--;
! 1410: }
! 1411: *--rd = evalsigne(sx) | evallgefint(lr);
! 1412: *--rd = evaltyp(t_INT) | evallg(lr);
! 1413: avma = (long)rd; return (GEN)rd;
! 1414: }
! 1415:
! 1416: lz = lx-j; lr = lz+2;
! 1417: if (lz)
! 1418: {
! 1419: xd = (ulong*)(x + lx);
! 1420: if (!sh)
! 1421: {
! 1422: rd = (ulong*)avma; r = new_chunk(lr);
! 1423: while (lz--) *--rd = *--xd;
! 1424: }
! 1425: else
! 1426: { /* shift remainder right by sh bits */
! 1427: const ulong shl = BITS_IN_LONG - sh;
! 1428: ulong l;
! 1429: rd = (ulong*)x; /* overwrite shifted y */
! 1430: xd--;
! 1431: while (--lz)
! 1432: {
! 1433: l = *xd >> sh;
! 1434: *--rd = l | (*--xd << shl);
! 1435: }
! 1436: l = *xd >> sh;
! 1437: if (l) *--rd = l; else lr--;
! 1438: }
! 1439: *--rd = evalsigne(sx) | evallgefint(lr);
! 1440: *--rd = evaltyp(t_INT) | evallg(lr);
! 1441: rd += lr;
! 1442: }
! 1443: qd = (ulong*)av;
! 1444: xd = (ulong*)(x + lq);
! 1445: if (x[1]) lq++;
! 1446: j = lq-2; while (j--) *--qd = *--xd;
! 1447: *--qd = evalsigne(sy) | evallgefint(lq);
! 1448: *--qd = evaltyp(t_INT) | evallg(lq);
! 1449: q = (GEN)qd;
! 1450: if (lr==2) *z = gzero;
! 1451: else
! 1452: {
! 1453: while (lr--) *--qd = *--rd;
! 1454: *z = (GEN)qd;
! 1455: }
! 1456: avma = (long)qd; return q;
! 1457: }
! 1458:
! 1459: /* assume y > x > 0. return y mod x */
! 1460: ulong
! 1461: mppgcd_resiu(GEN y, ulong x)
! 1462: {
! 1463: long i, ly = lgefint(y);
! 1464: LOCAL_HIREMAINDER;
! 1465:
! 1466: hiremainder = 0;
! 1467: for (i=2; i<ly; i++) (void)divll(y[i],x);
! 1468: return hiremainder;
! 1469: }
! 1470:
! 1471: /* Assume x>y>0, both of them odd. return x-y if x=y mod 4, x+y otherwise */
! 1472: void
! 1473: mppgcd_plus_minus(GEN x, GEN y, GEN res)
! 1474: {
! 1475: long av = avma;
! 1476: long lx = lgefint(x)-1;
! 1477: long ly = lgefint(y)-1, lt,m,i;
! 1478: GEN t;
! 1479:
! 1480: if ((x[lx]^y[ly]) & 3) /* x != y mod 4*/
! 1481: t = addiispec(x+2,y+2,lx-1,ly-1);
! 1482: else
! 1483: t = subiispec(x+2,y+2,lx-1,ly-1);
! 1484:
! 1485: lt = lgefint(t)-1; while (!t[lt]) lt--;
! 1486: m = vals(t[lt]); lt++;
! 1487: if (m == 0) /* 2^32 | t */
! 1488: {
! 1489: for (i = 2; i < lt; i++) res[i] = t[i];
! 1490: }
! 1491: else if (t[2] >> m)
! 1492: {
! 1493: shift_right(res,t, 2,lt, 0,m);
! 1494: }
! 1495: else
! 1496: {
! 1497: lt--; t++;
! 1498: shift_right(res,t, 2,lt, t[1],m);
! 1499: }
! 1500: res[1] = evalsigne(1)|evallgefint(lt);
! 1501: avma = av;
! 1502: }
! 1503:
! 1504: /* private version of mulss */
! 1505: ulong
! 1506: smulss(ulong x, ulong y, ulong *rem)
! 1507: {
! 1508: LOCAL_HIREMAINDER;
! 1509: x=mulll(x,y); *rem=hiremainder; return x;
! 1510: }
! 1511:
! 1512: #ifdef LONG_IS_64BIT
! 1513: # define DIVCONVI 7
! 1514: #else
! 1515: # define DIVCONVI 14
! 1516: #endif
! 1517:
! 1518: /* Conversion entier --> base 10^9. Assume x > 0 */
! 1519: GEN
! 1520: convi(GEN x)
! 1521: {
! 1522: ulong av=avma, lim;
! 1523: long lz;
! 1524: GEN z,p1;
! 1525:
! 1526: lz = 3 + ((lgefint(x)-2)*15)/DIVCONVI;
! 1527: z=new_chunk(lz); z[1] = -1; p1 = z+2;
! 1528: lim = stack_lim(av,1);
! 1529: for(;;)
! 1530: {
! 1531: x = divis(x,1000000000); *p1++ = hiremainder;
! 1532: if (!signe(x)) { avma=av; return p1; }
! 1533: if (low_stack(lim, stack_lim(av,1))) x = gerepileuptoint((long)z,x);
! 1534: }
! 1535: }
! 1536: #undef DIVCONVI
! 1537:
! 1538: /* Conversion partie fractionnaire --> base 10^9 (expo(x)<0) */
! 1539: GEN
! 1540: confrac(GEN x)
! 1541: {
! 1542: long lx=lg(x), ex = -expo(x)-1, d,m,ly,ey,lr,nbdec,i,j;
! 1543: GEN y,y1;
! 1544:
! 1545: ey = ex + bit_accuracy(lx);
! 1546: ly = 1 + ((ey-1) >> TWOPOTBITS_IN_LONG);
! 1547: d = ex >> TWOPOTBITS_IN_LONG;
! 1548: m = ex & (BITS_IN_LONG-1);
! 1549: y = new_chunk(ly); y1 = y + (d-2);
! 1550: while (d--) y[d]=0;
! 1551: if (!m)
! 1552: for (i=2; i<lx; i++) y1[i] = x[i];
! 1553: else
! 1554: { /* shift x left sh bits */
! 1555: const ulong sh=BITS_IN_LONG-m;
! 1556: ulong k=0, l;
! 1557: for (i=2; i<lx; i++)
! 1558: {
! 1559: l = x[i];
! 1560: y1[i] = (l>>m)|k;
! 1561: k = l<<sh;
! 1562: }
! 1563: y1[i] = k;
! 1564: }
! 1565: nbdec = (long) (ey*L2SL10)+1; lr=(nbdec+17)/9;
! 1566: y1=new_chunk(lr); *y1=nbdec;
! 1567: for (j=1; j<lr; j++)
! 1568: {
! 1569: LOCAL_HIREMAINDER;
! 1570: hiremainder=0;
! 1571: for (i=ly-1; i>=0; i--) y[i]=addmul(y[i],1000000000);
! 1572: y1[j]=hiremainder;
! 1573: }
! 1574: return y1;
! 1575: }
! 1576:
! 1577: GEN
! 1578: truedvmdii(GEN x, GEN y, GEN *z)
! 1579: {
! 1580: long av=avma,tetpil;
! 1581: GEN res, qu = dvmdii(x,y,&res);
! 1582: GEN *gptr[2];
! 1583:
! 1584: if (signe(res)>=0)
! 1585: {
! 1586: if (z == ONLY_REM) return gerepileuptoint(av,res);
! 1587: if (z) *z = res; else cgiv(res);
! 1588: return qu;
! 1589: }
! 1590:
! 1591: tetpil=avma;
! 1592: if (z == ONLY_REM)
! 1593: {
! 1594: res = subiispec(y+2,res+2, lgefint(y)-2,lgefint(res)-2);
! 1595: return gerepile(av,tetpil,res);
! 1596: }
! 1597: qu = addsi(-signe(y),qu);
! 1598: if (!z) return gerepile(av,tetpil,qu);
! 1599:
! 1600: *z = subiispec(y+2,res+2, lgefint(y)-2,lgefint(res)-2);
! 1601: gptr[0]=&qu; gptr[1]=z; gerepilemanysp(av,tetpil,gptr,2);
! 1602: return qu;
! 1603: }
! 1604:
! 1605: #if 0
! 1606: /* Exact integer division */
! 1607:
! 1608: static ulong
! 1609: invrev(ulong b)
! 1610: /* Find c such that 1=c*b mod B (where B = 2^32 or 2^64), assuming b odd,
! 1611: which is not checked */
! 1612: {
! 1613: int r;
! 1614: ulong x;
! 1615:
! 1616: r=b&7; x=(r==3 || r==5)? b+8: b; /* x=b^(-1) mod 2^4 */
! 1617: x=x*(2-b*x); x=x*(2-b*x); x=x*(2-b*x); /* x=b^(-1) mod 2^32 */
! 1618: #ifdef LONG_IS_64BIT
! 1619: x=x*(2-b*x); /* x=b^(-1) mod 2^64 */
! 1620: #endif
! 1621: return x;
! 1622: }
! 1623:
! 1624: #define divllrev(a,b) (((ulong)a)*invrev(b))
! 1625:
! 1626: /* 2-adic division */
! 1627: GEN
! 1628: diviirev(GEN x, GEN y, long a)
! 1629: /* Find z such that |x|=|y|*z mod B^a, where a<=lgefint(x)-2 */
! 1630: {
! 1631: long lx,lgx,ly,vy,av=avma,tetpil,i,j,ii;
! 1632: ulong binv,q;
! 1633: GEN z;
! 1634:
! 1635: if (!signe(y)) err(dvmer1);
! 1636: if (!signe(x)) return gzero;
! 1637: /* make y odd */
! 1638: vy=vali(y);
! 1639: if (vy)
! 1640: {
! 1641: if (vali(x)<vy) err(talker,"impossible division in diviirev");
! 1642: y=shifti(y,-vy); x=shifti(x,-vy); a-=(vy>>TWOPOTBITS_IN_LONG);
! 1643: }
! 1644: else x=icopy(x); /* necessary because we destroy x */
! 1645: /* improve the above by touching only a words */
! 1646: if (a<=0) {avma=a;return gzero;}
! 1647: /* now y is odd */
! 1648: lx=a+2; ly=lgefint(y); lgx=lgefint(x);
! 1649: if (lx>lgx) err(talker,"3rd parameter too large in diviirev");
! 1650: binv=invrev(y[ly-1]);
! 1651: z=cgeti(lx);
! 1652: for (ii=lgx-1,i=lx-1; i>=2; i--,ii--)
! 1653: {
! 1654: long limj;
! 1655: LOCAL_HIREMAINDER;
! 1656: LOCAL_OVERFLOW;
! 1657:
! 1658: z[i]=q=binv*((ulong)x[ii]); /* this is the i-th quotient */
! 1659: limj=max(lgx-a,3+ii-ly);
! 1660: mulll(q,y[ly-1]); overflow=0;
! 1661: for (j=ii-1; j>=limj; j--)
! 1662: x[j]=subllx(x[j],addmul(q,y[j+ly-ii-1]));
! 1663: }
! 1664: tetpil=avma; i=2; while((i<lx)&&(!z[i])) i++;
! 1665: if (i==lx) {avma=av; return gzero;}
! 1666: y=cgeti(lx-i+2); y[1]=evalsigne(1)|evallgefint(lx-i+2); j=2;
! 1667: for ( ; i<lx; i++) y[j++]=z[i];
! 1668: return gerepile(av,tetpil,y);
! 1669: }
! 1670:
! 1671: GEN
! 1672: diviiexactfullrev(GEN x, GEN y)
! 1673: /* Find z such that x=y*z knowing that y divides x */
! 1674: {
! 1675: long lx,lz,ly,vy,av=avma,tetpil,i,j,ii,sx=signe(x),sy=signe(y);
! 1676: ulong binv,q;
! 1677: GEN z;
! 1678:
! 1679: if (!sy) err(dvmer1);
! 1680: if (!sx) return gzero;
! 1681: /* make y odd */
! 1682: vy=vali(y);
! 1683: if (vy)
! 1684: {
! 1685: if (vali(x)<vy) err(talker,"impossible division in diviirev");
! 1686: y=shifti(y,-vy); x=shifti(x,-vy);
! 1687: }
! 1688: else x=icopy(x); /* necessary because we destroy x */
! 1689: /* now y is odd */
! 1690: ly=lgefint(y); lx=lgefint(x);
! 1691: if (ly>lx) err(talker,"impossible division in diviirev");
! 1692: binv=invrev(y[ly-1]);
! 1693: i=2; while (i<ly && y[i]==x[i]) i++;
! 1694: lz=(i==ly || y[i]<x[i]) ? lx-ly+3 : lx-ly+2;
! 1695: z=cgeti(lz);
! 1696: for (ii=lx-1,i=lz-1; i>=2; i--,ii--)
! 1697: {
! 1698: long limj;
! 1699: LOCAL_HIREMAINDER;
! 1700: LOCAL_OVERFLOW;
! 1701:
! 1702: z[i]=q=binv*((ulong)x[ii]); /* this is the i-th quotient */
! 1703: limj=max(lx-lz+2,3+ii-ly);
! 1704: mulll(q,y[ly-1]); overflow=0;
! 1705: for (j=ii-1; j>=limj; j--)
! 1706: x[j]=subllx(x[j],addmul(q,y[j+ly-ii-1]));
! 1707: }
! 1708: tetpil=avma; i=2; while((i<lz)&&(!z[i])) i++;
! 1709: if (i==lz) err(talker,"bug in diviiexact");
! 1710: y=cgeti(lz-i+2); y[1]=evalsigne(sx*sy) | evallgefint(lz-i+2); j=2;
! 1711: for ( ; i<lz; i++) y[j++]=z[i];
! 1712: return gerepile(av,tetpil,y);
! 1713: }
! 1714:
! 1715: GEN
! 1716: diviiexact2(GEN x, GEN y)
! 1717: /* Find z such that x=y*z assuming y divides x (which is not checked) */
! 1718: {
! 1719: long sx=signe(x),sy=signe(y),av=avma,tetpil,lyinv,lpr,a,lx,ly,lz,lzs,lp1;
! 1720: long i,j,vy;
! 1721: ulong aux;
! 1722: GEN yinv,p1,z,xinit,yinit;
! 1723:
! 1724: if (!sy) err(dvmer1);
! 1725: if (!sx) return gzero;
! 1726: xinit=x; yinit=y;
! 1727: setsigne(y,1);setsigne(x,1);
! 1728: /* make y odd */
! 1729: vy=vali(y);
! 1730: if (vy)
! 1731: {
! 1732: if (vali(x)<vy) err(talker,"impossible division in diviirev");
! 1733: y=shifti(y,-vy); x=shifti(x,-vy);
! 1734: }
! 1735: /* now y is odd */
! 1736: ly=lgefint(y); lx=lgefint(x);
! 1737: if (lx<ly) err(talker,"not an exact division in diviiexact");
! 1738: a=lx-ly+1; lz=a+2;
! 1739: aux=invrev(y[ly-1]);
! 1740: if (aux & HIGHBIT) {yinv=stoi(aux^HIGHBIT);yinv[2]|=HIGHBIT;}
! 1741: else yinv=stoi(aux); /* inverse of y mod 2^32 (or 2^64) */
! 1742: lpr=1; /* current accuracy */
! 1743: while(lpr<a)
! 1744: {
! 1745: long lycut;
! 1746: GEN ycut;
! 1747:
! 1748: lyinv=lgefint(yinv);
! 1749: lycut=min(2*lpr+2,ly);
! 1750: ycut=cgeti(lycut); ycut[1]=evalsigne(1) | evallgefint(lycut);
! 1751: for(i=2; i<lycut; i++) ycut[i]=y[ly+i-lycut];
! 1752: p1=mulii(yinv,ycut); lp1=lgefint(p1);
! 1753: if (lp1>lpr+2)
! 1754: {
! 1755: long lp1cut,lynew,lp2;
! 1756: GEN p1cut,p2,ynew;
! 1757: LOCAL_OVERFLOW;
! 1758:
! 1759: lp1cut=min(lp1-lpr,lpr+2);
! 1760: p1cut=cgeti(lp1cut); p1cut[1]=evalsigne(1) | evallgefint(lp1cut);
! 1761: overflow=0;
! 1762: for (i=lp1cut-1; i>=2; i--) p1cut[i]=subllx(0,p1[i+lp1-lpr-lp1cut]);
! 1763: p2=mulii(p1cut,yinv); lp2=lgefint(p2);
! 1764: lynew=(lp2<=lpr+2) ? lpr+lp2 : 2*lpr+2;
! 1765: ynew=cgeti(lynew); ynew[1]=evalsigne(1) | evallgefint(lynew);
! 1766: for (i=lynew-1; i>=lynew-lpr; i--) ynew[i]=yinv[i+lyinv-lynew];
! 1767: for (i=lynew-lpr-1; i>=2; i--) ynew[i]=p2[i+lp2+lpr-lynew];
! 1768: yinv=ynew;
! 1769: }
! 1770: lpr<<=1;
! 1771: }
! 1772: lyinv=lgefint(yinv); lzs=min(lz,lyinv);
! 1773: z=cgeti(lzs); z[1]=evalsigne(1) | evallgefint(lzs);
! 1774: for(i=2; i<lzs; i++) z[i]=yinv[i+lyinv-lzs];
! 1775: p1=mulii(z,x); lp1=lgefint(p1); lzs=min(lz,lp1);
! 1776: z=cgeti(lzs);
! 1777: for (i=2; i<lzs; i++) z[i]=p1[i+lp1-lzs];
! 1778: i=2; while (i<lzs && !z[i]) i++;
! 1779: if (i==lzs) err(talker,"bug in diviiexact");
! 1780: tetpil=avma; p1=cgeti(lzs-i+2);
! 1781: p1[1]=evalsigne(sx*sy) | evallgefint(lzs-i+2);
! 1782: for (j=2; j<=lzs-i+1; j++) p1[j]=z[j+i-2];
! 1783: setsigne(xinit,sx);setsigne(yinit,sy);
! 1784: return gerepile(av,tetpil,p1);
! 1785: }
! 1786: #endif
! 1787:
! 1788: long
! 1789: smodsi(long x, GEN y)
! 1790: {
! 1791: if (x<0) err(talker,"negative small integer in smodsi");
! 1792: (void)divsi(x,y); return hiremainder;
! 1793: }
! 1794:
! 1795: /* x and y are integers. Return 1 if |x| == |y|, 0 otherwise */
! 1796: int
! 1797: absi_equal(GEN x, GEN y)
! 1798: {
! 1799: long lx,i;
! 1800:
! 1801: if (!signe(x)) return !signe(y);
! 1802: if (!signe(y)) return 0;
! 1803:
! 1804: lx=lgefint(x); if (lx != lgefint(y)) return 0;
! 1805: i=2; while (i<lx && x[i]==y[i]) i++;
! 1806: return (i==lx);
! 1807: }
! 1808:
! 1809: /* x and y are integers. Return sign(|x| - |y|) */
! 1810: int
! 1811: absi_cmp(GEN x, GEN y)
! 1812: {
! 1813: long lx,ly,i;
! 1814:
! 1815: if (!signe(x)) return signe(y)? -1: 0;
! 1816: if (!signe(y)) return 1;
! 1817:
! 1818: lx=lgefint(x); ly=lgefint(y);
! 1819: if (lx>ly) return 1;
! 1820: if (lx<ly) return -1;
! 1821: i=2; while (i<lx && x[i]==y[i]) i++;
! 1822: if (i==lx) return 0;
! 1823: return ((ulong)x[i] > (ulong)y[i])? 1: -1;
! 1824: }
! 1825:
! 1826: /* x and y are reals. Return sign(|x| - |y|) */
! 1827: int
! 1828: absr_cmp(GEN x, GEN y)
! 1829: {
! 1830: long ex,ey,lx,ly,lz,i;
! 1831:
! 1832: if (!signe(x)) return signe(y)? -1: 0;
! 1833: if (!signe(y)) return 1;
! 1834:
! 1835: ex=expo(x); ey=expo(y);
! 1836: if (ex>ey) return 1;
! 1837: if (ex<ey) return -1;
! 1838:
! 1839: lx=lg(x); ly=lg(y); lz = (lx<ly)?lx:ly;
! 1840: i=2; while (i<lz && x[i]==y[i]) i++;
! 1841: if (i<lz) return ((ulong)x[i] > (ulong)y[i])? 1: -1;
! 1842: if (lx>=ly)
! 1843: {
! 1844: while (i<lx && !x[i]) i++;
! 1845: return (i==lx)? 0: 1;
! 1846: }
! 1847: while (i<ly && !y[i]) i++;
! 1848: return (i==ly)? 0: -1;
! 1849: }
! 1850:
! 1851: /********************************************************************/
! 1852: /** **/
! 1853: /** INTEGER MULTIPLICATION (KARATSUBA) **/
! 1854: /** **/
! 1855: /********************************************************************/
! 1856:
! 1857: #if 0 /* for tunings */
! 1858: long SQRI_LIMIT = 12;
! 1859: long MULII_LIMIT = 16;
! 1860:
! 1861: void
! 1862: setsqi(long a) { SQRI_LIMIT=a; }
! 1863: void
! 1864: setmuli(long a) { MULII_LIMIT=a; }
! 1865:
! 1866: GEN
! 1867: speci(GEN x, long nx)
! 1868: {
! 1869: GEN z = cgeti(nx+2);
! 1870: long i;
! 1871: for (i=0; i<nx; i++) z[i+2] = x[i];
! 1872: z[1]=evalsigne(1)|evallgefint(nx+2);
! 1873: return z;
! 1874: }
! 1875: #else
! 1876: # define SQRI_LIMIT 12
! 1877: # define MULII_LIMIT 16
! 1878: #endif
! 1879:
! 1880: /* nx >= ny = num. of digits of x, y (not GEN, see mulii) */
! 1881: #ifdef INLINE
! 1882: INLINE
! 1883: #endif
! 1884: GEN
! 1885: muliispec(GEN x, GEN y, long nx, long ny)
! 1886: {
! 1887: GEN z2e,z2d,yd,xd,ye,zd;
! 1888: long p1,lz;
! 1889: LOCAL_HIREMAINDER;
! 1890:
! 1891: if (!ny) return gzero;
! 1892: zd = (GEN)avma; lz = nx+ny+2;
! 1893: (void)new_chunk(lz);
! 1894: xd = x + nx;
! 1895: yd = y + ny;
! 1896: ye = yd; p1 = *--xd;
! 1897:
! 1898: *--zd = mulll(p1, *--yd); z2e = zd;
! 1899: while (yd > y) *--zd = addmul(p1, *--yd);
! 1900: *--zd = hiremainder;
! 1901:
! 1902: while (xd > x)
! 1903: {
! 1904: LOCAL_OVERFLOW;
! 1905: yd = ye; p1 = *--xd;
! 1906:
! 1907: z2d = --z2e;
! 1908: *z2d = addll(mulll(p1, *--yd), *z2d); z2d--;
! 1909: while (yd > y)
! 1910: {
! 1911: hiremainder += overflow;
! 1912: *z2d = addll(addmul(p1, *--yd), *z2d); z2d--;
! 1913: }
! 1914: *--zd = hiremainder + overflow;
! 1915: }
! 1916: if (*zd == 0) { zd++; lz--; } /* normalize */
! 1917: *--zd = evalsigne(1) | evallgefint(lz);
! 1918: *--zd = evaltyp(t_INT) | evallg(lz);
! 1919: avma=(long)zd; return zd;
! 1920: }
! 1921:
! 1922: /* return (x shifted left d words) + y. Assume d > 0, x > 0 and y >= 0 */
! 1923: static GEN
! 1924: addshiftw(GEN x, GEN y, long d)
! 1925: {
! 1926: GEN z,z0,y0,yd, zd = (GEN)avma;
! 1927: long a,lz,ly = lgefint(y);
! 1928:
! 1929: z0 = new_chunk(d);
! 1930: a = ly-2; yd = y+ly;
! 1931: if (a >= d)
! 1932: {
! 1933: y0 = yd-d; while (yd > y0) *--zd = *--yd; /* copy last d words of y */
! 1934: a -= d;
! 1935: if (a)
! 1936: z = addiispec(x+2, y+2, lgefint(x)-2, a);
! 1937: else
! 1938: z = icopy(x);
! 1939: }
! 1940: else
! 1941: {
! 1942: y0 = yd-a; while (yd > y0) *--zd = *--yd; /* copy last a words of y */
! 1943: while (zd >= z0) *--zd = 0; /* complete with 0s */
! 1944: z = icopy(x);
! 1945: }
! 1946: lz = lgefint(z)+d;
! 1947: z[1] = evalsigne(1) | evallgefint(lz);
! 1948: z[0] = evaltyp(t_INT) | evallg(lz); return z;
! 1949: }
! 1950:
! 1951: /* Fast product (Karatsuba) of integers a,b. These are not real GENs, a+2
! 1952: * and b+2 were sent instead. na, nb = number of digits of a, b.
! 1953: * c,c0,c1,c2 are genuine GENs.
! 1954: */
! 1955: static GEN
! 1956: quickmulii(GEN a, GEN b, long na, long nb)
! 1957: {
! 1958: GEN a0,c,c0;
! 1959: long av,n0,n0a,i;
! 1960:
! 1961: if (na < nb) swapspec(a,b, na,nb);
! 1962: if (nb == 1) return mulsispec(*b, a, na);
! 1963: if (nb == 0) return gzero;
! 1964: if (nb < MULII_LIMIT) return muliispec(a,b,na,nb);
! 1965: i=(na>>1); n0=na-i; na=i;
! 1966: av=avma; a0=a+na; n0a=n0;
! 1967: while (!*a0 && n0a) { a0++; n0a--; }
! 1968:
! 1969: if (n0a && nb > n0)
! 1970: { /* nb <= na <= n0 */
! 1971: GEN b0,c1,c2;
! 1972: long n0b;
! 1973:
! 1974: nb -= n0;
! 1975: c = quickmulii(a,b,na,nb);
! 1976: b0 = b+nb; n0b = n0;
! 1977: while (!*b0 && n0b) { b0++; n0b--; }
! 1978: if (n0b)
! 1979: {
! 1980: c0 = quickmulii(a0,b0, n0a,n0b);
! 1981:
! 1982: c2 = addiispec(a0,a, n0a,na);
! 1983: c1 = addiispec(b0,b, n0b,nb);
! 1984: c1 = quickmulii(c1+2,c2+2, lgefint(c1)-2,lgefint(c2)-2);
! 1985: c2 = addiispec(c0+2, c+2, lgefint(c0)-2,lgefint(c) -2);
! 1986:
! 1987: c1 = subiispec(c1+2,c2+2, lgefint(c1)-2,lgefint(c2)-2);
! 1988: }
! 1989: else
! 1990: {
! 1991: c0 = gzero;
! 1992: c1 = quickmulii(a0,b, n0a,nb);
! 1993: }
! 1994: c = addshiftw(c,c1, n0);
! 1995: }
! 1996: else
! 1997: {
! 1998: c = quickmulii(a,b,na,nb);
! 1999: c0 = quickmulii(a0,b,n0a,nb);
! 2000: }
! 2001: return gerepileuptoint(av, addshiftw(c,c0, n0));
! 2002: }
! 2003:
! 2004: /* actual operations will take place on a+2 and b+2: we strip the codewords */
! 2005: GEN
! 2006: mulii(GEN a,GEN b)
! 2007: {
! 2008: long sa,sb;
! 2009: GEN z;
! 2010:
! 2011: sa=signe(a); if (!sa) return gzero;
! 2012: sb=signe(b); if (!sb) return gzero;
! 2013: if (sb<0) sa = -sa;
! 2014: z = quickmulii(a+2,b+2, lgefint(a)-2,lgefint(b)-2);
! 2015: setsigne(z,sa); return z;
! 2016: }
! 2017:
! 2018: static GEN
! 2019: quicksqri(GEN a, long na)
! 2020: {
! 2021: GEN a0,c,c0,c1;
! 2022: long av,n0,n0a,i;
! 2023:
! 2024: if (na<SQRI_LIMIT) return muliispec(a,a,na,na);
! 2025: i=(na>>1); n0=na-i; na=i;
! 2026: av=avma; a0=a+na; n0a=n0;
! 2027: while (!*a0 && n0a) { a0++; n0a--; }
! 2028: c = quicksqri(a,na);
! 2029: if (n0a)
! 2030: {
! 2031: c0 = quicksqri(a0,n0a);
! 2032: c1 = shifti(quickmulii(a0,a, n0a,na),1);
! 2033: c = addshiftw(c,c1, n0);
! 2034: c = addshiftw(c,c0, n0);
! 2035: }
! 2036: else
! 2037: c = addshiftw(c,gzero,n0<<1);
! 2038: return gerepileuptoint(av, c);
! 2039: }
! 2040:
! 2041: GEN
! 2042: sqri(GEN a) { return quicksqri(a+2, lgefint(a)-2); }
! 2043:
! 2044: #define MULRR_LIMIT 32
! 2045: #define MULRR2_LIMIT 32
! 2046:
! 2047: #if 0
! 2048: GEN
! 2049: karamulrr1(GEN x, GEN y)
! 2050: {
! 2051: long sx,sy;
! 2052: long i,i1,i2,lx=lg(x),ly=lg(y),e,flag,garde;
! 2053: long lz2,lz3,lz4;
! 2054: GEN z,lo1,lo2,hi;
! 2055:
! 2056: /* ensure that lg(y) >= lg(x) */
! 2057: if (lx>ly) { lx=ly; z=x; x=y; y=z; flag=1; } else flag = (lx!=ly);
! 2058: if (lx < MULRR_LIMIT) return mulrr(x,y);
! 2059: sx=signe(x); sy=signe(y);
! 2060: e = evalexpo(expo(x)+expo(y)); if (e & ~EXPOBITS) err(muler4);
! 2061: if (!sx || !sy) { z=cgetr(3); z[2]=0; z[1]=e; return z; }
! 2062: if (sy<0) sx = -sx;
! 2063: ly=lx+flag; z=cgetr(lx);
! 2064: lz2 = (lx>>1); lz3 = lx-lz2;
! 2065: x += 2; lx -= 2;
! 2066: y += 2; ly -= 2;
! 2067: hi = quickmulii(x,y,lz2,lz2);
! 2068: i1=lz2; while (i1<lx && !x[i1]) i1++;
! 2069: lo1 = quickmulii(y,x+i1,lz2,lx-i1);
! 2070: i2=lz2; while (i2<ly && !y[i2]) i2++;
! 2071: lo2 = quickmulii(x,y+i2,lz2,ly-i2);
! 2072:
! 2073: if (signe(lo1))
! 2074: {
! 2075: if (flag) { lo2 = addshiftw(lo1,lo2,1); lz3++; } else lo2=addii(lo1,lo2);
! 2076: }
! 2077: lz4=lgefint(lo2)-lz3;
! 2078: if (lz4>0) hi = addiispec(hi+2,lo2+2, lgefint(hi)-2,lz4);
! 2079: if (hi[2] < 0)
! 2080: {
! 2081: e++; garde=hi[lx+2];
! 2082: for (i=lx+1; i>=2 ; i--) z[i]=hi[i];
! 2083: }
! 2084: else
! 2085: {
! 2086: garde = (hi[lx+2] << 1);
! 2087: shift_left(z,hi,2,lx+1, garde, 1);
! 2088: }
! 2089: z[1]=evalsigne(sx) | e;
! 2090: if (garde < 0)
! 2091: { /* round to nearest */
! 2092: i=lx+2; do z[--i]++; while (z[i]==0);
! 2093: if (i==1) z[2]=HIGHBIT;
! 2094: }
! 2095: avma=(long)z; return z;
! 2096: }
! 2097:
! 2098: GEN
! 2099: karamulrr2(GEN x, GEN y)
! 2100: {
! 2101: long sx,sy,i,lx=lg(x),ly=lg(y),e,flag,garde;
! 2102: GEN z,hi;
! 2103:
! 2104: if (lx>ly) { lx=ly; z=x; x=y; y=z; flag=1; } else flag = (lx!=ly);
! 2105: if (lx < MULRR2_LIMIT) return mulrr(x,y);
! 2106: ly=lx+flag; sx=signe(x); sy=signe(y);
! 2107: e = evalexpo(expo(x)+expo(y)); if (e & ~EXPOBITS) err(muler4);
! 2108: if (!sx || !sy) { z=cgetr(3); z[2]=0; z[1]=e; return z; }
! 2109: if (sy<0) sx = -sx;
! 2110: z=cgetr(lx);
! 2111: hi=quickmulii(y+2,x+2,ly-2,lx-2);
! 2112: if (hi[2] < 0)
! 2113: {
! 2114: e++; garde=hi[lx];
! 2115: for (i=2; i<lx ; i++) z[i]=hi[i];
! 2116: }
! 2117: else
! 2118: {
! 2119: garde = (hi[lx] << 1);
! 2120: shift_left(z,hi,2,lx-1, garde, 1);
! 2121: }
! 2122: z[1]=evalsigne(sx) | e;
! 2123: if (garde < 0)
! 2124: { /* round to nearest */
! 2125: i=lx; do z[--i]++; while (z[i]==0);
! 2126: if (i==1) z[2]=HIGHBIT;
! 2127: }
! 2128: avma=(long)z; return z;
! 2129: }
! 2130:
! 2131: GEN
! 2132: karamulrr(GEN x, GEN y, long flag)
! 2133: {
! 2134: switch(flag)
! 2135: {
! 2136: case 1: return karamulrr1(x,y);
! 2137: case 2: return karamulrr2(x,y);
! 2138: default: err(flagerr,"karamulrr");
! 2139: }
! 2140: return NULL; /* not reached */
! 2141: }
! 2142:
! 2143: GEN
! 2144: karamulir(GEN x, GEN y, long flag)
! 2145: {
! 2146: long sx=signe(x),lz,i;
! 2147: GEN z,temp,z1;
! 2148:
! 2149: if (!sx) return gzero;
! 2150: lz=lg(y); z=cgetr(lz);
! 2151: temp=cgetr(lz+1); affir(x,temp);
! 2152: z1=karamulrr(temp,y,flag);
! 2153: for (i=1; i<lz; i++) z[i]=z1[i];
! 2154: avma=(long)z; return z;
! 2155: }
! 2156: #endif
! 2157:
! 2158: #ifdef LONG_IS_64BIT
! 2159:
! 2160: #if PARI_BYTE_ORDER == LITTLE_ENDIAN_64 || PARI_BYTE_ORDER == BIG_ENDIAN_64
! 2161: #else
! 2162: error... unknown machine
! 2163: #endif
! 2164:
! 2165: GEN
! 2166: dbltor(double x)
! 2167: {
! 2168: GEN z = cgetr(3);
! 2169: long e;
! 2170: union { double f; ulong i; } fi;
! 2171: const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
! 2172: const int exp_mid = 0x3ff;/* exponent bias */
! 2173: const int expo_len = 11; /* number of bits of exponent */
! 2174: LOCAL_HIREMAINDER;
! 2175:
! 2176: if (x==0) { z[1]=evalexpo(-308); z[2]=0; return z; }
! 2177: fi.f = x;
! 2178: e = evalexpo(((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid);
! 2179: z[1] = e | evalsigne(x<0? -1: 1);
! 2180: z[2] = (fi.i << expo_len) | HIGHBIT;
! 2181: return z;
! 2182: }
! 2183:
! 2184: double
! 2185: rtodbl(GEN x)
! 2186: {
! 2187: long ex,s=signe(x);
! 2188: ulong a;
! 2189: union { double f; ulong i; } fi;
! 2190: const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
! 2191: const int exp_mid = 0x3ff;/* exponent bias */
! 2192: const int expo_len = 11; /* number of bits of exponent */
! 2193: LOCAL_HIREMAINDER;
! 2194:
! 2195: if (typ(x)==t_INT && !s) return 0.0;
! 2196: if (typ(x)!=t_REAL) err(typeer,"rtodbl");
! 2197: if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
! 2198:
! 2199: /* start by rounding to closest */
! 2200: a = (x[2] & (HIGHBIT-1)) + 0x400;
! 2201: if (a & HIGHBIT) { ex++; a=0; }
! 2202: if (ex >= exp_mid) err(rtodber);
! 2203: fi.i = ((ex + exp_mid) << mant_len) | (a >> expo_len);
! 2204: if (s<0) fi.i |= HIGHBIT;
! 2205: return fi.f;
! 2206: }
! 2207:
! 2208: #else
! 2209:
! 2210: #if PARI_BYTE_ORDER == LITTLE_ENDIAN
! 2211: # define INDEX0 1
! 2212: # define INDEX1 0
! 2213: #elif PARI_BYTE_ORDER == BIG_ENDIAN
! 2214: # define INDEX0 0
! 2215: # define INDEX1 1
! 2216: #else
! 2217: error... unknown machine
! 2218: #endif
! 2219:
! 2220: GEN
! 2221: dbltor(double x)
! 2222: {
! 2223: GEN z;
! 2224: long e;
! 2225: union { double f; ulong i[2]; } fi;
! 2226: const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
! 2227: const int exp_mid = 0x3ff;/* exponent bias */
! 2228: const int shift = mant_len-32;
! 2229: const int expo_len = 11; /* number of bits of exponent */
! 2230:
! 2231: if (x==0) { z=cgetr(3); z[1]=evalexpo(-308); z[2]=0; return z; }
! 2232: fi.f = x; z=cgetr(4);
! 2233: {
! 2234: const ulong a = fi.i[INDEX0];
! 2235: const ulong b = fi.i[INDEX1];
! 2236: e = evalexpo(((a & (HIGHBIT-1)) >> shift) - exp_mid);
! 2237: z[1] = e | evalsigne(x<0? -1: 1);
! 2238: z[3] = b << expo_len;
! 2239: z[2] = HIGHBIT | b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
! 2240: }
! 2241: return z;
! 2242: }
! 2243:
! 2244: double
! 2245: rtodbl(GEN x)
! 2246: {
! 2247: long ex,s=signe(x),lx=lg(x);
! 2248: ulong a,b,k;
! 2249: union { double f; ulong i[2]; } fi;
! 2250: const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
! 2251: const int exp_mid = 0x3ff;/* exponent bias */
! 2252: const int shift = mant_len-32;
! 2253: const int expo_len = 11; /* number of bits of exponent */
! 2254:
! 2255: if (typ(x)==t_INT && !s) return 0.0;
! 2256: if (typ(x)!=t_REAL) err(typeer,"rtodbl");
! 2257: if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
! 2258:
! 2259: /* start by rounding to closest */
! 2260: a = x[2] & (HIGHBIT-1);
! 2261: if (lx > 3)
! 2262: {
! 2263: b = x[3] + 0x400UL; if (b < 0x400UL) a++;
! 2264: if (a & HIGHBIT) { ex++; a=0; }
! 2265: }
! 2266: else b = 0;
! 2267: if (ex > exp_mid) err(rtodber);
! 2268: ex += exp_mid;
! 2269: k = (a >> expo_len) | (ex << shift);
! 2270: if (s<0) k |= HIGHBIT;
! 2271: fi.i[INDEX0] = k;
! 2272: fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len);
! 2273: return fi.f;
! 2274: }
! 2275: #endif
! 2276:
! 2277: /* Old cgiv without reference count (which was not used anyway)
! 2278: * Should be a macro.
! 2279: */
! 2280: void
! 2281: cgiv(GEN x)
! 2282: {
! 2283: if (x == (GEN) avma)
! 2284: avma = (long) (x+lg(x));
! 2285: }
! 2286:
! 2287: /********************************************************************/
! 2288: /** **/
! 2289: /** INTEGER EXTENDED GCD (AND INVMOD) **/
! 2290: /** **/
! 2291: /********************************************************************/
! 2292:
! 2293: /* GN 1998Oct25, originally developed in January 1998 under 2.0.4.alpha,
! 2294: * in the context of trying to improve elliptic curve cryptosystem attacking
! 2295: * algorithms.
! 2296: *
! 2297: * Two basic ideas - (1) avoid many integer divisions, especially when the
! 2298: * quotient is 1 (which happens more than 40% of the time). (2) Use Lehmer's
! 2299: * trick as modified by Jebelean of extracting a couple of words' worth of
! 2300: * leading bits from both operands, and compute partial quotients from them
! 2301: * as long as we can be sure of their values. The Jebelean modifications
! 2302: * consist in reliable inequalities from which we can decide fast whether
! 2303: * to carry on or to return to the outer loop, and in re-shifting after the
! 2304: * first word's worth of bits has been used up. All of this is described
! 2305: * in R. Lercier's these [pp148-153 & 163f.], except his outer loop isn't
! 2306: * quite right (the catch-up divisions needed when one partial quotient is
! 2307: * larger than a word are missing).
! 2308: *
! 2309: * The API is mpxgcd() below; the single-word routines xgcduu and xxgcduu
! 2310: * may be called directly if desired; lgcdii() probably doesn't make much
! 2311: * sense out of context.
! 2312: *
! 2313: * The whole lot is a factor 6 .. 8 faster on word-sized operands, and asym-
! 2314: * ptotically about a factor 2.5 .. 3, depending on processor architecture,
! 2315: * than the naive continued-division code. Unfortunately, thanks to the
! 2316: * unrolled loops and all, the code is a bit lengthy.
! 2317: */
! 2318:
! 2319: /*==================================
! 2320: * xgcduu(d,d1,f,v,v1,s)
! 2321: * xxgcduu(d,d1,f,u,u1,v,v1,s)
! 2322: *==================================*/
! 2323: /*
! 2324: * Fast `final' extended gcd algorithm, acting on two ulongs. Ideally this
! 2325: * should be replaced with assembler versions wherever possible. The present
! 2326: * code essentially does `subtract, compare, and possibly divide' at each step,
! 2327: * which is reasonable when hardware division (a) exists, (b) is a bit slowish
! 2328: * and (c) does not depend a lot on the operand values (as on i486). When
! 2329: * wordsize division is in fact an assembler routine based on subtraction,
! 2330: * this strategy may not be the most efficient one.
! 2331: *
! 2332: * xxgcduu() should be called with d > d1 > 0, returns gcd(d,d1), and assigns
! 2333: * the usual signless cont.frac. recurrence matrix to [u, u1; v, v1] (i.e.,
! 2334: * the product of all the [0, 1; 1 q_j] where the leftmost factor arises from
! 2335: * the quotient of the first division step), and the information about the
! 2336: * implied signs to s (-1 when an odd number of divisions has been done,
! 2337: * 1 otherwise). xgcduu() is exactly the same except that u,u1 are not com-
! 2338: * puted (and not returned, of course).
! 2339: *
! 2340: * The input flag f should be set to 1 if we know in advance that gcd(d,d1)==1
! 2341: * (so we can stop the chain division one step early: as soon as the remainder
! 2342: * equals 1). Use this when you intend to use only what would be v, or only
! 2343: * what would be u and v, after that final division step, but not u1 and v1.
! 2344: * With the flag in force and thus without that final step, the interesting
! 2345: * quantity/ies will still sit in [u1 and] v1, of course.
! 2346: *
! 2347: * For computing the inverse of a single-word INTMOD known to exist, pass f=1
! 2348: * to xgcduu(), and obtain the result from s and v1. (The routine does the
! 2349: * right thing when d1==1 already.) For finishing a multiword modinv known
! 2350: * to exist, pass f=1 to xxgcduu(), and multiply the returned matrix (with
! 2351: * properly adjusted signs) onto the values v' and v1' previously obtained
! 2352: * from the multiword division steps. Actually, just take the scalar product
! 2353: * of [v',v1'] with [u1,-v1], and change the sign if s==-1. (If the final
! 2354: * step had been carried out, it would be [-u,v], and s would also change.)
! 2355: * For reducing a rational number to lowest terms, pass f=0 to xgcduu().
! 2356: * Finally, f=0 with xxgcduu() is useful for Bezout computations.
! 2357: * [Harrumph. In the above prescription, the sign turns out to be precisely
! 2358: * wrong.]
! 2359: * (It is safe for inversemodulo() to call xgcduu() with f=1, because f&1
! 2360: * doesn't make a difference when gcd(d,d1)>1. The speedup is negligible.)
! 2361: *
! 2362: * In principle, when gcd(d,d1) is known to be 1, it is straightforward to
! 2363: * recover the final u,u1 given only v,v1 and s. However, it probably isn't
! 2364: * worthwhile, as it trades a few multiplications for a division.
! 2365: *
! 2366: * Note that these routines do not know and do not need to know about the
! 2367: * PARI stack.
! 2368: */
! 2369:
! 2370: ulong
! 2371: xgcduu(ulong d, ulong d1, int f, ulong* v, ulong* v1, long *s)
! 2372: {
! 2373: ulong xv,xv1, xs, q,res;
! 2374: LOCAL_HIREMAINDER;
! 2375:
! 2376: /* The above blurb contained a lie. The main loop always stops when d1
! 2377: * has become equal to 1. If (d1 == 1 && !(f&1)) after the loop, we do
! 2378: * the final `division' of d by 1 `by hand' as it were.
! 2379: *
! 2380: * The loop has already been unrolled once. Aggressive optimization could
! 2381: * well lead to a totally unrolled assembler version...
! 2382: *
! 2383: * On modern x86 architectures, this loop is a pig anyway. The division
! 2384: * instruction always puts its result into the same pair of registers, and
! 2385: * we always want to use one of them straight away, so pipeline performance
! 2386: * will suck big time. An assembler version should probably do a first loop
! 2387: * computing and storing all the quotients -- their number is bounded in
! 2388: * advance -- and then assembling the matrix in a second pass. On other
! 2389: * architectures where we can cycle through four or so groups of registers
! 2390: * and exploit a fast ALU result-to-operand feedback path, this is much less
! 2391: * of an issue. (Intel sucks. See http://www.x86.org/ ...)
! 2392: */
! 2393: xs = res = 0;
! 2394: xv = 0UL; xv1 = 1UL;
! 2395: while (d1 > 1UL)
! 2396: {
! 2397: d -= d1; /* no need to use subll */
! 2398: if (d >= d1)
! 2399: {
! 2400: hiremainder = 0; q = 1 + divll(d,d1); d = hiremainder;
! 2401: xv += q * xv1;
! 2402: }
! 2403: else
! 2404: xv += xv1;
! 2405: /* possible loop exit */
! 2406: if (d <= 1UL) { xs=1; break; }
! 2407: /* repeat with inverted roles */
! 2408: d1 -= d;
! 2409: if (d1 >= d)
! 2410: {
! 2411: hiremainder = 0; q = 1 + divll(d1,d); d1 = hiremainder;
! 2412: xv1 += q * xv;
! 2413: }
! 2414: else
! 2415: xv1 += xv;
! 2416: } /* while */
! 2417:
! 2418: if (!(f&1)) /* division by 1 postprocessing if needed */
! 2419: {
! 2420: if (xs && d==1)
! 2421: { xv1 += d1 * xv; xs = 0; res = 1UL; }
! 2422: else if (!xs && d1==1)
! 2423: { xv += d * xv1; xs = 1; res = 1UL; }
! 2424: }
! 2425:
! 2426: if (xs)
! 2427: {
! 2428: *s = -1; *v = xv1; *v1 = xv;
! 2429: return (res ? res : (d==1 ? 1UL : d1));
! 2430: }
! 2431: else
! 2432: {
! 2433: *s = 1; *v = xv; *v1 = xv1;
! 2434: return (res ? res : (d1==1 ? 1UL : d));
! 2435: }
! 2436: }
! 2437:
! 2438:
! 2439: ulong
! 2440: xxgcduu(ulong d, ulong d1, int f,
! 2441: ulong* u, ulong* u1, ulong* v, ulong* v1, long *s)
! 2442: {
! 2443: ulong xu,xu1, xv,xv1, xs, q,res;
! 2444: LOCAL_HIREMAINDER;
! 2445:
! 2446: xs = res = 0;
! 2447: xu = xv1 = 1UL;
! 2448: xu1 = xv = 0UL;
! 2449: while (d1 > 1UL)
! 2450: {
! 2451: d -= d1; /* no need to use subll */
! 2452: if (d >= d1)
! 2453: {
! 2454: hiremainder = 0; q = 1 + divll(d,d1); d = hiremainder;
! 2455: xv += q * xv1;
! 2456: xu += q * xu1;
! 2457: }
! 2458: else
! 2459: { xv += xv1; xu += xu1; }
! 2460: /* possible loop exit */
! 2461: if (d <= 1UL) { xs=1; break; }
! 2462: /* repeat with inverted roles */
! 2463: d1 -= d;
! 2464: if (d1 >= d)
! 2465: {
! 2466: hiremainder = 0; q = 1 + divll(d1,d); d1 = hiremainder;
! 2467: xv1 += q * xv;
! 2468: xu1 += q * xu;
! 2469: }
! 2470: else
! 2471: { xv1 += xv; xu1 += xu; }
! 2472: } /* while */
! 2473:
! 2474: if (!(f&1)) /* division by 1 postprocessing if needed */
! 2475: {
! 2476: if (xs && d==1)
! 2477: {
! 2478: xv1 += d1 * xv;
! 2479: xu1 += d1 * xu;
! 2480: xs = 0; res = 1UL;
! 2481: }
! 2482: else if (!xs && d1==1)
! 2483: {
! 2484: xv += d * xv1;
! 2485: xu += d * xu1;
! 2486: xs = 1; res = 1UL;
! 2487: }
! 2488: }
! 2489:
! 2490: if (xs)
! 2491: {
! 2492: *s = -1; *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;
! 2493: return (res ? res : (d==1 ? 1UL : d1));
! 2494: }
! 2495: else
! 2496: {
! 2497: *s = 1; *u = xu; *u1 = xu1; *v = xv; *v1 = xv1;
! 2498: return (res ? res : (d1==1 ? 1UL : d));
! 2499: }
! 2500: }
! 2501:
! 2502: /*==================================
! 2503: * lgcdii(d,d1,u,u1,v,v1)
! 2504: *==================================*/
! 2505: /* Lehmer's partial extended gcd algorithm, acting on two t_INT GENs.
! 2506: *
! 2507: * Tries to determine, using the leading 2*BITS_IN_LONG significant bits of d
! 2508: * and a quantity of bits from d1 obtained by a shift of the same displacement,
! 2509: * as many partial quotients of d/d1 as possible, and assigns to [u,u1;v,v1]
! 2510: * the product of all the [0, 1; 1, q_j] thus obtained, where the leftmost
! 2511: * factor arises from the quotient of the first division step.
! 2512: *
! 2513: * MUST be called with d > d1 > 0, and with d occupying more than one
! 2514: * significant word (if it doesn't, the caller has no business with us;
! 2515: * he/she/it should use xgcduu() instead). Returns the number of reduction/
! 2516: * swap steps carried out, possibly zero, or under certain conditions minus
! 2517: * that number. When the return value is nonzero, the caller should use the
! 2518: * returned recurrence matrix to update its own copies of d,d1. When the
! 2519: * return value is non-positive, and the latest remainder after updating
! 2520: * turns out to be nonzero, the caller should at once attempt a full division,
! 2521: * rather than first trying lgcdii() again -- this typically happens when we
! 2522: * are about to encounter a quotient larger than half a word. (This is not
! 2523: * detected infallibly -- after a positive return value, it is perfectly
! 2524: * possible that the next stage will end up needing a full division. After
! 2525: * a negative return value, however, this is certain, and should be acted
! 2526: * upon.)
! 2527: *
! 2528: * (The sign information, for which xgcduu() has its return argument s, is now
! 2529: * implicit in the LSB of our return value, and the caller may take advantage
! 2530: * of the fact that a return value of +-1 implies u==0,u1==v==1 [only v1 pro-
! 2531: * vides interesting information in this case]. One might also use the fact
! 2532: * that if the return value is +-2, then u==1, but this is rather marginal.)
! 2533: *
! 2534: * If it was not possible to determine even the first quotient, either because
! 2535: * we're too close to an integer quotient or because the quotient would be
! 2536: * larger than one word (if the `leading digit' of d1 after shifting is all
! 2537: * zeros), we return 0 and do not bother to assign anything to the last four
! 2538: * args.
! 2539: *
! 2540: * The division chain might (sometimes) even run to completion. It will be
! 2541: * up to the caller to detect this case.
! 2542: *
! 2543: * This routine does _not_ change d or d1; this will also be up to the caller.
! 2544: *
! 2545: * Note that this routine does not know and does not need to know about the
! 2546: * PARI stack.
! 2547: */
! 2548: /*#define DEBUG_LEHMER 1 */
! 2549: int
! 2550: lgcdii(ulong* d, ulong* d1,
! 2551: ulong* u, ulong* u1, ulong* v, ulong* v1)
! 2552: {
! 2553: /* Strategy: (1) Extract/shift most significant bits. We assume that d
! 2554: * has at least two significant words, but we can cope with a one-word d1.
! 2555: * Let dd,dd1 be the most significant dividend word and matching part of the
! 2556: * divisor.
! 2557: * (2) Check for overflow on the first division. For our purposes, this
! 2558: * happens when the upper half of dd1 is zero. (Actually this is detected
! 2559: * during extraction.)
! 2560: * (3) Get a fix on the first quotient. We compute q = floor(dd/dd1), which
! 2561: * is an upper bound for floor(d/d1), and which gives the true value of the
! 2562: * latter if (and-almost-only-if) the remainder dd' = dd-q*dd1 is >= q.
! 2563: * (If it isn't, we give up. This is annoying because the subsequent full
! 2564: * division will repeat some work already done, but it happens very infre-
! 2565: * quently. Doing the extra-bit-fetch in this case would be awkward.)
! 2566: * (4) Finish initializations.
! 2567: *
! 2568: * The remainder of the action is comparatively boring... The main loop has
! 2569: * been unrolled once (so we don't swap things and we can apply Jebelean's
! 2570: * termination conditions which alternatingly take two different forms during
! 2571: * successive iterations). When we first run out of sufficient bits to form
! 2572: * a quotient, and have an extra word of each operand, we pull out two whole
! 2573: * word's worth of dividend bits, and divisor bits of matching significance;
! 2574: * to these we apply our partial matrix (disregarding overflow because the
! 2575: * result mod 2^(2*BITS_IN_LONG) will in fact give the correct values), and
! 2576: * re-extract one word's worth of the current dividend and a matching amount
! 2577: * of divisor bits. The affair will normally terminate with matrix entries
! 2578: * just short of a whole word. (We terminate the inner loop before these can
! 2579: * possibly overflow.)
! 2580: */
! 2581: ulong dd,dd1,ddlo,dd1lo, sh,shc; /* `digits', shift count */
! 2582: ulong xu,xu1, xv,xv1, q,res; /* recurrences, partial quotient, count */
! 2583: ulong tmp0,tmp1,tmp2,tmpd,tmpu,tmpv; /* temps */
! 2584: long ld, ld1, lz; /* t_INT effective lengths */
! 2585: int skip = 0; /* a boolean flag */
! 2586: LOCAL_OVERFLOW;
! 2587: LOCAL_HIREMAINDER;
! 2588:
! 2589: #ifdef DEBUG_LEHMER
! 2590: voir(d, -1); voir(d1, -1);
! 2591: #endif
! 2592: ld = lgefint(d); ld1 = lgefint(d1); lz = ld - ld1; /* >= 0 */
! 2593: if (lz > 1) return 0; /* rare, quick and desperate exit */
! 2594:
! 2595: d += 2; d1 += 2; /* point at the leading `digits' */
! 2596: dd1lo = 0; /* unless we find something better */
! 2597: sh = bfffo(*d); /* obtain dividend left shift count */
! 2598:
! 2599: if (sh)
! 2600: { /* do the shifting */
! 2601: shc = BITS_IN_LONG - sh;
! 2602: if (lz)
! 2603: { /* dividend longer than divisor */
! 2604: dd1 = (*d1 >> shc);
! 2605: if (!(HIGHMASK & dd1)) return 0; /* overflow detected */
! 2606: if (ld1 > 3)
! 2607: dd1lo = (*d1 << sh) + (d1[1] >> shc);
! 2608: else
! 2609: dd1lo = (*d1 << sh);
! 2610: }
! 2611: else
! 2612: { /* dividend and divisor have the same length */
! 2613: /* dd1 = shiftl(d1,sh) would have left hiremainder==0, and dd1 != 0. */
! 2614: dd1 = (*d1 << sh);
! 2615: if (!(HIGHMASK & dd1)) return 0;
! 2616: if (ld1 > 3)
! 2617: {
! 2618: dd1 += (d1[1] >> shc);
! 2619: if (ld1 > 4)
! 2620: dd1lo = (d1[1] << sh) + (d1[2] >> shc);
! 2621: else
! 2622: dd1lo = (d1[1] << sh);
! 2623: }
! 2624: }
! 2625: /* following lines assume d to have 2 or more significant words */
! 2626: dd = (*d << sh) + (d[1] >> shc);
! 2627: if (ld > 4)
! 2628: ddlo = (d[1] << sh) + (d[2] >> shc);
! 2629: else
! 2630: ddlo = (d[1] << sh);
! 2631: }
! 2632: else
! 2633: { /* no shift needed */
! 2634: if (lz) return 0; /* div'd longer than div'r: o'flow automatic */
! 2635: dd1 = *d1;
! 2636: if (!(HIGHMASK & dd1)) return 0;
! 2637: if(ld1 > 3) dd1lo = d1[1];
! 2638: /* assume again that d has another significant word */
! 2639: dd = *d; ddlo = d[1];
! 2640: }
! 2641: #ifdef DEBUG_LEHMER
! 2642: fprintf(stderr, " %lx:%lx, %lx:%lx\n", dd, ddlo, dd1, dd1lo);
! 2643: #endif
! 2644:
! 2645: /* First subtraction/division stage. (If a subtraction initially suffices,
! 2646: * we don't divide at all.) If a Jebelean condition is violated, and we
! 2647: * can't fix it even by looking at the low-order bits in ddlo,dd1lo, we
! 2648: * give up and ask for a full division. Otherwise we commit the result,
! 2649: * possibly deciding to re-shift immediately afterwards.
! 2650: */
! 2651: dd -= dd1;
! 2652: if (dd < dd1)
! 2653: { /* first quotient known to be == 1 */
! 2654: xv1 = 1UL;
! 2655: if (!dd) /* !(Jebelean condition), extraspecial case */
! 2656: { /* note this can actually happen... Now
! 2657: * q==1 is known, but we underflow already.
! 2658: * OTOH we've just shortened d by a whole word.
! 2659: * Thus we feel pleased with ourselves and
! 2660: * return. (The re-shift code below would
! 2661: * do so anyway.) */
! 2662: *u = 0; *v = *u1 = *v1 = 1UL;
! 2663: return -1; /* Next step will be a full division. */
! 2664: } /* if !(Jebelean) then */
! 2665: }
! 2666: else
! 2667: { /* division indicated */
! 2668: hiremainder = 0;
! 2669: xv1 = 1 + divll(dd, dd1); /* xv1: alternative spelling of `q', here ;) */
! 2670: dd = hiremainder;
! 2671: if (dd < xv1) /* !(Jebelean cond'), non-extra special case */
! 2672: { /* Attempt to complete the division using the
! 2673: * less significant bits, before skipping right
! 2674: * past the 1st loop to the reshift stage. */
! 2675: ddlo = subll(ddlo, mulll(xv1, dd1lo));
! 2676: dd = subllx(dd, hiremainder);
! 2677:
! 2678: /* If we now have an overflow, q was _certainly_ too large. Thanks to
! 2679: * our decision not to get here unless the original dd1 had bits set in
! 2680: * the upper half of the word, however, we now do know that the correct
! 2681: * quotient is in fact q-1. Adjust our data accordingly. */
! 2682: if (overflow)
! 2683: {
! 2684: xv1--;
! 2685: ddlo = addll(ddlo,dd1lo);
! 2686: dd = addllx(dd,dd1); /* overflows again which cancels the borrow */
! 2687: /* ...and fall through to skip=1 below */
! 2688: }
! 2689: else
! 2690: /* Test Jebelean condition anew, at this point using _all_ the extracted
! 2691: * bits we have. This is clutching at straws; we have a more or less
! 2692: * even chance of succeeding this time. Note that if we fail, we really
! 2693: * do not know whether the correct quotient would have been q or some
! 2694: * smaller value. */
! 2695: if (!dd && ddlo < xv1) return 0;
! 2696:
! 2697: /* Otherwise, we now know that q is correct, but we cannot go into the
! 2698: * 1st loop. Raise a flag so we'll remember to skip past the loop.
! 2699: * Get here also after the q-1 adjustment case. */
! 2700: skip = 1;
! 2701: } /* if !(Jebelean) then */
! 2702: }
! 2703: xu = 0; xv = xu1 = 1UL; res = 1;
! 2704: #ifdef DEBUG_LEHMER
! 2705: fprintf(stderr, " q = %ld, %lx, %lx\n", xv1, dd1, dd);
! 2706: #endif
! 2707:
! 2708: /* Some invariants from here across the first loop:
! 2709: *
! 2710: * At this point, and again after we are finished with the first loop and
! 2711: * subsequent conditional, a division and the associated update of the
! 2712: * recurrence matrix have just been carried out completely. The matrix
! 2713: * xu,xu1;xv,xv1 has been initialized (or updated, possibly with permuted
! 2714: * columns), and the current remainder == next divisor (dd at the moment)
! 2715: * is nonzero (it might be zero here, but then skip will have been set).
! 2716: *
! 2717: * After the first loop, or when skip is set already, it will also be the
! 2718: * case that there aren't sufficiently many bits to continue without re-
! 2719: * shifting. If the divisor after reshifting is zero, or indeed if it
! 2720: * doesn't have more than half a word's worth of bits, we will have to
! 2721: * return at that point. Otherwise, we proceed into the second loop.
! 2722: *
! 2723: * Furthermore, when we reach the re-shift stage, dd:ddlo and dd1:dd1lo will
! 2724: * already reflect the result of applying the current matrix to the old
! 2725: * ddorig:ddlo and dd1orig:dd1lo. (For the first iteration above, this
! 2726: * was easy to achieve, and we didn't even need to peek into the (now
! 2727: * no longer existent!) saved words. After the loop, we'll stop for a
! 2728: * moment to merge in the ddlo,dd1lo contributions.)
! 2729: *
! 2730: * Note that after the first division, even an a priori quotient of 1 cannot
! 2731: * be trusted until we've checked Jebelean's condition -- it cannot be too
! 2732: * large, of course, but it might be too small.
! 2733: */
! 2734:
! 2735: if (!skip)
! 2736: {
! 2737: while (1)
! 2738: {
! 2739: /* First half of loop divides dd into dd1, and leaves the recurrence
! 2740: * matrix xu,...,xv1 groomed the wrong way round (xu,xv will be the newer
! 2741: * entries) when successful. */
! 2742: tmpd = dd1 - dd;
! 2743: if (tmpd < dd)
! 2744: { /* quotient suspected to be 1 */
! 2745: #ifdef DEBUG_LEHMER
! 2746: q = 1;
! 2747: #endif
! 2748: tmpu = xu + xu1; /* cannot overflow -- everything bounded by
! 2749: * the original dd during first loop */
! 2750: tmpv = xv + xv1;
! 2751: }
! 2752: else
! 2753: { /* division indicated */
! 2754: hiremainder = 0;
! 2755: q = 1 + divll(tmpd, dd);
! 2756: tmpd = hiremainder;
! 2757: tmpu = xu + q*xu1; /* can't overflow, but may need to be undone */
! 2758: tmpv = xv + q*xv1;
! 2759: }
! 2760:
! 2761: tmp0 = addll(tmpv, xv1);
! 2762: if ((tmpd < tmpu) || overflow ||
! 2763: (dd - tmpd < tmp0)) /* !(Jebelean cond.) */
! 2764: break; /* skip ahead to reshift stage */
! 2765: else
! 2766: { /* commit dd1, xu, xv */
! 2767: res++;
! 2768: dd1 = tmpd; xu = tmpu; xv = tmpv;
! 2769: #ifdef DEBUG_LEHMER
! 2770: fprintf(stderr, " q = %ld, %lx, %lx [%lu,%lu;%lu,%lu]\n",
! 2771: q, dd, dd1, xu1, xu, xv1, xv);
! 2772: #endif
! 2773: }
! 2774:
! 2775: /* Second half of loop divides dd1 into dd, and the matrix returns to its
! 2776: * normal arrangement. */
! 2777: tmpd = dd - dd1;
! 2778: if (tmpd < dd1)
! 2779: { /* quotient suspected to be 1 */
! 2780: #ifdef DEBUG_LEHMER
! 2781: q = 1;
! 2782: #endif
! 2783: tmpu = xu1 + xu; /* cannot overflow */
! 2784: tmpv = xv1 + xv;
! 2785: }
! 2786: else
! 2787: { /* division indicated */
! 2788: hiremainder = 0;
! 2789: q = 1 + divll(tmpd, dd1);
! 2790: tmpd = hiremainder;
! 2791: tmpu = xu1 + q*xu;
! 2792: tmpv = xv1 + q*xv;
! 2793: }
! 2794:
! 2795: tmp0 = addll(tmpu, xu);
! 2796: if ((tmpd < tmpv) || overflow ||
! 2797: (dd1 - tmpd < tmp0)) /* !(Jebelean cond.) */
! 2798: break; /* skip ahead to reshift stage */
! 2799: else
! 2800: { /* commit dd, xu1, xv1 */
! 2801: res++;
! 2802: dd = tmpd; xu1 = tmpu; xv1 = tmpv;
! 2803: #ifdef DEBUG_LEHMER
! 2804: fprintf(stderr, " q = %ld, %lx, %lx [%lu,%lu;%lu,%lu]\n",
! 2805: q, dd1, dd, xu, xu1, xv, xv1);
! 2806: #endif
! 2807: }
! 2808:
! 2809: } /* end of first loop */
! 2810:
! 2811: /* Intermezzo: update dd:ddlo, dd1:dd1lo. (But not if skip is set.) */
! 2812:
! 2813: if (res&1)
! 2814: {
! 2815: /* after failed division in 1st half of loop:
! 2816: * [dd1:dd1lo,dd:ddlo] = [ddorig:ddlo,dd1orig:dd1lo]
! 2817: * * [ -xu, xu1 ; xv, -xv1 ]
! 2818: * (Actually, we only multiply [ddlo,dd1lo] onto the matrix and
! 2819: * add the high-order remainders + overflows onto [dd1,dd].)
! 2820: */
! 2821: tmp1 = mulll(ddlo, xu); tmp0 = hiremainder;
! 2822: tmp1 = subll(mulll(dd1lo,xv), tmp1);
! 2823: dd1 += subllx(hiremainder, tmp0);
! 2824: tmp2 = mulll(ddlo, xu1); tmp0 = hiremainder;
! 2825: ddlo = subll(tmp2, mulll(dd1lo,xv1));
! 2826: dd += subllx(tmp0, hiremainder);
! 2827: dd1lo = tmp1;
! 2828: }
! 2829: else
! 2830: {
! 2831: /* after failed division in 2nd half of loop:
! 2832: * [dd:ddlo,dd1:dd1lo] = [ddorig:ddlo,dd1orig:dd1lo]
! 2833: * * [ xu1, -xu ; -xv1, xv ]
! 2834: * (Actually, we only multiply [ddlo,dd1lo] onto the matrix and
! 2835: * add the high-order remainders + overflows onto [dd,dd1].)
! 2836: */
! 2837: tmp1 = mulll(ddlo, xu1); tmp0 = hiremainder;
! 2838: tmp1 = subll(tmp1, mulll(dd1lo,xv1));
! 2839: dd += subllx(tmp0, hiremainder);
! 2840: tmp2 = mulll(ddlo, xu); tmp0 = hiremainder;
! 2841: dd1lo = subll(mulll(dd1lo,xv), tmp2);
! 2842: dd1 += subllx(hiremainder, tmp0);
! 2843: ddlo = tmp1;
! 2844: }
! 2845: #ifdef DEBUG_LEHMER
! 2846: fprintf(stderr, " %lx:%lx, %lx:%lx\n", dd, ddlo, dd1, dd1lo);
! 2847: #endif
! 2848: } /* end of skip-pable section: get here also, with res==1, when there
! 2849: * was a problem immediately after the very first division. */
! 2850:
! 2851: /* Re-shift. Note: the shift count _can_ be zero, viz. under the following
! 2852: * precise conditions: The original dd1 had its topmost bit set, so the 1st
! 2853: * q was 1, and after subtraction, dd had its topmost bit unset. If now
! 2854: * dd==0, we'd have taken the return exit already, so we couldn't have got
! 2855: * here. If not, then it must have been the second division which has gone
! 2856: * amiss (because dd1 was very close to an exact multiple of the remainder
! 2857: * dd value, so this will be very rare). At this point, we'd have a fairly
! 2858: * slim chance of fixing things by re-examining dd1:dd1lo vs. dd:ddlo, but
! 2859: * this is not guaranteed to work. Instead of trying, we return at once.
! 2860: * The caller will see to it that the initial subtraction is re-done using
! 2861: * _all_ the bits of both operands, which already helps, and the next round
! 2862: * will either be a full division (if dd occupied a halfword or less), or
! 2863: * another llgcdii() first step. In the latter case, since we try a little
! 2864: * harder during our first step, we may actually be able to fix the problem,
! 2865: * and get here again with improved low-order bits and with another step
! 2866: * under our belt. Otherwise we'll have given up above and forced a full-
! 2867: * blown division.
! 2868: *
! 2869: * If res is even, the shift count _cannot_ be zero. (The first step forces
! 2870: * a zero into the remainder's MSB, and all subsequent remainders will have
! 2871: * inherited it.)
! 2872: *
! 2873: * The re-shift stage exits if the next divisor has at most half a word's
! 2874: * worth of bits.
! 2875: *
! 2876: * For didactic reasons, the second loop will be arranged in the same way
! 2877: * as the first -- beginning with the division of dd into dd1, as if res
! 2878: * was odd. To cater for this, if res is actually even, we swap things
! 2879: * around during reshifting. (During the second loop, the parity of res
! 2880: * does not matter; we know in which half of the loop we are when we decide
! 2881: * to return.)
! 2882: */
! 2883: #ifdef DEBUG_LEHMER
! 2884: fprintf(stderr, "(sh)");
! 2885: #endif
! 2886:
! 2887: if (res&1)
! 2888: { /* after odd number of division(s) */
! 2889: if (dd1 && (sh = bfffo(dd1)))
! 2890: {
! 2891: shc = BITS_IN_LONG - sh;
! 2892: dd = (ddlo >> shc) + (dd << sh);
! 2893: if (!(HIGHMASK & dd))
! 2894: {
! 2895: *u = xu; *u1 = xu1; *v = xv; *v1 = xv1;
! 2896: return -res; /* full division asked for */
! 2897: }
! 2898: dd1 = (dd1lo >> shc) + (dd1 << sh);
! 2899: }
! 2900: else
! 2901: { /* time to return: <= 1 word left, or sh==0 */
! 2902: *u = xu; *u1 = xu1; *v = xv; *v1 = xv1;
! 2903: return res;
! 2904: }
! 2905: }
! 2906: else
! 2907: { /* after even number of divisions */
! 2908: if (dd)
! 2909: {
! 2910: sh = bfffo(dd); /* Known to be positive. */
! 2911: shc = BITS_IN_LONG - sh;
! 2912: /* dd:ddlo will become the new dd1, and v.v. */
! 2913: tmpd = (ddlo >> shc) + (dd << sh);
! 2914: dd = (dd1lo >> shc) + (dd1 << sh);
! 2915: dd1 = tmpd;
! 2916: /* This has completed the swap; now dd is again the current divisor.
! 2917: * The following test originally inspected dd1 -- a most subtle and
! 2918: * most annoying bug. The Management. */
! 2919: if (HIGHMASK & dd)
! 2920: {
! 2921: /* recurrence matrix is the wrong way round; swap it. */
! 2922: tmp0 = xu; xu = xu1; xu1 = tmp0;
! 2923: tmp0 = xv; xv = xv1; xv1 = tmp0;
! 2924: }
! 2925: else
! 2926: {
! 2927: /* recurrence matrix is the wrong way round; fix this. */
! 2928: *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;
! 2929: return -res; /* full division asked for */
! 2930: }
! 2931: }
! 2932: else
! 2933: { /* time to return: <= 1 word left */
! 2934: *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;
! 2935: return res;
! 2936: }
! 2937: } /* end reshift */
! 2938:
! 2939: #ifdef DEBUG_LEHMER
! 2940: fprintf(stderr, " %lx:%lx, %lx:%lx\n", dd, ddlo, dd1, dd1lo);
! 2941: #endif
! 2942:
! 2943: /* The Second Loop. Rip-off of the first, but we now check for overflow
! 2944: * in the recurrences. Returns instead of breaking when we cannot fix the
! 2945: * quotient any longer.
! 2946: */
! 2947:
! 2948: for(;;)
! 2949: {
! 2950: /* First half of loop divides dd into dd1, and leaves the recurrence
! 2951: * matrix xu,...,xv1 groomed the wrong way round (xu,xv will be the newer
! 2952: * entries) when successful. */
! 2953: tmpd = dd1 - dd;
! 2954: if (tmpd < dd)
! 2955: { /* quotient suspected to be 1 */
! 2956: #ifdef DEBUG_LEHMER
! 2957: q = 1;
! 2958: #endif
! 2959: tmpu = xu + xu1;
! 2960: tmpv = addll(xv, xv1); /* xv,xv1 will overflow first */
! 2961: tmp1 = overflow;
! 2962: }
! 2963: else
! 2964: { /* division indicated */
! 2965: hiremainder = 0;
! 2966: q = 1 + divll(tmpd, dd);
! 2967: tmpd = hiremainder;
! 2968: tmpu = xu + q*xu1;
! 2969: tmpv = addll(xv, mulll(q,xv1));
! 2970: tmp1 = overflow | hiremainder;
! 2971: }
! 2972:
! 2973: tmp0 = addll(tmpv, xv1);
! 2974: if ((tmpd < tmpu) || overflow || tmp1 ||
! 2975: (dd - tmpd < tmp0)) /* !(Jebelean cond.) */
! 2976: {
! 2977: /* The recurrence matrix has not yet been warped... */
! 2978: *u = xu; *u1 = xu1; *v = xv; *v1 = xv1;
! 2979: return res;
! 2980: }
! 2981: else
! 2982: { /* commit dd1, xu, xv */
! 2983: res++;
! 2984: dd1 = tmpd; xu = tmpu; xv = tmpv;
! 2985: #ifdef DEBUG_LEHMER
! 2986: fprintf(stderr, " q = %ld, %lx, %lx\n", q, dd, dd1);
! 2987: #endif
! 2988: }
! 2989:
! 2990: /* Second half of loop divides dd1 into dd, and the matrix returns to its
! 2991: * normal arrangement. */
! 2992: tmpd = dd - dd1;
! 2993: if (tmpd < dd1)
! 2994: { /* quotient suspected to be 1 */
! 2995: #ifdef DEBUG_LEHMER
! 2996: q = 1;
! 2997: #endif
! 2998: tmpu = xu1 + xu;
! 2999: tmpv = addll(xv1, xv);
! 3000: tmp1 = overflow;
! 3001: }
! 3002: else
! 3003: { /* division indicated */
! 3004: hiremainder = 0;
! 3005: q = 1 + divll(tmpd, dd1);
! 3006: tmpd = hiremainder;
! 3007: tmpu = xu1 + q*xu;
! 3008: tmpv = addll(xv1, mulll(q, xv));
! 3009: tmp1 = overflow | hiremainder;
! 3010: }
! 3011:
! 3012: tmp0 = addll(tmpu, xu);
! 3013: if ((tmpd < tmpv) || overflow || tmp1 ||
! 3014: (dd1 - tmpd < tmp0)) /* !(Jebelean cond.) */
! 3015: {
! 3016: /* The recurrence matrix has not yet been unwarped, so it is
! 3017: * the wrong way round; fix this. */
! 3018: *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;
! 3019: return res;
! 3020: }
! 3021: else
! 3022: { /* commit dd, xu1, xv1 */
! 3023: res++;
! 3024: dd = tmpd; xu1 = tmpu; xv1 = tmpv;
! 3025: #ifdef DEBUG_LEHMER
! 3026: fprintf(stderr, " q = %ld, %lx, %lx\n", q, dd1, dd);
! 3027: #endif
! 3028: }
! 3029:
! 3030: } /* end of second loop */
! 3031:
! 3032: return(0); /* never reached */
! 3033: }
! 3034:
! 3035: /*==================================
! 3036: * invmod(a,b,res)
! 3037: *==================================
! 3038: * If a is invertible, return 1, and set res = a^{ -1 }
! 3039: * Otherwise, return 0, and set res = gcd(a,b)
! 3040:
! 3041: * Temporary - to be replaced with mpxgcd()
! 3042: * I think we won't need to change much ... just the parameter-passing and
! 3043: * return value conventions, and computing a final u in addition to a final
! 3044: * v value before returning (and doing so even when the gcd isn't 1)
! 3045: */
! 3046:
! 3047: #define mysetsigne(d,s) { fprintferr("%ld %ld\n",s,signe(d)); if (signe(d)==-(s)) setsigne(d,s); }
! 3048:
! 3049: int
! 3050: invmod(GEN a, GEN b, GEN *res)
! 3051: {
! 3052: GEN v,v1,d,d1,q,r;
! 3053: long av,av1,lim;
! 3054: long s;
! 3055: ulong g;
! 3056: ulong xu,xu1,xv,xv1; /* Lehmer stage recurrence matrix */
! 3057: int lhmres; /* Lehmer stage return value */
! 3058:
! 3059: if (typ(a) != t_INT || typ(b) != t_INT) err(arither1);
! 3060: if (!signe(b)) { *res=absi(a); return 0; }
! 3061: av = avma;
! 3062: if (lgefint(b) == 3) /* single-word affair */
! 3063: {
! 3064: d1 = modiu(a, (ulong)b[2]);
! 3065: if (d1 == gzero)
! 3066: {
! 3067: if (b[2] == 1L)
! 3068: { *res = gzero; return 1; }
! 3069: else
! 3070: { *res = absi(b); return 0; }
! 3071: }
! 3072: g = xgcduu((ulong)b[2], (ulong)d1[2], 1, &xv, &xv1, &s);
! 3073: #ifdef DEBUG_LEHMER
! 3074: fprintferr(" <- %lu,%lu\n", (ulong)b[2], (ulong)d1[2]);
! 3075: fprintferr(" -> %lu,%ld,%lu; %lx\n", g,s,xv1,avma);
! 3076: #endif
! 3077: avma = av;
! 3078: if (g != 1UL) { *res = utoi(g); return 0; }
! 3079: xv = xv1 % (ulong)b[2]; if (s < 0) xv = ((ulong)b[2]) - xv;
! 3080: *res = utoi(xv); return 1;
! 3081: }
! 3082:
! 3083: (void)new_chunk(lgefint(b));
! 3084: d = absi(b); d1 = modii(a,d);
! 3085:
! 3086: v=gzero; v1=gun; /* general case */
! 3087: #ifdef DEBUG_LEHMER
! 3088: fprintferr("INVERT: -------------------------\n");
! 3089: output(d1);
! 3090: #endif
! 3091: av1 = avma; lim = stack_lim(av,1);
! 3092:
! 3093: while (lgefint(d) > 3 && signe(d1))
! 3094: {
! 3095: #ifdef DEBUG_LEHMER
! 3096: fprintferr("Calling Lehmer:\n");
! 3097: #endif
! 3098: lhmres = lgcdii((ulong*)d, (ulong*)d1, &xu, &xu1, &xv, &xv1);
! 3099: if (lhmres != 0) /* check progress */
! 3100: { /* apply matrix */
! 3101: #ifdef DEBUG_LEHMER
! 3102: fprintferr("Lehmer returned %d [%lu,%lu;%lu,%lu].\n",
! 3103: lhmres, xu, xu1, xv, xv1);
! 3104: #endif
! 3105: if ((lhmres == 1) || (lhmres == -1))
! 3106: {
! 3107: if (xv1 == 1)
! 3108: {
! 3109: r = subii(d,d1); d=d1; d1=r;
! 3110: a = subii(v,v1); v=v1; v1=a;
! 3111: }
! 3112: else
! 3113: {
! 3114: r = subii(d, mului(xv1,d1)); d=d1; d1=r;
! 3115: a = subii(v, mului(xv1,v1)); v=v1; v1=a;
! 3116: }
! 3117: }
! 3118: else
! 3119: {
! 3120: r = subii(muliu(d,xu), muliu(d1,xv));
! 3121: a = subii(muliu(v,xu), muliu(v1,xv));
! 3122: d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
! 3123: v1 = subii(muliu(v,xu1), muliu(v1,xv1)); v = a;
! 3124: if (lhmres&1) {
! 3125: setsigne(d,-signe(d));
! 3126: setsigne(v,-signe(v));
! 3127: }
! 3128: else {
! 3129: setsigne(d1,-signe(d1));
! 3130: setsigne(v1,-signe(v1));
! 3131: }
! 3132: }
! 3133: }
! 3134: #ifdef DEBUG_LEHMER
! 3135: else
! 3136: fprintferr("Lehmer returned 0.\n");
! 3137: output(d); output(d1); output(v); output(v1);
! 3138: sleep(1);
! 3139: #endif
! 3140:
! 3141: if (lhmres <= 0 && signe(d1))
! 3142: {
! 3143: q = dvmdii(d,d1,&r);
! 3144: #ifdef DEBUG_LEHMER
! 3145: fprintferr("Full division:\n");
! 3146: printf(" q = "); output(q); sleep (1);
! 3147: #endif
! 3148: a = subii(v,mulii(q,v1));
! 3149: v=v1; v1=a;
! 3150: d=d1; d1=r;
! 3151: }
! 3152: if (low_stack(lim, stack_lim(av,1)))
! 3153: {
! 3154: GEN *gptr[4]; gptr[0]=&d; gptr[1]=&d1; gptr[2]=&v; gptr[3]=&v1;
! 3155: if(DEBUGMEM>1) err(warnmem,"invmod");
! 3156: gerepilemany(av1,gptr,4);
! 3157: }
! 3158: } /* end while */
! 3159:
! 3160: /* Postprocessing - final sprint */
! 3161: if (signe(d1))
! 3162: {
! 3163: /* Assertions: lgefint(d)==lgefint(d1)==3, and
! 3164: * gcd(d,d1) is nonzero and fits into one word
! 3165: */
! 3166: g = xxgcduu(d[2], d1[2], 1, &xu, &xu1, &xv, &xv1, &s);
! 3167: #ifdef DEBUG_LEHMER
! 3168: output(d);output(d1);output(v);output(v1);
! 3169: fprintferr(" <- %lu,%lu\n", (ulong)d[2], (ulong)d1[2]);
! 3170: fprintferr(" -> %lu,%ld,%lu; %lx\n", g,s,xv1,avma);
! 3171: #endif
! 3172: if (g != 1UL) { avma = av; *res = utoi(g); return 0; }
! 3173: /* (From the xgcduu() blurb:)
! 3174: * For finishing the multiword modinv, we now have to multiply the
! 3175: * returned matrix (with properly adjusted signs) onto the values
! 3176: * v' and v1' previously obtained from the multiword division steps.
! 3177: * Actually, it is sufficient to take the scalar product of [v',v1']
! 3178: * with [u1,-v1], and change the sign if s==1.
! 3179: */
! 3180: v = subii(muliu(v,xu1),muliu(v1,xv1));
! 3181: if (s > 0) setsigne(v,-signe(v));
! 3182: avma = av; *res = modii(v,b);
! 3183: #ifdef DEBUG_LEHMER
! 3184: output(*res); fprintfderr("============================Done.\n");
! 3185: sleep(1);
! 3186: #endif
! 3187: return 1;
! 3188: }
! 3189: /* get here when the final sprint was skipped (d1 was zero already) */
! 3190: avma = av;
! 3191: if (!egalii(d,gun)) { *res = icopy(d); return 0; }
! 3192: *res = modii(v,b);
! 3193: #ifdef DEBUG_LEHMER
! 3194: output(*res); fprintferr("============================Done.\n");
! 3195: sleep(1);
! 3196: #endif
! 3197: return 1;
! 3198: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mp.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / kernel / none