00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 static char *rcsid="@(#)$Id$";
00019 #include "eus.h"
00020
00021 #define ltz(x) ((x)<0)
00022 #define gez(x) ((x)>=0)
00023
00024 extern pointer makebig();
00025
00026 #if (WORD_SIZE == 64)
00027
00028 void extended_mul(d, q, r, hp, lp)
00029 eusinteger_t d, q, r;
00030 eusinteger_t *hp, *lp;
00031 {
00032 eusinteger_t ld, hd, lq, hq, m, lm, hm, lr;
00033 eusinteger_t hz, lz;
00034
00035
00036 ld = (d & 0x00000000ffffffff);
00037 hd = (d >> 32);
00038 hd = (hd & 0x000000007fffffff);
00039 lq = (q & 0x00000000ffffffff);
00040 hq = (q >> 32);
00041 hq = (hq & 0x000000007fffffff);
00042 lr = r;
00043
00044 m = ld * lq;
00045 hm = (m & 0x8000000000000000);
00046 hm = (hm >> 32);
00047 hm = (hm & 0x0000000080000000);
00048 m = (m & 0x7fffffffffffffff);
00049
00050 m = m + lr;
00051
00052 lz = (m & 0x00000000ffffffff);
00053
00054 m = (m >> 32);
00055 m = (m & 0x00000000ffffffff);
00056 m = m + hm;
00057 m = m + (lq * hd);
00058 hm = (m & 0x8000000000000000);
00059 hm = (hm >> 31);
00060 hm = (hm & 0x00000001ffffffff);
00061 m = (m & 0x7fffffffffffffff);
00062
00063 m = m + (hq * ld);
00064
00065 lm = (m & 0x000000007fffffff);
00066
00067 lz = lz + (lm << 32);
00068
00069 m = (m >> 31);
00070 m = (m & 0x00000001ffffffff);
00071 m = m + hm;
00072
00073 hz = m + (hq * hd * 2);
00074
00075
00076 *hp = hz;
00077 *lp = lz;
00078
00079 }
00080
00081 void extended_div(d, h, l, qp, rp)
00082 eusinteger_t d, h, l;
00083 eusinteger_t *qp, *rp;
00084 {
00085 eusinteger_t lh, ld, ll, ans;
00086 int i;
00087
00088 ld = d;
00089 lh = h;
00090 ll = l;
00091
00092
00093 ans = 0;
00094 for(i=63;i>0;i--) {
00095 ans = (ans << 1);
00096 lh = (lh << 1);
00097 ll = (ll << 1);
00098 if (ll < 0) lh += 1;
00099 if (lh >= d) {
00100 lh -= ld;
00101 ans += 1;
00102 }
00103 if (lh < 0) {
00104 lh = (0x7fffffffffffffff - ld) +
00105 (lh & 0x7fffffffffffffff) + 1;
00106 ans += 1;
00107 }
00108 }
00109
00110
00111 *qp = ans;
00112 *rp = lh;
00113 }
00114
00115 #else
00116 #if (!Solaris2 || !sun4)
00117
00118 void extended_mul(d, q, r, hp, lp)
00119 eusinteger_t d, q, r;
00120 eusinteger_t *hp, *lp;
00121 {
00122 long int ld, hd, lq, hq, m, lm, hm, lr;
00123 int hz, lz;
00124
00125
00126 ld = (d & 0x0000ffff);
00127 hd = (d >> 16);
00128 hd = (hd & 0x00007fff);
00129 lq = (q & 0x0000ffff);
00130 hq = (q >> 16);
00131 hq = (hq & 0x00007fff);
00132 lr = r;
00133
00134 m = ld * lq;
00135 hm = (m & 0x80000000);
00136 hm = (hm >> 16);
00137 hm = (hm & 0x00008000);
00138 m = (m & 0x7fffffff);
00139
00140 m = m + lr;
00141
00142 lz = (m & 0x0000ffff);
00143
00144 m = (m >> 16);
00145 m = (m & 0x0000ffff);
00146 m = m + hm;
00147 m = m + (lq * hd);
00148 hm = (m & 0x80000000);
00149 hm = (hm >> 15);
00150 hm = (hm & 0x0001ffff);
00151 m = (m & 0x7fffffff);
00152
00153 m = m + (hq * ld);
00154
00155 lm = (m & 0x00007fff);
00156
00157 lz = lz + (lm << 16);
00158
00159 m = (m >> 15);
00160 m = (m & 0x0001ffff);
00161 m = m + hm;
00162
00163 hz = m + (hq * hd * 2);
00164
00165
00166 *hp = hz;
00167 *lp = lz;
00168
00169 }
00170
00171 void extended_div(d, h, l, qp, rp)
00172 eusinteger_t d, h, l;
00173 eusinteger_t *qp, *rp;
00174 {
00175 long int lh, ld, ll;
00176 int i,ans;
00177
00178 ld = d;
00179 lh = h;
00180 ll = l;
00181
00182
00183 ans = 0;
00184 for(i=31;i>0;i--) {
00185 ans = (ans << 1);
00186 lh = (lh << 1);
00187 ll = (ll << 1);
00188 if (ll < 0) lh += 1;
00189 if (lh >= d) {
00190 lh -= ld;
00191 ans += 1;
00192 }
00193 if (lh < 0) {
00194 lh = (0x7fffffff - ld) + (lh & 0x7fffffff) + 1;
00195 ans += 1;
00196 }
00197 }
00198
00199
00200 *qp = ans;
00201 *rp = lh;
00202 }
00203
00204 #else
00205
00206
00207 #define u_int unsigned int
00208
00209 #define DEBUG
00210 #undef DEBUG
00211 #ifdef DEBUG
00212 extended_mul(d,q,r,hp,lp)
00213 u_int d,q,r,*hp,*lp;
00214 {
00215 printf("ext.mul (%x %x %x)",d,q,r);fflush(stdout);
00216 low_extended_mul(d,q,r,hp,lp);
00217 printf(" =%x:%08x\n",*hp,*lp);fflush(stdout);
00218 }
00219 static low_extended_mul(d, q, r, hp, lp)
00220 #else
00221 extended_mul(d,q,r,hp,lp)
00222 #endif
00223 register u_int d, q, r;
00224 register u_int *hp, *lp;
00225 {
00226 asm(" save %sp, -0x60, %sp");
00227 asm(" mov %i0,%o0");
00228 asm(" call .umul,2");
00229 asm(" mov %i1,%o1");
00230 asm(" addcc %o0,%i2,%o0");
00231 asm(" addx %o1, 0,%o1");
00232
00233 asm(" sll %o1, 2,%o1");
00234 asm(" srl %o1, 1,%o1");
00235 asm(" srl %o0, 31,%i2");
00236 asm(" or %o1,%i2,%o1");
00237
00238 asm(" sll %o0, 1,%o0");
00239 asm(" srl %o0, 1,%o0");
00240
00241 asm(" st %o0,[%i4]");
00242 asm(" st %o1,[%i3]");
00243 asm(" ret");
00244 asm(" restore");
00245 }
00246
00247 extended_div(d, h, l, qp, rp)
00248 register u_int d, h, l;
00249 u_int *qp, *rp;
00250 {
00251 register u_int q,one;
00252 #ifdef DEBUG
00253 printf("ext.div (%x %x:%x)",d,h,l);fflush(stdout);
00254 #endif
00255 l=l<<1; one=1<<31; q=0;
00256 if (h>=d) { q+= one;h-=d; }
00257 h<<=1; h|=l>>31; l<<=1; one>>=1;
00258 if (h>=d) { q+= one;h-=d; }
00259 h<<=1; h|=l>>31; l<<=1; one>>=1;
00260 if (h>=d) { q+= one;h-=d; }
00261 h<<=1; h|=l>>31; l<<=1; one>>=1;
00262 if (h>=d) { q+= one;h-=d; }
00263 h<<=1; h|=l>>31; l<<=1; one>>=1;
00264 if (h>=d) { q+= one;h-=d; }
00265 h<<=1; h|=l>>31; l<<=1; one>>=1;
00266 if (h>=d) { q+= one;h-=d; }
00267 h<<=1; h|=l>>31; l<<=1; one>>=1;
00268 if (h>=d) { q+= one;h-=d; }
00269 h<<=1; h|=l>>31; l<<=1; one>>=1;
00270 if (h>=d) { q+= one;h-=d; }
00271 h<<=1; h|=l>>31; l<<=1; one>>=1;
00272
00273 if (h>=d) { q+= one;h-=d; }
00274 h<<=1; h|=l>>31; l<<=1; one>>=1;
00275 if (h>=d) { q+= one;h-=d; }
00276 h<<=1; h|=l>>31; l<<=1; one>>=1;
00277 if (h>=d) { q+= one;h-=d; }
00278 h<<=1; h|=l>>31; l<<=1; one>>=1;
00279 if (h>=d) { q+= one;h-=d; }
00280 h<<=1; h|=l>>31; l<<=1; one>>=1;
00281 if (h>=d) { q+= one;h-=d; }
00282 h<<=1; h|=l>>31; l<<=1; one>>=1;
00283 if (h>=d) { q+= one;h-=d; }
00284 h<<=1; h|=l>>31; l<<=1; one>>=1;
00285 if (h>=d) { q+= one;h-=d; }
00286 h<<=1; h|=l>>31; l<<=1; one>>=1;
00287 if (h>=d) { q+= one;h-=d; }
00288 h<<=1; h|=l>>31; l<<=1; one>>=1;
00289
00290 if (h>=d) { q+= one;h-=d; }
00291 h<<=1; h|=l>>31; l<<=1; one>>=1;
00292 if (h>=d) { q+= one;h-=d; }
00293 h<<=1; h|=l>>31; l<<=1; one>>=1;
00294 if (h>=d) { q+= one;h-=d; }
00295 h<<=1; h|=l>>31; l<<=1; one>>=1;
00296 if (h>=d) { q+= one;h-=d; }
00297 h<<=1; h|=l>>31; l<<=1; one>>=1;
00298 if (h>=d) { q+= one;h-=d; }
00299 h<<=1; h|=l>>31; l<<=1; one>>=1;
00300 if (h>=d) { q+= one;h-=d; }
00301 h<<=1; h|=l>>31; l<<=1; one>>=1;
00302 if (h>=d) { q+= one;h-=d; }
00303 h<<=1; h|=l>>31; l<<=1; one>>=1;
00304 if (h>=d) { q+= one;h-=d; }
00305 h<<=1; h|=l>>31; l<<=1; one>>=1;
00306
00307 if (h>=d) { q+= one;h-=d; }
00308 h<<=1; h|=l>>31; l<<=1; one>>=1;
00309 if (h>=d) { q+= one;h-=d; }
00310 h<<=1; h|=l>>31; l<<=1; one>>=1;
00311 if (h>=d) { q+= one;h-=d; }
00312 h<<=1; h|=l>>31; l<<=1; one>>=1;
00313 if (h>=d) { q+= one;h-=d; }
00314 h<<=1; h|=l>>31; l<<=1; one>>=1;
00315 if (h>=d) { q+= one;h-=d; }
00316 h<<=1; h|=l>>31; l<<=1; one>>=1;
00317 if (h>=d) { q+= one;h-=d; }
00318 h<<=1; h|=l>>31; l<<=1; one>>=1;
00319 if (h>=d) { q+= one;h-=d; }
00320 h<<=1; h|=l>>31; l<<=1; one>>=1;
00321 if (h>=d) { q+= one;h-=d; }
00322 h<<=1; h|=l>>31; l<<=1; one>>=1;
00323
00324 *qp=q;
00325 *rp=h>>1;
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338 #ifdef DEBUG
00339 printf(" =%x:%x\n",*qp,*rp);fflush(stdout);
00340 #endif
00341 }
00342
00343 #endif
00344 #endif
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355 inline pointer stretch_big(x, i)
00356 pointer x;
00357 eusinteger_t i;
00358 { pointer bn=x->c.bgnm.bv;
00359 int vlen=vecsize(bn);
00360 pointer newv, oldv;
00361 register int j;
00362
00363 newv=makevector(C_INTVECTOR, vlen+1);
00364 for (j=0; j<vlen; j++) newv->c.ivec.iv[j]=bn->c.ivec.iv[j];
00365 newv->c.ivec.iv[vlen]=i;
00366 pointer_update(x->c.bgnm.bv, newv);
00367 return(newv);
00368 }
00369
00370 pointer copy_big(x)
00371 pointer x;
00372 { pointer y;
00373 register int size, i;
00374 register eusinteger_t *xv, *yv;
00375
00376 size=bigsize(x);
00377 y=makebig(size);
00378 xv=bigvec(x);
00379 yv=bigvec(y);
00380 for (i=0; i<size; i++) yv[i]=xv[i];
00381 #ifdef SAFETY
00382 take_care(y);
00383 #endif
00384 return(y);
00385 }
00386
00387
00388
00389
00390
00391 pointer extend_big(pointer b, int newsize)
00392 { pointer nb;
00393 eusinteger_t *bv, *nbv;
00394 int i,bsize;
00395 nb=makebig(newsize);
00396 bv=bigvec(b); nbv=bigvec(nb); bsize=bigsize(b);
00397 for (i=0; i<bsize; i++) nbv[i]=bv[i];
00398 if (big_sign(b)<0) {
00399 for (i=bsize; i<newsize; i++) nbv[i]= MASK;
00400 nbv[newsize-1] = -1;}
00401 return(nb); }
00402
00403
00404
00405
00406
00407 eusinteger_t big_zerop(x)
00408 pointer x;
00409 { register eusinteger_t *xv;
00410 register int i, size;
00411 xv=bigvec(x); size=bigsize(x);
00412 for (i=0; i<size; i++) if (xv[i] != 0) return(0);
00413 return(1);
00414 }
00415
00416
00417
00418
00419
00420
00421
00422 eusinteger_t big_sign(x)
00423 pointer x;
00424 { pointer y;
00425 y=x->c.bgnm.bv;
00426 return(y->c.ivec.iv[vecsize(y)-1]);
00427 }
00428
00429
00430
00431
00432
00433
00434
00435
00436 int big_compare(x, y)
00437 register pointer x, y;
00438 { register eusinteger_t *xv, *yv;
00439 register int xsize,ysize,i,j=0;
00440 eusinteger_t xsign, ysign;
00441 xv=bigvec(x); yv=bigvec(y);
00442 xsize=bigsize(x); ysize=bigsize(y);
00443
00444 xsign=xv[xsize-1]; ysign=yv[ysize-1];
00445 if (xsign<0 && ysign>=0) return(-1);
00446 else if (xsign>=0 && ysign<0) return(1);
00447 if (xsign<0) {
00448
00449 if (xsize>ysize) return(-1);
00450 else if (xsize<ysize) return(1);
00451 else {
00452 if (xv[xsize-1] < yv[xsize-1]) return(-1);
00453 else if (xv[xsize-1]>yv[xsize-1]) return(1);
00454 for (i=xsize-2; i>=0; i--) {
00455 if (xv[i] < yv[i]) return(1);
00456 else if (xv[i] > yv[i]) return(-1); }
00457 return(0);} }
00458 else {
00459 if (xsize>ysize) return(1);
00460 else if (xsize<ysize) return(-1);
00461 else {
00462 for (i=xsize-1; i>=0; i--) {
00463 if (xv[i] > yv[i]) return(1);
00464 else if (xv[i] < yv[i]) return(-1); }
00465 return(0);} }
00466 }
00467
00468
00469
00470
00471
00472
00473
00474 void complement_big(x)
00475 pointer x;
00476 {
00477 int size, i=0;
00478 eusinteger_t *xv;
00479 size=bigsize(x); xv=bigvec(x);
00480
00481 while (i != size-1) {
00482 if (xv[i] != 0) {
00483 xv[i] = (-xv[i]) & MASK;
00484 goto ONE;}
00485 i++;}
00486 if (xv[i] == ~MASK) {
00487 xv[i] = 0;
00488 stretch_big(x, 1);
00489 xv=bigvec(x); }
00490 else xv[i] = -(xv[i]);
00491 return;
00492
00493 ONE:
00494 for (;;) {
00495 i++;
00496 if (i==size-1) break;
00497 xv[i] = (~xv[i]) & MASK;}
00498 xv[i] = ~xv[i];
00499 return;
00500 }
00501
00502
00503
00504
00505
00506 pointer big_minus(x)
00507 pointer x;
00508 { int size, i;
00509 eusinteger_t *xv, *yv;
00510 pointer y;
00511
00512 size=bigsize(x); xv=bigvec(x);
00513 y = makebig(size);
00514 yv=bigvec(y);
00515 i=0;
00516 while (i<size-1) {
00517 if (xv[i] != 0) {
00518 yv[i]= (-xv[i]) & MASK; goto ONE;}
00519 i++;}
00520 if (xv[i] == MSB) {
00521 yv[i]= 0;
00522 stretch_big(y, 1);
00523 yv=bigvec(y);
00524 }
00525 else yv[i] = -(xv[i]);
00526 #ifdef SAFETY
00527 take_care(y);
00528 #endif
00529 return(y);
00530
00531 ONE:
00532 for (;;) {
00533 i++;
00534 if (i==size-1) break;
00535 yv[i] = (~xv[i]) & MASK;
00536 }
00537 yv[i] = ~xv[i];
00538 #ifdef SAFETY
00539 take_care(y);
00540 #endif
00541 return(y);
00542 }
00543
00544
00545
00546
00547
00548
00549
00550
00551 void add_int_big(c, x)
00552 eusinteger_t c;
00553 pointer x;
00554 { register int size, i=0;
00555 register eusinteger_t *xv;
00556
00557 size=bigsize(x); xv=bigvec(x);
00558 while (i<size-1) {
00559 xv[i] += c;
00560 if (xv[i] < 0) {
00561 c = 1;
00562 xv[i] &= MASK;
00563 i++; }
00564 else return;}
00565
00566 if (xv[i] >= 0) {
00567 xv[i] += c;
00568 if (xv[i] < 0) {
00569 xv[i] &= MASK;
00570 stretch_big(x, 1);
00571 xv=bigvec(x);} }
00572 else xv[i] += c;
00573 return; }
00574
00575
00576
00577
00578
00579
00580
00581 void sub_int_big(c, x)
00582 eusinteger_t c;
00583 pointer x;
00584 { register int size, i=0;
00585 register eusinteger_t *xv;
00586 size=bigsize(x); xv=bigvec(x);
00587 while (i<size-1) {
00588 xv[i] -= c;
00589 if (xv[i] < 0) {
00590 c = 1; xv[i] &= MASK; i++;}
00591 else return;}
00592 if (xv[i] < 0) {
00593 xv[i] -= c;
00594 if (xv[i] >= 0) {
00595 xv[i] &= MASK;
00596 stretch_big(x, -2); } }
00597 else xv[i] -= c;
00598 return;}
00599
00600
00601
00602
00603
00604
00605
00606 void mul_int_big(c, x)
00607 eusinteger_t c;
00608 pointer x;
00609 { int size, i=0;
00610 eusinteger_t *xv, h=0;
00611
00612 size=bigsize(x); xv=bigvec(x);
00613 while (i<size) {
00614 extended_mul(c, xv[i], h, &h, &xv[i]);
00615 i++;}
00616 if (h > 0) stretch_big(x, h);
00617 return;}
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627 eusinteger_t div_int_big(c, x)
00628 eusinteger_t c;
00629 pointer x;
00630 { int i, size;
00631 eusinteger_t *xv, r;
00632 if (c == 0) error(E_USER,(pointer)"divide by zero in bignum div");
00633 size=bigsize(x); xv=bigvec(x);
00634
00635 r = xv[size-1] % c;
00636 xv[size-1] /= c;
00637 i=size-2;
00638 while (i>=0) {
00639 extended_div(c, r, xv[i], &(xv[i]), &r);
00640 i--;}
00641 return(r); }
00642
00643
00644
00645
00646
00647 pointer big_plus(x, y)
00648 pointer x, y;
00649 { pointer z;
00650 int i, xsize, ysize;
00651 eusinteger_t c, zlast, *xv, *yv, *zv;
00652
00653 xsize=bigsize(x); ysize=bigsize(y);
00654 if (xsize<ysize) {
00655 z=x; x=y; y=z;
00656 i=xsize; xsize=ysize; ysize=i;}
00657
00658
00659
00660 xv=bigvec(x); yv=bigvec(y);
00661 z=copy_big(x); zv=bigvec(z);
00662 c=0;
00663 for (i=0; i<ysize; i++) {
00664 zv[i]=zlast=zv[i]+yv[i]+c;
00665 if (zlast<0) { c=1; zv[i] &= MASK;}
00666 else c=0;}
00667 i= ysize - 1;
00668 if (ysize==xsize) {
00669 if (xv[i]>=0 && yv[i]>=0 && zlast<0) stretch_big(z, 1);
00670 else if (xv[i]<0 && yv[i]<0 && zlast>=0) {
00671 stretch_big(z, -2);}
00672 else if (zlast<0) zv[i] |= MSB;
00673 return(z);}
00674 else {
00675 if (yv[ysize-1]>=0) c=1; else c= -1;
00676 while (i<xsize-1) {
00677 if (zlast<0) zv[i] &= MASK;
00678 else return(z);
00679 i++;
00680 zv[i] += c;
00681 zlast=zv[i];}
00682 if (c>=0 && xv[i]>=0 && zv[i]<0) {
00683 zv[i] &= MASK; stretch_big(z, 1);}
00684 else if (c<0 && xv[i]<0 && zv[i]>=0) {
00685 stretch_big(z, -2);}
00686 return(z); }
00687 }
00688
00689
00690
00691
00692
00693
00694 pointer big_times(x, y)
00695 pointer x, y;
00696 { int xsize, ysize, zsize;
00697 eusinteger_t *xv, *yv, *zv;
00698 int i, j, k, m;
00699 pointer z;
00700 eusinteger_t hi, lo;
00701
00702 xsize=bigsize(x); ysize=bigsize(y);
00703 zsize=xsize+ysize;
00704 z=makebig(zsize);
00705 xv=bigvec(x); yv=bigvec(y); zv=bigvec(z);
00706
00707 zv[0]=0;
00708
00709 for (i=0; i<xsize; i++) {
00710 m=xv[i];
00711 hi=0;
00712 for (j=0; j<ysize; j++) {
00713 k=i+j;
00714 extended_mul(m, yv[j], hi, &hi, &lo);
00715 zv[k]+=lo;
00716 if (zv[k] & MSB) { zv[k] &= MASK; hi++; }
00717 if (j==ysize-1) {
00718 while (hi>0) {
00719 k++;
00720 if (k>=xsize+ysize) error(E_USER,(pointer)"bignum mult overflow");
00721 zv[k] += hi;
00722 if (zv[k] & MSB) { zv[k] &= MASK; hi=1; }
00723 else hi=0;}
00724 }
00725 }
00726 }
00727 #ifdef SAFETY
00728 take_care(z);
00729 #endif
00730 return(z);}
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740 void sub_int_big_big(c, x, y)
00741 eusinteger_t c;
00742 pointer x, y;
00743 { int i, j;
00744 int xsize, ysize;
00745 eusinteger_t *xv, *yv;
00746 eusinteger_t hi, lo;
00747
00748 xsize=bigsize(x); ysize=bigsize(y);
00749 xv=bigvec(x); yv=bigvec(y);
00750
00751 hi = 0;
00752 for (i=0;i<xsize;i++) {
00753 extended_mul(c, xv[i], hi, &hi, &lo);
00754 yv[i] -= lo;
00755 if (yv[i] & MSB) {
00756 yv[i] &= MASK;
00757 hi++;
00758 }
00759 if (i==xsize-1) {
00760 while (hi > 0) {
00761 i++;
00762 yv[i] -= hi;
00763 if (yv[i] & MSB) {
00764 yv[i] &= MASK;
00765 hi = 1;}
00766 else break;}
00767 break;}
00768 }
00769 y=normalize_big(y);
00770 }
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782 inline eusinteger_t get_standardizing_factor_and_normalize(x)
00783 pointer x;
00784 { int size, s, newsize;
00785 eusinteger_t *xv, i, j;
00786
00787 size=bigsize(x); xv=bigvec(x);
00788 s=size-1;
00789 newsize=size;
00790 while (xv[newsize-1]==0 && newsize>=0) newsize--;
00791 x->c.bgnm.bv->c.ivec.length=makeint(newsize);
00792 if (newsize==0) return(-1);
00793 for (j = 1, i = xv[newsize-1]; (i <<= 1) >= 0; j <<= 1) ;
00794 return(j);
00795 }
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808 eusinteger_t div_big_big(x, y)
00809 pointer x, y;
00810 { int xsize, ysize, zsize;
00811 eusinteger_t *xv, *yv, *zv;
00812 eusinteger_t q,r;
00813
00814
00815 xsize=bigsize(x); ysize=bigsize(y);
00816 xv=bigvec(x); yv=bigvec(y);
00817
00818 if (xsize != ysize) {
00819 if (yv[xsize] >= xv[xsize-1])
00820 q = MASK-1;
00821
00822
00823
00824 else {
00825 extended_div(xv[xsize-1], yv[xsize], yv[xsize-1], &q, &r);
00826 q -= 2;
00827
00828
00829
00830 } }
00831 else
00832 q = yv[xsize-1] / xv[xsize-1] - 2;
00833 if (q <= 0) q = 0;
00834 else sub_int_big_big(q, x, y);
00835 while (big_compare(x, y) <= 0) {
00836 q++;
00837 sub_int_big_big(1, x, y); }
00838 return(q); }
00839
00840 inline pointer big_quotient_remainder_auxiliary(x, y, i)
00841 pointer x, y;
00842 int i;
00843 { pointer q, qq;
00844 int xsize, ysize;
00845 eusinteger_t *xv, *yv;
00846
00847 xsize=bigsize(x); ysize=bigsize(y);
00848
00849 if (i < 0) return(NULL);
00850 if (i == 0) {
00851 i = div_big_big(y, x);
00852 if (i == 0) return(NULL);
00853 else {
00854 qq = makebig(1);
00855 bigvec(qq)[0] = i;
00856 #ifdef SAFETY
00857 take_care(qq);
00858 #endif
00859 return(qq);
00860 }
00861 }
00862
00863 qq=makebig(i);
00864
00865
00866 #ifdef SAFETY
00867 take_care(qq);
00868 #endif
00869
00870 return(qq);}
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885 void big_quotient_remainder(x0, y0, qp, rp)
00886 pointer x0, y0, *qp, *rp;
00887 {
00888 context *ctx=euscontexts[thr_self()];
00889 pointer x, y;
00890 eusinteger_t i;
00891 int l, m;
00892
00893 x = copy_big(x0);
00894 vpush(x);
00895 y = y0;
00896 i = get_standardizing_factor_and_normalize(y);
00897 mul_int_big(i, x);
00898 mul_int_big(i, y);
00899 l = bigsize(x);
00900 m = bigsize(y);
00901 *qp = big_quotient_remainder_auxiliary(x, y, l - m);
00902 if (*qp == NULL) {
00903 *qp = makebig1(0);
00904 }
00905 div_int_big(i, x);
00906 div_int_big(i, y);
00907 vpop();
00908 *rp = x;
00909 }
00910
00911
00912
00913 int big_length(x)
00914 pointer x;
00915 {
00916 return(bigsize(x));
00917 }
00918
00919 pointer normalize_big(pointer x)
00920 { int size, i, newsize, positive;
00921 eusinteger_t *xv;
00922
00923 size=bigsize(x); xv=bigvec(x);
00924 newsize=size;
00925 positive = gez(xv[size-1]);
00926 for (i=size-1; i>=0; i--) {
00927 if (positive && xv[i]==0) newsize--;
00928 else if (!positive && xv[i]==-1 && size!=1) newsize--;
00929 else break;}
00930 if (newsize != size) {
00931 x->c.bgnm.bv->c.ivec.length=makeint(newsize);
00932 x->c.bgnm.size=makeint(newsize);
00933 if (!positive) xv[newsize-1] |= MSB;}
00934 return(x);
00935 }
00936
00937 pointer normalize_bignum(x)
00938 pointer x;
00939 { eusinteger_t y, msb3;
00940 int bs;
00941
00942 if (!isbignum(x)) return(x);
00943 normalize_big(x);
00944
00945 bs=bigsize(x);
00946 if (bs==0) return(makeint(0));
00947 else if (bs==1) {
00948 y=x->c.bgnm.bv->c.ivec.iv[0];
00949 msb3= (y >> (WORD_SIZE-3)) & 0x7;
00950 if (msb3==0 || msb3==7) return(makeint(y));}
00951 return(x);}
00952
00953 eusfloat_t big_to_float(pointer x)
00954 { int size, i;
00955 eusinteger_t *xv;
00956 double d, e;
00957 size=bigsize(x); xv=bigvec(x);
00958
00959 for (i=0, d = 0.0, e = 1.0; i<size-1; i++) {
00960 d += e * (double)(xv[i]);
00961 e *= (double)((unsigned long)1 << (WORD_SIZE-1));
00962 }
00963 d += e * (double)(xv[size-1]);
00964 return(d);
00965 }
00966
00967 pointer eusfloat_to_big(float f)
00968 { int i, j, exp, bsize, iz;
00969 pointer b;
00970 eusinteger_t *bv;
00971 double sign, z;
00972 extern double frexp(double, int *);
00973 extern double ldexp(double, int);
00974
00975 if (f<0.0) {
00976 f= -f;
00977 if (f<=MAXPOSFIXNUM+1) {
00978 i=f;
00979 return(makeint(-i)); }
00980 else sign= -1;}
00981 else if (f<=MAXPOSFIXNUM) {
00982 i=f;
00983 return(makeint(i));}
00984 else sign= 1;
00985
00986
00987 z=frexp(f, &exp);
00988
00989 z= ldexp(z, 24);
00990 bsize=(exp + WORD_SIZE - 2)/(WORD_SIZE - 1);
00991 b=makebig(bsize);
00992 exp -= 24;
00993 iz=z;
00994 bv=b->c.bgnm.bv->c.ivec.iv;
00995 i=exp / (WORD_SIZE -1);
00996 j=exp % (WORD_SIZE -1);
00997 bv[i]= (iz << j) & MASK;
00998 if (j+24 > 31) bv[i+1]= (iz >> (WORD_SIZE -1 - j)) & MASK;
00999 if (sign<0) complement_big(b);
01000 return(b);
01001 }
01002
01003