Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | 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 < DBL_DIG + 1) {
  304.                         if (nd < 9)
  305.                                 y = 10*y + c - '0';
  306.                         else
  307.                                 z = 10*z + c - '0';
  308.                 }
  309.         }
  310.         nd0 = nd;
  311.         if (strncmp (s, _localeconv_r (ptr)->decimal_point,
  312.                      strlen (_localeconv_r (ptr)->decimal_point)) == 0) {
  313.                 decpt = 1;
  314.                 c = *(s += strlen (_localeconv_r (ptr)->decimal_point));
  315.                 if (!nd) {
  316.                         for(; c == '0'; c = *++s)
  317.                                 nz++;
  318.                         if (c > '0' && c <= '9') {
  319.                                 s0 = s;
  320.                                 nf += nz;
  321.                                 nz = 0;
  322.                                 goto have_dig;
  323.                                 }
  324.                         goto dig_done;
  325.                         }
  326.                 for(; c >= '0' && c <= '9'; c = *++s) {
  327.  have_dig:
  328.                         nz++;
  329.                         if (c -= '0') {
  330.                                 for(i = 1; i < nz; i++) {
  331.                                         if (nd <= DBL_DIG + 1) {
  332.                                                 if (nd + i < 10)
  333.                                                         y *= 10;
  334.                                                 else
  335.                                                         z *= 10;
  336.                                         }
  337.                                 }
  338.                                 if (nd <= DBL_DIG + 1) {
  339.                                         if (nd + i < 10)
  340.                                                 y = 10*y + c;
  341.                                         else
  342.                                                 z = 10*z + c;
  343.                                 }
  344.                                 if (nd <= DBL_DIG + 1) {
  345.                                         nf += nz;
  346.                                         nd += nz;
  347.                                 }
  348.                                 nz = 0;
  349.                                 }
  350.                         }
  351.                 }
  352.  dig_done:
  353.         e = 0;
  354.         if (c == 'e' || c == 'E') {
  355.                 if (!nd && !nz && !nz0) {
  356.                         goto ret0;
  357.                         }
  358.                 s00 = s;
  359.                 esign = 0;
  360.                 switch(c = *++s) {
  361.                         case '-':
  362.                                 esign = 1;
  363.                         case '+':
  364.                                 c = *++s;
  365.                         }
  366.                 if (c >= '0' && c <= '9') {
  367.                         while(c == '0')
  368.                                 c = *++s;
  369.                         if (c > '0' && c <= '9') {
  370.                                 L = c - '0';
  371.                                 s1 = s;
  372.                                 while((c = *++s) >= '0' && c <= '9')
  373.                                         L = 10*L + c - '0';
  374.                                 if (s - s1 > 8 || L > 19999)
  375.                                         /* Avoid confusion from exponents
  376.                                          * so large that e might overflow.
  377.                                          */
  378.                                         e = 19999; /* safe for 16 bit ints */
  379.                                 else
  380.                                         e = (int)L;
  381.                                 if (esign)
  382.                                         e = -e;
  383.                                 }
  384.                         else
  385.                                 e = 0;
  386.                         }
  387.                 else
  388.                         s = s00;
  389.                 }
  390.         if (!nd) {
  391.                 if (!nz && !nz0) {
  392. #ifdef INFNAN_CHECK
  393.                         /* Check for Nan and Infinity */
  394.                         __ULong bits[2];
  395.                         static FPI fpinan =     /* only 52 explicit bits */
  396.                                 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
  397.                         if (!decpt)
  398.                          switch(c) {
  399.                           case 'i':
  400.                           case 'I':
  401.                                 if (match(&s,"nf")) {
  402.                                         --s;
  403.                                         if (!match(&s,"inity"))
  404.                                                 ++s;
  405.                                         dword0(rv) = 0x7ff00000;
  406. #ifndef _DOUBLE_IS_32BITS
  407.                                         dword1(rv) = 0;
  408. #endif /*!_DOUBLE_IS_32BITS*/
  409.                                         goto ret;
  410.                                         }
  411.                                 break;
  412.                           case 'n':
  413.                           case 'N':
  414.                                 if (match(&s, "an")) {
  415. #ifndef No_Hex_NaN
  416.                                         if (*s == '(' /*)*/
  417.                                          && hexnan(&s, &fpinan, bits)
  418.                                                         == STRTOG_NaNbits) {
  419.                                                 dword0(rv) = 0x7ff00000 | bits[1];
  420. #ifndef _DOUBLE_IS_32BITS
  421.                                                 dword1(rv) = bits[0];
  422. #endif /*!_DOUBLE_IS_32BITS*/
  423.                                                 }
  424.                                         else {
  425. #endif
  426.                                                 dword0(rv) = NAN_WORD0;
  427. #ifndef _DOUBLE_IS_32BITS
  428.                                                 dword1(rv) = NAN_WORD1;
  429. #endif /*!_DOUBLE_IS_32BITS*/
  430. #ifndef No_Hex_NaN
  431.                                                 }
  432. #endif
  433.                                         goto ret;
  434.                                         }
  435.                           }
  436. #endif /* INFNAN_CHECK */
  437.  ret0:
  438.                         s = s00;
  439.                         sign = 0;
  440.                         }
  441.                 goto ret;
  442.                 }
  443.         e1 = e -= nf;
  444.  
  445.         /* Now we have nd0 digits, starting at s0, followed by a
  446.          * decimal point, followed by nd-nd0 digits.  The number we're
  447.          * after is the integer represented by those digits times
  448.          * 10**e */
  449.  
  450.         if (!nd0)
  451.                 nd0 = nd;
  452.         k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
  453.         dval(rv) = y;
  454.         if (k > 9) {
  455. #ifdef SET_INEXACT
  456.                 if (k > DBL_DIG)
  457.                         oldinexact = get_inexact();
  458. #endif
  459.                 dval(rv) = tens[k - 9] * dval(rv) + z;
  460.                 }
  461.         bd0 = 0;
  462.         if (nd <= DBL_DIG
  463. #ifndef RND_PRODQUOT
  464. #ifndef Honor_FLT_ROUNDS
  465.                 && Flt_Rounds == 1
  466. #endif
  467. #endif
  468.                         ) {
  469.                 if (!e)
  470.                         goto ret;
  471.                 if (e > 0) {
  472.                         if (e <= Ten_pmax) {
  473. #ifdef VAX
  474.                                 goto vax_ovfl_check;
  475. #else
  476. #ifdef Honor_FLT_ROUNDS
  477.                                 /* round correctly FLT_ROUNDS = 2 or 3 */
  478.                                 if (sign) {
  479.                                         dval(rv) = -dval(rv);
  480.                                         sign = 0;
  481.                                         }
  482. #endif
  483.                                 /* rv = */ rounded_product(dval(rv), tens[e]);
  484.                                 goto ret;
  485. #endif
  486.                                 }
  487.                         i = DBL_DIG - nd;
  488.                         if (e <= Ten_pmax + i) {
  489.                                 /* A fancier test would sometimes let us do
  490.                                  * this for larger i values.
  491.                                  */
  492. #ifdef Honor_FLT_ROUNDS
  493.                                 /* round correctly FLT_ROUNDS = 2 or 3 */
  494.                                 if (sign) {
  495.                                         dval(rv) = -dval(rv);
  496.                                         sign = 0;
  497.                                         }
  498. #endif
  499.                                 e -= i;
  500.                                 dval(rv) *= tens[i];
  501. #ifdef VAX
  502.                                 /* VAX exponent range is so narrow we must
  503.                                  * worry about overflow here...
  504.                                  */
  505.  vax_ovfl_check:
  506.                                 dword0(rv) -= P*Exp_msk1;
  507.                                 /* rv = */ rounded_product(dval(rv), tens[e]);
  508.                                 if ((dword0(rv) & Exp_mask)
  509.                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
  510.                                         goto ovfl;
  511.                                 dword0(rv) += P*Exp_msk1;
  512. #else
  513.                                 /* rv = */ rounded_product(dval(rv), tens[e]);
  514. #endif
  515.                                 goto ret;
  516.                                 }
  517.                         }
  518. #ifndef Inaccurate_Divide
  519.                 else if (e >= -Ten_pmax) {
  520. #ifdef Honor_FLT_ROUNDS
  521.                         /* round correctly FLT_ROUNDS = 2 or 3 */
  522.                         if (sign) {
  523.                                 dval(rv) = -dval(rv);
  524.                                 sign = 0;
  525.                                 }
  526. #endif
  527.                         /* rv = */ rounded_quotient(dval(rv), tens[-e]);
  528.                         goto ret;
  529.                         }
  530. #endif
  531.                 }
  532.         e1 += nd - k;
  533.  
  534. #ifdef IEEE_Arith
  535. #ifdef SET_INEXACT
  536.         inexact = 1;
  537.         if (k <= DBL_DIG)
  538.                 oldinexact = get_inexact();
  539. #endif
  540. #ifdef Avoid_Underflow
  541.         scale = 0;
  542. #endif
  543. #ifdef Honor_FLT_ROUNDS
  544.         if ((rounding = Flt_Rounds) >= 2) {
  545.                 if (sign)
  546.                         rounding = rounding == 2 ? 0 : 2;
  547.                 else
  548.                         if (rounding != 2)
  549.                                 rounding = 0;
  550.                 }
  551. #endif
  552. #endif /*IEEE_Arith*/
  553.  
  554.         /* Get starting approximation = rv * 10**e1 */
  555.  
  556.         if (e1 > 0) {
  557.                 if ( (i = e1 & 15) !=0)
  558.                         dval(rv) *= tens[i];
  559.                 if (e1 &= ~15) {
  560.                         if (e1 > DBL_MAX_10_EXP) {
  561.  ovfl:
  562. #ifndef NO_ERRNO
  563.                                 ptr->_errno = ERANGE;
  564. #endif
  565.                                 /* Can't trust HUGE_VAL */
  566. #ifdef IEEE_Arith
  567. #ifdef Honor_FLT_ROUNDS
  568.                                 switch(rounding) {
  569.                                   case 0: /* toward 0 */
  570.                                   case 3: /* toward -infinity */
  571.                                         dword0(rv) = Big0;
  572. #ifndef _DOUBLE_IS_32BITS
  573.                                         dword1(rv) = Big1;
  574. #endif /*!_DOUBLE_IS_32BITS*/
  575.                                         break;
  576.                                   default:
  577.                                         dword0(rv) = Exp_mask;
  578. #ifndef _DOUBLE_IS_32BITS
  579.                                         dword1(rv) = 0;
  580. #endif /*!_DOUBLE_IS_32BITS*/
  581.                                   }
  582. #else /*Honor_FLT_ROUNDS*/
  583.                                 dword0(rv) = Exp_mask;
  584. #ifndef _DOUBLE_IS_32BITS
  585.                                 dword1(rv) = 0;
  586. #endif /*!_DOUBLE_IS_32BITS*/
  587. #endif /*Honor_FLT_ROUNDS*/
  588. #ifdef SET_INEXACT
  589.                                 /* set overflow bit */
  590.                                 dval(rv0) = 1e300;
  591.                                 dval(rv0) *= dval(rv0);
  592. #endif
  593. #else /*IEEE_Arith*/
  594.                                 dword0(rv) = Big0;
  595. #ifndef _DOUBLE_IS_32BITS
  596.                                 dword1(rv) = Big1;
  597. #endif /*!_DOUBLE_IS_32BITS*/
  598. #endif /*IEEE_Arith*/
  599.                                 if (bd0)
  600.                                         goto retfree;
  601.                                 goto ret;
  602.                                 }
  603.                         e1 >>= 4;
  604.                         for(j = 0; e1 > 1; j++, e1 >>= 1)
  605.                                 if (e1 & 1)
  606.                                         dval(rv) *= bigtens[j];
  607.                 /* The last multiplication could overflow. */
  608.                         dword0(rv) -= P*Exp_msk1;
  609.                         dval(rv) *= bigtens[j];
  610.                         if ((z = dword0(rv) & Exp_mask)
  611.                          > Exp_msk1*(DBL_MAX_EXP+Bias-P))
  612.                                 goto ovfl;
  613.                         if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
  614.                                 /* set to largest number */
  615.                                 /* (Can't trust DBL_MAX) */
  616.                                 dword0(rv) = Big0;
  617. #ifndef _DOUBLE_IS_32BITS
  618.                                 dword1(rv) = Big1;
  619. #endif /*!_DOUBLE_IS_32BITS*/
  620.                                 }
  621.                         else
  622.                                 dword0(rv) += P*Exp_msk1;
  623.                         }
  624.                 }
  625.         else if (e1 < 0) {
  626.                 e1 = -e1;
  627.                 if ( (i = e1 & 15) !=0)
  628.                         dval(rv) /= tens[i];
  629.                 if (e1 >>= 4) {
  630.                         if (e1 >= 1 << n_bigtens)
  631.                                 goto undfl;
  632. #ifdef Avoid_Underflow
  633.                         if (e1 & Scale_Bit)
  634.                                 scale = 2*P;
  635.                         for(j = 0; e1 > 0; j++, e1 >>= 1)
  636.                                 if (e1 & 1)
  637.                                         dval(rv) *= tinytens[j];
  638.                         if (scale && (j = 2*P + 1 - ((dword0(rv) & Exp_mask)
  639.                                                 >> Exp_shift)) > 0) {
  640.                                 /* scaled rv is denormal; zap j low bits */
  641.                                 if (j >= 32) {
  642. #ifndef _DOUBLE_IS_32BITS
  643.                                         dword1(rv) = 0;
  644. #endif /*!_DOUBLE_IS_32BITS*/
  645.                                         if (j >= 53)
  646.                                          dword0(rv) = (P+2)*Exp_msk1;
  647.                                         else
  648.                                          dword0(rv) &= 0xffffffff << (j-32);
  649.                                         }
  650. #ifndef _DOUBLE_IS_32BITS
  651.                                 else
  652.                                         dword1(rv) &= 0xffffffff << j;
  653. #endif /*!_DOUBLE_IS_32BITS*/
  654.                                 }
  655. #else
  656.                         for(j = 0; e1 > 1; j++, e1 >>= 1)
  657.                                 if (e1 & 1)
  658.                                         dval(rv) *= tinytens[j];
  659.                         /* The last multiplication could underflow. */
  660.                         dval(rv0) = dval(rv);
  661.                         dval(rv) *= tinytens[j];
  662.                         if (!dval(rv)) {
  663.                                 dval(rv) = 2.*dval(rv0);
  664.                                 dval(rv) *= tinytens[j];
  665. #endif
  666.                                 if (!dval(rv)) {
  667.  undfl:
  668.                                         dval(rv) = 0.;
  669. #ifndef NO_ERRNO
  670.                                         ptr->_errno = ERANGE;
  671. #endif
  672.                                         if (bd0)
  673.                                                 goto retfree;
  674.                                         goto ret;
  675.                                         }
  676. #ifndef Avoid_Underflow
  677. #ifndef _DOUBLE_IS_32BITS
  678.                                 dword0(rv) = Tiny0;
  679.                                 dword1(rv) = Tiny1;
  680. #else
  681.                                 dword0(rv) = Tiny1;
  682. #endif /*_DOUBLE_IS_32BITS*/
  683.                                 /* The refinement below will clean
  684.                                  * this approximation up.
  685.                                  */
  686.                                 }
  687. #endif
  688.                         }
  689.                 }
  690.  
  691.         /* Now the hard part -- adjusting rv to the correct value.*/
  692.  
  693.         /* Put digits into bd: true value = bd * 10^e */
  694.  
  695.         bd0 = s2b(ptr, s0, nd0, nd, y);
  696.  
  697.         for(;;) {
  698.                 bd = Balloc(ptr,bd0->_k);
  699.                 Bcopy(bd, bd0);
  700.                 bb = d2b(ptr,dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
  701.                 bs = i2b(ptr,1);
  702.  
  703.                 if (e >= 0) {
  704.                         bb2 = bb5 = 0;
  705.                         bd2 = bd5 = e;
  706.                         }
  707.                 else {
  708.                         bb2 = bb5 = -e;
  709.                         bd2 = bd5 = 0;
  710.                         }
  711.                 if (bbe >= 0)
  712.                         bb2 += bbe;
  713.                 else
  714.                         bd2 -= bbe;
  715.                 bs2 = bb2;
  716. #ifdef Honor_FLT_ROUNDS
  717.                 if (rounding != 1)
  718.                         bs2++;
  719. #endif
  720. #ifdef Avoid_Underflow
  721.                 j = bbe - scale;
  722.                 i = j + bbbits - 1;     /* logb(rv) */
  723.                 if (i < Emin)   /* denormal */
  724.                         j += P - Emin;
  725.                 else
  726.                         j = P + 1 - bbbits;
  727. #else /*Avoid_Underflow*/
  728. #ifdef Sudden_Underflow
  729. #ifdef IBM
  730.                 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
  731. #else
  732.                 j = P + 1 - bbbits;
  733. #endif
  734. #else /*Sudden_Underflow*/
  735.                 j = bbe;
  736.                 i = j + bbbits - 1;     /* logb(rv) */
  737.                 if (i < Emin)   /* denormal */
  738.                         j += P - Emin;
  739.                 else
  740.                         j = P + 1 - bbbits;
  741. #endif /*Sudden_Underflow*/
  742. #endif /*Avoid_Underflow*/
  743.                 bb2 += j;
  744.                 bd2 += j;
  745. #ifdef Avoid_Underflow
  746.                 bd2 += scale;
  747. #endif
  748.                 i = bb2 < bd2 ? bb2 : bd2;
  749.                 if (i > bs2)
  750.                         i = bs2;
  751.                 if (i > 0) {
  752.                         bb2 -= i;
  753.                         bd2 -= i;
  754.                         bs2 -= i;
  755.                         }
  756.                 if (bb5 > 0) {
  757.                         bs = pow5mult(ptr, bs, bb5);
  758.                         bb1 = mult(ptr, bs, bb);
  759.                         Bfree(ptr, bb);
  760.                         bb = bb1;
  761.                         }
  762.                 if (bb2 > 0)
  763.                         bb = lshift(ptr, bb, bb2);
  764.                 if (bd5 > 0)
  765.                         bd = pow5mult(ptr, bd, bd5);
  766.                 if (bd2 > 0)
  767.                         bd = lshift(ptr, bd, bd2);
  768.                 if (bs2 > 0)
  769.                         bs = lshift(ptr, bs, bs2);
  770.                 delta = diff(ptr, bb, bd);
  771.                 dsign = delta->_sign;
  772.                 delta->_sign = 0;
  773.                 i = cmp(delta, bs);
  774. #ifdef Honor_FLT_ROUNDS
  775.                 if (rounding != 1) {
  776.                         if (i < 0) {
  777.                                 /* Error is less than an ulp */
  778.                                 if (!delta->_x[0] && delta->_wds <= 1) {
  779.                                         /* exact */
  780. #ifdef SET_INEXACT
  781.                                         inexact = 0;
  782. #endif
  783.                                         break;
  784.                                         }
  785.                                 if (rounding) {
  786.                                         if (dsign) {
  787.                                                 adj = 1.;
  788.                                                 goto apply_adj;
  789.                                                 }
  790.                                         }
  791.                                 else if (!dsign) {
  792.                                         adj = -1.;
  793.                                         if (!dword1(rv)
  794.                                          && !(dword0(rv) & Frac_mask)) {
  795.                                                 y = dword0(rv) & Exp_mask;
  796. #ifdef Avoid_Underflow
  797.                                                 if (!scale || y > 2*P*Exp_msk1)
  798. #else
  799.                                                 if (y)
  800. #endif
  801.                                                   {
  802.                                                   delta = lshift(ptr, delta,Log2P);
  803.                                                   if (cmp(delta, bs) <= 0)
  804.                                                         adj = -0.5;
  805.                                                   }
  806.                                                 }
  807.  apply_adj:
  808. #ifdef Avoid_Underflow
  809.                                         if (scale && (y = dword0(rv) & Exp_mask)
  810.                                                 <= 2*P*Exp_msk1)
  811.                                           dword0(adj) += (2*P+1)*Exp_msk1 - y;
  812. #else
  813. #ifdef Sudden_Underflow
  814.                                         if ((dword0(rv) & Exp_mask) <=
  815.                                                         P*Exp_msk1) {
  816.                                                 dword0(rv) += P*Exp_msk1;
  817.                                                 dval(rv) += adj*ulp(dval(rv));
  818.                                                 dword0(rv) -= P*Exp_msk1;
  819.                                                 }
  820.                                         else
  821. #endif /*Sudden_Underflow*/
  822. #endif /*Avoid_Underflow*/
  823.                                         dval(rv) += adj*ulp(dval(rv));
  824.                                         }
  825.                                 break;
  826.                                 }
  827.                         adj = ratio(delta, bs);
  828.                         if (adj < 1.)
  829.                                 adj = 1.;
  830.                         if (adj <= 0x7ffffffe) {
  831.                                 /* adj = rounding ? ceil(adj) : floor(adj); */
  832.                                 y = adj;
  833.                                 if (y != adj) {
  834.                                         if (!((rounding>>1) ^ dsign))
  835.                                                 y++;
  836.                                         adj = y;
  837.                                         }
  838.                                 }
  839. #ifdef Avoid_Underflow
  840.                         if (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1)
  841.                                 dword0(adj) += (2*P+1)*Exp_msk1 - y;
  842. #else
  843. #ifdef Sudden_Underflow
  844.                         if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
  845.                                 dword0(rv) += P*Exp_msk1;
  846.                                 adj *= ulp(dval(rv));
  847.                                 if (dsign)
  848.                                         dval(rv) += adj;
  849.                                 else
  850.                                         dval(rv) -= adj;
  851.                                 dword0(rv) -= P*Exp_msk1;
  852.                                 goto cont;
  853.                                 }
  854. #endif /*Sudden_Underflow*/
  855. #endif /*Avoid_Underflow*/
  856.                         adj *= ulp(dval(rv));
  857.                         if (dsign)
  858.                                 dval(rv) += adj;
  859.                         else
  860.                                 dval(rv) -= adj;
  861.                         goto cont;
  862.                         }
  863. #endif /*Honor_FLT_ROUNDS*/
  864.  
  865.                 if (i < 0) {
  866.                         /* Error is less than half an ulp -- check for
  867.                          * special case of mantissa a power of two.
  868.                          */
  869.                         if (dsign || dword1(rv) || dword0(rv) & Bndry_mask
  870. #ifdef IEEE_Arith
  871. #ifdef Avoid_Underflow
  872.                          || (dword0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
  873. #else
  874.                          || (dword0(rv) & Exp_mask) <= Exp_msk1
  875. #endif
  876. #endif
  877.                                 ) {
  878. #ifdef SET_INEXACT
  879.                                 if (!delta->x[0] && delta->wds <= 1)
  880.                                         inexact = 0;
  881. #endif
  882.                                 break;
  883.                                 }
  884.                         if (!delta->_x[0] && delta->_wds <= 1) {
  885.                                 /* exact result */
  886. #ifdef SET_INEXACT
  887.                                 inexact = 0;
  888. #endif
  889.                                 break;
  890.                                 }
  891.                         delta = lshift(ptr,delta,Log2P);
  892.                         if (cmp(delta, bs) > 0)
  893.                                 goto drop_down;
  894.                         break;
  895.                         }
  896.                 if (i == 0) {
  897.                         /* exactly half-way between */
  898.                         if (dsign) {
  899.                                 if ((dword0(rv) & Bndry_mask1) == Bndry_mask1
  900.                                  &&  dword1(rv) == (
  901. #ifdef Avoid_Underflow
  902.                         (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1)
  903.                 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
  904. #endif
  905.                                                    0xffffffff)) {
  906.                                         /*boundary case -- increment exponent*/
  907.                                         dword0(rv) = (dword0(rv) & Exp_mask)
  908.                                                 + Exp_msk1
  909. #ifdef IBM
  910.                                                 | Exp_msk1 >> 4
  911. #endif
  912.                                                 ;
  913. #ifndef _DOUBLE_IS_32BITS
  914.                                         dword1(rv) = 0;
  915. #endif /*!_DOUBLE_IS_32BITS*/
  916. #ifdef Avoid_Underflow
  917.                                         dsign = 0;
  918. #endif
  919.                                         break;
  920.                                         }
  921.                                 }
  922.                         else if (!(dword0(rv) & Bndry_mask) && !dword1(rv)) {
  923.  drop_down:
  924.                                 /* boundary case -- decrement exponent */
  925. #ifdef Sudden_Underflow /*{{*/
  926.                                 L = dword0(rv) & Exp_mask;
  927. #ifdef IBM
  928.                                 if (L <  Exp_msk1)
  929. #else
  930. #ifdef Avoid_Underflow
  931.                                 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
  932. #else
  933.                                 if (L <= Exp_msk1)
  934. #endif /*Avoid_Underflow*/
  935. #endif /*IBM*/
  936.                                         goto undfl;
  937.                                 L -= Exp_msk1;
  938. #else /*Sudden_Underflow}{*/
  939. #ifdef Avoid_Underflow
  940.                                 if (scale) {
  941.                                         L = dword0(rv) & Exp_mask;
  942.                                         if (L <= (2*P+1)*Exp_msk1) {
  943.                                                 if (L > (P+2)*Exp_msk1)
  944.                                                         /* round even ==> */
  945.                                                         /* accept rv */
  946.                                                         break;
  947.                                                 /* rv = smallest denormal */
  948.                                                 goto undfl;
  949.                                                 }
  950.                                         }
  951. #endif /*Avoid_Underflow*/
  952.                                 L = (dword0(rv) & Exp_mask) - Exp_msk1;
  953. #endif /*Sudden_Underflow}*/
  954.                                 dword0(rv) = L | Bndry_mask1;
  955. #ifndef _DOUBLE_IS_32BITS
  956.                                 dword1(rv) = 0xffffffff;
  957. #endif /*!_DOUBLE_IS_32BITS*/
  958. #ifdef IBM
  959.                                 goto cont;
  960. #else
  961.                                 break;
  962. #endif
  963.                                 }
  964. #ifndef ROUND_BIASED
  965.                         if (!(dword1(rv) & LSB))
  966.                                 break;
  967. #endif
  968.                         if (dsign)
  969.                                 dval(rv) += ulp(dval(rv));
  970. #ifndef ROUND_BIASED
  971.                         else {
  972.                                 dval(rv) -= ulp(dval(rv));
  973. #ifndef Sudden_Underflow
  974.                                 if (!dval(rv))
  975.                                         goto undfl;
  976. #endif
  977.                                 }
  978. #ifdef Avoid_Underflow
  979.                         dsign = 1 - dsign;
  980. #endif
  981. #endif
  982.                         break;
  983.                         }
  984.                 if ((aadj = ratio(delta, bs)) <= 2.) {
  985.                         if (dsign)
  986.                                 aadj = dval(aadj1) = 1.;
  987.                         else if (dword1(rv) || dword0(rv) & Bndry_mask) {
  988. #ifndef Sudden_Underflow
  989.                                 if (dword1(rv) == Tiny1 && !dword0(rv))
  990.                                         goto undfl;
  991. #endif
  992.                                 aadj = 1.;
  993.                                 dval(aadj1) = -1.;
  994.                                 }
  995.                         else {
  996.                                 /* special case -- power of FLT_RADIX to be */
  997.                                 /* rounded down... */
  998.  
  999.                                 if (aadj < 2./FLT_RADIX)
  1000.                                         aadj = 1./FLT_RADIX;
  1001.                                 else
  1002.                                         aadj *= 0.5;
  1003.                                 dval(aadj1) = -aadj;
  1004.                                 }
  1005.                         }
  1006.                 else {
  1007.                         aadj *= 0.5;
  1008.                         dval(aadj1) = dsign ? aadj : -aadj;
  1009. #ifdef Check_FLT_ROUNDS
  1010.                         switch(Rounding) {
  1011.                                 case 2: /* towards +infinity */
  1012.                                         dval(aadj1) -= 0.5;
  1013.                                         break;
  1014.                                 case 0: /* towards 0 */
  1015.                                 case 3: /* towards -infinity */
  1016.                                         dval(aadj1) += 0.5;
  1017.                                 }
  1018. #else
  1019.                         if (Flt_Rounds == 0)
  1020.                                 dval(aadj1) += 0.5;
  1021. #endif /*Check_FLT_ROUNDS*/
  1022.                         }
  1023.                 y = dword0(rv) & Exp_mask;
  1024.  
  1025.                 /* Check for overflow */
  1026.  
  1027.                 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
  1028.                         dval(rv0) = dval(rv);
  1029.                         dword0(rv) -= P*Exp_msk1;
  1030.                         adj = dval(aadj1) * ulp(dval(rv));
  1031.                         dval(rv) += adj;
  1032.                         if ((dword0(rv) & Exp_mask) >=
  1033.                                         Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
  1034.                                 if (dword0(rv0) == Big0 && dword1(rv0) == Big1)
  1035.                                         goto ovfl;
  1036.                                 dword0(rv) = Big0;
  1037. #ifndef _DOUBLE_IS_32BITS
  1038.                                 dword1(rv) = Big1;
  1039. #endif /*!_DOUBLE_IS_32BITS*/
  1040.                                 goto cont;
  1041.                                 }
  1042.                         else
  1043.                                 dword0(rv) += P*Exp_msk1;
  1044.                         }
  1045.                 else {
  1046. #ifdef Avoid_Underflow
  1047.                         if (scale && y <= 2*P*Exp_msk1) {
  1048.                                 if (aadj <= 0x7fffffff) {
  1049.                                         if ((z = aadj) <= 0)
  1050.                                                 z = 1;
  1051.                                         aadj = z;
  1052.                                         dval(aadj1) = dsign ? aadj : -aadj;
  1053.                                         }
  1054.                                 dword0(aadj1) += (2*P+1)*Exp_msk1 - y;
  1055.                                 }
  1056.                         adj = dval(aadj1) * ulp(dval(rv));
  1057.                         dval(rv) += adj;
  1058. #else
  1059. #ifdef Sudden_Underflow
  1060.                         if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
  1061.                                 dval(rv0) = dval(rv);
  1062.                                 dword0(rv) += P*Exp_msk1;
  1063.                                 adj = dval(aadj1) * ulp(dval(rv));
  1064.                                 dval(rv) += adj;
  1065. #ifdef IBM
  1066.                                 if ((dword0(rv) & Exp_mask) <  P*Exp_msk1)
  1067. #else
  1068.                                 if ((dword0(rv) & Exp_mask) <= P*Exp_msk1)
  1069. #endif
  1070.                                         {
  1071.                                         if (dword0(rv0) == Tiny0
  1072.                                          && dword1(rv0) == Tiny1)
  1073.                                                 goto undfl;
  1074. #ifndef _DOUBLE_IS_32BITS
  1075.                                         dword0(rv) = Tiny0;
  1076.                                         dword1(rv) = Tiny1;
  1077. #else
  1078.                                         dword0(rv) = Tiny1;
  1079. #endif /*_DOUBLE_IS_32BITS*/
  1080.                                         goto cont;
  1081.                                         }
  1082.                                 else
  1083.                                         dword0(rv) -= P*Exp_msk1;
  1084.                                 }
  1085.                         else {
  1086.                                 adj = dval(aadj1) * ulp(dval(rv));
  1087.                                 dval(rv) += adj;
  1088.                                 }
  1089. #else /*Sudden_Underflow*/
  1090.                         /* Compute adj so that the IEEE rounding rules will
  1091.                          * correctly round rv + adj in some half-way cases.
  1092.                          * If rv * ulp(rv) is denormalized (i.e.,
  1093.                          * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
  1094.                          * trouble from bits lost to denormalization;
  1095.                          * example: 1.2e-307 .
  1096.                          */
  1097.                         if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
  1098.                                 dval(aadj1) = (double)(int)(aadj + 0.5);
  1099.                                 if (!dsign)
  1100.                                         dval(aadj1) = -dval(aadj1);
  1101.                                 }
  1102.                         adj = dval(aadj1) * ulp(dval(rv));
  1103.                         dval(rv) += adj;
  1104. #endif /*Sudden_Underflow*/
  1105. #endif /*Avoid_Underflow*/
  1106.                         }
  1107.                 z = dword0(rv) & Exp_mask;
  1108. #ifndef SET_INEXACT
  1109. #ifdef Avoid_Underflow
  1110.                 if (!scale)
  1111. #endif
  1112.                 if (y == z) {
  1113.                         /* Can we stop now? */
  1114.                         L = (Long)aadj;
  1115.                         aadj -= L;
  1116.                         /* The tolerances below are conservative. */
  1117.                         if (dsign || dword1(rv) || dword0(rv) & Bndry_mask) {
  1118.                                 if (aadj < .4999999 || aadj > .5000001)
  1119.                                         break;
  1120.                                 }
  1121.                         else if (aadj < .4999999/FLT_RADIX)
  1122.                                 break;
  1123.                         }
  1124. #endif
  1125.  cont:
  1126.                 Bfree(ptr,bb);
  1127.                 Bfree(ptr,bd);
  1128.                 Bfree(ptr,bs);
  1129.                 Bfree(ptr,delta);
  1130.                 }
  1131. #ifdef SET_INEXACT
  1132.         if (inexact) {
  1133.                 if (!oldinexact) {
  1134.                         dword0(rv0) = Exp_1 + (70 << Exp_shift);
  1135. #ifndef _DOUBLE_IS_32BITS
  1136.                         dword1(rv0) = 0;
  1137. #endif /*!_DOUBLE_IS_32BITS*/
  1138.                         dval(rv0) += 1.;
  1139.                         }
  1140.                 }
  1141.         else if (!oldinexact)
  1142.                 clear_inexact();
  1143. #endif
  1144. #ifdef Avoid_Underflow
  1145.         if (scale) {
  1146.                 dword0(rv0) = Exp_1 - 2*P*Exp_msk1;
  1147. #ifndef _DOUBLE_IS_32BITS
  1148.                 dword1(rv0) = 0;
  1149. #endif /*!_DOUBLE_IS_32BITS*/
  1150.                 dval(rv) *= dval(rv0);
  1151. #ifndef NO_ERRNO
  1152.                 /* try to avoid the bug of testing an 8087 register value */
  1153.                 if (dword0(rv) == 0 && dword1(rv) == 0)
  1154.                         ptr->_errno = ERANGE;
  1155. #endif
  1156.                 }
  1157. #endif /* Avoid_Underflow */
  1158. #ifdef SET_INEXACT
  1159.         if (inexact && !(dword0(rv) & Exp_mask)) {
  1160.                 /* set underflow bit */
  1161.                 dval(rv0) = 1e-300;
  1162.                 dval(rv0) *= dval(rv0);
  1163.                 }
  1164. #endif
  1165.  retfree:
  1166.         Bfree(ptr,bb);
  1167.         Bfree(ptr,bd);
  1168.         Bfree(ptr,bs);
  1169.         Bfree(ptr,bd0);
  1170.         Bfree(ptr,delta);
  1171.  ret:
  1172.         if (se)
  1173.                 *se = (char *)s;
  1174.         return sign ? -dval(rv) : dval(rv);
  1175. }
  1176.  
  1177. #ifndef _REENT_ONLY
  1178.  
  1179. double
  1180. _DEFUN (strtod, (s00, se),
  1181.         _CONST char *s00 _AND char **se)
  1182. {
  1183.   return _strtod_r (_REENT, s00, se);
  1184. }
  1185.  
  1186. float
  1187. _DEFUN (strtof, (s00, se),
  1188.         _CONST char *s00 _AND
  1189.         char **se)
  1190. {
  1191.   double retval = _strtod_r (_REENT, s00, se);
  1192.   if (isnan (retval))
  1193.     return nanf (NULL);
  1194.   return (float)retval;
  1195. }
  1196.  
  1197. #endif
  1198.