Contiki 2.6
|
00001 /**************************************************************** 00002 * 00003 * The author of this software is David M. Gay. 00004 * 00005 * Copyright (c) 1991 by AT&T. 00006 * 00007 * Permission to use, copy, modify, and distribute this software for any 00008 * purpose without fee is hereby granted, provided that this entire notice 00009 * is included in all copies of any software which is or includes a copy 00010 * or modification of this software and in all copies of the supporting 00011 * documentation for such software. 00012 * 00013 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 00014 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY 00015 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 00016 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 00017 * 00018 ***************************************************************/ 00019 00020 /* Please send bug reports to 00021 David M. Gay 00022 AT&T Bell Laboratories, Room 2C-463 00023 600 Mountain Avenue 00024 Murray Hill, NJ 07974-2070 00025 U.S.A. 00026 dmg@research.att.com or research!dmg 00027 */ 00028 00029 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines. 00030 * 00031 * This strtod returns a nearest machine number to the input decimal 00032 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 00033 * broken by the IEEE round-even rule. Otherwise ties are broken by 00034 * biased rounding (add half and chop). 00035 * 00036 * Inspired loosely by William D. Clinger's paper "How to Read Floating 00037 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 00038 * 00039 * Modifications: 00040 * 00041 * 1. We only require IEEE, IBM, or VAX double-precision 00042 * arithmetic (not IEEE double-extended). 00043 * 2. We get by with floating-point arithmetic in a case that 00044 * Clinger missed -- when we're computing d * 10^n 00045 * for a small integer d and the integer n is not too 00046 * much larger than 22 (the maximum integer k for which 00047 * we can represent 10^k exactly), we may be able to 00048 * compute (d*10^k) * 10^(e-k) with just one roundoff. 00049 * 3. Rather than a bit-at-a-time adjustment of the binary 00050 * result in the hard case, we use floating-point 00051 * arithmetic to determine the adjustment to within 00052 * one bit; only in really hard cases do we need to 00053 * compute a second residual. 00054 * 4. Because of 3., we don't need a large table of powers of 10 00055 * for ten-to-e (just some small tables, e.g. of 10^k 00056 * for 0 <= k <= 22). 00057 */ 00058 00059 /* 00060 * #define IEEE_8087 for IEEE-arithmetic machines where the least 00061 * significant byte has the lowest address. 00062 * #define IEEE_MC68k for IEEE-arithmetic machines where the most 00063 * significant byte has the lowest address. 00064 * #define Sudden_Underflow for IEEE-format machines without gradual 00065 * underflow (i.e., that flush to zero on underflow). 00066 * #define IBM for IBM mainframe-style floating-point arithmetic. 00067 * #define VAX for VAX-style floating-point arithmetic. 00068 * #define Unsigned_Shifts if >> does treats its left operand as unsigned. 00069 * #define No_leftright to omit left-right logic in fast floating-point 00070 * computation of dtoa. 00071 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3. 00072 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines 00073 * that use extended-precision instructions to compute rounded 00074 * products and quotients) with IBM. 00075 * #define ROUND_BIASED for IEEE-format with biased rounding. 00076 * #define Inaccurate_Divide for IEEE-format with correctly rounded 00077 * products but inaccurate quotients, e.g., for Intel i860. 00078 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision 00079 * integer arithmetic. Whether this speeds things up or slows things 00080 * down depends on the machine and the number being converted. 00081 */ 00082 00083 00084 00085 #include <_ansi.h> 00086 #include <stdlib.h> 00087 #include <string.h> 00088 #include <reent.h> 00089 #include "small_mprec.h" 00090 00091 /* reent.c knows this value */ 00092 #define _Kmax 15 00093 00094 #if defined (_SMALL_PRINTF) || defined(SMALL_SCANF) 00095 #define SMALL_LIB 00096 #else 00097 #define FULL_LIB 00098 #endif 00099 00100 /* SMALL_LIB is only defined if _SMALL_PRINTF or SMALL_SCANF have been defined 00101 * which means that if you do only specified _SMALL_PRINTF and not SMALL_SCANF 00102 * optimization about call of balloc in small_dtoa_r.c or smal_strtod.c will be 00103 * performed for both printf and scanf 00104 */ 00105 00106 00107 00108 #ifdef FULL_LIB 00109 00110 00111 00112 _Bigint * 00113 _DEFUN (Balloc, (ptr, k), struct _reent *ptr _AND int k) 00114 { 00115 int x; 00116 _Bigint *rv ; 00117 00118 _REENT_CHECK_MP(ptr); 00119 if (_REENT_MP_FREELIST(ptr) == NULL) 00120 { 00121 /* Allocate a list of pointers to the mprec objects */ 00122 _REENT_MP_FREELIST(ptr) = (struct _Bigint **) _calloc_r (ptr, 00123 sizeof (struct _Bigint *), 00124 _Kmax + 1); 00125 if (_REENT_MP_FREELIST(ptr) == NULL) 00126 { 00127 return NULL; 00128 } 00129 } 00130 00131 if ((rv = _REENT_MP_FREELIST(ptr)[k]) != 0) 00132 { 00133 _REENT_MP_FREELIST(ptr)[k] = rv->_next; 00134 } 00135 else 00136 { 00137 x = 1 << k; 00138 /* Allocate an mprec Bigint and stick in in the freelist */ 00139 rv = (_Bigint *) _calloc_r (ptr, 00140 1, 00141 sizeof (_Bigint) + 00142 (x-1) * sizeof(rv->_x)); 00143 if (rv == NULL) return NULL; 00144 rv->_k = k; 00145 rv->_maxwds = x; 00146 } 00147 rv->_sign = rv->_wds = 0; 00148 return rv; 00149 } 00150 00151 void 00152 _DEFUN (Bfree, (ptr, v), struct _reent *ptr _AND _Bigint * v) 00153 { 00154 _REENT_CHECK_MP(ptr); 00155 if (v) 00156 { 00157 v->_next = _REENT_MP_FREELIST(ptr)[v->_k]; 00158 _REENT_MP_FREELIST(ptr)[v->_k] = v; 00159 } 00160 } 00161 00162 00163 _Bigint * 00164 _DEFUN (multadd, (ptr, b, m, a), 00165 struct _reent *ptr _AND 00166 _Bigint * b _AND 00167 int m _AND 00168 int a) 00169 { 00170 int i, wds; 00171 __ULong *x, y; 00172 #ifdef Pack_32 00173 __ULong xi, z; 00174 #endif 00175 _Bigint *b1; 00176 00177 wds = b->_wds; 00178 x = b->_x; 00179 i = 0; 00180 do 00181 { 00182 #ifdef Pack_32 00183 xi = *x; 00184 y = (xi & 0xffff) * m + a; 00185 z = (xi >> 16) * m + (y >> 16); 00186 a = (int) (z >> 16); 00187 *x++ = (z << 16) + (y & 0xffff); 00188 #else 00189 y = *x * m + a; 00190 a = (int) (y >> 16); 00191 *x++ = y & 0xffff; 00192 #endif 00193 } 00194 while (++i < wds); 00195 if (a) 00196 { 00197 if (wds >= b->_maxwds) 00198 { 00199 b1 = Balloc (ptr, b->_k + 1); 00200 Bcopy (b1, b); 00201 Bfree (ptr, b); 00202 b = b1; 00203 } 00204 b->_x[wds++] = a; 00205 b->_wds = wds; 00206 } 00207 return b; 00208 } 00209 00210 _Bigint * 00211 _DEFUN (s2b, (ptr, s, nd0, nd, y9), 00212 struct _reent * ptr _AND 00213 _CONST char *s _AND 00214 int nd0 _AND 00215 int nd _AND 00216 __ULong y9) 00217 { 00218 _Bigint *b; 00219 int i, k; 00220 __Long x, y; 00221 00222 x = (nd + 8) / 9; 00223 for (k = 0, y = 1; x > y; y <<= 1, k++); 00224 #ifdef Pack_32 00225 b = Balloc (ptr, k); 00226 b->_x[0] = y9; 00227 b->_wds = 1; 00228 #else 00229 b = Balloc (ptr, k + 1); 00230 b->_x[0] = y9 & 0xffff; 00231 b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1; 00232 #endif 00233 00234 i = 9; 00235 if (9 < nd0) 00236 { 00237 s += 9; 00238 do 00239 b = multadd (ptr, b, 10, *s++ - '0'); 00240 while (++i < nd0); 00241 s++; 00242 } 00243 else 00244 s += 10; 00245 for (; i < nd; i++) 00246 b = multadd (ptr, b, 10, *s++ - '0'); 00247 return b; 00248 } 00249 00250 int 00251 _DEFUN (hi0bits, 00252 (x), register __ULong x) 00253 { 00254 register int k = 0; 00255 00256 if (!(x & 0xffff0000)) 00257 { 00258 k = 16; 00259 x <<= 16; 00260 } 00261 if (!(x & 0xff000000)) 00262 { 00263 k += 8; 00264 x <<= 8; 00265 } 00266 if (!(x & 0xf0000000)) 00267 { 00268 k += 4; 00269 x <<= 4; 00270 } 00271 if (!(x & 0xc0000000)) 00272 { 00273 k += 2; 00274 x <<= 2; 00275 } 00276 if (!(x & 0x80000000)) 00277 { 00278 k++; 00279 if (!(x & 0x40000000)) 00280 return 32; 00281 } 00282 return k; 00283 } 00284 00285 int 00286 _DEFUN (lo0bits, (y), __ULong *y) 00287 { 00288 register int k; 00289 register __ULong x = *y; 00290 00291 if (x & 7) 00292 { 00293 if (x & 1) 00294 return 0; 00295 if (x & 2) 00296 { 00297 *y = x >> 1; 00298 return 1; 00299 } 00300 *y = x >> 2; 00301 return 2; 00302 } 00303 k = 0; 00304 if (!(x & 0xffff)) 00305 { 00306 k = 16; 00307 x >>= 16; 00308 } 00309 if (!(x & 0xff)) 00310 { 00311 k += 8; 00312 x >>= 8; 00313 } 00314 if (!(x & 0xf)) 00315 { 00316 k += 4; 00317 x >>= 4; 00318 } 00319 if (!(x & 0x3)) 00320 { 00321 k += 2; 00322 x >>= 2; 00323 } 00324 if (!(x & 1)) 00325 { 00326 k++; 00327 x >>= 1; 00328 if (!x & 1) 00329 return 32; 00330 } 00331 *y = x; 00332 return k; 00333 } 00334 00335 _Bigint * 00336 _DEFUN (i2b, (ptr, i), struct _reent * ptr _AND int i) 00337 { 00338 _Bigint *b; 00339 00340 b = Balloc (ptr, 1); 00341 b->_x[0] = i; 00342 b->_wds = 1; 00343 return b; 00344 } 00345 00346 _Bigint * 00347 _DEFUN (mult, (ptr, a, b), struct _reent * ptr _AND _Bigint * a _AND _Bigint * b) 00348 { 00349 _Bigint *c; 00350 int k, wa, wb, wc; 00351 __ULong carry, y, z; 00352 __ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 00353 #ifdef Pack_32 00354 __ULong z2; 00355 #endif 00356 00357 if (a->_wds < b->_wds) 00358 { 00359 c = a; 00360 a = b; 00361 b = c; 00362 } 00363 k = a->_k; 00364 wa = a->_wds; 00365 wb = b->_wds; 00366 wc = wa + wb; 00367 if (wc > a->_maxwds) 00368 k++; 00369 c = Balloc (ptr, k); 00370 for (x = c->_x, xa = x + wc; x < xa; x++) 00371 *x = 0; 00372 xa = a->_x; 00373 xae = xa + wa; 00374 xb = b->_x; 00375 xbe = xb + wb; 00376 xc0 = c->_x; 00377 #ifdef Pack_32 00378 for (; xb < xbe; xb++, xc0++) 00379 { 00380 if ((y = *xb & 0xffff) != 0) 00381 { 00382 x = xa; 00383 xc = xc0; 00384 carry = 0; 00385 do 00386 { 00387 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 00388 carry = z >> 16; 00389 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 00390 carry = z2 >> 16; 00391 Storeinc (xc, z2, z); 00392 } 00393 while (x < xae); 00394 *xc = carry; 00395 } 00396 if ((y = *xb >> 16) != 0) 00397 { 00398 x = xa; 00399 xc = xc0; 00400 carry = 0; 00401 z2 = *xc; 00402 do 00403 { 00404 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 00405 carry = z >> 16; 00406 Storeinc (xc, z, z2); 00407 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 00408 carry = z2 >> 16; 00409 } 00410 while (x < xae); 00411 *xc = z2; 00412 } 00413 } 00414 #else 00415 for (; xb < xbe; xc0++) 00416 { 00417 if (y = *xb++) 00418 { 00419 x = xa; 00420 xc = xc0; 00421 carry = 0; 00422 do 00423 { 00424 z = *x++ * y + *xc + carry; 00425 carry = z >> 16; 00426 *xc++ = z & 0xffff; 00427 } 00428 while (x < xae); 00429 *xc = carry; 00430 } 00431 } 00432 #endif 00433 for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc); 00434 c->_wds = wc; 00435 return c; 00436 } 00437 00438 _Bigint * 00439 _DEFUN (pow5mult, 00440 (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k) 00441 { 00442 _Bigint *b1, *p5, *p51; 00443 int i; 00444 static _CONST int p05[3] = {5, 25, 125}; 00445 00446 if ((i = k & 3) != 0) 00447 b = multadd (ptr, b, p05[i - 1], 0); 00448 00449 if (!(k >>= 2)) 00450 return b; 00451 _REENT_CHECK_MP(ptr); 00452 if (!(p5 = _REENT_MP_P5S(ptr))) 00453 { 00454 /* first time */ 00455 p5 = _REENT_MP_P5S(ptr) = i2b (ptr, 625); 00456 p5->_next = 0; 00457 } 00458 for (;;) 00459 { 00460 if (k & 1) 00461 { 00462 b1 = mult (ptr, b, p5); 00463 Bfree (ptr, b); 00464 b = b1; 00465 } 00466 if (!(k >>= 1)) 00467 break; 00468 if (!(p51 = p5->_next)) 00469 { 00470 p51 = p5->_next = mult (ptr, p5, p5); 00471 p51->_next = 0; 00472 } 00473 p5 = p51; 00474 } 00475 return b; 00476 } 00477 00478 _Bigint * 00479 _DEFUN (lshift, (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k) 00480 { 00481 int i, k1, n, n1; 00482 _Bigint *b1; 00483 __ULong *x, *x1, *xe, z; 00484 00485 #ifdef Pack_32 00486 n = k >> 5; 00487 #else 00488 n = k >> 4; 00489 #endif 00490 k1 = b->_k; 00491 n1 = n + b->_wds + 1; 00492 for (i = b->_maxwds; n1 > i; i <<= 1) 00493 k1++; 00494 b1 = Balloc (ptr, k1); 00495 x1 = b1->_x; 00496 for (i = 0; i < n; i++) 00497 *x1++ = 0; 00498 x = b->_x; 00499 xe = x + b->_wds; 00500 #ifdef Pack_32 00501 if (k &= 0x1f) 00502 { 00503 k1 = 32 - k; 00504 z = 0; 00505 do 00506 { 00507 *x1++ = *x << k | z; 00508 z = *x++ >> k1; 00509 } 00510 while (x < xe); 00511 if ((*x1 = z) != 0) 00512 ++n1; 00513 } 00514 #else 00515 if (k &= 0xf) 00516 { 00517 k1 = 16 - k; 00518 z = 0; 00519 do 00520 { 00521 *x1++ = *x << k & 0xffff | z; 00522 z = *x++ >> k1; 00523 } 00524 while (x < xe); 00525 if (*x1 = z) 00526 ++n1; 00527 } 00528 #endif 00529 else 00530 do 00531 *x1++ = *x++; 00532 while (x < xe); 00533 b1->_wds = n1 - 1; 00534 Bfree (ptr, b); 00535 return b1; 00536 } 00537 00538 int 00539 _DEFUN (cmp, (a, b), _Bigint * a _AND _Bigint * b) 00540 { 00541 __ULong *xa, *xa0, *xb, *xb0; 00542 int i, j; 00543 00544 i = a->_wds; 00545 j = b->_wds; 00546 #ifdef DEBUG 00547 if (i > 1 && !a->_x[i - 1]) 00548 Bug ("cmp called with a->_x[a->_wds-1] == 0"); 00549 if (j > 1 && !b->_x[j - 1]) 00550 Bug ("cmp called with b->_x[b->_wds-1] == 0"); 00551 #endif 00552 if (i -= j) 00553 return i; 00554 xa0 = a->_x; 00555 xa = xa0 + j; 00556 xb0 = b->_x; 00557 xb = xb0 + j; 00558 for (;;) 00559 { 00560 if (*--xa != *--xb) 00561 return *xa < *xb ? -1 : 1; 00562 if (xa <= xa0) 00563 break; 00564 } 00565 return 0; 00566 } 00567 00568 _Bigint * 00569 _DEFUN (diff, (ptr, a, b), struct _reent * ptr _AND 00570 _Bigint * a _AND _Bigint * b) 00571 { 00572 _Bigint *c; 00573 int i, wa, wb; 00574 __Long borrow, y; /* We need signed shifts here. */ 00575 __ULong *xa, *xae, *xb, *xbe, *xc; 00576 #ifdef Pack_32 00577 __Long z; 00578 #endif 00579 00580 i = cmp (a, b); 00581 if (!i) 00582 { 00583 c = Balloc (ptr, 0); 00584 c->_wds = 1; 00585 c->_x[0] = 0; 00586 return c; 00587 } 00588 if (i < 0) 00589 { 00590 c = a; 00591 a = b; 00592 b = c; 00593 i = 1; 00594 } 00595 else 00596 i = 0; 00597 c = Balloc (ptr, a->_k); 00598 c->_sign = i; 00599 wa = a->_wds; 00600 xa = a->_x; 00601 xae = xa + wa; 00602 wb = b->_wds; 00603 xb = b->_x; 00604 xbe = xb + wb; 00605 xc = c->_x; 00606 borrow = 0; 00607 #ifdef Pack_32 00608 do 00609 { 00610 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow; 00611 borrow = y >> 16; 00612 Sign_Extend (borrow, y); 00613 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow; 00614 borrow = z >> 16; 00615 Sign_Extend (borrow, z); 00616 Storeinc (xc, z, y); 00617 } 00618 while (xb < xbe); 00619 while (xa < xae) 00620 { 00621 y = (*xa & 0xffff) + borrow; 00622 borrow = y >> 16; 00623 Sign_Extend (borrow, y); 00624 z = (*xa++ >> 16) + borrow; 00625 borrow = z >> 16; 00626 Sign_Extend (borrow, z); 00627 Storeinc (xc, z, y); 00628 } 00629 #else 00630 do 00631 { 00632 y = *xa++ - *xb++ + borrow; 00633 borrow = y >> 16; 00634 Sign_Extend (borrow, y); 00635 *xc++ = y & 0xffff; 00636 } 00637 while (xb < xbe); 00638 while (xa < xae) 00639 { 00640 y = *xa++ + borrow; 00641 borrow = y >> 16; 00642 Sign_Extend (borrow, y); 00643 *xc++ = y & 0xffff; 00644 } 00645 #endif 00646 while (!*--xc) 00647 wa--; 00648 c->_wds = wa; 00649 return c; 00650 } 00651 00652 double 00653 _DEFUN (ulp, (_x), double _x) 00654 { 00655 union double_union x, a; 00656 register __Long L; 00657 00658 x.d = _x; 00659 00660 L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1; 00661 #ifndef Sudden_Underflow 00662 if (L > 0) 00663 { 00664 #endif 00665 #ifdef IBM 00666 L |= Exp_msk1 >> 4; 00667 #endif 00668 word0 (a) = L; 00669 #ifndef _DOUBLE_IS_32BITS 00670 word1 (a) = 0; 00671 #endif 00672 00673 #ifndef Sudden_Underflow 00674 } 00675 else 00676 { 00677 L = -L >> Exp_shift; 00678 if (L < Exp_shift) 00679 { 00680 word0 (a) = 0x80000 >> L; 00681 #ifndef _DOUBLE_IS_32BITS 00682 word1 (a) = 0; 00683 #endif 00684 } 00685 else 00686 { 00687 word0 (a) = 0; 00688 L -= Exp_shift; 00689 #ifndef _DOUBLE_IS_32BITS 00690 word1 (a) = L >= 31 ? 1 : 1 << (31 - L); 00691 #endif 00692 } 00693 } 00694 #endif 00695 return a.d; 00696 } 00697 00698 double 00699 _DEFUN (b2d, (a, e), 00700 _Bigint * a _AND int *e) 00701 { 00702 __ULong *xa, *xa0, w, y, z; 00703 int k; 00704 union double_union d; 00705 #ifdef VAX 00706 __ULong d0, d1; 00707 #else 00708 #define d0 word0(d) 00709 #define d1 word1(d) 00710 #endif 00711 00712 xa0 = a->_x; 00713 xa = xa0 + a->_wds; 00714 y = *--xa; 00715 #ifdef DEBUG 00716 if (!y) 00717 Bug ("zero y in b2d"); 00718 #endif 00719 k = hi0bits (y); 00720 *e = 32 - k; 00721 #ifdef Pack_32 00722 if (k < Ebits) 00723 { 00724 d0 = Exp_1 | y >> (Ebits - k); 00725 w = xa > xa0 ? *--xa : 0; 00726 #ifndef _DOUBLE_IS_32BITS 00727 d1 = y << ((32 - Ebits) + k) | w >> (Ebits - k); 00728 #endif 00729 goto ret_d; 00730 } 00731 z = xa > xa0 ? *--xa : 0; 00732 if (k -= Ebits) 00733 { 00734 d0 = Exp_1 | y << k | z >> (32 - k); 00735 y = xa > xa0 ? *--xa : 0; 00736 #ifndef _DOUBLE_IS_32BITS 00737 d1 = z << k | y >> (32 - k); 00738 #endif 00739 } 00740 else 00741 { 00742 d0 = Exp_1 | y; 00743 #ifndef _DOUBLE_IS_32BITS 00744 d1 = z; 00745 #endif 00746 } 00747 #else 00748 if (k < Ebits + 16) 00749 { 00750 z = xa > xa0 ? *--xa : 0; 00751 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 00752 w = xa > xa0 ? *--xa : 0; 00753 y = xa > xa0 ? *--xa : 0; 00754 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 00755 goto ret_d; 00756 } 00757 z = xa > xa0 ? *--xa : 0; 00758 w = xa > xa0 ? *--xa : 0; 00759 k -= Ebits + 16; 00760 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 00761 y = xa > xa0 ? *--xa : 0; 00762 d1 = w << k + 16 | y << k; 00763 #endif 00764 ret_d: 00765 #ifdef VAX 00766 word0 (d) = d0 >> 16 | d0 << 16; 00767 word1 (d) = d1 >> 16 | d1 << 16; 00768 #else 00769 #undef d0 00770 #undef d1 00771 #endif 00772 return d.d; 00773 } 00774 00775 _Bigint * 00776 _DEFUN (d2b, 00777 (ptr, _d, e, bits), 00778 struct _reent * ptr _AND 00779 double _d _AND 00780 int *e _AND 00781 int *bits) 00782 00783 { 00784 union double_union d; 00785 _Bigint *b; 00786 int de, i, k; 00787 __ULong *x, y, z; 00788 #ifdef VAX 00789 __ULong d0, d1; 00790 #endif 00791 d.d = _d; 00792 #ifdef VAX 00793 d0 = word0 (d) >> 16 | word0 (d) << 16; 00794 d1 = word1 (d) >> 16 | word1 (d) << 16; 00795 #else 00796 #define d0 word0(d) 00797 #define d1 word1(d) 00798 d.d = _d; 00799 #endif 00800 00801 #ifdef Pack_32 00802 b = Balloc (ptr, 1); 00803 #else 00804 b = Balloc (ptr, 2); 00805 #endif 00806 x = b->_x; 00807 00808 z = d0 & Frac_mask; 00809 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 00810 #ifdef Sudden_Underflow 00811 de = (int) (d0 >> Exp_shift); 00812 #ifndef IBM 00813 z |= Exp_msk11; 00814 #endif 00815 #else 00816 if ((de = (int) (d0 >> Exp_shift)) != 0) 00817 z |= Exp_msk1; 00818 #endif 00819 #ifdef Pack_32 00820 #ifndef _DOUBLE_IS_32BITS 00821 if (d1) 00822 { 00823 y = d1; 00824 k = lo0bits (&y); 00825 if (k) 00826 { 00827 x[0] = y | z << (32 - k); 00828 z >>= k; 00829 } 00830 else 00831 x[0] = y; 00832 i = b->_wds = (x[1] = z) ? 2 : 1; 00833 } 00834 else 00835 #endif 00836 { 00837 #ifdef DEBUG 00838 if (!z) 00839 Bug ("Zero passed to d2b"); 00840 #endif 00841 k = lo0bits (&z); 00842 x[0] = z; 00843 i = b->_wds = 1; 00844 #ifndef _DOUBLE_IS_32BITS 00845 k += 32; 00846 #endif 00847 } 00848 #else 00849 if (d1) 00850 { 00851 y = d1; 00852 k = lo0bits (&y); 00853 if (k) 00854 if (k >= 16) 00855 { 00856 x[0] = y | z << 32 - k & 0xffff; 00857 x[1] = z >> k - 16 & 0xffff; 00858 x[2] = z >> k; 00859 i = 2; 00860 } 00861 else 00862 { 00863 x[0] = y & 0xffff; 00864 x[1] = y >> 16 | z << 16 - k & 0xffff; 00865 x[2] = z >> k & 0xffff; 00866 x[3] = z >> k + 16; 00867 i = 3; 00868 } 00869 else 00870 { 00871 x[0] = y & 0xffff; 00872 x[1] = y >> 16; 00873 x[2] = z & 0xffff; 00874 x[3] = z >> 16; 00875 i = 3; 00876 } 00877 } 00878 else 00879 { 00880 #ifdef DEBUG 00881 if (!z) 00882 Bug ("Zero passed to d2b"); 00883 #endif 00884 k = lo0bits (&z); 00885 if (k >= 16) 00886 { 00887 x[0] = z; 00888 i = 0; 00889 } 00890 else 00891 { 00892 x[0] = z & 0xffff; 00893 x[1] = z >> 16; 00894 i = 1; 00895 } 00896 k += 32; 00897 } 00898 while (!x[i]) 00899 --i; 00900 b->_wds = i + 1; 00901 #endif 00902 #ifndef Sudden_Underflow 00903 if (de) 00904 { 00905 #endif 00906 #ifdef IBM 00907 *e = (de - Bias - (P - 1) << 2) + k; 00908 *bits = 4 * P + 8 - k - hi0bits (word0 (d) & Frac_mask); 00909 #else 00910 *e = de - Bias - (P - 1) + k; 00911 *bits = P - k; 00912 #endif 00913 #ifndef Sudden_Underflow 00914 } 00915 else 00916 { 00917 *e = de - Bias - (P - 1) + 1 + k; 00918 #ifdef Pack_32 00919 *bits = 32 * i - hi0bits (x[i - 1]); 00920 #else 00921 *bits = (i + 2) * 16 - hi0bits (x[i]); 00922 #endif 00923 } 00924 #endif 00925 return b; 00926 } 00927 #undef d0 00928 #undef d1 00929 00930 double 00931 _DEFUN (ratio, (a, b), _Bigint * a _AND _Bigint * b) 00932 00933 { 00934 union double_union da, db; 00935 int k, ka, kb; 00936 00937 da.d = b2d (a, &ka); 00938 db.d = b2d (b, &kb); 00939 #ifdef Pack_32 00940 k = ka - kb + 32 * (a->_wds - b->_wds); 00941 #else 00942 k = ka - kb + 16 * (a->_wds - b->_wds); 00943 #endif 00944 #ifdef IBM 00945 if (k > 0) 00946 { 00947 word0 (da) += (k >> 2) * Exp_msk1; 00948 if (k &= 3) 00949 da.d *= 1 << k; 00950 } 00951 else 00952 { 00953 k = -k; 00954 word0 (db) += (k >> 2) * Exp_msk1; 00955 if (k &= 3) 00956 db.d *= 1 << k; 00957 } 00958 #else 00959 if (k > 0) 00960 word0 (da) += k * Exp_msk1; 00961 else 00962 { 00963 k = -k; 00964 word0 (db) += k * Exp_msk1; 00965 } 00966 #endif 00967 return da.d / db.d; 00968 } 00969 00970 00971 _CONST double 00972 tens[] = 00973 { 00974 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 00975 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 00976 1e20, 1e21, 1e22, 1e23, 1e24 00977 00978 }; 00979 00980 #if !defined(_DOUBLE_IS_32BITS) && !defined(__v800) 00981 _CONST double bigtens[] = 00982 {1e16, 1e32, 1e64, 1e128, 1e256}; 00983 00984 _CONST double tinytens[] = 00985 {1e-16, 1e-32, 1e-64, 1e-128, 1e-256}; 00986 #else 00987 _CONST double bigtens[] = 00988 {1e16, 1e32}; 00989 00990 _CONST double tinytens[] = 00991 {1e-16, 1e-32}; 00992 #endif 00993 00994 00995 double 00996 _DEFUN (_mprec_log10, (dig), 00997 int dig) 00998 { 00999 double v = 1.0; 01000 if (dig < 24) 01001 return tens[dig]; 01002 while (dig > 0) 01003 { 01004 v *= 10; 01005 dig--; 01006 } 01007 return v; 01008 } 01009 01010 #endif // SMALL_LIB 01011 01012 01013 /*****************************************************************************************/ 01014 /* FOR SMALL_LIB */ 01015 /*****************************************************************************************/ 01016 01017 01018 #ifdef SMALL_LIB 01019 01020 /* 01021 * For SMALL_LIB, all references to balloc or bfree have been taken off 01022 * so as to avoided call of malloc libraries 01023 * For each function returning a _Bigint * we added to the function a parameter 01024 * of array type. The purpose of this array is to iniatialized pointers of _Bigint 01025 * with the address of the array instead of the adress provided by balloc. 01026 * 01027 */ 01028 01029 01030 01031 _Bigint * 01032 _DEFUN (small_multadd, (ptr, b, m, a, tab), 01033 struct _reent *ptr _AND 01034 _Bigint * b _AND 01035 int m _AND 01036 int a _AND 01037 _Bigint tab[]) 01038 { 01039 int i, wds; 01040 __ULong *x, y; 01041 #ifdef Pack_32 01042 __ULong xi, z; 01043 #endif 01044 _Bigint *b1; 01045 01046 wds = b->_wds; 01047 x = b->_x; 01048 i = 0; 01049 do 01050 { 01051 #ifdef Pack_32 01052 xi = *x; 01053 y = (xi & 0xffff) * m + a; 01054 z = (xi >> 16) * m + (y >> 16); 01055 a = (int) (z >> 16); 01056 *x++ = (z << 16) + (y & 0xffff); 01057 #else 01058 y = *x * m + a; 01059 a = (int) (y >> 16); 01060 *x++ = y & 0xffff; 01061 #endif 01062 } 01063 while (++i < wds); 01064 if (a) 01065 { 01066 if (wds >= b->_maxwds) 01067 { 01068 b1=&tab[0]; 01069 b1->_k=b->_k+1; 01070 b1->_maxwds = (1 <<(b->_k+1)); 01071 b1->_sign=b1->_wds=0; 01072 Bcopy(b1,b); 01073 b=b1; 01074 } 01075 b->_x[wds++] = a; 01076 b->_wds = wds; 01077 } 01078 return b; 01079 } 01080 01081 01082 _Bigint * 01083 _DEFUN (small_s2b, (ptr, s, nd0, nd, y9,tab), 01084 struct _reent * ptr _AND 01085 _CONST char *s _AND 01086 int nd0 _AND 01087 int nd _AND 01088 __ULong y9 _AND 01089 _Bigint tab[]) 01090 { 01091 _Bigint *b; 01092 _Bigint tab_b[50]; 01093 int i, k; 01094 __Long x, y; 01095 01096 x = (nd + 8) / 9; 01097 for (k = 0, y = 1; x > y; y <<= 1, k++); 01098 #ifdef Pack_32 01099 b = &tab[0]; 01100 b->_k=k; 01101 b->_maxwds = 1 <<k; 01102 b->_sign=0; 01103 b->_x[0] = y9; 01104 b->_wds = 1; 01105 #else 01106 b = &tab[0]; 01107 b = &tab[0]; 01108 b->_k=k+1; 01109 b->_mawxds = 1 <<(k+1); 01110 b->_sign=0; 01111 b->_x[0] = y9 & 0xffff; 01112 b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1; 01113 #endif 01114 01115 i = 9; 01116 if (9 < nd0) 01117 { 01118 s += 9; 01119 01120 while (++i < nd0); 01121 s++; 01122 } 01123 else 01124 s += 10; 01125 for (; i < nd; i++) 01126 b = small_multadd (ptr, b, 10, *s++ - '0',tab); 01127 return b; 01128 } 01129 01130 int 01131 _DEFUN (small_hi0bits, 01132 (x), register __ULong x) 01133 { 01134 register int k = 0; 01135 01136 if (!(x & 0xffff0000)) 01137 { 01138 k = 16; 01139 x <<= 16; 01140 } 01141 if (!(x & 0xff000000)) 01142 { 01143 k += 8; 01144 x <<= 8; 01145 } 01146 if (!(x & 0xf0000000)) 01147 { 01148 k += 4; 01149 x <<= 4; 01150 } 01151 if (!(x & 0xc0000000)) 01152 { 01153 k += 2; 01154 x <<= 2; 01155 } 01156 if (!(x & 0x80000000)) 01157 { 01158 k++; 01159 if (!(x & 0x40000000)) 01160 return 32; 01161 } 01162 return k; 01163 } 01164 01165 int 01166 _DEFUN (small_lo0bits, (y), __ULong *y) 01167 { 01168 register int k; 01169 register __ULong x = *y; 01170 01171 if (x & 7) 01172 { 01173 if (x & 1) 01174 return 0; 01175 if (x & 2) 01176 { 01177 *y = x >> 1; 01178 return 1; 01179 } 01180 *y = x >> 2; 01181 return 2; 01182 } 01183 k = 0; 01184 if (!(x & 0xffff)) 01185 { 01186 k = 16; 01187 x >>= 16; 01188 } 01189 if (!(x & 0xff)) 01190 { 01191 k += 8; 01192 x >>= 8; 01193 } 01194 if (!(x & 0xf)) 01195 { 01196 k += 4; 01197 x >>= 4; 01198 } 01199 if (!(x & 0x3)) 01200 { 01201 k += 2; 01202 x >>= 2; 01203 } 01204 if (!(x & 1)) 01205 { 01206 k++; 01207 x >>= 1; 01208 if (!x & 1) 01209 return 32; 01210 } 01211 *y = x; 01212 return k; 01213 } 01214 01215 01216 _Bigint * 01217 _DEFUN (small_i2b, (ptr, i,tab), struct _reent * ptr _AND int i _AND _Bigint tab[]) 01218 { 01219 _Bigint *b; 01220 01221 b=&tab[0]; 01222 b->_k=1; 01223 b->_maxwds = 1 << 1; 01224 b->_sign =0; 01225 b->_x[0] = i; 01226 b->_wds = 1; 01227 return b; 01228 } 01229 01230 01231 01232 _Bigint * 01233 _DEFUN (small_mult, (ptr, a, b,tab), struct _reent * ptr _AND _Bigint * a _AND _Bigint * b _AND _Bigint tab[]) 01234 { 01235 _Bigint *c; 01236 int k, wa, wb, wc; 01237 __ULong carry, y, z; 01238 __ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 01239 #ifdef Pack_32 01240 __ULong z2; 01241 #endif 01242 01243 if (a->_wds < b->_wds) 01244 { 01245 c = a; 01246 a = b; 01247 b = c; 01248 } 01249 k = a->_k; 01250 wa = a->_wds; 01251 wb = b->_wds; 01252 wc = wa + wb; 01253 if (wc > a->_maxwds) 01254 k++; 01255 c=&tab[0]; 01256 c->_k=k; 01257 c->_maxwds = 1 << k; 01258 c->_sign= c->_wds =0; 01259 01260 for (x = c->_x, xa = x + wc; x < xa; x++) 01261 *x = 0; 01262 xa = a->_x; 01263 xae = xa + wa; 01264 xb = b->_x; 01265 xbe = xb + wb; 01266 xc0 = c->_x; 01267 #ifdef Pack_32 01268 for (; xb < xbe; xb++, xc0++) 01269 { 01270 if ((y = *xb & 0xffff) != 0) 01271 { 01272 x = xa; 01273 xc = xc0; 01274 carry = 0; 01275 do 01276 { 01277 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 01278 carry = z >> 16; 01279 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 01280 carry = z2 >> 16; 01281 Storeinc (xc, z2, z); 01282 } 01283 while (x < xae); 01284 *xc = carry; 01285 } 01286 if ((y = *xb >> 16) != 0) 01287 { 01288 x = xa; 01289 xc = xc0; 01290 carry = 0; 01291 z2 = *xc; 01292 do 01293 { 01294 z = (*x & 0xffff) * y + (*xc >> 16) + carry; 01295 carry = z >> 16; 01296 Storeinc (xc, z, z2); 01297 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 01298 carry = z2 >> 16; 01299 } 01300 while (x < xae); 01301 *xc = z2; 01302 } 01303 } 01304 #else 01305 for (; xb < xbe; xc0++) 01306 { 01307 if (y = *xb++) 01308 { 01309 x = xa; 01310 xc = xc0; 01311 carry = 0; 01312 do 01313 { 01314 z = *x++ * y + *xc + carry; 01315 carry = z >> 16; 01316 *xc++ = z & 0xffff; 01317 } 01318 while (x < xae); 01319 *xc = carry; 01320 } 01321 } 01322 #endif 01323 for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc); 01324 c->_wds = wc; 01325 return c; 01326 } 01327 01328 01329 01330 _Bigint * 01331 _DEFUN (small_pow5mult, 01332 (ptr, b, k, tab), struct _reent * ptr _AND _Bigint * b _AND int k _AND _Bigint tab[]) 01333 { 01334 _Bigint *b1, *p5, *p51; 01335 _Bigint tab_p5[50], tab_p5mult[50]; 01336 _Bigint tab_b2[40], tab_b1[40]; 01337 int i; 01338 int k_sauv; 01339 static _CONST int p05[3] = {5, 25, 125}; 01340 if ((i = k & 3) != 0){ 01341 k_sauv=k; 01342 k_sauv >>= 2; 01343 if (!(k_sauv)){ // we anticipate the case !(k>>=2) with return so as 01344 //to provide the good address &tab[0] 01345 b = small_multadd (ptr, b, p05[i - 1], 0,&tab[0]); 01346 return b; 01347 } 01348 else{ 01349 b = small_multadd (ptr, b, p05[i - 1], 0,&tab_b1[0]); //&tab_b1[0] is an temporary address that 01350 // that we provide to b 01351 } 01352 } 01353 if (!(k>>= 2)){ 01354 return b; 01355 } 01356 if (p5 != tab_p5 ){ /* first time */ 01357 p5 =small_i2b (ptr, 625,tab_p5); 01358 p5->_next = 0; 01359 } 01360 01361 for (;;) {// We are in a loop for each variable and call of mult we must provide a 01362 // an address which differs from the current address of the variable 01363 // that is why we add a test on the current address 01364 01365 if (k & 1){ 01366 01367 k_sauv = k ; 01368 if (!(k_sauv >>= 1) ){ // it is the last passage in the loop 01369 // we must provide the good address 01370 b1 = small_mult (ptr, b, p5,&tab[0]); 01371 b = b1; 01372 break; 01373 } 01374 else { 01375 if (b == &tab_b1[0]){ 01376 b1 = small_mult (ptr, b, p5,&tab_b2[0]); 01377 } 01378 else { 01379 b1 = small_mult (ptr, b, p5,&tab_b1[0]); 01380 } 01381 } 01382 b = b1; 01383 } 01384 if (!(k >>= 1)) 01385 break; 01386 if (!(p51 = p5->_next)) 01387 { 01388 if ( p5 == &tab_p5[0] ){ 01389 p51 = p5->_next = small_mult (ptr, p5, p5,&tab_p5mult[0]); 01390 } 01391 else { 01392 p51 = p5->_next = small_mult (ptr, p5, p5,&tab_p5[0]); 01393 } 01394 p51->_next = 0; 01395 } 01396 p5 = p51; 01397 } 01398 return b; 01399 } 01400 01401 01402 _Bigint * 01403 _DEFUN (small_lshift, (ptr, b, k,tab), struct _reent * ptr _AND _Bigint * b _AND int k _AND _Bigint tab[] ) 01404 { 01405 int i, k1, n, n1; 01406 _Bigint *b1; 01407 __ULong *x, *x1, *xe, z; 01408 01409 #ifdef Pack_32 01410 n = k >> 5; 01411 #else 01412 n = k >> 4; 01413 #endif 01414 k1 = b->_k; 01415 n1 = n + b->_wds + 1; 01416 for (i = b->_maxwds; n1 > i; i <<= 1) 01417 k1++; 01418 //b1 = Balloc (ptr, k1); 01419 b1=&tab[0]; 01420 b1->_k=k1; 01421 b1->_maxwds = 1 << k1; 01422 b1->_sign= b1->_wds =0; 01423 x1 = b1->_x; 01424 for (i = 0; i < n; i++) 01425 *x1++ = 0; 01426 x = b->_x; 01427 xe = x + b->_wds; 01428 #ifdef Pack_32 01429 if (k &= 0x1f) 01430 { 01431 k1 = 32 - k; 01432 z = 0; 01433 do 01434 { 01435 *x1++ = *x << k | z; 01436 z = *x++ >> k1; 01437 } 01438 while (x < xe); 01439 if ((*x1 = z) != 0) 01440 ++n1; 01441 } 01442 #else 01443 if (k &= 0xf) 01444 { 01445 k1 = 16 - k; 01446 z = 0; 01447 do 01448 { 01449 *x1++ = *x << k & 0xffff | z; 01450 z = *x++ >> k1; 01451 } 01452 while (x < xe); 01453 if (*x1 = z) 01454 ++n1; 01455 } 01456 #endif 01457 else 01458 do 01459 *x1++ = *x++; 01460 while (x < xe); 01461 b1->_wds = n1 - 1; 01462 //Bfree (ptr, b); 01463 return b1; 01464 } 01465 01466 int 01467 _DEFUN (small_cmp, (a, b), _Bigint * a _AND _Bigint * b) 01468 { 01469 __ULong *xa, *xa0, *xb, *xb0; 01470 int i, j; 01471 01472 i = a->_wds; 01473 j = b->_wds; 01474 #ifdef DEBUG 01475 if (i > 1 && !a->_x[i - 1]) 01476 Bug ("cmp called with a->_x[a->_wds-1] == 0"); 01477 if (j > 1 && !b->_x[j - 1]) 01478 Bug ("cmp called with b->_x[b->_wds-1] == 0"); 01479 #endif 01480 if (i -= j) 01481 return i; 01482 xa0 = a->_x; 01483 xa = xa0 + j; 01484 xb0 = b->_x; 01485 xb = xb0 + j; 01486 for (;;) 01487 { 01488 if (*--xa != *--xb) 01489 return *xa < *xb ? -1 : 1; 01490 if (xa <= xa0) 01491 break; 01492 } 01493 return 0; 01494 } 01495 01496 01497 01498 _Bigint * 01499 _DEFUN (small_diff, (ptr, a, b,tab), struct _reent * ptr _AND 01500 _Bigint * a _AND _Bigint * b _AND _Bigint tab[]) 01501 { 01502 _Bigint *c; 01503 int i, wa, wb; 01504 __Long borrow, y; /* We need signed shifts here. */ 01505 __ULong *xa, *xae, *xb, *xbe, *xc; 01506 #ifdef Pack_32 01507 __Long z; 01508 #endif 01509 01510 i = small_cmp (a, b); 01511 if (!i) 01512 { 01513 c=&tab[0]; 01514 c->_k = 0; 01515 c->_maxwds = (1 << 0) ; 01516 c->_sign = 0 ; 01517 c->_wds = 1; 01518 c->_x[0] = 0; 01519 return c; 01520 } 01521 if (i < 0) 01522 { 01523 c = a; 01524 a = b; 01525 b = c; 01526 i = 1; 01527 } 01528 else 01529 i = 0; 01530 c=&tab[0]; 01531 c->_k = a->_k; 01532 c->_maxwds = (1 << (a->_k)) ; 01533 c->_wds = 0 ; 01534 c->_sign = i; 01535 wa = a->_wds; 01536 xa = a->_x; 01537 xae = xa + wa; 01538 wb = b->_wds; 01539 xb = b->_x; 01540 xbe = xb + wb; 01541 xc = c->_x; 01542 borrow = 0; 01543 #ifdef Pack_32 01544 do 01545 { 01546 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow; 01547 borrow = y >> 16; 01548 Sign_Extend (borrow, y); 01549 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow; 01550 borrow = z >> 16; 01551 Sign_Extend (borrow, z); 01552 Storeinc (xc, z, y); 01553 } 01554 while (xb < xbe); 01555 while (xa < xae) 01556 { 01557 y = (*xa & 0xffff) + borrow; 01558 borrow = y >> 16; 01559 Sign_Extend (borrow, y); 01560 z = (*xa++ >> 16) + borrow; 01561 borrow = z >> 16; 01562 Sign_Extend (borrow, z); 01563 Storeinc (xc, z, y); 01564 } 01565 #else 01566 do 01567 { 01568 y = *xa++ - *xb++ + borrow; 01569 borrow = y >> 16; 01570 Sign_Extend (borrow, y); 01571 *xc++ = y & 0xffff; 01572 } 01573 while (xb < xbe); 01574 while (xa < xae) 01575 { 01576 y = *xa++ + borrow; 01577 borrow = y >> 16; 01578 Sign_Extend (borrow, y); 01579 *xc++ = y & 0xffff; 01580 } 01581 #endif 01582 while (!*--xc) 01583 wa--; 01584 c->_wds = wa; 01585 return c; 01586 } 01587 01588 double 01589 _DEFUN (small_ulp, (_x), double _x) 01590 { 01591 union double_union x, a; 01592 register __Long L; 01593 01594 x.d = _x; 01595 01596 L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1; 01597 #ifndef Sudden_Underflow 01598 if (L > 0) 01599 { 01600 #endif 01601 #ifdef IBM 01602 L |= Exp_msk1 >> 4; 01603 #endif 01604 word0 (a) = L; 01605 #ifndef _DOUBLE_IS_32BITS 01606 word1 (a) = 0; 01607 #endif 01608 01609 #ifndef Sudden_Underflow 01610 } 01611 else 01612 { 01613 L = -L >> Exp_shift; 01614 if (L < Exp_shift) 01615 { 01616 word0 (a) = 0x80000 >> L; 01617 #ifndef _DOUBLE_IS_32BITS 01618 word1 (a) = 0; 01619 #endif 01620 } 01621 else 01622 { 01623 word0 (a) = 0; 01624 L -= Exp_shift; 01625 #ifndef _DOUBLE_IS_32BITS 01626 word1 (a) = L >= 31 ? 1 : 1 << (31 - L); 01627 #endif 01628 } 01629 } 01630 #endif 01631 return a.d; 01632 } 01633 01634 double 01635 _DEFUN (small_b2d, (a, e), 01636 _Bigint * a _AND int *e) 01637 { 01638 __ULong *xa, *xa0, w, y, z; 01639 int k; 01640 union double_union d; 01641 #ifdef VAX 01642 __ULong d0, d1; 01643 #else 01644 #define d0 word0(d) 01645 #define d1 word1(d) 01646 #endif 01647 01648 xa0 = a->_x; 01649 xa = xa0 + a->_wds; 01650 y = *--xa; 01651 #ifdef DEBUG 01652 if (!y) 01653 Bug ("zero y in b2d"); 01654 #endif 01655 k = small_hi0bits (y); 01656 *e = 32 - k; 01657 #ifdef Pack_32 01658 if (k < Ebits) 01659 { 01660 d0 = Exp_1 | y >> (Ebits - k); 01661 w = xa > xa0 ? *--xa : 0; 01662 #ifndef _DOUBLE_IS_32BITS 01663 d1 = y << ((32 - Ebits) + k) | w >> (Ebits - k); 01664 #endif 01665 goto ret_d; 01666 } 01667 z = xa > xa0 ? *--xa : 0; 01668 if (k -= Ebits) 01669 { 01670 d0 = Exp_1 | y << k | z >> (32 - k); 01671 y = xa > xa0 ? *--xa : 0; 01672 #ifndef _DOUBLE_IS_32BITS 01673 d1 = z << k | y >> (32 - k); 01674 #endif 01675 } 01676 else 01677 { 01678 d0 = Exp_1 | y; 01679 #ifndef _DOUBLE_IS_32BITS 01680 d1 = z; 01681 #endif 01682 } 01683 #else 01684 if (k < Ebits + 16) 01685 { 01686 z = xa > xa0 ? *--xa : 0; 01687 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 01688 w = xa > xa0 ? *--xa : 0; 01689 y = xa > xa0 ? *--xa : 0; 01690 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 01691 goto ret_d; 01692 } 01693 z = xa > xa0 ? *--xa : 0; 01694 w = xa > xa0 ? *--xa : 0; 01695 k -= Ebits + 16; 01696 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 01697 y = xa > xa0 ? *--xa : 0; 01698 d1 = w << k + 16 | y << k; 01699 #endif 01700 ret_d: 01701 #ifdef VAX 01702 word0 (d) = d0 >> 16 | d0 << 16; 01703 word1 (d) = d1 >> 16 | d1 << 16; 01704 #else 01705 #undef d0 01706 #undef d1 01707 #endif 01708 return d.d; 01709 } 01710 01711 01712 _Bigint * 01713 _DEFUN (small_d2b, 01714 (ptr, _d, e, bits,tab), 01715 struct _reent * ptr _AND 01716 double _d _AND 01717 int *e _AND 01718 int *bits _AND 01719 _Bigint tab[] 01720 ) 01721 01722 { 01723 union double_union d; 01724 _Bigint *b; 01725 01726 int de, i, k; 01727 __ULong *x, y, z; 01728 #ifdef VAX 01729 __ULong d0, d1; 01730 #endif 01731 d.d = _d; 01732 #ifdef VAX 01733 d0 = word0 (d) >> 16 | word0 (d) << 16; 01734 d1 = word1 (d) >> 16 | word1 (d) << 16; 01735 #else 01736 #define d0 word0(d) 01737 #define d1 word1(d) 01738 d.d = _d; 01739 #endif 01740 01741 #ifdef Pack_32 01742 b=&tab[0]; 01743 b->_k = 1; 01744 b->_maxwds = (1 << 1) ; 01745 b->_sign = b->_wds = 0 ; 01746 01747 #else 01748 b=&tab[0]; 01749 b->_k = 2; 01750 b->_maxwds = (1 << 2) ; 01751 b->_sign = b->_wds = 0 ; 01752 #endif 01753 01754 01755 x = b->_x; 01756 01757 z = d0 & Frac_mask; 01758 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 01759 #ifdef Sudden_Underflow 01760 de = (int) (d0 >> Exp_shift); 01761 #ifndef IBM 01762 z |= Exp_msk11; 01763 #endif 01764 #else 01765 if ((de = (int) (d0 >> Exp_shift)) != 0) 01766 z |= Exp_msk1; 01767 #endif 01768 #ifdef Pack_32 01769 #ifndef _DOUBLE_IS_32BITS 01770 if (d1) 01771 { 01772 y = d1; 01773 k = small_lo0bits (&y); 01774 if (k) 01775 { 01776 x[0] = y | z << (32 - k); 01777 z >>= k; 01778 } 01779 else 01780 x[0] = y; 01781 i = b->_wds = (x[1] = z) ? 2 : 1; 01782 } 01783 else 01784 #endif 01785 { 01786 #ifdef DEBUG 01787 if (!z) 01788 Bug ("Zero passed to d2b"); 01789 #endif 01790 k = small_lo0bits (&z); 01791 x[0] = z; 01792 i = b->_wds = 1; 01793 #ifndef _DOUBLE_IS_32BITS 01794 k += 32; 01795 #endif 01796 } 01797 #else 01798 if (d1) 01799 { 01800 y = d1; 01801 k = small_lo0bits (&y); 01802 if (k) 01803 if (k >= 16) 01804 { 01805 x[0] = y | z << 32 - k & 0xffff; 01806 x[1] = z >> k - 16 & 0xffff; 01807 x[2] = z >> k; 01808 i = 2; 01809 } 01810 else 01811 { 01812 x[0] = y & 0xffff; 01813 x[1] = y >> 16 | z << 16 - k & 0xffff; 01814 x[2] = z >> k & 0xffff; 01815 x[3] = z >> k + 16; 01816 i = 3; 01817 } 01818 else 01819 { 01820 x[0] = y & 0xffff; 01821 x[1] = y >> 16; 01822 x[2] = z & 0xffff; 01823 x[3] = z >> 16; 01824 i = 3; 01825 } 01826 } 01827 else 01828 { 01829 #ifdef DEBUG 01830 if (!z) 01831 Bug ("Zero passed to d2b"); 01832 #endif 01833 k = lo0bits (&z); 01834 if (k >= 16) 01835 { 01836 x[0] = z; 01837 i = 0; 01838 } 01839 else 01840 { 01841 x[0] = z & 0xffff; 01842 x[1] = z >> 16; 01843 i = 1; 01844 } 01845 k += 32; 01846 } 01847 while (!x[i]) 01848 --i; 01849 b->_wds = i + 1; 01850 #endif 01851 #ifndef Sudden_Underflow 01852 if (de) 01853 { 01854 #endif 01855 #ifdef IBM 01856 *e = (de - Bias - (P - 1) << 2) + k; 01857 *bits = 4 * P + 8 - k - small_hi0bits (word0 (d) & Frac_mask); 01858 #else 01859 *e = de - Bias - (P - 1) + k; 01860 *bits = P - k; 01861 #endif 01862 #ifndef Sudden_Underflow 01863 } 01864 else 01865 { 01866 *e = de - Bias - (P - 1) + 1 + k; 01867 #ifdef Pack_32 01868 *bits = 32 * i - small_hi0bits (x[i - 1]); 01869 #else 01870 *bits = (i + 2) * 16 - small_hi0bits (x[i]); 01871 #endif 01872 } 01873 #endif 01874 return b; 01875 } 01876 #undef d0 01877 #undef d1 01878 01879 double 01880 _DEFUN (small_ratio, (a, b), _Bigint * a _AND _Bigint * b) 01881 01882 { 01883 union double_union da, db; 01884 int k, ka, kb; 01885 01886 da.d = small_b2d (a, &ka); 01887 db.d = small_b2d (b, &kb); 01888 #ifdef Pack_32 01889 k = ka - kb + 32 * (a->_wds - b->_wds); 01890 #else 01891 k = ka - kb + 16 * (a->_wds - b->_wds); 01892 #endif 01893 #ifdef IBM 01894 if (k > 0) 01895 { 01896 word0 (da) += (k >> 2) * Exp_msk1; 01897 if (k &= 3) 01898 da.d *= 1 << k; 01899 } 01900 else 01901 { 01902 k = -k; 01903 word0 (db) += (k >> 2) * Exp_msk1; 01904 if (k &= 3) 01905 db.d *= 1 << k; 01906 } 01907 #else 01908 if (k > 0) 01909 word0 (da) += k * Exp_msk1; 01910 else 01911 { 01912 k = -k; 01913 word0 (db) += k * Exp_msk1; 01914 } 01915 #endif 01916 return da.d / db.d; 01917 } 01918 01919 01920 _CONST double 01921 small_tens[] = 01922 { 01923 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 01924 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 01925 1e20, 1e21, 1e22, 1e23, 1e24 01926 01927 }; 01928 01929 #if !defined(_DOUBLE_IS_32BITS) && !defined(__v800) 01930 _CONST double small_bigtens[] = 01931 {1e16, 1e32, 1e64, 1e128, 1e256}; 01932 01933 _CONST double small_tinytens[] = 01934 {1e-16, 1e-32, 1e-64, 1e-128, 1e-256}; 01935 #else 01936 _CONST double small_bigtens[] = 01937 {1e16, 1e32}; 01938 01939 _CONST double small_tinytens[] = 01940 {1e-16, 1e-32}; 01941 #endif 01942 01943 01944 01945 01946 double 01947 _DEFUN (small__mprec_log10, (dig), 01948 int dig) 01949 { 01950 double v = 1.0; 01951 if (dig < 24) 01952 return small_tens[dig]; 01953 while (dig > 0) 01954 { 01955 v *= 10; 01956 dig--; 01957 } 01958 return v; 01959 } 01960 01961 #endif // SMALL_PRINTF 01962 01963