Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1. /*
  2. FUNCTION
  3.         <<strtod>>, <<strtof>>---string to double or float
  4.  
  5. INDEX
  6.         strtod
  7. INDEX
  8.         _strtod_r
  9. INDEX
  10.         strtof
  11.  
  12. ANSI_SYNOPSIS
  13.         #include <stdlib.h>
  14.         double strtod(const char *<[str]>, char **<[tail]>);
  15.         float strtof(const char *<[str]>, char **<[tail]>);
  16.  
  17.         double _strtod_r(void *<[reent]>,
  18.                          const char *<[str]>, char **<[tail]>);
  19.  
  20. TRAD_SYNOPSIS
  21.         #include <stdlib.h>
  22.         double strtod(<[str]>,<[tail]>)
  23.         char *<[str]>;
  24.         char **<[tail]>;
  25.  
  26.         float strtof(<[str]>,<[tail]>)
  27.         char *<[str]>;
  28.         char **<[tail]>;
  29.  
  30.         double _strtod_r(<[reent]>,<[str]>,<[tail]>)
  31.         char *<[reent]>;
  32.         char *<[str]>;
  33.         char **<[tail]>;
  34.  
  35. DESCRIPTION
  36.         The function <<strtod>> parses the character string <[str]>,
  37.         producing a substring which can be converted to a double
  38.         value.  The substring converted is the longest initial
  39.         subsequence of <[str]>, beginning with the first
  40.         non-whitespace character, that has one of these formats:
  41.         .[+|-]<[digits]>[.[<[digits]>]][(e|E)[+|-]<[digits]>]
  42.         .[+|-].<[digits]>[(e|E)[+|-]<[digits]>]
  43.         .[+|-](i|I)(n|N)(f|F)[(i|I)(n|N)(i|I)(t|T)(y|Y)]
  44.         .[+|-](n|N)(a|A)(n|N)[<(>[<[hexdigits]>]<)>]
  45.         .[+|-]0(x|X)<[hexdigits]>[.[<[hexdigits]>]][(p|P)[+|-]<[digits]>]
  46.         .[+|-]0(x|X).<[hexdigits]>[(p|P)[+|-]<[digits]>]
  47.         The substring contains no characters if <[str]> is empty, consists
  48.         entirely of whitespace, or if the first non-whitespace
  49.         character is something other than <<+>>, <<->>, <<.>>, or a
  50.         digit, and cannot be parsed as infinity or NaN. If the platform
  51.         does not support NaN, then NaN is treated as an empty substring.
  52.         If the substring is empty, no conversion is done, and
  53.         the value of <[str]> is stored in <<*<[tail]>>>.  Otherwise,
  54.         the substring is converted, and a pointer to the final string
  55.         (which will contain at least the terminating null character of
  56.         <[str]>) is stored in <<*<[tail]>>>.  If you want no
  57.         assignment to <<*<[tail]>>>, pass a null pointer as <[tail]>.
  58.         <<strtof>> is identical to <<strtod>> except for its return type.
  59.  
  60.         This implementation returns the nearest machine number to the
  61.         input decimal string.  Ties are broken by using the IEEE
  62.         round-even rule.  However, <<strtof>> is currently subject to
  63.         double rounding errors.
  64.  
  65.         The alternate function <<_strtod_r>> is a reentrant version.
  66.         The extra argument <[reent]> is a pointer to a reentrancy structure.
  67.  
  68. RETURNS
  69.         <<strtod>> returns the converted substring value, if any.  If
  70.         no conversion could be performed, 0 is returned.  If the
  71.         correct value is out of the range of representable values,
  72.         plus or minus <<HUGE_VAL>> is returned, and <<ERANGE>> is
  73.         stored in errno. If the correct value would cause underflow, 0
  74.         is returned and <<ERANGE>> is stored in errno.
  75.  
  76. Supporting OS subroutines required: <<close>>, <<fstat>>, <<isatty>>,
  77. <<lseek>>, <<read>>, <<sbrk>>, <<write>>.
  78. */
  79.  
  80. /****************************************************************
  81.  
  82. The author of this software is David M. Gay.
  83.  
  84. Copyright (C) 1998-2001 by Lucent Technologies
  85. All Rights Reserved
  86.  
  87. Permission to use, copy, modify, and distribute this software and
  88. its documentation for any purpose and without fee is hereby
  89. granted, provided that the above copyright notice appear in all
  90. copies and that both that the copyright notice and this
  91. permission notice and warranty disclaimer appear in supporting
  92. documentation, and that the name of Lucent or any of its entities
  93. not be used in advertising or publicity pertaining to
  94. distribution of the software without specific, written prior
  95. permission.
  96.  
  97. LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
  98. INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
  99. IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
  100. SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  101. WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
  102. IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
  103. ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
  104. THIS SOFTWARE.
  105.  
  106. ****************************************************************/
  107.  
  108. /* Please send bug reports to David M. Gay (dmg at acm dot org,
  109.  * with " at " changed at "@" and " dot " changed to ".").      */
  110.  
  111. /* Original file gdtoa-strtod.c Modified 06-21-2006 by Jeff Johnston to work within newlib.  */
  112.  
  113. #include <_ansi.h>
  114. #include <errno.h>
  115. #include <stdlib.h>
  116. #include <string.h>
  117. #include "mprec.h"
  118. #include "gdtoa.h"
  119. #include "gd_qnan.h"
  120.  
  121. /* #ifndef NO_FENV_H */
  122. /* #include <fenv.h> */
  123. /* #endif */
  124.  
  125. #include "locale.h"
  126.  
  127. #ifdef IEEE_Arith
  128. #ifndef NO_IEEE_Scale
  129. #define Avoid_Underflow
  130. #undef tinytens
  131. /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
  132. /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
  133. static _CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
  134.                 9007199254740992.e-256
  135.                 };
  136. #endif
  137. #endif
  138.  
  139. #ifdef Honor_FLT_ROUNDS
  140. #define Rounding rounding
  141. #undef Check_FLT_ROUNDS
  142. #define Check_FLT_ROUNDS
  143. #else
  144. #define Rounding Flt_Rounds
  145. #endif
  146.  
  147. #ifndef NO_HEX_FP
  148.  
  149. static void
  150. _DEFUN (ULtod, (L, bits, exp, k),
  151.         __ULong *L _AND
  152.         __ULong *bits _AND
  153.         Long exp _AND
  154.         int k)
  155. {
  156.         switch(k & STRTOG_Retmask) {
  157.           case STRTOG_NoNumber:
  158.           case STRTOG_Zero:
  159.                 L[0] = L[1] = 0;
  160.                 break;
  161.  
  162.           case STRTOG_Denormal:
  163.                 L[_1] = bits[0];
  164.                 L[_0] = bits[1];
  165.                 break;
  166.  
  167.           case STRTOG_Normal:
  168.           case STRTOG_NaNbits:
  169.                 L[_1] = bits[0];
  170.                 L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
  171.                 break;
  172.  
  173.           case STRTOG_Infinite:
  174.                 L[_0] = 0x7ff00000;
  175.                 L[_1] = 0;
  176.                 break;
  177.  
  178.           case STRTOG_NaN:
  179.                 L[_0] = 0x7fffffff;
  180.                 L[_1] = (__ULong)-1;
  181.           }
  182.         if (k & STRTOG_Neg)
  183.                 L[_0] |= 0x80000000L;
  184. }
  185. #endif /* !NO_HEX_FP */
  186.  
  187. #ifdef INFNAN_CHECK
  188. static int
  189. _DEFUN (match, (sp, t),
  190.         _CONST char **sp _AND
  191.         char *t)
  192. {
  193.         int c, d;
  194.         _CONST char *s = *sp;
  195.  
  196.         while( (d = *t++) !=0) {
  197.                 if ((c = *++s) >= 'A' && c <= 'Z')
  198.                         c += 'a' - 'A';
  199.                 if (c != d)
  200.                         return 0;
  201.                 }
  202.         *sp = s + 1;
  203.         return 1;
  204. }
  205. #endif /* INFNAN_CHECK */
  206.  
  207.  
  208. double
  209. _DEFUN (_strtod_r, (ptr, s00, se),
  210.         struct _reent *ptr _AND
  211.         _CONST char *s00 _AND
  212.         char **se)
  213. {
  214. #ifdef Avoid_Underflow
  215.         int scale;
  216. #endif
  217.         int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
  218.                  e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
  219.         _CONST char *s, *s0, *s1;
  220.         double aadj, adj;
  221.         U aadj1, rv, rv0;
  222.         Long L;
  223.         __ULong y, z;
  224.         _Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
  225. #ifdef SET_INEXACT
  226.         int inexact, oldinexact;
  227. #endif
  228. #ifdef Honor_FLT_ROUNDS
  229.         int rounding;
  230. #endif
  231.  
  232.         delta = bs = bd = NULL;
  233.         sign = nz0 = nz = decpt = 0;
  234.         dval(rv) = 0.;
  235.         for(s = s00;;s++) switch(*s) {
  236.                 case '-':
  237.                         sign = 1;
  238.                         /* no break */
  239.                 case '+':
  240.                         if (*++s)
  241.                                 goto break2;
  242.                         /* no break */
  243.                 case 0:
  244.                         goto ret0;
  245.                 case '\t':
  246.                 case '\n':
  247.                 case '\v':
  248.                 case '\f':
  249.                 case '\r':
  250.                 case ' ':
  251.                         continue;
  252.                 default:
  253.                         goto break2;
  254.                 }
  255.  break2:
  256.         if (*s == '0') {
  257. #ifndef NO_HEX_FP
  258.                 {
  259.                 static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
  260.                 Long exp;
  261.                 __ULong bits[2];
  262.                 switch(s[1]) {
  263.                   case 'x':
  264.                   case 'X':
  265.                         /* If the number is not hex, then the parse of
  266.                            0 is still valid.  */
  267.                         s00 = s + 1;
  268.                         {
  269. #if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD)
  270.                         FPI fpi1 = fpi;
  271.                         switch(fegetround()) {
  272.                           case FE_TOWARDZERO:   fpi1.rounding = 0; break;
  273.                           case FE_UPWARD:       fpi1.rounding = 2; break;
  274.                           case FE_DOWNWARD:     fpi1.rounding = 3;
  275.                           }
  276. #else
  277. #define fpi1 fpi
  278. #endif
  279.                         switch((i = gethex(ptr, &s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
  280.                           case STRTOG_NoNumber:
  281.                                 s = s00;
  282.                           case STRTOG_Zero:
  283.                                 break;
  284.                           default:
  285.                                 if (bb) {
  286.                                         copybits(bits, fpi.nbits, bb);
  287.                                         Bfree(ptr,bb);
  288.                                         }
  289.                                 ULtod(rv.i, bits, exp, i);
  290.                           }}
  291.                         goto ret;
  292.                   }
  293.                 }
  294. #endif
  295.                 nz0 = 1;
  296.                 while(*++s == '0') ;
  297.                 if (!*s)
  298.                         goto ret;
  299.                 }
  300.         s0 = s;
  301.         y = z = 0;
  302.         for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
  303.                 if (nd < 9)
  304.                         y = 10*y + c - '0';
  305.                 else if (nd < 16)
  306.                         z = 10*z + c - '0';
  307.         nd0 = nd;
  308.         if (strncmp (s, _localeconv_r (ptr)->decimal_point,
  309.                      strlen (_localeconv_r (ptr)->decimal_point)) == 0)
  310.                 {
  311.                 decpt = 1;
  312.                 c = *(s += strlen (_localeconv_r (ptr)->decimal_point));
  313.                 if (!nd) {
  314.                         for(; c == '0'; c = *++s)
  315.                                 nz++;
  316.                         if (c > '0' && c <= '9') {
  317.                                 s0 = s;
  318.                                 nf += nz;
  319.                                 nz = 0;
  320.                                 goto have_dig;
  321.                                 }
  322.                         goto dig_done;
  323.                         }
  324.                 for(; c >= '0' && c <= '9'; c = *++s) {
  325.  have_dig:
  326.                         nz++;
  327.                         if (c -= '0') {
  328.                                 nf += nz;
  329.                                 for(i = 1; i < nz; i++)
  330.                                         if (nd++ < 9)
  331.                                                 y *= 10;
  332.                                         else if (nd <= DBL_DIG + 1)
  333.                                                 z *= 10;
  334.                                 if (nd++ < 9)
  335.                                         y = 10*y + c;
  336.                                 else if (nd <= DBL_DIG + 1)
  337.                                         z = 10*z + c;
  338.                                 nz = 0;
  339.                                 }
  340.                         }
  341.                 }
  342.  dig_done:
  343.         e = 0;
  344.         if (c == 'e' || c == 'E') {
  345.                 if (!nd && !nz && !nz0) {
  346.                         goto ret0;
  347.                         }
  348.                 s00 = s;
  349.                 esign = 0;
  350.                 switch(c = *++s) {
  351.                         case '-':
  352.                                 esign = 1;
  353.                         case '+':
  354.                                 c = *++s;
  355.                         }
  356.                 if (c >= '0' && c <= '9') {
  357.                         while(c == '0')
  358.                                 c = *++s;
  359.                         if (c > '0' && c <= '9') {
  360.                                 L = c - '0';
  361.                                 s1 = s;
  362.                                 while((c = *++s) >= '0' && c <= '9')
  363.                                         L = 10*L + c - '0';
  364.                                 if (s - s1 > 8 || L > 19999)
  365.                                         /* Avoid confusion from exponents
  366.                                          * so large that e might overflow.
  367.                                          */
  368.                                         e = 19999; /* safe for 16 bit ints */
  369.                                 else
  370.                                         e = (int)L;
  371.                                 if (esign)
  372.                                         e = -e;
  373.                                 }
  374.                         else
  375.                                 e = 0;
  376.                         }
  377.                 else
  378.                         s = s00;
  379.                 }
  380.         if (!nd) {
  381.                 if (!nz && !nz0) {
  382. #ifdef INFNAN_CHECK
  383.                         /* Check for Nan and Infinity */
  384.                         __ULong bits[2];
  385.                         static FPI fpinan =     /* only 52 explicit bits */
  386.                                 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
  387.                         if (!decpt)
  388.                          switch(c) {
  389.                           case 'i':
  390.                           case 'I':
  391.                                 if (match(&s,"nf")) {
  392.                                         --s;
  393.                                         if (!match(&s,"inity"))
  394.                                                 ++s;
  395.                                         dword0(rv) = 0x7ff00000;
  396. #ifndef _DOUBLE_IS_32BITS
  397.                                         dword1(rv) = 0;
  398. #endif /*!_DOUBLE_IS_32BITS*/
  399.                                         goto ret;
  400.                                         }
  401.                                 break;
  402.                           case 'n':
  403.                           case 'N':
  404.                                 if (match(&s, "an")) {
  405. #ifndef No_Hex_NaN
  406.                                         if (*s == '(' /*)*/
  407.                                          && hexnan(&s, &fpinan, bits)
  408.                                                         == STRTOG_NaNbits) {
  409.                                                 dword0(rv) = 0x7ff00000 | bits[1];
  410. #ifndef _DOUBLE_IS_32BITS
  411.                                                 dword1(rv) = bits[0];
  412. #endif /*!_DOUBLE_IS_32BITS*/
  413.                                                 }
  414.                                         else {
  415. #endif
  416.                                                 dword0(rv) = NAN_WORD0;
  417. #ifndef _DOUBLE_IS_32BITS
  418.                                                 dword1(rv) = NAN_WORD1;
  419. #endif /*!_DOUBLE_IS_32BITS*/
  420. #ifndef No_Hex_NaN
  421.                                                 }
  422. #endif
  423.                                         goto ret;
  424.                                         }
  425.                           }
  426. #endif /* INFNAN_CHECK */
  427.  ret0:
  428.                         s = s00;
  429.                         sign = 0;
  430.                         }
  431.                 goto ret;
  432.                 }
  433.         e1 = e -= nf;
  434.  
  435.         /* Now we have nd0 digits, starting at s0, followed by a
  436.          * decimal point, followed by nd-nd0 digits.  The number we're
  437.          * after is the integer represented by those digits times
  438.          * 10**e */
  439.  
  440.         if (!nd0)
  441.                 nd0 = nd;
  442.         k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
  443.         dval(rv) = y;
  444.         if (k > 9) {
  445. #ifdef SET_INEXACT
  446.                 if (k > DBL_DIG)
  447.                         oldinexact = get_inexact();
  448. #endif
  449.                 dval(rv) = tens[k - 9] * dval(rv) + z;
  450.                 }
  451.         bd0 = 0;
  452.         if (nd <= DBL_DIG
  453. #ifndef RND_PRODQUOT
  454. #ifndef Honor_FLT_ROUNDS
  455.                 && Flt_Rounds == 1
  456. #endif
  457. #endif
  458.                         ) {
  459.                 if (!e)
  460.                         goto ret;
  461.                 if (e > 0) {
  462.                         if (e <= Ten_pmax) {
  463. #ifdef VAX
  464.                                 goto vax_ovfl_check;
  465. #else
  466. #ifdef Honor_FLT_ROUNDS
  467.                                 /* round correctly FLT_ROUNDS = 2 or 3 */
  468.                                 if (sign) {
  469.                                         dval(rv) = -dval(rv);
  470.                                         sign = 0;
  471.                                         }
  472. #endif
  473.                                 /* rv = */ rounded_product(dval(rv), tens[e]);
  474.                                 goto ret;
  475. #endif
  476.                                 }
  477.                         i = DBL_DIG - nd;
  478.                         if (e <= Ten_pmax + i) {
  479.                                 /* A fancier test would sometimes let us do
  480.                                  * this for larger i values.
  481.                                  */
  482. #ifdef Honor_FLT_ROUNDS
  483.                                 /* round correctly FLT_ROUNDS = 2 or 3 */
  484.                                 if (sign) {
  485.                                         dval(rv) = -dval(rv);
  486.                                         sign = 0;
  487.                                         }
  488. #endif
  489.                                 e -= i;
  490.                                 dval(rv) *= tens[i];
  491. #ifdef VAX
  492.                                 /* VAX exponent range is so narrow we must
  493.                                  * worry about overflow here...
  494.                                  */
  495.  vax_ovfl_check:
  496.                                 dword0(rv) -= P*Exp_msk1;
  497.                                 /* rv = */ rounded_product(dval(rv), tens[e]);
  498.                                 if ((dword0(rv) & Exp_mask)
  499.                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
  500.                                         goto ovfl;
  501.                                 dword0(rv) += P*Exp_msk1;
  502. #else
  503.                                 /* rv = */ rounded_product(dval(rv), tens[e]);
  504. #endif
  505.                                 goto ret;
  506.                                 }
  507.                         }
  508. #ifndef Inaccurate_Divide
  509.                 else if (e >= -Ten_pmax) {
  510. #ifdef Honor_FLT_ROUNDS
  511.                         /* round correctly FLT_ROUNDS = 2 or 3 */
  512.                         if (sign) {
  513.                                 dval(rv) = -dval(rv);
  514.                                 sign = 0;
  515.                                 }
  516. #endif
  517.                         /* rv = */ rounded_quotient(dval(rv), tens[-e]);
  518.                         goto ret;
  519.                         }
  520. #endif
  521.                 }
  522.         e1 += nd - k;
  523.  
  524. #ifdef IEEE_Arith
  525. #ifdef SET_INEXACT
  526.         inexact = 1;
  527.         if (k <= DBL_DIG)
  528.                 oldinexact = get_inexact();
  529. #endif
  530. #ifdef Avoid_Underflow
  531.         scale = 0;
  532. #endif
  533. #ifdef Honor_FLT_ROUNDS
  534.         if ((rounding = Flt_Rounds) >= 2) {
  535.                 if (sign)
  536.                         rounding = rounding == 2 ? 0 : 2;
  537.                 else
  538.                         if (rounding != 2)
  539.                                 rounding = 0;
  540.                 }
  541. #endif
  542. #endif /*IEEE_Arith*/
  543.  
  544.         /* Get starting approximation = rv * 10**e1 */
  545.  
  546.         if (e1 > 0) {
  547.                 if ( (i = e1 & 15) !=0)
  548.                         dval(rv) *= tens[i];
  549.                 if (e1 &= ~15) {
  550.                         if (e1 > DBL_MAX_10_EXP) {
  551.  ovfl:
  552. #ifndef NO_ERRNO
  553.                                 ptr->_errno = ERANGE;
  554. #endif
  555.                                 /* Can't trust HUGE_VAL */
  556. #ifdef IEEE_Arith
  557. #ifdef Honor_FLT_ROUNDS
  558.                                 switch(rounding) {
  559.                                   case 0: /* toward 0 */
  560.                                   case 3: /* toward -infinity */
  561.                                         dword0(rv) = Big0;
  562. #ifndef _DOUBLE_IS_32BITS
  563.                                         dword1(rv) = Big1;
  564. #endif /*!_DOUBLE_IS_32BITS*/
  565.                                         break;
  566.                                   default:
  567.                                         dword0(rv) = Exp_mask;
  568. #ifndef _DOUBLE_IS_32BITS
  569.                                         dword1(rv) = 0;
  570. #endif /*!_DOUBLE_IS_32BITS*/
  571.                                   }
  572. #else /*Honor_FLT_ROUNDS*/
  573.                                 dword0(rv) = Exp_mask;
  574. #ifndef _DOUBLE_IS_32BITS
  575.                                 dword1(rv) = 0;
  576. #endif /*!_DOUBLE_IS_32BITS*/
  577. #endif /*Honor_FLT_ROUNDS*/
  578. #ifdef SET_INEXACT
  579.                                 /* set overflow bit */
  580.                                 dval(rv0) = 1e300;
  581.                                 dval(rv0) *= dval(rv0);
  582. #endif
  583. #else /*IEEE_Arith*/
  584.                                 dword0(rv) = Big0;
  585. #ifndef _DOUBLE_IS_32BITS
  586.                                 dword1(rv) = Big1;
  587. #endif /*!_DOUBLE_IS_32BITS*/
  588. #endif /*IEEE_Arith*/
  589.                                 if (bd0)
  590.                                         goto retfree;
  591.                                 goto ret;
  592.                                 }
  593.                         e1 >>= 4;
  594.                         for(j = 0; e1 > 1; j++, e1 >>= 1)
  595.                                 if (e1 & 1)
  596.                                         dval(rv) *= bigtens[j];
  597.                 /* The last multiplication could overflow. */
  598.                         dword0(rv) -= P*Exp_msk1;
  599.                         dval(rv) *= bigtens[j];
  600.                         if ((z = dword0(rv) & Exp_mask)
  601.                          > Exp_msk1*(DBL_MAX_EXP+Bias-P))
  602.                                 goto ovfl;
  603.                         if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
  604.                                 /* set to largest number */
  605.                                 /* (Can't trust DBL_MAX) */
  606.                                 dword0(rv) = Big0;
  607. #ifndef _DOUBLE_IS_32BITS
  608.                                 dword1(rv) = Big1;
  609. #endif /*!_DOUBLE_IS_32BITS*/
  610.                                 }
  611.                         else
  612.                                 dword0(rv) += P*Exp_msk1;
  613.                         }
  614.                 }
  615.         else if (e1 < 0) {
  616.                 e1 = -e1;
  617.                 if ( (i = e1 & 15) !=0)
  618.                         dval(rv) /= tens[i];
  619.                 if (e1 >>= 4) {
  620.                         if (e1 >= 1 << n_bigtens)
  621.                                 goto undfl;
  622. #ifdef Avoid_Underflow
  623.                         if (e1 & Scale_Bit)
  624.                                 scale = 2*P;
  625.                         for(j = 0; e1 > 0; j++, e1 >>= 1)
  626.                                 if (e1 & 1)
  627.                                         dval(rv) *= tinytens[j];
  628.                         if (scale && (j = 2*P + 1 - ((dword0(rv) & Exp_mask)
  629.                                                 >> Exp_shift)) > 0) {
  630.                                 /* scaled rv is denormal; zap j low bits */
  631.                                 if (j >= 32) {
  632. #ifndef _DOUBLE_IS_32BITS
  633.                                         dword1(rv) = 0;
  634. #endif /*!_DOUBLE_IS_32BITS*/
  635.                                         if (j >= 53)
  636.                                          dword0(rv) = (P+2)*Exp_msk1;
  637.                                         else
  638.                                          dword0(rv) &= 0xffffffff << (j-32);
  639.                                         }
  640. #ifndef _DOUBLE_IS_32BITS
  641.                                 else
  642.                                         dword1(rv) &= 0xffffffff << j;
  643. #endif /*!_DOUBLE_IS_32BITS*/
  644.                                 }
  645. #else
  646.                         for(j = 0; e1 > 1; j++, e1 >>= 1)
  647.                                 if (e1 & 1)
  648.                                         dval(rv) *= tinytens[j];
  649.                         /* The last multiplication could underflow. */
  650.                         dval(rv0) = dval(rv);
  651.                         dval(rv) *= tinytens[j];
  652.                         if (!dval(rv)) {
  653.                                 dval(rv) = 2.*dval(rv0);
  654.                                 dval(rv) *= tinytens[j];
  655. #endif
  656.                                 if (!dval(rv)) {
  657.  undfl:
  658.                                         dval(rv) = 0.;
  659. #ifndef NO_ERRNO
  660.                                         ptr->_errno = ERANGE;
  661. #endif
  662.                                         if (bd0)
  663.                                                 goto retfree;
  664.                                         goto ret;
  665.                                         }
  666. #ifndef Avoid_Underflow
  667. #ifndef _DOUBLE_IS_32BITS
  668.                                 dword0(rv) = Tiny0;
  669.                                 dword1(rv) = Tiny1;
  670. #else
  671.                                 dword0(rv) = Tiny1;
  672. #endif /*_DOUBLE_IS_32BITS*/
  673.                                 /* The refinement below will clean
  674.                                  * this approximation up.
  675.                                  */
  676.                                 }
  677. #endif
  678.                         }
  679.                 }
  680.  
  681.         /* Now the hard part -- adjusting rv to the correct value.*/
  682.  
  683.         /* Put digits into bd: true value = bd * 10^e */
  684.  
  685.         bd0 = s2b(ptr, s0, nd0, nd, y);
  686.  
  687.         for(;;) {
  688.                 bd = Balloc(ptr,bd0->_k);
  689.                 Bcopy(bd, bd0);
  690.                 bb = d2b(ptr,dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
  691.                 bs = i2b(ptr,1);
  692.  
  693.                 if (e >= 0) {
  694.                         bb2 = bb5 = 0;
  695.                         bd2 = bd5 = e;
  696.                         }
  697.                 else {
  698.                         bb2 = bb5 = -e;
  699.                         bd2 = bd5 = 0;
  700.                         }
  701.                 if (bbe >= 0)
  702.                         bb2 += bbe;
  703.                 else
  704.                         bd2 -= bbe;
  705.                 bs2 = bb2;
  706. #ifdef Honor_FLT_ROUNDS
  707.                 if (rounding != 1)
  708.                         bs2++;
  709. #endif
  710. #ifdef Avoid_Underflow
  711.                 j = bbe - scale;
  712.                 i = j + bbbits - 1;     /* logb(rv) */
  713.                 if (i < Emin)   /* denormal */
  714.                         j += P - Emin;
  715.                 else
  716.                         j = P + 1 - bbbits;
  717. #else /*Avoid_Underflow*/
  718. #ifdef Sudden_Underflow
  719. #ifdef IBM
  720.                 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
  721. #else
  722.                 j = P + 1 - bbbits;
  723. #endif
  724. #else /*Sudden_Underflow*/
  725.                 j = bbe;
  726.                 i = j + bbbits - 1;     /* logb(rv) */
  727.                 if (i < Emin)   /* denormal */
  728.                         j += P - Emin;
  729.                 else
  730.                         j = P + 1 - bbbits;
  731. #endif /*Sudden_Underflow*/
  732. #endif /*Avoid_Underflow*/
  733.                 bb2 += j;
  734.                 bd2 += j;
  735. #ifdef Avoid_Underflow
  736.                 bd2 += scale;
  737. #endif
  738.                 i = bb2 < bd2 ? bb2 : bd2;
  739.                 if (i > bs2)
  740.                         i = bs2;
  741.                 if (i > 0) {
  742.                         bb2 -= i;
  743.                         bd2 -= i;
  744.                         bs2 -= i;
  745.                         }
  746.                 if (bb5 > 0) {
  747.                         bs = pow5mult(ptr, bs, bb5);
  748.                         bb1 = mult(ptr, bs, bb);
  749.                         Bfree(ptr, bb);
  750.                         bb = bb1;
  751.                         }
  752.                 if (bb2 > 0)
  753.                         bb = lshift(ptr, bb, bb2);
  754.                 if (bd5 > 0)
  755.                         bd = pow5mult(ptr, bd, bd5);
  756.                 if (bd2 > 0)
  757.                         bd = lshift(ptr, bd, bd2);
  758.                 if (bs2 > 0)
  759.                         bs = lshift(ptr, bs, bs2);
  760.                 delta = diff(ptr, bb, bd);
  761.                 dsign = delta->_sign;
  762.                 delta->_sign = 0;
  763.                 i = cmp(delta, bs);
  764. #ifdef Honor_FLT_ROUNDS
  765.                 if (rounding != 1) {
  766.                         if (i < 0) {
  767.                                 /* Error is less than an ulp */
  768.                                 if (!delta->_x[0] && delta->_wds <= 1) {
  769.                                         /* exact */
  770. #ifdef SET_INEXACT
  771.                                         inexact = 0;
  772. #endif
  773.                                         break;
  774.                                         }
  775.                                 if (rounding) {
  776.                                         if (dsign) {
  777.                                                 adj = 1.;
  778.                                                 goto apply_adj;
  779.                                                 }
  780.                                         }
  781.                                 else if (!dsign) {
  782.                                         adj = -1.;
  783.                                         if (!dword1(rv)
  784.                                          && !(dword0(rv) & Frac_mask)) {
  785.                                                 y = dword0(rv) & Exp_mask;
  786. #ifdef Avoid_Underflow
  787.                                                 if (!scale || y > 2*P*Exp_msk1)
  788. #else
  789.                                                 if (y)
  790. #endif
  791.                                                   {
  792.                                                   delta = lshift(ptr, delta,Log2P);
  793.                                                   if (cmp(delta, bs) <= 0)
  794.                                                         adj = -0.5;
  795.                                                   }
  796.                                                 }
  797.  apply_adj:
  798. #ifdef Avoid_Underflow
  799.                                         if (scale && (y = dword0(rv) & Exp_mask)
  800.                                                 <= 2*P*Exp_msk1)
  801.                                           dword0(adj) += (2*P+1)*Exp_msk1 - y;
  802. #else
  803. #ifdef Sudden_Underflow
  804.                                         if ((dword0(rv) & Exp_mask) <=
  805.                                                         P*Exp_msk1) {
  806.                                                 dword0(rv) += P*Exp_msk1;
  807.                                                 dval(rv) += adj*ulp(dval(rv));
  808.                                                 dword0(rv) -= P*Exp_msk1;
  809.                                                 }
  810.                                         else
  811. #endif /*Sudden_Underflow*/
  812. #endif /*Avoid_Underflow*/
  813.                                         dval(rv) += adj*ulp(dval(rv));
  814.                                         }
  815.                                 break;
  816.                                 }
  817.                         adj = ratio(delta, bs);
  818.                         if (adj < 1.)
  819.                                 adj = 1.;
  820.                         if (adj <= 0x7ffffffe) {
  821.                                 /* adj = rounding ? ceil(adj) : floor(adj); */
  822.                                 y = adj;
  823.                                 if (y != adj) {
  824.                                         if (!((rounding>>1) ^ dsign))
  825.                                                 y++;
  826.                                         adj = y;
  827.                                         }
  828.                                 }
  829. #ifdef Avoid_Underflow
  830.                         if (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1)
  831.                                 dword0(adj) += (2*P+1)*Exp_msk1 - y;
  832. #else
  833. #ifdef Sudden_Underflow
  834.                         if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
  835.                                 dword0(rv) += P*Exp_msk1;
  836.                                 adj *= ulp(dval(rv));
  837.                                 if (dsign)
  838.                                         dval(rv) += adj;
  839.                                 else
  840.                                         dval(rv) -= adj;
  841.                                 dword0(rv) -= P*Exp_msk1;
  842.                                 goto cont;
  843.                                 }
  844. #endif /*Sudden_Underflow*/
  845. #endif /*Avoid_Underflow*/
  846.                         adj *= ulp(dval(rv));
  847.                         if (dsign)
  848.                                 dval(rv) += adj;
  849.                         else
  850.                                 dval(rv) -= adj;
  851.                         goto cont;
  852.                         }
  853. #endif /*Honor_FLT_ROUNDS*/
  854.  
  855.                 if (i < 0) {
  856.                         /* Error is less than half an ulp -- check for
  857.                          * special case of mantissa a power of two.
  858.                          */
  859.                         if (dsign || dword1(rv) || dword0(rv) & Bndry_mask
  860. #ifdef IEEE_Arith
  861. #ifdef Avoid_Underflow
  862.                          || (dword0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
  863. #else
  864.                          || (dword0(rv) & Exp_mask) <= Exp_msk1
  865. #endif
  866. #endif
  867.                                 ) {
  868. #ifdef SET_INEXACT
  869.                                 if (!delta->x[0] && delta->wds <= 1)
  870.                                         inexact = 0;
  871. #endif
  872.                                 break;
  873.                                 }
  874.                         if (!delta->_x[0] && delta->_wds <= 1) {
  875.                                 /* exact result */
  876. #ifdef SET_INEXACT
  877.                                 inexact = 0;
  878. #endif
  879.                                 break;
  880.                                 }
  881.                         delta = lshift(ptr,delta,Log2P);
  882.                         if (cmp(delta, bs) > 0)
  883.                                 goto drop_down;
  884.                         break;
  885.                         }
  886.                 if (i == 0) {
  887.                         /* exactly half-way between */
  888.                         if (dsign) {
  889.                                 if ((dword0(rv) & Bndry_mask1) == Bndry_mask1
  890.                                  &&  dword1(rv) == (
  891. #ifdef Avoid_Underflow
  892.                         (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1)
  893.                 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
  894. #endif
  895.                                                    0xffffffff)) {
  896.                                         /*boundary case -- increment exponent*/
  897.                                         dword0(rv) = (dword0(rv) & Exp_mask)
  898.                                                 + Exp_msk1
  899. #ifdef IBM
  900.                                                 | Exp_msk1 >> 4
  901. #endif
  902.                                                 ;
  903. #ifndef _DOUBLE_IS_32BITS
  904.                                         dword1(rv) = 0;
  905. #endif /*!_DOUBLE_IS_32BITS*/
  906. #ifdef Avoid_Underflow
  907.                                         dsign = 0;
  908. #endif
  909.                                         break;
  910.                                         }
  911.                                 }
  912.                         else if (!(dword0(rv) & Bndry_mask) && !dword1(rv)) {
  913.  drop_down:
  914.                                 /* boundary case -- decrement exponent */
  915. #ifdef Sudden_Underflow /*{{*/
  916.                                 L = dword0(rv) & Exp_mask;
  917. #ifdef IBM
  918.                                 if (L <  Exp_msk1)
  919. #else
  920. #ifdef Avoid_Underflow
  921.                                 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
  922. #else
  923.                                 if (L <= Exp_msk1)
  924. #endif /*Avoid_Underflow*/
  925. #endif /*IBM*/
  926.                                         goto undfl;
  927.                                 L -= Exp_msk1;
  928. #else /*Sudden_Underflow}{*/
  929. #ifdef Avoid_Underflow
  930.                                 if (scale) {
  931.                                         L = dword0(rv) & Exp_mask;
  932.                                         if (L <= (2*P+1)*Exp_msk1) {
  933.                                                 if (L > (P+2)*Exp_msk1)
  934.                                                         /* round even ==> */
  935.                                                         /* accept rv */
  936.                                                         break;
  937.                                                 /* rv = smallest denormal */
  938.                                                 goto undfl;
  939.                                                 }
  940.                                         }
  941. #endif /*Avoid_Underflow*/
  942.                                 L = (dword0(rv) & Exp_mask) - Exp_msk1;
  943. #endif /*Sudden_Underflow}*/
  944.                                 dword0(rv) = L | Bndry_mask1;
  945. #ifndef _DOUBLE_IS_32BITS
  946.                                 dword1(rv) = 0xffffffff;
  947. #endif /*!_DOUBLE_IS_32BITS*/
  948. #ifdef IBM
  949.                                 goto cont;
  950. #else
  951.                                 break;
  952. #endif
  953.                                 }
  954. #ifndef ROUND_BIASED
  955.                         if (!(dword1(rv) & LSB))
  956.                                 break;
  957. #endif
  958.                         if (dsign)
  959.                                 dval(rv) += ulp(dval(rv));
  960. #ifndef ROUND_BIASED
  961.                         else {
  962.                                 dval(rv) -= ulp(dval(rv));
  963. #ifndef Sudden_Underflow
  964.                                 if (!dval(rv))
  965.                                         goto undfl;
  966. #endif
  967.                                 }
  968. #ifdef Avoid_Underflow
  969.                         dsign = 1 - dsign;
  970. #endif
  971. #endif
  972.                         break;
  973.                         }
  974.                 if ((aadj = ratio(delta, bs)) <= 2.) {
  975.                         if (dsign)
  976.                                 aadj = dval(aadj1) = 1.;
  977.                         else if (dword1(rv) || dword0(rv) & Bndry_mask) {
  978. #ifndef Sudden_Underflow
  979.                                 if (dword1(rv) == Tiny1 && !dword0(rv))
  980.                                         goto undfl;
  981. #endif
  982.                                 aadj = 1.;
  983.                                 dval(aadj1) = -1.;
  984.                                 }
  985.                         else {
  986.                                 /* special case -- power of FLT_RADIX to be */
  987.                                 /* rounded down... */
  988.  
  989.                                 if (aadj < 2./FLT_RADIX)
  990.                                         aadj = 1./FLT_RADIX;
  991.                                 else
  992.                                         aadj *= 0.5;
  993.                                 dval(aadj1) = -aadj;
  994.                                 }
  995.                         }
  996.                 else {
  997.                         aadj *= 0.5;
  998.                         dval(aadj1) = dsign ? aadj : -aadj;
  999. #ifdef Check_FLT_ROUNDS
  1000.                         switch(Rounding) {
  1001.                                 case 2: /* towards +infinity */
  1002.                                         dval(aadj1) -= 0.5;
  1003.                                         break;
  1004.                                 case 0: /* towards 0 */
  1005.                                 case 3: /* towards -infinity */
  1006.                                         dval(aadj1) += 0.5;
  1007.                                 }
  1008. #else
  1009.                         if (Flt_Rounds == 0)
  1010.                                 dval(aadj1) += 0.5;
  1011. #endif /*Check_FLT_ROUNDS*/
  1012.                         }
  1013.                 y = dword0(rv) & Exp_mask;
  1014.  
  1015.                 /* Check for overflow */
  1016.  
  1017.                 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
  1018.                         dval(rv0) = dval(rv);
  1019.                         dword0(rv) -= P*Exp_msk1;
  1020.                         adj = dval(aadj1) * ulp(dval(rv));
  1021.                         dval(rv) += adj;
  1022.                         if ((dword0(rv) & Exp_mask) >=
  1023.                                         Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
  1024.                                 if (dword0(rv0) == Big0 && dword1(rv0) == Big1)
  1025.                                         goto ovfl;
  1026.                                 dword0(rv) = Big0;
  1027. #ifndef _DOUBLE_IS_32BITS
  1028.                                 dword1(rv) = Big1;
  1029. #endif /*!_DOUBLE_IS_32BITS*/
  1030.                                 goto cont;
  1031.                                 }
  1032.                         else
  1033.                                 dword0(rv) += P*Exp_msk1;
  1034.                         }
  1035.                 else {
  1036. #ifdef Avoid_Underflow
  1037.                         if (scale && y <= 2*P*Exp_msk1) {
  1038.                                 if (aadj <= 0x7fffffff) {
  1039.                                         if ((z = aadj) <= 0)
  1040.                                                 z = 1;
  1041.                                         aadj = z;
  1042.                                         dval(aadj1) = dsign ? aadj : -aadj;
  1043.                                         }
  1044.                                 dword0(aadj1) += (2*P+1)*Exp_msk1 - y;
  1045.                                 }
  1046.                         adj = dval(aadj1) * ulp(dval(rv));
  1047.                         dval(rv) += adj;
  1048. #else
  1049. #ifdef Sudden_Underflow
  1050.                         if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
  1051.                                 dval(rv0) = dval(rv);
  1052.                                 dword0(rv) += P*Exp_msk1;
  1053.                                 adj = dval(aadj1) * ulp(dval(rv));
  1054.                                 dval(rv) += adj;
  1055. #ifdef IBM
  1056.                                 if ((dword0(rv) & Exp_mask) <  P*Exp_msk1)
  1057. #else
  1058.                                 if ((dword0(rv) & Exp_mask) <= P*Exp_msk1)
  1059. #endif
  1060.                                         {
  1061.                                         if (dword0(rv0) == Tiny0
  1062.                                          && dword1(rv0) == Tiny1)
  1063.                                                 goto undfl;
  1064. #ifndef _DOUBLE_IS_32BITS
  1065.                                         dword0(rv) = Tiny0;
  1066.                                         dword1(rv) = Tiny1;
  1067. #else
  1068.                                         dword0(rv) = Tiny1;
  1069. #endif /*_DOUBLE_IS_32BITS*/
  1070.                                         goto cont;
  1071.                                         }
  1072.                                 else
  1073.                                         dword0(rv) -= P*Exp_msk1;
  1074.                                 }
  1075.                         else {
  1076.                                 adj = dval(aadj1) * ulp(dval(rv));
  1077.                                 dval(rv) += adj;
  1078.                                 }
  1079. #else /*Sudden_Underflow*/
  1080.                         /* Compute adj so that the IEEE rounding rules will
  1081.                          * correctly round rv + adj in some half-way cases.
  1082.                          * If rv * ulp(rv) is denormalized (i.e.,
  1083.                          * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
  1084.                          * trouble from bits lost to denormalization;
  1085.                          * example: 1.2e-307 .
  1086.                          */
  1087.                         if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
  1088.                                 dval(aadj1) = (double)(int)(aadj + 0.5);
  1089.                                 if (!dsign)
  1090.                                         dval(aadj1) = -dval(aadj1);
  1091.                                 }
  1092.                         adj = dval(aadj1) * ulp(dval(rv));
  1093.                         dval(rv) += adj;
  1094. #endif /*Sudden_Underflow*/
  1095. #endif /*Avoid_Underflow*/
  1096.                         }
  1097.                 z = dword0(rv) & Exp_mask;
  1098. #ifndef SET_INEXACT
  1099. #ifdef Avoid_Underflow
  1100.                 if (!scale)
  1101. #endif
  1102.                 if (y == z) {
  1103.                         /* Can we stop now? */
  1104.                         L = (Long)aadj;
  1105.                         aadj -= L;
  1106.                         /* The tolerances below are conservative. */
  1107.                         if (dsign || dword1(rv) || dword0(rv) & Bndry_mask) {
  1108.                                 if (aadj < .4999999 || aadj > .5000001)
  1109.                                         break;
  1110.                                 }
  1111.                         else if (aadj < .4999999/FLT_RADIX)
  1112.                                 break;
  1113.                         }
  1114. #endif
  1115.  cont:
  1116.                 Bfree(ptr,bb);
  1117.                 Bfree(ptr,bd);
  1118.                 Bfree(ptr,bs);
  1119.                 Bfree(ptr,delta);
  1120.                 }
  1121. #ifdef SET_INEXACT
  1122.         if (inexact) {
  1123.                 if (!oldinexact) {
  1124.                         dword0(rv0) = Exp_1 + (70 << Exp_shift);
  1125. #ifndef _DOUBLE_IS_32BITS
  1126.                         dword1(rv0) = 0;
  1127. #endif /*!_DOUBLE_IS_32BITS*/
  1128.                         dval(rv0) += 1.;
  1129.                         }
  1130.                 }
  1131.         else if (!oldinexact)
  1132.                 clear_inexact();
  1133. #endif
  1134. #ifdef Avoid_Underflow
  1135.         if (scale) {
  1136.                 dword0(rv0) = Exp_1 - 2*P*Exp_msk1;
  1137. #ifndef _DOUBLE_IS_32BITS
  1138.                 dword1(rv0) = 0;
  1139. #endif /*!_DOUBLE_IS_32BITS*/
  1140.                 dval(rv) *= dval(rv0);
  1141. #ifndef NO_ERRNO
  1142.                 /* try to avoid the bug of testing an 8087 register value */
  1143.                 if (dword0(rv) == 0 && dword1(rv) == 0)
  1144.                         ptr->_errno = ERANGE;
  1145. #endif
  1146.                 }
  1147. #endif /* Avoid_Underflow */
  1148. #ifdef SET_INEXACT
  1149.         if (inexact && !(dword0(rv) & Exp_mask)) {
  1150.                 /* set underflow bit */
  1151.                 dval(rv0) = 1e-300;
  1152.                 dval(rv0) *= dval(rv0);
  1153.                 }
  1154. #endif
  1155.  retfree:
  1156.         Bfree(ptr,bb);
  1157.         Bfree(ptr,bd);
  1158.         Bfree(ptr,bs);
  1159.         Bfree(ptr,bd0);
  1160.         Bfree(ptr,delta);
  1161.  ret:
  1162.         if (se)
  1163.                 *se = (char *)s;
  1164.         return sign ? -dval(rv) : dval(rv);
  1165. }
  1166.  
  1167. #ifndef _REENT_ONLY
  1168.  
  1169. double
  1170. _DEFUN (strtod, (s00, se),
  1171.         _CONST char *s00 _AND char **se)
  1172. {
  1173.   return _strtod_r (_REENT, s00, se);
  1174. }
  1175.  
  1176. float
  1177. _DEFUN (strtof, (s00, se),
  1178.         _CONST char *s00 _AND
  1179.         char **se)
  1180. {
  1181.   double retval = _strtod_r (_REENT, s00, se);
  1182.   if (isnan (retval))
  1183.     return nanf (NULL);
  1184.   return (float)retval;
  1185. }
  1186.  
  1187. #endif
  1188.