big.c
Go to the documentation of this file.
00001 /*
00002         big.c
00003 
00004         bignum routines for EusLisp
00005         1995 (c) Toshihiro Matsui, Electrotechnical Laboraotory
00006  Skeltal codes are borrowed from KCL's big.c and earith.c.
00007 
00008  EusLisp's bignum is represented by an integer-vector.
00009  A value is placed in a vector from LSB to MSB.
00010  Sign bit is not used in the first to (MSB -1)th words.
00011  Sign is held in the sign bit of the MSB word. 
00012  example: #x1122334455667788 
00013   bv[0]= #x55667788
00014   bv[1]= #x22446688
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   /* WORD_SIZE==32 */
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 /* !(sun4 && Solaris2) */
00205 
00206 /* Compile with -O option */
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");      /* d*q+r */
00228         asm("   call    .umul,2");      /* delayed inst. */
00229         asm("   mov     %i1,%o1"); 
00230         asm("   addcc   %o0,%i2,%o0");
00231         asm("   addx    %o1,  0,%o1");
00232 
00233         asm("   sll     %o1,  2,%o1");  /* ((H<<1) & 0x3fffffff) */
00234         asm("   srl     %o1,  1,%o1");
00235         asm("   srl     %o0, 31,%i2");  /*      | (L>>31)  -> %o1 */
00236         asm("   or      %o1,%i2,%o1");
00237 
00238         asm("   sll     %o0,  1,%o0");  /* L & 0x7ffffffff */
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 /*                      l=l<<1; one=1<<31;
00328                         for(q=0;one;one>>=1)
00329                           { 
00330                             if (h>=d) { q+= one;h-=d; }
00331                             h<<=1; 
00332                             h|=l>>31;
00333                             l<<=1;
00334                           }
00335                         *qp=q;
00336                         *rp=h>>1;
00337 */
00338 #ifdef DEBUG
00339 printf(" =%x:%x\n",*qp,*rp);fflush(stdout);
00340 #endif
00341 }
00342 
00343 #endif /* !Solaris2 */
00344 #endif
00345 
00346 
00347 
00348 /****************************************************************/
00349 /* stretch_big(x, i)
00350 /*  'x' is a bignum and 'i' is integer.
00351 /*  Allocate bigvector whose size is one word bigger than x,
00352 /*  and put 'i' in the MSB word.
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 /* extend_big(b,newsize) allocates a new bignum which has the same value
00388    as 'b' in the 'newsize' word. Sign is copied in the upper extra words.
00389    Newsize should be greater than the original size of b.
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         Big_zerop(x) answers if bignum x is zero or not.
00405         X may be any bignum.
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         Big_sign(x) returns 
00418                 something < 0   if x < -1
00419                 something >= 0  if x > 0.
00420         X may be any bignum.
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 /*  Big_compare(x, y) returns
00430                 -1      if x < y
00431                 0       if x = y
00432                 1       if x > y.
00433     X and y may be any bignum.
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) { /* ysign is also negative*/
00448     /*assume both x and y are normalized*/
00449     if (xsize>ysize) return(-1);
00450     else if (xsize<ysize) return(1);
00451     else { /* xsize=ysize*/
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 {  /* both x and y are positive */
00459     if (xsize>ysize) return(1);
00460     else if (xsize<ysize) return(-1);
00461     else { /* xsize=ysize*/
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         Complement_big(x) destructively takes
00471         the 2's complement of bignum x.
00472         X may be any bignum.
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); /*reload*/ }
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         big_minus(x) returns the 2's complement of bignum x.
00504         X may be any bignum.
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) { /* 0x80000000 */
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         Add_int_big(c, x) destructively adds non-negative int c
00547         to bignum x.
00548         I should be non-negative.
00549         X may be any bignum.
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) { /*  carry  */
00561       c = 1;
00562       xv[i] &= MASK;
00563       i++; }
00564     else return;}       /* no carry propagation to upper bits */
00565 
00566   if (xv[i] >= 0) {     /* MSB */
00567     xv[i] += c;
00568     if (xv[i] < 0) {            /*  overflow  */
00569         xv[i] &= MASK;
00570         stretch_big(x, 1);
00571         xv=bigvec(x);}   }
00572   else xv[i] += c;
00573   return; }
00574 
00575 /*
00576         Sub_int_big(c, x) destructively subtracts non-negative int c
00577         from bignum x.
00578         c should be non-negative.
00579         X may be any bignum.
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) {            /*  borrow  */
00590         c = 1; xv[i] &= MASK; i++;}
00591     else return;}
00592   if (xv[i] < 0) {      /* MSB */
00593     xv[i] -= c;
00594     if (xv[i] >= 0) { /*  if sign changes, underflow  */
00595       xv[i] &= MASK;
00596       stretch_big(x, -2); } }
00597   else xv[i] -= c;
00598   return;}
00599 
00600 /*
00601         Mul_int_big(i, x) destructively multiplies non-negative bignum x
00602         by non-negative int i.
00603         I should be non-negative.
00604         X should be non-negative.
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         Div_int_big(c, x) destructively divides non-negative bignum x
00621         by positive int c.
00622         X will hold the remainder of the division.
00623         Div_int_big(c, x) returns the remainder of the division.
00624         c should be positive.
00625         X should be non-negative.
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   /* divide from MSB */ 
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         Big_plus(x, y) returns the sum of bignum x and bignum y.
00645         X and y may be any bignum.
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) { /*exchange x and y so that x is bigger than y */
00655     z=x; x=y; y=z;
00656     i=xsize; xsize=ysize; ysize=i;}
00657 
00658   /*fprintf(stderr, "big_plus xsize=%d ysize=%d\n", xsize,ysize);*/
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 {  /* xsize>ysize */
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         Big_times(x0, y0) returns the product
00691         of non-negative bignum x0 and non-negative bignum y0.
00692         X0 and y0 should be non-negative.
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         Sub_int_big_big(c, x, y) destructively subtracts c*x from y.
00734         C should be a non-negative int.
00735         X should be a normalized non-negative bignum.
00736         Y should be non-negative bignum and should be such that
00737                 y <= c*x.
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         Get_standardizing_factor_and_normalize(x)
00774         returns the standardizing factor of x.
00775         As a side-effect, x will be normalized.
00776         X should be a positive bignum.
00777         If x is multiplied by the standardizing factor,
00778         the most significant digit of x will be not less than 2**30,
00779         i.e. the most significant bit of that digit will be set.
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         Div_big_big(x0, y0) divides y0 by x0
00799         and destructively places the remainder in y0.
00800         X0 should be a normalized positive bignum,
00801         whose most significant digit is not less than 2**30.
00802         Y0 should be a non-negative bignum,
00803         whose length must be equal to the length of x0
00804         or one bigger than that.
00805         Div_big_big(x0, y0) returns the quotient of the division,
00806         which must be less than 2**31.
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; /*quotient, remainder*/
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         /*  This is the most critical point.  */
00822         /*  The long division will fail,  */
00823         /*  if the quotient can't lie in 31 bits.  */
00824     else {
00825         extended_div(xv[xsize-1], yv[xsize], yv[xsize-1], &q, &r);
00826         q -= 2;
00827                 /*  You have to prove that 2 is enough,  */
00828                 /*  by doing elementary arithmetic,  */
00829                 /*  or refer to Knuth's dictionary.  */
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) { /* x and y are of the same size */
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         Big_quotient_remainder(x0, y0, qp, rp)
00878         sets the quotient and the remainder of the division of x0 by y0
00879         to *qp and *rp respectively.
00880         qp and rp should be place holders on the value stack to protect
00881         from GC.
00882         X0 should be a positive bignum.
00883         Y0 should be a non-negative bignum.
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) { /* 0x1fffffff on 32bit machines */
00982     i=f;
00983     return(makeint(i));}
00984   else sign= 1;
00985 
00986   /* here, f>MAXPOSFIXNUM */
00987   z=frexp(f, &exp);   /* f=z * 2^exp; 0.5<z<=1.0; exp>=30 */
00988                       /* z has at most 22bit or 52bit precision */
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 


euslisp
Author(s): Toshihiro Matsui
autogenerated on Thu Sep 3 2015 10:36:19