+++ /dev/null
-/** @file\r
-\r
- Copyright (c) 2012, Intel Corporation. All rights reserved.<BR>\r
- This program and the accompanying materials are licensed and made available under\r
- the terms and conditions of the BSD License that accompanies this distribution.\r
- The full text of the license may be found at\r
- http://opensource.org/licenses/bsd-license.\r
-\r
- THE PROGRAM IS DISTRIBUTED UNDER THE BSD LICENSE ON AN "AS IS" BASIS,\r
- WITHOUT WARRANTIES OR REPRESENTATIONS OF ANY KIND, EITHER EXPRESS OR IMPLIED.\r
-\r
- *****************************************************************\r
-\r
- The author of this software is David M. Gay.\r
-\r
- Copyright (C) 1998-2001 by Lucent Technologies\r
- All Rights Reserved\r
-\r
- Permission to use, copy, modify, and distribute this software and\r
- its documentation for any purpose and without fee is hereby\r
- granted, provided that the above copyright notice appear in all\r
- copies and that both that the copyright notice and this\r
- permission notice and warranty disclaimer appear in supporting\r
- documentation, and that the name of Lucent or any of its entities\r
- not be used in advertising or publicity pertaining to\r
- distribution of the software without specific, written prior\r
- permission.\r
-\r
- LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,\r
- INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.\r
- IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY\r
- SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES\r
- WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER\r
- IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,\r
- ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF\r
- THIS SOFTWARE.\r
-\r
-\r
- Please send bug reports to David M. Gay (dmg at acm dot org,\r
- with " at " changed at "@" and " dot " changed to ".").\r
-\r
- *****************************************************************\r
-\r
- NetBSD: strtod.c,v 1.4.14.1 2008/04/08 21:10:55 jdc Exp\r
-**/\r
-#include <LibConfig.h>\r
-\r
-#include "gdtoaimp.h"\r
-#ifndef NO_FENV_H\r
-#include <fenv.h>\r
-#endif\r
-\r
-#ifdef USE_LOCALE\r
-#include "locale.h"\r
-#endif\r
-\r
-#ifdef IEEE_Arith\r
-#ifndef NO_IEEE_Scale\r
-#define Avoid_Underflow\r
-#undef tinytens\r
-/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */\r
-/* flag unnecessarily. It leads to a song and dance at the end of strtod. */\r
-static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,\r
- 9007199254740992.e-256\r
- };\r
-#endif\r
-#endif\r
-\r
-#ifdef Honor_FLT_ROUNDS\r
-#define Rounding rounding\r
-#undef Check_FLT_ROUNDS\r
-#define Check_FLT_ROUNDS\r
-#else\r
-#define Rounding Flt_Rounds\r
-#endif\r
-\r
-//#ifndef __HAVE_LONG_DOUBLE\r
-//__strong_alias(_strtold, strtod)\r
-//__weak_alias(strtold, _strtold)\r
-//#endif\r
-\r
-#if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */\r
-// Disable: warning C4700: uninitialized local variable 'xx' used\r
-#pragma warning ( disable : 4700 )\r
-#endif /* defined(_MSC_VER) */\r
-\r
-double\r
-strtod(CONST char *s00, char **se)\r
-{\r
-#ifdef Avoid_Underflow\r
- int scale;\r
- #endif\r
- int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,\r
- e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;\r
- CONST char *s, *s0, *s1;\r
- double aadj, aadj1, adj, rv, rv0;\r
- Long L;\r
- ULong y, z;\r
- Bigint *bb = NULL, *bb1, *bd0;\r
- Bigint *bd = NULL, *bs = NULL, *delta = NULL; /* pacify gcc */\r
-#ifdef SET_INEXACT\r
- int inexact, oldinexact;\r
-#endif\r
-#ifdef Honor_FLT_ROUNDS\r
- int rounding;\r
-#endif\r
-\r
- sign = nz0 = nz = decpt = 0;\r
- dval(rv) = 0.;\r
- for(s = s00;;s++) {\r
- switch(*s) {\r
- case '-':\r
- sign = 1;\r
- /* FALLTHROUGH */\r
- case '+':\r
- if (*++s)\r
- goto break2;\r
- /* FALLTHROUGH */\r
- case 0:\r
- goto ret0;\r
- case '\t':\r
- case '\n':\r
- case '\v':\r
- case '\f':\r
- case '\r':\r
- case ' ':\r
- continue;\r
- default:\r
- goto break2;\r
- }\r
- }\r
- break2:\r
- if (*s == '0') {\r
-#ifndef NO_HEX_FP\r
- {\r
- static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };\r
- Long expt;\r
- ULong bits[2];\r
- switch(s[1]) {\r
- case 'x':\r
- case 'X':\r
- {\r
-#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD)\r
- FPI fpi1 = fpi;\r
- switch(fegetround()) {\r
- case FE_TOWARDZERO: fpi1.rounding = 0; break;\r
- case FE_UPWARD: fpi1.rounding = 2; break;\r
- case FE_DOWNWARD: fpi1.rounding = 3;\r
- }\r
-#else\r
-#endif\r
- switch((i = gethex(&s, &fpi, &expt, &bb, sign)) & STRTOG_Retmask) {\r
- case STRTOG_NoNumber:\r
- s = s00;\r
- sign = 0;\r
- /* FALLTHROUGH */\r
- case STRTOG_Zero:\r
- break;\r
- default:\r
- if (bb) {\r
- copybits(bits, fpi.nbits, bb);\r
- Bfree(bb);\r
- }\r
- ULtod((/* LINTED */(U*)&rv)->L, bits, expt, i);\r
- }}\r
- goto ret;\r
- }\r
- }\r
-#endif\r
- nz0 = 1;\r
- while(*++s == '0') ;\r
- if (!*s)\r
- goto ret;\r
- }\r
- s0 = s;\r
- y = z = 0;\r
- for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)\r
- if (nd < 9)\r
- y = 10*y + c - '0';\r
- else if (nd < 16)\r
- z = 10*z + c - '0';\r
- nd0 = nd;\r
-#ifdef USE_LOCALE\r
- if (c == *localeconv()->decimal_point)\r
-#else\r
- if (c == '.')\r
-#endif\r
- {\r
- decpt = 1;\r
- c = *++s;\r
- if (!nd) {\r
- for(; c == '0'; c = *++s)\r
- nz++;\r
- if (c > '0' && c <= '9') {\r
- s0 = s;\r
- nf += nz;\r
- nz = 0;\r
- goto have_dig;\r
- }\r
- goto dig_done;\r
- }\r
- for(; c >= '0' && c <= '9'; c = *++s) {\r
- have_dig:\r
- nz++;\r
- if (c -= '0') {\r
- nf += nz;\r
- for(i = 1; i < nz; i++)\r
- if (nd++ < 9)\r
- y *= 10;\r
- else if (nd <= DBL_DIG + 1)\r
- z *= 10;\r
- if (nd++ < 9)\r
- y = 10*y + c;\r
- else if (nd <= DBL_DIG + 1)\r
- z = 10*z + c;\r
- nz = 0;\r
- }\r
- }\r
- }\r
- dig_done:\r
- e = 0;\r
- if (c == 'e' || c == 'E') {\r
- if (!nd && !nz && !nz0) {\r
- goto ret0;\r
- }\r
- s00 = s;\r
- esign = 0;\r
- switch(c = *++s) {\r
- case '-':\r
- esign = 1;\r
- /* FALLTHROUGH */\r
- case '+':\r
- c = *++s;\r
- }\r
- if (c >= '0' && c <= '9') {\r
- while(c == '0')\r
- c = *++s;\r
- if (c > '0' && c <= '9') {\r
- L = c - '0';\r
- s1 = s;\r
- while((c = *++s) >= '0' && c <= '9')\r
- L = 10*L + c - '0';\r
- if (s - s1 > 8 || L > 19999)\r
- /* Avoid confusion from exponents\r
- * so large that e might overflow.\r
- */\r
- e = 19999; /* safe for 16 bit ints */\r
- else\r
- e = (int)L;\r
- if (esign)\r
- e = -e;\r
- }\r
- else\r
- e = 0;\r
- }\r
- else\r
- s = s00;\r
- }\r
- if (!nd) {\r
- if (!nz && !nz0) {\r
-#ifdef INFNAN_CHECK\r
- /* Check for Nan and Infinity */\r
-#ifndef No_Hex_NaN\r
- ULong bits[2];\r
- static FPI fpinan = /* only 52 explicit bits */\r
- { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };\r
-#endif // No_Hex_NaN\r
- if (!decpt)\r
- switch(c) {\r
- case 'i':\r
- case 'I':\r
- if (match(&s,"nf")) {\r
- --s;\r
- if (!match(&s,"inity"))\r
- ++s;\r
- word0(rv) = 0x7ff00000;\r
- word1(rv) = 0;\r
- goto ret;\r
- }\r
- break;\r
- case 'n':\r
- case 'N':\r
- if (match(&s, "an")) {\r
-#ifndef No_Hex_NaN\r
- if (*s == '(' /*)*/\r
- && hexnan(&s, &fpinan, bits)\r
- == STRTOG_NaNbits) {\r
- word0(rv) = (UINT32)(0x7ff00000U | bits[1]);\r
- word1(rv) = (UINT32)bits[0];\r
- }\r
- else {\r
-#endif\r
- word0(rv) = NAN_WORD0;\r
- word1(rv) = NAN_WORD1;\r
-#ifndef No_Hex_NaN\r
- }\r
-#endif\r
- goto ret;\r
- }\r
- }\r
-#endif /* INFNAN_CHECK */\r
- ret0:\r
- s = s00;\r
- sign = 0;\r
- }\r
- goto ret;\r
- }\r
- e1 = e -= nf;\r
-\r
- /* Now we have nd0 digits, starting at s0, followed by a\r
- * decimal point, followed by nd-nd0 digits. The number we're\r
- * after is the integer represented by those digits times\r
- * 10**e */\r
-\r
- if (!nd0)\r
- nd0 = nd;\r
- k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;\r
- dval(rv) = (double)y;\r
- if (k > 9) {\r
-#ifdef SET_INEXACT\r
- if (k > DBL_DIG)\r
- oldinexact = get_inexact();\r
-#endif\r
- dval(rv) = tens[k - 9] * dval(rv) + z;\r
- }\r
- bd0 = 0;\r
- if (nd <= DBL_DIG\r
-#ifndef RND_PRODQUOT\r
-#ifndef Honor_FLT_ROUNDS\r
- && Flt_Rounds == 1\r
-#endif\r
-#endif\r
- ) {\r
- if (!e)\r
- goto ret;\r
- if (e > 0) {\r
- if (e <= Ten_pmax) {\r
-#ifdef VAX\r
- goto vax_ovfl_check;\r
-#else\r
-#ifdef Honor_FLT_ROUNDS\r
- /* round correctly FLT_ROUNDS = 2 or 3 */\r
- if (sign) {\r
- rv = -rv;\r
- sign = 0;\r
- }\r
-#endif\r
- /* rv = */ rounded_product(dval(rv), tens[e]);\r
- goto ret;\r
-#endif\r
- }\r
- i = DBL_DIG - nd;\r
- if (e <= Ten_pmax + i) {\r
- /* A fancier test would sometimes let us do\r
- * this for larger i values.\r
- */\r
-#ifdef Honor_FLT_ROUNDS\r
- /* round correctly FLT_ROUNDS = 2 or 3 */\r
- if (sign) {\r
- rv = -rv;\r
- sign = 0;\r
- }\r
-#endif\r
- e -= i;\r
- dval(rv) *= tens[i];\r
-#ifdef VAX\r
- /* VAX exponent range is so narrow we must\r
- * worry about overflow here...\r
- */\r
- vax_ovfl_check:\r
- word0(rv) -= P*Exp_msk1;\r
- /* rv = */ rounded_product(dval(rv), tens[e]);\r
- if ((word0(rv) & Exp_mask)\r
- > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))\r
- goto ovfl;\r
- word0(rv) += P*Exp_msk1;\r
-#else\r
- /* rv = */ rounded_product(dval(rv), tens[e]);\r
-#endif\r
- goto ret;\r
- }\r
- }\r
-#ifndef Inaccurate_Divide\r
- else if (e >= -Ten_pmax) {\r
-#ifdef Honor_FLT_ROUNDS\r
- /* round correctly FLT_ROUNDS = 2 or 3 */\r
- if (sign) {\r
- rv = -rv;\r
- sign = 0;\r
- }\r
-#endif\r
- /* rv = */ rounded_quotient(dval(rv), tens[-e]);\r
- goto ret;\r
- }\r
-#endif\r
- }\r
- e1 += nd - k;\r
-\r
-#ifdef IEEE_Arith\r
-#ifdef SET_INEXACT\r
- inexact = 1;\r
- if (k <= DBL_DIG)\r
- oldinexact = get_inexact();\r
-#endif\r
-#ifdef Avoid_Underflow\r
- scale = 0;\r
-#endif\r
-#ifdef Honor_FLT_ROUNDS\r
- if ((rounding = Flt_Rounds) >= 2) {\r
- if (sign)\r
- rounding = rounding == 2 ? 0 : 2;\r
- else\r
- if (rounding != 2)\r
- rounding = 0;\r
- }\r
-#endif\r
-#endif /*IEEE_Arith*/\r
-\r
- /* Get starting approximation = rv * 10**e1 */\r
-\r
- if (e1 > 0) {\r
- if ( (i = e1 & 15) !=0)\r
- dval(rv) *= tens[i];\r
- if (e1 &= ~15) {\r
- if (e1 > DBL_MAX_10_EXP) {\r
- ovfl:\r
-#ifndef NO_ERRNO\r
- errno = ERANGE;\r
-#endif\r
- /* Can't trust HUGE_VAL */\r
-#ifdef IEEE_Arith\r
-#ifdef Honor_FLT_ROUNDS\r
- switch(rounding) {\r
- case 0: /* toward 0 */\r
- case 3: /* toward -infinity */\r
- word0(rv) = Big0;\r
- word1(rv) = Big1;\r
- break;\r
- default:\r
- word0(rv) = Exp_mask;\r
- word1(rv) = 0;\r
- }\r
-#else /*Honor_FLT_ROUNDS*/\r
- word0(rv) = Exp_mask;\r
- word1(rv) = 0;\r
-#endif /*Honor_FLT_ROUNDS*/\r
-#ifdef SET_INEXACT\r
- /* set overflow bit */\r
- dval(rv0) = 1e300;\r
- dval(rv0) *= dval(rv0);\r
-#endif\r
-#else /*IEEE_Arith*/\r
- word0(rv) = Big0;\r
- word1(rv) = Big1;\r
-#endif /*IEEE_Arith*/\r
- if (bd0)\r
- goto retfree;\r
- goto ret;\r
- }\r
- e1 = (unsigned int)e1 >> 4;\r
- for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1)\r
- if (e1 & 1)\r
- dval(rv) *= bigtens[j];\r
- /* The last multiplication could overflow. */\r
- word0(rv) -= P*Exp_msk1;\r
- dval(rv) *= bigtens[j];\r
- if ((z = word0(rv) & Exp_mask)\r
- > Exp_msk1*(DBL_MAX_EXP+Bias-P))\r
- goto ovfl;\r
- if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {\r
- /* set to largest number */\r
- /* (Can't trust DBL_MAX) */\r
- word0(rv) = Big0;\r
- word1(rv) = Big1;\r
- }\r
- else\r
- word0(rv) += P*Exp_msk1;\r
- }\r
- }\r
- else if (e1 < 0) {\r
- e1 = -e1;\r
- if ( (i = e1 & 15) !=0)\r
- dval(rv) /= tens[i];\r
- if (e1 >>= 4) {\r
- if (e1 >= 1 << n_bigtens)\r
- goto undfl;\r
-#ifdef Avoid_Underflow\r
- if (e1 & Scale_Bit)\r
- scale = 2*P;\r
- for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1)\r
- if (e1 & 1)\r
- dval(rv) *= tinytens[j];\r
- if (scale && (j = 2*P + 1 - (unsigned int)((word0(rv) & Exp_mask)\r
- >> Exp_shift)) > 0) {\r
- /* scaled rv is denormal; zap j low bits */\r
- if (j >= 32) {\r
- word1(rv) = 0;\r
- if (j >= 53)\r
- word0(rv) = (P+2)*Exp_msk1;\r
- else\r
- word0(rv) &= 0xffffffff << (j-32);\r
- }\r
- else\r
- word1(rv) &= 0xffffffff << j;\r
- }\r
-#else\r
- for(j = 0; e1 > 1; j++, e1 >>= 1)\r
- if (e1 & 1)\r
- dval(rv) *= tinytens[j];\r
- /* The last multiplication could underflow. */\r
- dval(rv0) = dval(rv);\r
- dval(rv) *= tinytens[j];\r
- if (!dval(rv)) {\r
- dval(rv) = 2.*dval(rv0);\r
- dval(rv) *= tinytens[j];\r
-#endif\r
- if (!dval(rv)) {\r
- undfl:\r
- dval(rv) = 0.;\r
-#ifndef NO_ERRNO\r
- errno = ERANGE;\r
-#endif\r
- if (bd0)\r
- goto retfree;\r
- goto ret;\r
- }\r
-#ifndef Avoid_Underflow\r
- word0(rv) = Tiny0;\r
- word1(rv) = Tiny1;\r
- /* The refinement below will clean\r
- * this approximation up.\r
- */\r
- }\r
-#endif\r
- }\r
- }\r
-\r
- /* Now the hard part -- adjusting rv to the correct value.*/\r
-\r
- /* Put digits into bd: true value = bd * 10^e */\r
-\r
- bd0 = s2b(s0, nd0, nd, y);\r
- if (bd0 == NULL)\r
- goto ovfl;\r
-\r
- for(;;) {\r
- bd = Balloc(bd0->k);\r
- if (bd == NULL)\r
- goto ovfl;\r
- Bcopy(bd, bd0);\r
- bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */\r
- if (bb == NULL)\r
- goto ovfl;\r
- bs = i2b(1);\r
- if (bs == NULL)\r
- goto ovfl;\r
-\r
- if (e >= 0) {\r
- bb2 = bb5 = 0;\r
- bd2 = bd5 = e;\r
- }\r
- else {\r
- bb2 = bb5 = -e;\r
- bd2 = bd5 = 0;\r
- }\r
- if (bbe >= 0)\r
- bb2 += bbe;\r
- else\r
- bd2 -= bbe;\r
- bs2 = bb2;\r
-#ifdef Honor_FLT_ROUNDS\r
- if (rounding != 1)\r
- bs2++;\r
-#endif\r
-#ifdef Avoid_Underflow\r
- j = bbe - scale;\r
- i = j + bbbits - 1; /* logb(rv) */\r
- if (i < Emin) /* denormal */\r
- j += P - Emin;\r
- else\r
- j = P + 1 - bbbits;\r
-#else /*Avoid_Underflow*/\r
-#ifdef Sudden_Underflow\r
-#ifdef IBM\r
- j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);\r
-#else\r
- j = P + 1 - bbbits;\r
-#endif\r
-#else /*Sudden_Underflow*/\r
- j = bbe;\r
- i = j + bbbits - 1; /* logb(rv) */\r
- if (i < Emin) /* denormal */\r
- j += P - Emin;\r
- else\r
- j = P + 1 - bbbits;\r
-#endif /*Sudden_Underflow*/\r
-#endif /*Avoid_Underflow*/\r
- bb2 += j;\r
- bd2 += j;\r
-#ifdef Avoid_Underflow\r
- bd2 += scale;\r
-#endif\r
- i = bb2 < bd2 ? bb2 : bd2;\r
- if (i > bs2)\r
- i = bs2;\r
- if (i > 0) {\r
- bb2 -= i;\r
- bd2 -= i;\r
- bs2 -= i;\r
- }\r
- if (bb5 > 0) {\r
- bs = pow5mult(bs, bb5);\r
- if (bs == NULL)\r
- goto ovfl;\r
- bb1 = mult(bs, bb);\r
- if (bb1 == NULL)\r
- goto ovfl;\r
- Bfree(bb);\r
- bb = bb1;\r
- }\r
- if (bb2 > 0) {\r
- bb = lshift(bb, bb2);\r
- if (bb == NULL)\r
- goto ovfl;\r
- }\r
- if (bd5 > 0) {\r
- bd = pow5mult(bd, bd5);\r
- if (bd == NULL)\r
- goto ovfl;\r
- }\r
- if (bd2 > 0) {\r
- bd = lshift(bd, bd2);\r
- if (bd == NULL)\r
- goto ovfl;\r
- }\r
- if (bs2 > 0) {\r
- bs = lshift(bs, bs2);\r
- if (bs == NULL)\r
- goto ovfl;\r
- }\r
- delta = diff(bb, bd);\r
- if (delta == NULL)\r
- goto ovfl;\r
- dsign = delta->sign;\r
- delta->sign = 0;\r
- i = cmp(delta, bs);\r
-#ifdef Honor_FLT_ROUNDS\r
- if (rounding != 1) {\r
- if (i < 0) {\r
- /* Error is less than an ulp */\r
- if (!delta->x[0] && delta->wds <= 1) {\r
- /* exact */\r
-#ifdef SET_INEXACT\r
- inexact = 0;\r
-#endif\r
- break;\r
- }\r
- if (rounding) {\r
- if (dsign) {\r
- adj = 1.;\r
- goto apply_adj;\r
- }\r
- }\r
- else if (!dsign) {\r
- adj = -1.;\r
- if (!word1(rv)\r
- && !(word0(rv) & Frac_mask)) {\r
- y = word0(rv) & Exp_mask;\r
-#ifdef Avoid_Underflow\r
- if (!scale || y > 2*P*Exp_msk1)\r
-#else\r
- if (y)\r
-#endif\r
- {\r
- delta = lshift(delta,Log2P);\r
- if (cmp(delta, bs) <= 0)\r
- adj = -0.5;\r
- }\r
- }\r
- apply_adj:\r
-#ifdef Avoid_Underflow\r
- if (scale && (y = word0(rv) & Exp_mask)\r
- <= 2*P*Exp_msk1)\r
- word0(adj) += (2*P+1)*Exp_msk1 - y;\r
-#else\r
-#ifdef Sudden_Underflow\r
- if ((word0(rv) & Exp_mask) <=\r
- P*Exp_msk1) {\r
- word0(rv) += P*Exp_msk1;\r
- dval(rv) += adj*ulp(dval(rv));\r
- word0(rv) -= P*Exp_msk1;\r
- }\r
- else\r
-#endif /*Sudden_Underflow*/\r
-#endif /*Avoid_Underflow*/\r
- dval(rv) += adj*ulp(dval(rv));\r
- }\r
- break;\r
- }\r
- adj = ratio(delta, bs);\r
- if (adj < 1.)\r
- adj = 1.;\r
- if (adj <= 0x7ffffffe) {\r
- /* adj = rounding ? ceil(adj) : floor(adj); */\r
- y = adj;\r
- if (y != adj) {\r
- if (!((rounding>>1) ^ dsign))\r
- y++;\r
- adj = y;\r
- }\r
- }\r
-#ifdef Avoid_Underflow\r
- if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)\r
- word0(adj) += (2*P+1)*Exp_msk1 - y;\r
-#else\r
-#ifdef Sudden_Underflow\r
- if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {\r
- word0(rv) += P*Exp_msk1;\r
- adj *= ulp(dval(rv));\r
- if (dsign)\r
- dval(rv) += adj;\r
- else\r
- dval(rv) -= adj;\r
- word0(rv) -= P*Exp_msk1;\r
- goto cont;\r
- }\r
-#endif /*Sudden_Underflow*/\r
-#endif /*Avoid_Underflow*/\r
- adj *= ulp(dval(rv));\r
- if (dsign)\r
- dval(rv) += adj;\r
- else\r
- dval(rv) -= adj;\r
- goto cont;\r
- }\r
-#endif /*Honor_FLT_ROUNDS*/\r
-\r
- if (i < 0) {\r
- /* Error is less than half an ulp -- check for\r
- * special case of mantissa a power of two.\r
- */\r
- if (dsign || word1(rv) || word0(rv) & Bndry_mask\r
-#ifdef IEEE_Arith\r
-#ifdef Avoid_Underflow\r
- || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1\r
-#else\r
- || (word0(rv) & Exp_mask) <= Exp_msk1\r
-#endif\r
-#endif\r
- ) {\r
-#ifdef SET_INEXACT\r
- if (!delta->x[0] && delta->wds <= 1)\r
- inexact = 0;\r
-#endif\r
- break;\r
- }\r
- if (!delta->x[0] && delta->wds <= 1) {\r
- /* exact result */\r
-#ifdef SET_INEXACT\r
- inexact = 0;\r
-#endif\r
- break;\r
- }\r
- delta = lshift(delta,Log2P);\r
- if (cmp(delta, bs) > 0)\r
- goto drop_down;\r
- break;\r
- }\r
- if (i == 0) {\r
- /* exactly half-way between */\r
- if (dsign) {\r
- if ((word0(rv) & Bndry_mask1) == Bndry_mask1\r
- && word1(rv) == (\r
-#ifdef Avoid_Underflow\r
- (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)\r
- ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :\r
-#endif\r
- 0xffffffff)) {\r
- /*boundary case -- increment exponent*/\r
- word0(rv) = (word0(rv) & Exp_mask)\r
- + Exp_msk1\r
-#ifdef IBM\r
- | Exp_msk1 >> 4\r
-#endif\r
- ;\r
- word1(rv) = 0;\r
-#ifdef Avoid_Underflow\r
- dsign = 0;\r
-#endif\r
- break;\r
- }\r
- }\r
- else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {\r
- drop_down:\r
- /* boundary case -- decrement exponent */\r
-#ifdef Sudden_Underflow /*{{*/\r
- L = word0(rv) & Exp_mask;\r
-#ifdef IBM\r
- if (L < Exp_msk1)\r
-#else\r
-#ifdef Avoid_Underflow\r
- if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))\r
-#else\r
- if (L <= Exp_msk1)\r
-#endif /*Avoid_Underflow*/\r
-#endif /*IBM*/\r
- goto undfl;\r
- L -= Exp_msk1;\r
-#else /*Sudden_Underflow}{*/\r
-#ifdef Avoid_Underflow\r
- if (scale) {\r
- L = word0(rv) & Exp_mask;\r
- if (L <= (2*P+1)*Exp_msk1) {\r
- if (L > (P+2)*Exp_msk1)\r
- /* round even ==> */\r
- /* accept rv */\r
- break;\r
- /* rv = smallest denormal */\r
- goto undfl;\r
- }\r
- }\r
-#endif /*Avoid_Underflow*/\r
- L = (word0(rv) & Exp_mask) - Exp_msk1;\r
-#endif /*Sudden_Underflow}*/\r
- word0(rv) = (UINT32)(L | Bndry_mask1);\r
- word1(rv) = 0xffffffffU;\r
-#ifdef IBM\r
- goto cont;\r
-#else\r
- break;\r
-#endif\r
- }\r
-#ifndef ROUND_BIASED\r
- if (!(word1(rv) & LSB))\r
- break;\r
-#endif\r
- if (dsign)\r
- dval(rv) += ulp(dval(rv));\r
-#ifndef ROUND_BIASED\r
- else {\r
- dval(rv) -= ulp(dval(rv));\r
-#ifndef Sudden_Underflow\r
- if (!dval(rv))\r
- goto undfl;\r
-#endif\r
- }\r
-#ifdef Avoid_Underflow\r
- dsign = 1 - dsign;\r
-#endif\r
-#endif\r
- break;\r
- }\r
- if ((aadj = ratio(delta, bs)) <= 2.) {\r
- if (dsign)\r
- aadj = aadj1 = 1.;\r
- else if (word1(rv) || word0(rv) & Bndry_mask) {\r
-#ifndef Sudden_Underflow\r
- if (word1(rv) == Tiny1 && !word0(rv))\r
- goto undfl;\r
-#endif\r
- aadj = 1.;\r
- aadj1 = -1.;\r
- }\r
- else {\r
- /* special case -- power of FLT_RADIX to be */\r
- /* rounded down... */\r
-\r
- if (aadj < 2./FLT_RADIX)\r
- aadj = 1./FLT_RADIX;\r
- else\r
- aadj *= 0.5;\r
- aadj1 = -aadj;\r
- }\r
- }\r
- else {\r
- aadj *= 0.5;\r
- aadj1 = dsign ? aadj : -aadj;\r
-#ifdef Check_FLT_ROUNDS\r
- switch(Rounding) {\r
- case 2: /* towards +infinity */\r
- aadj1 -= 0.5;\r
- break;\r
- case 0: /* towards 0 */\r
- case 3: /* towards -infinity */\r
- aadj1 += 0.5;\r
- }\r
-#else\r
- if (Flt_Rounds == 0)\r
- aadj1 += 0.5;\r
-#endif /*Check_FLT_ROUNDS*/\r
- }\r
- y = word0(rv) & Exp_mask;\r
-\r
- /* Check for overflow */\r
-\r
- if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {\r
- dval(rv0) = dval(rv);\r
- word0(rv) -= P*Exp_msk1;\r
- adj = aadj1 * ulp(dval(rv));\r
- dval(rv) += adj;\r
- if ((word0(rv) & Exp_mask) >=\r
- Exp_msk1*(DBL_MAX_EXP+Bias-P)) {\r
- if (word0(rv0) == Big0 && word1(rv0) == Big1)\r
- goto ovfl;\r
- word0(rv) = Big0;\r
- word1(rv) = Big1;\r
- goto cont;\r
- }\r
- else\r
- word0(rv) += P*Exp_msk1;\r
- }\r
- else {\r
-#ifdef Avoid_Underflow\r
- if (scale && y <= 2*P*Exp_msk1) {\r
- if (aadj <= 0x7fffffff) {\r
- if ((z = (uint32_t)aadj) == 0)\r
- z = 1;\r
- aadj = (double)z;\r
- aadj1 = dsign ? aadj : -aadj;\r
- }\r
- word0(aadj1) += (UINT32)((2*P+1)*Exp_msk1 - y);\r
- }\r
- adj = aadj1 * ulp(dval(rv));\r
- dval(rv) += adj;\r
-#else\r
-#ifdef Sudden_Underflow\r
- if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {\r
- dval(rv0) = dval(rv);\r
- word0(rv) += P*Exp_msk1;\r
- adj = aadj1 * ulp(dval(rv));\r
- dval(rv) += adj;\r
-#ifdef IBM\r
- if ((word0(rv) & Exp_mask) < P*Exp_msk1)\r
-#else\r
- if ((word0(rv) & Exp_mask) <= P*Exp_msk1)\r
-#endif\r
- {\r
- if (word0(rv0) == Tiny0\r
- && word1(rv0) == Tiny1)\r
- goto undfl;\r
- word0(rv) = Tiny0;\r
- word1(rv) = Tiny1;\r
- goto cont;\r
- }\r
- else\r
- word0(rv) -= P*Exp_msk1;\r
- }\r
- else {\r
- adj = aadj1 * ulp(dval(rv));\r
- dval(rv) += adj;\r
- }\r
-#else /*Sudden_Underflow*/\r
- /* Compute adj so that the IEEE rounding rules will\r
- * correctly round rv + adj in some half-way cases.\r
- * If rv * ulp(rv) is denormalized (i.e.,\r
- * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid\r
- * trouble from bits lost to denormalization;\r
- * example: 1.2e-307 .\r
- */\r
- if (y <= (P-1)*Exp_msk1 && aadj > 1.) {\r
- aadj1 = (double)(int)(aadj + 0.5);\r
- if (!dsign)\r
- aadj1 = -aadj1;\r
- }\r
- adj = aadj1 * ulp(dval(rv));\r
- dval(rv) += adj;\r
-#endif /*Sudden_Underflow*/\r
-#endif /*Avoid_Underflow*/\r
- }\r
- z = word0(rv) & Exp_mask;\r
-#ifndef SET_INEXACT\r
-#ifdef Avoid_Underflow\r
- if (!scale)\r
-#endif\r
- if (y == z) {\r
- /* Can we stop now? */\r
- L = (Long)aadj;\r
- aadj -= L;\r
- /* The tolerances below are conservative. */\r
- if (dsign || word1(rv) || word0(rv) & Bndry_mask) {\r
- if (aadj < .4999999 || aadj > .5000001)\r
- break;\r
- }\r
- else if (aadj < .4999999/FLT_RADIX)\r
- break;\r
- }\r
-#endif\r
- cont:\r
- Bfree(bb);\r
- Bfree(bd);\r
- Bfree(bs);\r
- Bfree(delta);\r
- }\r
-#ifdef SET_INEXACT\r
- if (inexact) {\r
- if (!oldinexact) {\r
- word0(rv0) = Exp_1 + (70 << Exp_shift);\r
- word1(rv0) = 0;\r
- dval(rv0) += 1.;\r
- }\r
- }\r
- else if (!oldinexact)\r
- clear_inexact();\r
-#endif\r
-#ifdef Avoid_Underflow\r
- if (scale) {\r
- word0(rv0) = Exp_1 - 2*P*Exp_msk1;\r
- word1(rv0) = 0;\r
- dval(rv) *= dval(rv0);\r
-#ifndef NO_ERRNO\r
- /* try to avoid the bug of testing an 8087 register value */\r
- if (word0(rv) == 0 && word1(rv) == 0)\r
- errno = ERANGE;\r
-#endif\r
- }\r
-#endif /* Avoid_Underflow */\r
-#ifdef SET_INEXACT\r
- if (inexact && !(word0(rv) & Exp_mask)) {\r
- /* set underflow bit */\r
- dval(rv0) = 1e-300;\r
- dval(rv0) *= dval(rv0);\r
- }\r
-#endif\r
- retfree:\r
- Bfree(bb);\r
- Bfree(bd);\r
- Bfree(bs);\r
- Bfree(bd0);\r
- Bfree(delta);\r
- ret:\r
- if (se)\r
- *se = __UNCONST(s);\r
- return sign ? -dval(rv) : dval(rv);\r
-}\r
-\r