Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

  1. /****************************************************************
  2.  
  3. The author of this software is David M. Gay.
  4.  
  5. Copyright (C) 1998-2001 by Lucent Technologies
  6. All Rights Reserved
  7.  
  8. Permission to use, copy, modify, and distribute this software and
  9. its documentation for any purpose and without fee is hereby
  10. granted, provided that the above copyright notice appear in all
  11. copies and that both that the copyright notice and this
  12. permission notice and warranty disclaimer appear in supporting
  13. documentation, and that the name of Lucent or any of its entities
  14. not be used in advertising or publicity pertaining to
  15. distribution of the software without specific, written prior
  16. permission.
  17.  
  18. LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
  19. INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
  20. IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
  21. SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  22. WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
  23. IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
  24. ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
  25. THIS SOFTWARE.
  26.  
  27. ****************************************************************/
  28.  
  29. /* Please send bug reports to David M. Gay (dmg at acm dot org,
  30.  * with " at " changed at "@" and " dot " changed to ".").      */
  31.  
  32. #include <_ansi.h>
  33. #include <errno.h>
  34. #include <stdlib.h>
  35. #include <string.h>
  36. #include "mprec.h"
  37. #include "gdtoa.h"
  38. #include "gd_qnan.h"
  39.  
  40. #include "locale.h"
  41.  
  42. #if defined (_HAVE_LONG_DOUBLE) && !defined (_LDBL_EQ_DBL)
  43.  
  44. #define USE_LOCALE
  45.  
  46.  static const int
  47. fivesbits[] = {  0,  3,  5,  7, 10, 12, 14, 17, 19, 21,
  48.                 24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
  49.                 47, 49, 52
  50. #ifdef VAX
  51.                 , 54, 56
  52. #endif
  53.                 };
  54.  
  55. static _Bigint *
  56. #ifdef KR_headers
  57. sum(p, a, b) struct _reent *p; _Bigint *a; _Bigint *b;
  58. #else
  59. sum(struct _reent *p, _Bigint *a, _Bigint *b)
  60. #endif
  61. {
  62.         _Bigint *c;
  63.         __ULong carry, *xc, *xa, *xb, *xe, y;
  64. #ifdef Pack_32
  65.         __ULong z;
  66. #endif
  67.  
  68.         if (a->_wds < b->_wds) {
  69.                 c = b; b = a; a = c;
  70.                 }
  71.         c = Balloc(p, a->_k);
  72.         c->_wds = a->_wds;
  73.         carry = 0;
  74.         xa = a->_x;
  75.         xb = b->_x;
  76.         xc = c->_x;
  77.         xe = xc + b->_wds;
  78. #ifdef Pack_32
  79.         do {
  80.                 y = (*xa & 0xffff) + (*xb & 0xffff) + carry;
  81.                 carry = (y & 0x10000) >> 16;
  82.                 z = (*xa++ >> 16) + (*xb++ >> 16) + carry;
  83.                 carry = (z & 0x10000) >> 16;
  84.                 Storeinc(xc, z, y);
  85.                 }
  86.                 while(xc < xe);
  87.         xe += a->_wds - b->_wds;
  88.         while(xc < xe) {
  89.                 y = (*xa & 0xffff) + carry;
  90.                 carry = (y & 0x10000) >> 16;
  91.                 z = (*xa++ >> 16) + carry;
  92.                 carry = (z & 0x10000) >> 16;
  93.                 Storeinc(xc, z, y);
  94.                 }
  95. #else
  96.         do {
  97.                 y = *xa++ + *xb++ + carry;
  98.                 carry = (y & 0x10000) >> 16;
  99.                 *xc++ = y & 0xffff;
  100.                 }
  101.                 while(xc < xe);
  102.         xe += a->_wds - b->_wds;
  103.         while(xc < xe) {
  104.                 y = *xa++ + carry;
  105.                 carry = (y & 0x10000) >> 16;
  106.                 *xc++ = y & 0xffff;
  107.                 }
  108. #endif
  109.         if (carry) {
  110.                 if (c->_wds == c->_maxwds) {
  111.                         b = Balloc(p, c->_k + 1);
  112.                         Bcopy(b, c);
  113.                         Bfree(p, c);
  114.                         c = b;
  115.                         }
  116.                 c->_x[c->_wds++] = 1;
  117.                 }
  118.         return c;
  119.         }
  120.  
  121. static void
  122. #ifdef KR_headers
  123. rshift(b, k) _Bigint *b; int k;
  124. #else
  125. rshift(_Bigint *b, int k)
  126. #endif
  127. {
  128.         __ULong *x, *x1, *xe, y;
  129.         int n;
  130.  
  131.         x = x1 = b->_x;
  132.         n = k >> kshift;
  133.         if (n < b->_wds) {
  134.                 xe = x + b->_wds;
  135.                 x += n;
  136.                 if (k &= kmask) {
  137.                         n = ULbits - k;
  138.                         y = *x++ >> k;
  139.                         while(x < xe) {
  140.                                 *x1++ = (y | (*x << n)) & ALL_ON;
  141.                                 y = *x++ >> k;
  142.                                 }
  143.                         if ((*x1 = y) !=0)
  144.                                 x1++;
  145.                         }
  146.                 else
  147.                         while(x < xe)
  148.                                 *x1++ = *x++;
  149.                 }
  150.         if ((b->_wds = x1 - b->_x) == 0)
  151.                 b->_x[0] = 0;
  152.         }
  153.  
  154. static int
  155. #ifdef KR_headers
  156. trailz(b) _Bigint *b;
  157. #else
  158. trailz(_Bigint *b)
  159. #endif
  160. {
  161.         __ULong L, *x, *xe;
  162.         int n = 0;
  163.  
  164.         x = b->_x;
  165.         xe = x + b->_wds;
  166.         for(n = 0; x < xe && !*x; x++)
  167.                 n += ULbits;
  168.         if (x < xe) {
  169.                 L = *x;
  170.                 n += lo0bits(&L);
  171.                 }
  172.         return n;
  173.         }
  174.  
  175.  _Bigint *
  176. #ifdef KR_headers
  177. increment(p, b) struct _reent *p; _Bigint *b;
  178. #else
  179. increment(struct _reent *p, _Bigint *b)
  180. #endif
  181. {
  182.         __ULong *x, *xe;
  183.         _Bigint *b1;
  184. #ifdef Pack_16
  185.         __ULong carry = 1, y;
  186. #endif
  187.  
  188.         x = b->_x;
  189.         xe = x + b->_wds;
  190. #ifdef Pack_32
  191.         do {
  192.                 if (*x < (__ULong)0xffffffffL) {
  193.                         ++*x;
  194.                         return b;
  195.                         }
  196.                 *x++ = 0;
  197.                 } while(x < xe);
  198. #else
  199.         do {
  200.                 y = *x + carry;
  201.                 carry = y >> 16;
  202.                 *x++ = y & 0xffff;
  203.                 if (!carry)
  204.                         return b;
  205.                 } while(x < xe);
  206.         if (carry)
  207. #endif
  208.         {
  209.                 if (b->_wds >= b->_maxwds) {
  210.                         b1 = Balloc(p,b->_k+1);
  211.                         Bcopy(b1,b);
  212.                         Bfree(p,b);
  213.                         b = b1;
  214.                         }
  215.                 b->_x[b->_wds++] = 1;
  216.                 }
  217.         return b;
  218.         }
  219.  
  220.  int
  221. #ifdef KR_headers
  222. decrement(b) _Bigint *b;
  223. #else
  224. decrement(_Bigint *b)
  225. #endif
  226. {
  227.         __ULong *x, *xe;
  228. #ifdef Pack_16
  229.         __ULong borrow = 1, y;
  230. #endif
  231.  
  232.         x = b->_x;
  233.         xe = x + b->_wds;
  234. #ifdef Pack_32
  235.         do {
  236.                 if (*x) {
  237.                         --*x;
  238.                         break;
  239.                         }
  240.                 *x++ = 0xffffffffL;
  241.                 }
  242.                 while(x < xe);
  243. #else
  244.         do {
  245.                 y = *x - borrow;
  246.                 borrow = (y & 0x10000) >> 16;
  247.                 *x++ = y & 0xffff;
  248.                 } while(borrow && x < xe);
  249. #endif
  250.         return STRTOG_Inexlo;
  251.         }
  252.  
  253.  static int
  254. #ifdef KR_headers
  255. all_on(b, n) _Bigint *b; int n;
  256. #else
  257. all_on(_Bigint *b, int n)
  258. #endif
  259. {
  260.         __ULong *x, *xe;
  261.  
  262.         x = b->_x;
  263.         xe = x + (n >> kshift);
  264.         while(x < xe)
  265.                 if ((*x++ & ALL_ON) != ALL_ON)
  266.                         return 0;
  267.         if (n &= kmask)
  268.                 return ((*x | (ALL_ON << n)) & ALL_ON) == ALL_ON;
  269.         return 1;
  270.         }
  271.  
  272.  _Bigint *
  273. #ifdef KR_headers
  274. set_ones(p, b, n) struct _reent *p; _Bigint *b; int n;
  275. #else
  276. set_ones(struct _reent *p, _Bigint *b, int n)
  277. #endif
  278. {
  279.         int k;
  280.         __ULong *x, *xe;
  281.  
  282.         k = (n + ((1 << kshift) - 1)) >> kshift;
  283.         if (b->_k < k) {
  284.                 Bfree(p,b);
  285.                 b = Balloc(p,k);
  286.                 }
  287.         k = n >> kshift;
  288.         if (n &= kmask)
  289.                 k++;
  290.         b->_wds = k;
  291.         x = b->_x;
  292.         xe = x + k;
  293.         while(x < xe)
  294.                 *x++ = ALL_ON;
  295.         if (n)
  296.                 x[-1] >>= ULbits - n;
  297.         return b;
  298.         }
  299.  
  300.  static int
  301. rvOK
  302. #ifdef KR_headers
  303.  (p, d, fpi, exp, bits, exact, rd, irv)
  304.  struct _reent *p; double d; FPI *fpi; Long *exp; __ULong *bits; int exact, rd, *irv;
  305. #else
  306.  (struct _reent *p, double d, FPI *fpi, Long *exp, __ULong *bits, int exact, int rd, int *irv)
  307. #endif
  308. {
  309.         _Bigint *b;
  310.         __ULong carry, inex, lostbits;
  311.         int bdif, e, j, k, k1, nb, rv;
  312.  
  313.         carry = rv = 0;
  314.         b = d2b(p, d, &e, &bdif);
  315.         bdif -= nb = fpi->nbits;
  316.         e += bdif;
  317.         if (bdif <= 0) {
  318.                 if (exact)
  319.                         goto trunc;
  320.                 goto ret;
  321.                 }
  322.         if (P == nb) {
  323.                 if (
  324. #ifndef IMPRECISE_INEXACT
  325.                         exact &&
  326. #endif
  327.                         fpi->rounding ==
  328. #ifdef RND_PRODQUOT
  329.                                         FPI_Round_near
  330. #else
  331.                                         Flt_Rounds
  332. #endif
  333.                         ) goto trunc;
  334.                 goto ret;
  335.                 }
  336.         switch(rd) {
  337.           case 1:
  338.                 goto trunc;
  339.           case 2:
  340.                 break;
  341.           default: /* round near */
  342.                 k = bdif - 1;
  343.                 if (k < 0)
  344.                         goto trunc;
  345.                 if (!k) {
  346.                         if (!exact)
  347.                                 goto ret;
  348.                         if (b->_x[0] & 2)
  349.                                 break;
  350.                         goto trunc;
  351.                         }
  352.                 if (b->_x[k>>kshift] & ((__ULong)1 << (k & kmask)))
  353.                         break;
  354.                 goto trunc;
  355.           }
  356.         /* "break" cases: round up 1 bit, then truncate; bdif > 0 */
  357.         carry = 1;
  358.  trunc:
  359.         inex = lostbits = 0;
  360.         if (bdif > 0) {
  361.                 if ( (lostbits = any_on(b, bdif)) !=0)
  362.                         inex = STRTOG_Inexlo;
  363.                 rshift(b, bdif);
  364.                 if (carry) {
  365.                         inex = STRTOG_Inexhi;
  366.                         b = increment(p, b);
  367.                         if ( (j = nb & kmask) !=0)
  368.                                 j = ULbits - j;
  369.                         if (hi0bits(b->_x[b->_wds - 1]) != j) {
  370.                                 if (!lostbits)
  371.                                         lostbits = b->_x[0] & 1;
  372.                                 rshift(b, 1);
  373.                                 e++;
  374.                                 }
  375.                         }
  376.                 }
  377.         else if (bdif < 0)
  378.                 b = lshift(p, b, -bdif);
  379.         if (e < fpi->emin) {
  380.                 k = fpi->emin - e;
  381.                 e = fpi->emin;
  382.                 if (k > nb || fpi->sudden_underflow) {
  383.                         b->_wds = inex = 0;
  384.                         *irv = STRTOG_Underflow | STRTOG_Inexlo;
  385.                         }
  386.                 else {
  387.                         k1 = k - 1;
  388.                         if (k1 > 0 && !lostbits)
  389.                                 lostbits = any_on(b, k1);
  390.                         if (!lostbits && !exact)
  391.                                 goto ret;
  392.                         lostbits |=
  393.                           carry = b->_x[k1>>kshift] & (1 << (k1 & kmask));
  394.                         rshift(b, k);
  395.                         *irv = STRTOG_Denormal;
  396.                         if (carry) {
  397.                                 b = increment(p, b);
  398.                                 inex = STRTOG_Inexhi | STRTOG_Underflow;
  399.                                 }
  400.                         else if (lostbits)
  401.                                 inex = STRTOG_Inexlo | STRTOG_Underflow;
  402.                         }
  403.                 }
  404.         else if (e > fpi->emax) {
  405.                 e = fpi->emax + 1;
  406.                 *irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
  407. #ifndef NO_ERRNO
  408.                 errno = ERANGE;
  409. #endif
  410.                 b->_wds = inex = 0;
  411.                 }
  412.         *exp = e;
  413.         copybits(bits, nb, b);
  414.         *irv |= inex;
  415.         rv = 1;
  416.  ret:
  417.         Bfree(p,b);
  418.         return rv;
  419.         }
  420.  
  421.  static int
  422. #ifdef KR_headers
  423. mantbits(d) double d;
  424. #else
  425. mantbits(U d)
  426. #endif
  427. {
  428.         __ULong L;
  429. #ifdef VAX
  430.         L = word1(d) << 16 | word1(d) >> 16;
  431.         if (L)
  432. #else
  433.         if ( (L = word1(d)) !=0)
  434. #endif
  435.                 return P - lo0bits(&L);
  436. #ifdef VAX
  437.         L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11;
  438. #else
  439.         L = word0(d) | Exp_msk1;
  440. #endif
  441.         return P - 32 - lo0bits(&L);
  442.         }
  443.  
  444.  int
  445. _strtodg_r
  446. #ifdef KR_headers
  447.         (p, s00, se, fpi, exp, bits)
  448.         struct _reent *p; const char *s00; char **se; FPI *fpi; Long *exp; __ULong *bits;
  449. #else
  450.         (struct _reent *p, const char *s00, char **se, FPI *fpi, Long *exp, __ULong *bits)
  451. #endif
  452. {
  453.         int abe, abits, asub;
  454.         int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm;
  455.         int dsign, e, e1, e2, emin, esign, finished, i, inex, irv;
  456.         int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
  457.         int sudden_underflow;
  458.         const char *s, *s0, *s1;
  459.         //double adj, adj0, rv, tol;
  460.         double adj0, tol;
  461.         U adj, rv;
  462.         Long L;
  463.         __ULong y, z;
  464.         _Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
  465.  
  466.         irv = STRTOG_Zero;
  467.         denorm = sign = nz0 = nz = 0;
  468.         dval(rv) = 0.;
  469.         rvb = 0;
  470.         nbits = fpi->nbits;
  471.         for(s = s00;;s++) switch(*s) {
  472.                 case '-':
  473.                         sign = 1;
  474.                         /* no break */
  475.                 case '+':
  476.                         if (*++s)
  477.                                 goto break2;
  478.                         /* no break */
  479.                 case 0:
  480.                         sign = 0;
  481.                         irv = STRTOG_NoNumber;
  482.                         s = s00;
  483.                         goto ret;
  484.                 case '\t':
  485.                 case '\n':
  486.                 case '\v':
  487.                 case '\f':
  488.                 case '\r':
  489.                 case ' ':
  490.                         continue;
  491.                 default:
  492.                         goto break2;
  493.                 }
  494.  break2:
  495.         if (*s == '0') {
  496. #ifndef NO_HEX_FP
  497.                 switch(s[1]) {
  498.                   case 'x':
  499.                   case 'X':
  500.                         irv = gethex(p, &s, fpi, exp, &rvb, sign);
  501.                         if (irv == STRTOG_NoNumber) {
  502.                                 s = s00;
  503.                                 sign = 0;
  504.                                 }
  505.                         goto ret;
  506.                   }
  507. #endif
  508.                 nz0 = 1;
  509.                 while(*++s == '0') ;
  510.                 if (!*s)
  511.                         goto ret;
  512.                 }
  513.         sudden_underflow = fpi->sudden_underflow;
  514.         s0 = s;
  515.         y = z = 0;
  516.         for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
  517.                 if (nd < 9)
  518.                         y = 10*y + c - '0';
  519.                 else if (nd < 16)
  520.                         z = 10*z + c - '0';
  521.         nd0 = nd;
  522. #ifdef USE_LOCALE
  523.         if (strncmp (s, _localeconv_r (p)->decimal_point,
  524.                      strlen (_localeconv_r (p)->decimal_point)) == 0)
  525. #else
  526.         if (c == '.')
  527. #endif
  528.                 {
  529.                 decpt = 1;
  530. #ifdef USE_LOCALE
  531.                 c = *(s += strlen (_localeconv_r (p)->decimal_point));
  532. #else
  533.                 c = *++s;
  534. #endif
  535.                 if (!nd) {
  536.                         for(; c == '0'; c = *++s)
  537.                                 nz++;
  538.                         if (c > '0' && c <= '9') {
  539.                                 s0 = s;
  540.                                 nf += nz;
  541.                                 nz = 0;
  542.                                 goto have_dig;
  543.                                 }
  544.                         goto dig_done;
  545.                         }
  546.                 for(; c >= '0' && c <= '9'; c = *++s) {
  547.  have_dig:
  548.                         nz++;
  549.                         if (c -= '0') {
  550.                                 nf += nz;
  551.                                 for(i = 1; i < nz; i++)
  552.                                         if (nd++ < 9)
  553.                                                 y *= 10;
  554.                                         else if (nd <= DBL_DIG + 1)
  555.                                                 z *= 10;
  556.                                 if (nd++ < 9)
  557.                                         y = 10*y + c;
  558.                                 else if (nd <= DBL_DIG + 1)
  559.                                         z = 10*z + c;
  560.                                 nz = 0;
  561.                                 }
  562.                         }
  563.                 }
  564.  dig_done:
  565.         e = 0;
  566.         if (c == 'e' || c == 'E') {
  567.                 if (!nd && !nz && !nz0) {
  568.                         irv = STRTOG_NoNumber;
  569.                         s = s00;
  570.                         goto ret;
  571.                         }
  572.                 s00 = s;
  573.                 esign = 0;
  574.                 switch(c = *++s) {
  575.                         case '-':
  576.                                 esign = 1;
  577.                         case '+':
  578.                                 c = *++s;
  579.                         }
  580.                 if (c >= '0' && c <= '9') {
  581.                         while(c == '0')
  582.                                 c = *++s;
  583.                         if (c > '0' && c <= '9') {
  584.                                 L = c - '0';
  585.                                 s1 = s;
  586.                                 while((c = *++s) >= '0' && c <= '9')
  587.                                         L = 10*L + c - '0';
  588.                                 if (s - s1 > 8 || L > 19999)
  589.                                         /* Avoid confusion from exponents
  590.                                          * so large that e might overflow.
  591.                                          */
  592.                                         e = 19999; /* safe for 16 bit ints */
  593.                                 else
  594.                                         e = (int)L;
  595.                                 if (esign)
  596.                                         e = -e;
  597.                                 }
  598.                         else
  599.                                 e = 0;
  600.                         }
  601.                 else
  602.                         s = s00;
  603.                 }
  604.         if (!nd) {
  605.                 if (!nz && !nz0) {
  606. #ifdef INFNAN_CHECK
  607.                         /* Check for Nan and Infinity */
  608.                         if (!decpt)
  609.                          switch(c) {
  610.                           case 'i':
  611.                           case 'I':
  612.                                 if (match(&s,"nf")) {
  613.                                         --s;
  614.                                         if (!match(&s,"inity"))
  615.                                                 ++s;
  616.                                         irv = STRTOG_Infinite;
  617.                                         goto infnanexp;
  618.                                         }
  619.                                 break;
  620.                           case 'n':
  621.                           case 'N':
  622.                                 if (match(&s, "an")) {
  623.                                         irv = STRTOG_NaN;
  624.                                         *exp = fpi->emax + 1;
  625. #ifndef No_Hex_NaN
  626.                                         if (*s == '(') /*)*/
  627.                                                 irv = hexnan(&s, fpi, bits);
  628. #endif
  629.                                         goto infnanexp;
  630.                                         }
  631.                           }
  632. #endif /* INFNAN_CHECK */
  633.                         irv = STRTOG_NoNumber;
  634.                         s = s00;
  635.                         }
  636.                 goto ret;
  637.                 }
  638.  
  639.         irv = STRTOG_Normal;
  640.         e1 = e -= nf;
  641.         rd = 0;
  642.         switch(fpi->rounding & 3) {
  643.           case FPI_Round_up:
  644.                 rd = 2 - sign;
  645.                 break;
  646.           case FPI_Round_zero:
  647.                 rd = 1;
  648.                 break;
  649.           case FPI_Round_down:
  650.                 rd = 1 + sign;
  651.           }
  652.  
  653.         /* Now we have nd0 digits, starting at s0, followed by a
  654.          * decimal point, followed by nd-nd0 digits.  The number we're
  655.          * after is the integer represented by those digits times
  656.          * 10**e */
  657.  
  658.         if (!nd0)
  659.                 nd0 = nd;
  660.         k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
  661.         dval(rv) = y;
  662.         if (k > 9)
  663.                 dval(rv) = tens[k - 9] * dval(rv) + z;
  664.         bd0 = 0;
  665.         if (nbits <= P && nd <= DBL_DIG) {
  666.                 if (!e) {
  667.                         if (rvOK(p, dval(rv), fpi, exp, bits, 1, rd, &irv))
  668.                                 goto ret;
  669.                         }
  670.                 else if (e > 0) {
  671.                         if (e <= Ten_pmax) {
  672. #ifdef VAX
  673.                                 goto vax_ovfl_check;
  674. #else
  675.                                 i = fivesbits[e] + mantbits(rv) <= P;
  676.                                 /* rv = */ rounded_product(dval(rv), tens[e]);
  677.                                 if (rvOK(p, dval(rv), fpi, exp, bits, i, rd, &irv))
  678.                                         goto ret;
  679.                                 e1 -= e;
  680.                                 goto rv_notOK;
  681. #endif
  682.                                 }
  683.                         i = DBL_DIG - nd;
  684.                         if (e <= Ten_pmax + i) {
  685.                                 /* A fancier test would sometimes let us do
  686.                                  * this for larger i values.
  687.                                  */
  688.                                 e2 = e - i;
  689.                                 e1 -= i;
  690.                                 dval(rv) *= tens[i];
  691. #ifdef VAX
  692.                                 /* VAX exponent range is so narrow we must
  693.                                  * worry about overflow here...
  694.                                  */
  695.  vax_ovfl_check:
  696.                                 dval(adj) = dval(rv);
  697.                                 word0(adj) -= P*Exp_msk1;
  698.                                 /* adj = */ rounded_product(dval(adj), tens[e2]);
  699.                                 if ((word0(adj) & Exp_mask)
  700.                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
  701.                                         goto rv_notOK;
  702.                                 word0(adj) += P*Exp_msk1;
  703.                                 dval(rv) = dval(adj);
  704. #else
  705.                                 /* rv = */ rounded_product(dval(rv), tens[e2]);
  706. #endif
  707.                                 if (rvOK(p, dval(rv), fpi, exp, bits, 0, rd, &irv))
  708.                                         goto ret;
  709.                                 e1 -= e2;
  710.                                 }
  711.                         }
  712. #ifndef Inaccurate_Divide
  713.                 else if (e >= -Ten_pmax) {
  714.                         /* rv = */ rounded_quotient(dval(rv), tens[-e]);
  715.                         if (rvOK(p, dval(rv), fpi, exp, bits, 0, rd, &irv))
  716.                                 goto ret;
  717.                         e1 -= e;
  718.                         }
  719. #endif
  720.                 }
  721.  rv_notOK:
  722.         e1 += nd - k;
  723.  
  724.         /* Get starting approximation = rv * 10**e1 */
  725.  
  726.         e2 = 0;
  727.         if (e1 > 0) {
  728.                 if ( (i = e1 & 15) !=0)
  729.                         dval(rv) *= tens[i];
  730.                 if (e1 &= ~15) {
  731.                         e1 >>= 4;
  732.                         while(e1 >= (1 << n_bigtens-1)) {
  733.                                 e2 += ((word0(rv) & Exp_mask)
  734.                                         >> Exp_shift1) - Bias;
  735.                                 word0(rv) &= ~Exp_mask;
  736.                                 word0(rv) |= Bias << Exp_shift1;
  737.                                 dval(rv) *= bigtens[n_bigtens-1];
  738.                                 e1 -= 1 << n_bigtens-1;
  739.                                 }
  740.                         e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
  741.                         word0(rv) &= ~Exp_mask;
  742.                         word0(rv) |= Bias << Exp_shift1;
  743.                         for(j = 0; e1 > 0; j++, e1 >>= 1)
  744.                                 if (e1 & 1)
  745.                                         dval(rv) *= bigtens[j];
  746.                         }
  747.                 }
  748.         else if (e1 < 0) {
  749.                 e1 = -e1;
  750.                 if ( (i = e1 & 15) !=0)
  751.                         dval(rv) /= tens[i];
  752.                 if (e1 &= ~15) {
  753.                         e1 >>= 4;
  754.                         while(e1 >= (1 << n_bigtens-1)) {
  755.                                 e2 += ((word0(rv) & Exp_mask)
  756.                                         >> Exp_shift1) - Bias;
  757.                                 word0(rv) &= ~Exp_mask;
  758.                                 word0(rv) |= Bias << Exp_shift1;
  759.                                 dval(rv) *= tinytens[n_bigtens-1];
  760.                                 e1 -= 1 << n_bigtens-1;
  761.                                 }
  762.                         e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
  763.                         word0(rv) &= ~Exp_mask;
  764.                         word0(rv) |= Bias << Exp_shift1;
  765.                         for(j = 0; e1 > 0; j++, e1 >>= 1)
  766.                                 if (e1 & 1)
  767.                                         dval(rv) *= tinytens[j];
  768.                         }
  769.                 }
  770. #ifdef IBM
  771.         /* e2 is a correction to the (base 2) exponent of the return
  772.          * value, reflecting adjustments above to avoid overflow in the
  773.          * native arithmetic.  For native IBM (base 16) arithmetic, we
  774.          * must multiply e2 by 4 to change from base 16 to 2.
  775.          */
  776.         e2 <<= 2;
  777. #endif
  778.         rvb = d2b(p, dval(rv), &rve, &rvbits);  /* rv = rvb * 2^rve */
  779.         rve += e2;
  780.         if ((j = rvbits - nbits) > 0) {
  781.                 rshift(rvb, j);
  782.                 rvbits = nbits;
  783.                 rve += j;
  784.                 }
  785.         bb0 = 0;        /* trailing zero bits in rvb */
  786.         e2 = rve + rvbits - nbits;
  787.         if (e2 > fpi->emax + 1)
  788.                 goto huge;
  789.         rve1 = rve + rvbits - nbits;
  790.         if (e2 < (emin = fpi->emin)) {
  791.                 denorm = 1;
  792.                 j = rve - emin;
  793.                 if (j > 0) {
  794.                         rvb = lshift(p, rvb, j);
  795.                         rvbits += j;
  796.                         }
  797.                 else if (j < 0) {
  798.                         rvbits += j;
  799.                         if (rvbits <= 0) {
  800.                                 if (rvbits < -1) {
  801.  ufl:
  802.                                         rvb->_wds = 0;
  803.                                         rvb->_x[0] = 0;
  804.                                         *exp = emin;
  805.                                         irv = STRTOG_Underflow | STRTOG_Inexlo;
  806.                                         goto ret;
  807.                                         }
  808.                                 rvb->_x[0] = rvb->_wds = rvbits = 1;
  809.                                 }
  810.                         else
  811.                                 rshift(rvb, -j);
  812.                         }
  813.                 rve = rve1 = emin;
  814.                 if (sudden_underflow && e2 + 1 < emin)
  815.                         goto ufl;
  816.                 }
  817.  
  818.         /* Now the hard part -- adjusting rv to the correct value.*/
  819.  
  820.         /* Put digits into bd: true value = bd * 10^e */
  821.  
  822.         bd0 = s2b(p, s0, nd0, nd, y);
  823.  
  824.         for(;;) {
  825.                 bd = Balloc(p,bd0->_k);
  826.                 Bcopy(bd, bd0);
  827.                 bb = Balloc(p,rvb->_k);
  828.                 Bcopy(bb, rvb);
  829.                 bbbits = rvbits - bb0;
  830.                 bbe = rve + bb0;
  831.                 bs = i2b(p, 1);
  832.  
  833.                 if (e >= 0) {
  834.                         bb2 = bb5 = 0;
  835.                         bd2 = bd5 = e;
  836.                         }
  837.                 else {
  838.                         bb2 = bb5 = -e;
  839.                         bd2 = bd5 = 0;
  840.                         }
  841.                 if (bbe >= 0)
  842.                         bb2 += bbe;
  843.                 else
  844.                         bd2 -= bbe;
  845.                 bs2 = bb2;
  846.                 j = nbits + 1 - bbbits;
  847.                 i = bbe + bbbits - nbits;
  848.                 if (i < emin)   /* denormal */
  849.                         j += i - emin;
  850.                 bb2 += j;
  851.                 bd2 += j;
  852.                 i = bb2 < bd2 ? bb2 : bd2;
  853.                 if (i > bs2)
  854.                         i = bs2;
  855.                 if (i > 0) {
  856.                         bb2 -= i;
  857.                         bd2 -= i;
  858.                         bs2 -= i;
  859.                         }
  860.                 if (bb5 > 0) {
  861.                         bs = pow5mult(p, bs, bb5);
  862.                         bb1 = mult(p, bs, bb);
  863.                         Bfree(p,bb);
  864.                         bb = bb1;
  865.                         }
  866.                 bb2 -= bb0;
  867.                 if (bb2 > 0)
  868.                         bb = lshift(p, bb, bb2);
  869.                 else if (bb2 < 0)
  870.                         rshift(bb, -bb2);
  871.                 if (bd5 > 0)
  872.                         bd = pow5mult(p, bd, bd5);
  873.                 if (bd2 > 0)
  874.                         bd = lshift(p, bd, bd2);
  875.                 if (bs2 > 0)
  876.                         bs = lshift(p, bs, bs2);
  877.                 asub = 1;
  878.                 inex = STRTOG_Inexhi;
  879.                 delta = diff(p, bb, bd);
  880.                 if (delta->_wds <= 1 && !delta->_x[0])
  881.                         break;
  882.                 dsign = delta->_sign;
  883.                 delta->_sign = finished = 0;
  884.                 L = 0;
  885.                 i = cmp(delta, bs);
  886.                 if (rd && i <= 0) {
  887.                         irv = STRTOG_Normal;
  888.                         if ( (finished = dsign ^ (rd&1)) !=0) {
  889.                                 if (dsign != 0) {
  890.                                         irv |= STRTOG_Inexhi;
  891.                                         goto adj1;
  892.                                         }
  893.                                 irv |= STRTOG_Inexlo;
  894.                                 if (rve1 == emin)
  895.                                         goto adj1;
  896.                                 for(i = 0, j = nbits; j >= ULbits;
  897.                                                 i++, j -= ULbits) {
  898.                                         if (rvb->_x[i] & ALL_ON)
  899.                                                 goto adj1;
  900.                                         }
  901.                                 if (j > 1 && lo0bits(rvb->_x + i) < j - 1)
  902.                                         goto adj1;
  903.                                 rve = rve1 - 1;
  904.                                 rvb = set_ones(p, rvb, rvbits = nbits);
  905.                                 break;
  906.                                 }
  907.                         irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
  908.                         break;
  909.                         }
  910.                 if (i < 0) {
  911.                         /* Error is less than half an ulp -- check for
  912.                          * special case of mantissa a power of two.
  913.                          */
  914.                         irv = dsign
  915.                                 ? STRTOG_Normal | STRTOG_Inexlo
  916.                                 : STRTOG_Normal | STRTOG_Inexhi;
  917.                         if (dsign || bbbits > 1 || denorm || rve1 == emin)
  918.                                 break;
  919.                         delta = lshift(p, delta,1);
  920.                         if (cmp(delta, bs) > 0) {
  921.                                 irv = STRTOG_Normal | STRTOG_Inexlo;
  922.                                 goto drop_down;
  923.                                 }
  924.                         break;
  925.                         }
  926.                 if (i == 0) {
  927.                         /* exactly half-way between */
  928.                         if (dsign) {
  929.                                 if (denorm && all_on(rvb, rvbits)) {
  930.                                         /*boundary case -- increment exponent*/
  931.                                         rvb->_wds = 1;
  932.                                         rvb->_x[0] = 1;
  933.                                         rve = emin + nbits - (rvbits = 1);
  934.                                         irv = STRTOG_Normal | STRTOG_Inexhi;
  935.                                         denorm = 0;
  936.                                         break;
  937.                                         }
  938.                                 irv = STRTOG_Normal | STRTOG_Inexlo;
  939.                                 }
  940.                         else if (bbbits == 1) {
  941.                                 irv = STRTOG_Normal;
  942.  drop_down:
  943.                                 /* boundary case -- decrement exponent */
  944.                                 if (rve1 == emin) {
  945.                                         irv = STRTOG_Normal | STRTOG_Inexhi;
  946.                                         if (rvb->_wds == 1 && rvb->_x[0] == 1)
  947.                                                 sudden_underflow = 1;
  948.                                         break;
  949.                                         }
  950.                                 rve -= nbits;
  951.                                 rvb = set_ones(p, rvb, rvbits = nbits);
  952.                                 break;
  953.                                 }
  954.                         else
  955.                                 irv = STRTOG_Normal | STRTOG_Inexhi;
  956.                         if (bbbits < nbits && !denorm || !(rvb->_x[0] & 1))
  957.                                 break;
  958.                         if (dsign) {
  959.                                 rvb = increment(p, rvb);
  960.                                 j = kmask & (ULbits - (rvbits & kmask));
  961.                                 if (hi0bits(rvb->_x[rvb->_wds - 1]) != j)
  962.                                         rvbits++;
  963.                                 irv = STRTOG_Normal | STRTOG_Inexhi;
  964.                                 }
  965.                         else {
  966.                                 if (bbbits == 1)
  967.                                         goto undfl;
  968.                                 decrement(rvb);
  969.                                 irv = STRTOG_Normal | STRTOG_Inexlo;
  970.                                 }
  971.                         break;
  972.                         }
  973.                 if ((dval(adj) = ratio(delta, bs)) <= 2.) {
  974.  adj1:
  975.                         inex = STRTOG_Inexlo;
  976.                         if (dsign) {
  977.                                 asub = 0;
  978.                                 inex = STRTOG_Inexhi;
  979.                                 }
  980.                         else if (denorm && bbbits <= 1) {
  981.  undfl:
  982.                                 rvb->_wds = 0;
  983.                                 rve = emin;
  984.                                 irv = STRTOG_Underflow | STRTOG_Inexlo;
  985.                                 break;
  986.                                 }
  987.                         adj0 = dval(adj) = 1.;
  988.                         }
  989.                 else {
  990.                         adj0 = dval(adj) *= 0.5;
  991.                         if (dsign) {
  992.                                 asub = 0;
  993.                                 inex = STRTOG_Inexlo;
  994.                                 }
  995.                         if (dval(adj) < 2147483647.) {
  996.                                 L = adj0;
  997.                                 adj0 -= L;
  998.                                 switch(rd) {
  999.                                   case 0:
  1000.                                         if (adj0 >= .5)
  1001.                                                 goto inc_L;
  1002.                                         break;
  1003.                                   case 1:
  1004.                                         if (asub && adj0 > 0.)
  1005.                                                 goto inc_L;
  1006.                                         break;
  1007.                                   case 2:
  1008.                                         if (!asub && adj0 > 0.) {
  1009.  inc_L:
  1010.                                                 L++;
  1011.                                                 inex = STRTOG_Inexact - inex;
  1012.                                                 }
  1013.                                   }
  1014.                                 dval(adj) = L;
  1015.                                 }
  1016.                         }
  1017.                 y = rve + rvbits;
  1018.  
  1019.                 /* adj *= ulp(dval(rv)); */
  1020.                 /* if (asub) rv -= adj; else rv += adj; */
  1021.  
  1022.                 if (!denorm && rvbits < nbits) {
  1023.                         rvb = lshift(p, rvb, j = nbits - rvbits);
  1024.                         rve -= j;
  1025.                         rvbits = nbits;
  1026.                         }
  1027.                 ab = d2b(p, dval(adj), &abe, &abits);
  1028.                 if (abe < 0)
  1029.                         rshift(ab, -abe);
  1030.                 else if (abe > 0)
  1031.                         ab = lshift(p, ab, abe);
  1032.                 rvb0 = rvb;
  1033.                 if (asub) {
  1034.                         /* rv -= adj; */
  1035.                         j = hi0bits(rvb->_x[rvb->_wds-1]);
  1036.                         rvb = diff(p, rvb, ab);
  1037.                         k = rvb0->_wds - 1;
  1038.                         if (denorm)
  1039.                                 /* do nothing */;
  1040.                         else if (rvb->_wds <= k
  1041.                                 || hi0bits( rvb->_x[k]) >
  1042.                                    hi0bits(rvb0->_x[k])) {
  1043.                                 /* unlikely; can only have lost 1 high bit */
  1044.                                 if (rve1 == emin) {
  1045.                                         --rvbits;
  1046.                                         denorm = 1;
  1047.                                         }
  1048.                                 else {
  1049.                                         rvb = lshift(p, rvb, 1);
  1050.                                         --rve;
  1051.                                         --rve1;
  1052.                                         L = finished = 0;
  1053.                                         }
  1054.                                 }
  1055.                         }
  1056.                 else {
  1057.                         rvb = sum(p, rvb, ab);
  1058.                         k = rvb->_wds - 1;
  1059.                         if (k >= rvb0->_wds
  1060.                          || hi0bits(rvb->_x[k]) < hi0bits(rvb0->_x[k])) {
  1061.                                 if (denorm) {
  1062.                                         if (++rvbits == nbits)
  1063.                                                 denorm = 0;
  1064.                                         }
  1065.                                 else {
  1066.                                         rshift(rvb, 1);
  1067.                                         rve++;
  1068.                                         rve1++;
  1069.                                         L = 0;
  1070.                                         }
  1071.                                 }
  1072.                         }
  1073.                 Bfree(p,ab);
  1074.                 Bfree(p,rvb0);
  1075.                 if (finished)
  1076.                         break;
  1077.  
  1078.                 z = rve + rvbits;
  1079.                 if (y == z && L) {
  1080.                         /* Can we stop now? */
  1081.                         tol = dval(adj) * 5e-16; /* > max rel error */
  1082.                         dval(adj) = adj0 - .5;
  1083.                         if (dval(adj) < -tol) {
  1084.                                 if (adj0 > tol) {
  1085.                                         irv |= inex;
  1086.                                         break;
  1087.                                         }
  1088.                                 }
  1089.                         else if (dval(adj) > tol && adj0 < 1. - tol) {
  1090.                                 irv |= inex;
  1091.                                 break;
  1092.                                 }
  1093.                         }
  1094.                 bb0 = denorm ? 0 : trailz(rvb);
  1095.                 Bfree(p,bb);
  1096.                 Bfree(p,bd);
  1097.                 Bfree(p,bs);
  1098.                 Bfree(p,delta);
  1099.                 }
  1100.         if (!denorm && (j = nbits - rvbits)) {
  1101.                 if (j > 0)
  1102.                         rvb = lshift(p, rvb, j);
  1103.                 else
  1104.                         rshift(rvb, -j);
  1105.                 rve -= j;
  1106.                 }
  1107.         *exp = rve;
  1108.         Bfree(p,bb);
  1109.         Bfree(p,bd);
  1110.         Bfree(p,bs);
  1111.         Bfree(p,bd0);
  1112.         Bfree(p,delta);
  1113.         if (rve > fpi->emax) {
  1114.  huge:
  1115.                 rvb->_wds = 0;
  1116.                 irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
  1117. #ifndef NO_ERRNO
  1118.                 errno = ERANGE;
  1119. #endif
  1120.  infnanexp:
  1121.                 *exp = fpi->emax + 1;
  1122.                 }
  1123.  ret:
  1124.         if (denorm) {
  1125.                 if (sudden_underflow) {
  1126.                         rvb->_wds = 0;
  1127.                         irv = STRTOG_Underflow | STRTOG_Inexlo;
  1128.                         }
  1129.                 else  {
  1130.                         irv = (irv & ~STRTOG_Retmask) |
  1131.                                 (rvb->_wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
  1132.                         if (irv & STRTOG_Inexact)
  1133.                                 irv |= STRTOG_Underflow;
  1134.                         }
  1135.                 }
  1136.         if (se)
  1137.                 *se = (char *)s;
  1138.         if (sign)
  1139.                 irv |= STRTOG_Neg;
  1140.         if (rvb) {
  1141.                 copybits(bits, nbits, rvb);
  1142.                 Bfree(p,rvb);
  1143.                 }
  1144.         return irv;
  1145.         }
  1146.  
  1147. #endif /* _HAVE_LONG_DOUBLE && !_LDBL_EQ_DBL */
  1148.