1/** @file
2
3  Copyright (c) 2012, Intel Corporation. All rights reserved.<BR>
4  This program and the accompanying materials are licensed and made available under
5  the terms and conditions of the BSD License that accompanies this distribution.
6  The full text of the license may be found at
7  http://opensource.org/licenses/bsd-license.
8
9  THE PROGRAM IS DISTRIBUTED UNDER THE BSD LICENSE ON AN "AS IS" BASIS,
10  WITHOUT WARRANTIES OR REPRESENTATIONS OF ANY KIND, EITHER EXPRESS OR IMPLIED.
11
12  *****************************************************************
13
14  The author of this software is David M. Gay.
15
16  Copyright (C) 1998-2001 by Lucent Technologies
17  All Rights Reserved
18
19  Permission to use, copy, modify, and distribute this software and
20  its documentation for any purpose and without fee is hereby
21  granted, provided that the above copyright notice appear in all
22  copies and that both that the copyright notice and this
23  permission notice and warranty disclaimer appear in supporting
24  documentation, and that the name of Lucent or any of its entities
25  not be used in advertising or publicity pertaining to
26  distribution of the software without specific, written prior
27  permission.
28
29  LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
30  INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
31  IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
32  SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
33  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
34  IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
35  ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
36  THIS SOFTWARE.
37
38
39  Please send bug reports to David M. Gay (dmg at acm dot org,
40  with " at " changed at "@" and " dot " changed to ".").
41
42  *****************************************************************
43
44  NetBSD: strtod.c,v 1.4.14.1 2008/04/08 21:10:55 jdc Exp
45**/
46#include  <LibConfig.h>
47
48#include "gdtoaimp.h"
49#ifndef NO_FENV_H
50#include <fenv.h>
51#endif
52
53#ifdef USE_LOCALE
54#include "locale.h"
55#endif
56
57#ifdef IEEE_Arith
58#ifndef NO_IEEE_Scale
59#define Avoid_Underflow
60#undef tinytens
61/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
62/* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
63static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
64    9007199254740992.e-256
65    };
66#endif
67#endif
68
69#ifdef Honor_FLT_ROUNDS
70#define Rounding rounding
71#undef Check_FLT_ROUNDS
72#define Check_FLT_ROUNDS
73#else
74#define Rounding Flt_Rounds
75#endif
76
77//#ifndef __HAVE_LONG_DOUBLE
78//__strong_alias(_strtold, strtod)
79//__weak_alias(strtold, _strtold)
80//#endif
81
82#if defined(_MSC_VER)           /* Handle Microsoft VC++ compiler specifics. */
83// Disable: warning C4700: uninitialized local variable 'xx' used
84#pragma warning ( disable : 4700 )
85#endif  /* defined(_MSC_VER) */
86
87double
88strtod(CONST char *s00, char **se)
89{
90#ifdef Avoid_Underflow
91  int scale;
92  #endif
93  int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
94      e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
95  CONST char *s, *s0, *s1;
96  double aadj, aadj1, adj, rv, rv0;
97  Long L;
98  ULong y, z;
99  Bigint *bb = NULL, *bb1, *bd0;
100  Bigint *bd = NULL, *bs = NULL, *delta = NULL; /* pacify gcc */
101#ifdef SET_INEXACT
102  int inexact, oldinexact;
103#endif
104#ifdef Honor_FLT_ROUNDS
105  int rounding;
106#endif
107
108  sign = nz0 = nz = decpt = 0;
109  dval(rv) = 0.;
110  for(s = s00;;s++) {
111    switch(*s) {
112      case '-':
113        sign = 1;
114        /* FALLTHROUGH */
115      case '+':
116        if (*++s)
117          goto break2;
118        /* FALLTHROUGH */
119      case 0:
120        goto ret0;
121      case '\t':
122      case '\n':
123      case '\v':
124      case '\f':
125      case '\r':
126      case ' ':
127        continue;
128      default:
129        goto break2;
130    }
131  }
132 break2:
133  if (*s == '0') {
134#ifndef NO_HEX_FP
135    {
136    static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
137    Long expt;
138    ULong bits[2];
139    switch(s[1]) {
140      case 'x':
141      case 'X':
142      {
143#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD)
144      FPI fpi1 = fpi;
145      switch(fegetround()) {
146        case FE_TOWARDZERO: fpi1.rounding = 0; break;
147        case FE_UPWARD: fpi1.rounding = 2; break;
148        case FE_DOWNWARD: fpi1.rounding = 3;
149        }
150#else
151#endif
152      switch((i = gethex(&s, &fpi, &expt, &bb, sign)) & STRTOG_Retmask) {
153        case STRTOG_NoNumber:
154        s = s00;
155        sign = 0;
156        /* FALLTHROUGH */
157        case STRTOG_Zero:
158        break;
159        default:
160        if (bb) {
161          copybits(bits, fpi.nbits, bb);
162          Bfree(bb);
163          }
164        ULtod((/* LINTED */(U*)&rv)->L, bits, expt, i);
165        }}
166      goto ret;
167      }
168    }
169#endif
170    nz0 = 1;
171    while(*++s == '0') ;
172    if (!*s)
173      goto ret;
174  }
175  s0 = s;
176  y = z = 0;
177  for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
178    if (nd < 9)
179      y = 10*y + c - '0';
180    else if (nd < 16)
181      z = 10*z + c - '0';
182  nd0 = nd;
183#ifdef USE_LOCALE
184  if (c == *localeconv()->decimal_point)
185#else
186  if (c == '.')
187#endif
188    {
189    decpt = 1;
190    c = *++s;
191    if (!nd) {
192      for(; c == '0'; c = *++s)
193        nz++;
194      if (c > '0' && c <= '9') {
195        s0 = s;
196        nf += nz;
197        nz = 0;
198        goto have_dig;
199        }
200      goto dig_done;
201      }
202    for(; c >= '0' && c <= '9'; c = *++s) {
203 have_dig:
204      nz++;
205      if (c -= '0') {
206        nf += nz;
207        for(i = 1; i < nz; i++)
208          if (nd++ < 9)
209            y *= 10;
210          else if (nd <= DBL_DIG + 1)
211            z *= 10;
212        if (nd++ < 9)
213          y = 10*y + c;
214        else if (nd <= DBL_DIG + 1)
215          z = 10*z + c;
216        nz = 0;
217        }
218      }
219    }
220 dig_done:
221  e = 0;
222  if (c == 'e' || c == 'E') {
223    if (!nd && !nz && !nz0) {
224      goto ret0;
225      }
226    s00 = s;
227    esign = 0;
228    switch(c = *++s) {
229      case '-':
230        esign = 1;
231        /* FALLTHROUGH */
232      case '+':
233        c = *++s;
234      }
235    if (c >= '0' && c <= '9') {
236      while(c == '0')
237        c = *++s;
238      if (c > '0' && c <= '9') {
239        L = c - '0';
240        s1 = s;
241        while((c = *++s) >= '0' && c <= '9')
242          L = 10*L + c - '0';
243        if (s - s1 > 8 || L > 19999)
244          /* Avoid confusion from exponents
245           * so large that e might overflow.
246           */
247          e = 19999; /* safe for 16 bit ints */
248        else
249          e = (int)L;
250        if (esign)
251          e = -e;
252        }
253      else
254        e = 0;
255      }
256    else
257      s = s00;
258    }
259  if (!nd) {
260    if (!nz && !nz0) {
261#ifdef INFNAN_CHECK
262      /* Check for Nan and Infinity */
263#ifndef No_Hex_NaN
264      ULong bits[2];
265      static FPI fpinan = /* only 52 explicit bits */
266        { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
267#endif  // No_Hex_NaN
268      if (!decpt)
269       switch(c) {
270        case 'i':
271        case 'I':
272        if (match(&s,"nf")) {
273          --s;
274          if (!match(&s,"inity"))
275            ++s;
276          word0(rv) = 0x7ff00000;
277          word1(rv) = 0;
278          goto ret;
279          }
280        break;
281        case 'n':
282        case 'N':
283        if (match(&s, "an")) {
284#ifndef No_Hex_NaN
285          if (*s == '(' /*)*/
286           && hexnan(&s, &fpinan, bits)
287              == STRTOG_NaNbits) {
288            word0(rv) = (UINT32)(0x7ff00000U | bits[1]);
289            word1(rv) = (UINT32)bits[0];
290            }
291          else {
292#endif
293            word0(rv) = NAN_WORD0;
294            word1(rv) = NAN_WORD1;
295#ifndef No_Hex_NaN
296            }
297#endif
298          goto ret;
299          }
300        }
301#endif /* INFNAN_CHECK */
302 ret0:
303      s = s00;
304      sign = 0;
305      }
306    goto ret;
307    }
308  e1 = e -= nf;
309
310  /* Now we have nd0 digits, starting at s0, followed by a
311   * decimal point, followed by nd-nd0 digits.  The number we're
312   * after is the integer represented by those digits times
313   * 10**e */
314
315  if (!nd0)
316    nd0 = nd;
317  k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
318  dval(rv) = (double)y;
319  if (k > 9) {
320#ifdef SET_INEXACT
321    if (k > DBL_DIG)
322      oldinexact = get_inexact();
323#endif
324    dval(rv) = tens[k - 9] * dval(rv) + z;
325    }
326  bd0 = 0;
327  if (nd <= DBL_DIG
328#ifndef RND_PRODQUOT
329#ifndef Honor_FLT_ROUNDS
330    && Flt_Rounds == 1
331#endif
332#endif
333      ) {
334    if (!e)
335      goto ret;
336    if (e > 0) {
337      if (e <= Ten_pmax) {
338#ifdef VAX
339        goto vax_ovfl_check;
340#else
341#ifdef Honor_FLT_ROUNDS
342        /* round correctly FLT_ROUNDS = 2 or 3 */
343        if (sign) {
344          rv = -rv;
345          sign = 0;
346          }
347#endif
348        /* rv = */ rounded_product(dval(rv), tens[e]);
349        goto ret;
350#endif
351        }
352      i = DBL_DIG - nd;
353      if (e <= Ten_pmax + i) {
354        /* A fancier test would sometimes let us do
355         * this for larger i values.
356         */
357#ifdef Honor_FLT_ROUNDS
358        /* round correctly FLT_ROUNDS = 2 or 3 */
359        if (sign) {
360          rv = -rv;
361          sign = 0;
362          }
363#endif
364        e -= i;
365        dval(rv) *= tens[i];
366#ifdef VAX
367        /* VAX exponent range is so narrow we must
368         * worry about overflow here...
369         */
370 vax_ovfl_check:
371        word0(rv) -= P*Exp_msk1;
372        /* rv = */ rounded_product(dval(rv), tens[e]);
373        if ((word0(rv) & Exp_mask)
374         > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
375          goto ovfl;
376        word0(rv) += P*Exp_msk1;
377#else
378        /* rv = */ rounded_product(dval(rv), tens[e]);
379#endif
380        goto ret;
381        }
382      }
383#ifndef Inaccurate_Divide
384    else if (e >= -Ten_pmax) {
385#ifdef Honor_FLT_ROUNDS
386      /* round correctly FLT_ROUNDS = 2 or 3 */
387      if (sign) {
388        rv = -rv;
389        sign = 0;
390        }
391#endif
392      /* rv = */ rounded_quotient(dval(rv), tens[-e]);
393      goto ret;
394      }
395#endif
396    }
397  e1 += nd - k;
398
399#ifdef IEEE_Arith
400#ifdef SET_INEXACT
401  inexact = 1;
402  if (k <= DBL_DIG)
403    oldinexact = get_inexact();
404#endif
405#ifdef Avoid_Underflow
406  scale = 0;
407#endif
408#ifdef Honor_FLT_ROUNDS
409  if ((rounding = Flt_Rounds) >= 2) {
410    if (sign)
411      rounding = rounding == 2 ? 0 : 2;
412    else
413      if (rounding != 2)
414        rounding = 0;
415    }
416#endif
417#endif /*IEEE_Arith*/
418
419  /* Get starting approximation = rv * 10**e1 */
420
421  if (e1 > 0) {
422    if ( (i = e1 & 15) !=0)
423      dval(rv) *= tens[i];
424    if (e1 &= ~15) {
425      if (e1 > DBL_MAX_10_EXP) {
426 ovfl:
427#ifndef NO_ERRNO
428        errno = ERANGE;
429#endif
430        /* Can't trust HUGE_VAL */
431#ifdef IEEE_Arith
432#ifdef Honor_FLT_ROUNDS
433        switch(rounding) {
434          case 0: /* toward 0 */
435          case 3: /* toward -infinity */
436          word0(rv) = Big0;
437          word1(rv) = Big1;
438          break;
439          default:
440          word0(rv) = Exp_mask;
441          word1(rv) = 0;
442          }
443#else /*Honor_FLT_ROUNDS*/
444        word0(rv) = Exp_mask;
445        word1(rv) = 0;
446#endif /*Honor_FLT_ROUNDS*/
447#ifdef SET_INEXACT
448        /* set overflow bit */
449        dval(rv0) = 1e300;
450        dval(rv0) *= dval(rv0);
451#endif
452#else /*IEEE_Arith*/
453        word0(rv) = Big0;
454        word1(rv) = Big1;
455#endif /*IEEE_Arith*/
456        if (bd0)
457          goto retfree;
458        goto ret;
459        }
460      e1 = (unsigned int)e1 >> 4;
461      for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1)
462        if (e1 & 1)
463          dval(rv) *= bigtens[j];
464    /* The last multiplication could overflow. */
465      word0(rv) -= P*Exp_msk1;
466      dval(rv) *= bigtens[j];
467      if ((z = word0(rv) & Exp_mask)
468       > Exp_msk1*(DBL_MAX_EXP+Bias-P))
469        goto ovfl;
470      if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
471        /* set to largest number */
472        /* (Can't trust DBL_MAX) */
473        word0(rv) = Big0;
474        word1(rv) = Big1;
475        }
476      else
477        word0(rv) += P*Exp_msk1;
478      }
479    }
480  else if (e1 < 0) {
481    e1 = -e1;
482    if ( (i = e1 & 15) !=0)
483      dval(rv) /= tens[i];
484    if (e1 >>= 4) {
485      if (e1 >= 1 << n_bigtens)
486        goto undfl;
487#ifdef Avoid_Underflow
488      if (e1 & Scale_Bit)
489        scale = 2*P;
490      for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1)
491        if (e1 & 1)
492          dval(rv) *= tinytens[j];
493      if (scale && (j = 2*P + 1 - (unsigned int)((word0(rv) & Exp_mask)
494            >> Exp_shift)) > 0) {
495        /* scaled rv is denormal; zap j low bits */
496        if (j >= 32) {
497          word1(rv) = 0;
498          if (j >= 53)
499           word0(rv) = (P+2)*Exp_msk1;
500          else
501           word0(rv) &= 0xffffffff << (j-32);
502          }
503        else
504          word1(rv) &= 0xffffffff << j;
505        }
506#else
507      for(j = 0; e1 > 1; j++, e1 >>= 1)
508        if (e1 & 1)
509          dval(rv) *= tinytens[j];
510      /* The last multiplication could underflow. */
511      dval(rv0) = dval(rv);
512      dval(rv) *= tinytens[j];
513      if (!dval(rv)) {
514        dval(rv) = 2.*dval(rv0);
515        dval(rv) *= tinytens[j];
516#endif
517        if (!dval(rv)) {
518 undfl:
519          dval(rv) = 0.;
520#ifndef NO_ERRNO
521          errno = ERANGE;
522#endif
523          if (bd0)
524            goto retfree;
525          goto ret;
526          }
527#ifndef Avoid_Underflow
528        word0(rv) = Tiny0;
529        word1(rv) = Tiny1;
530        /* The refinement below will clean
531         * this approximation up.
532         */
533        }
534#endif
535      }
536    }
537
538  /* Now the hard part -- adjusting rv to the correct value.*/
539
540  /* Put digits into bd: true value = bd * 10^e */
541
542  bd0 = s2b(s0, nd0, nd, y);
543  if (bd0 == NULL)
544    goto ovfl;
545
546  for(;;) {
547    bd = Balloc(bd0->k);
548    if (bd == NULL)
549      goto ovfl;
550    Bcopy(bd, bd0);
551    bb = d2b(dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
552    if (bb == NULL)
553      goto ovfl;
554    bs = i2b(1);
555    if (bs == NULL)
556      goto ovfl;
557
558    if (e >= 0) {
559      bb2 = bb5 = 0;
560      bd2 = bd5 = e;
561      }
562    else {
563      bb2 = bb5 = -e;
564      bd2 = bd5 = 0;
565      }
566    if (bbe >= 0)
567      bb2 += bbe;
568    else
569      bd2 -= bbe;
570    bs2 = bb2;
571#ifdef Honor_FLT_ROUNDS
572    if (rounding != 1)
573      bs2++;
574#endif
575#ifdef Avoid_Underflow
576    j = bbe - scale;
577    i = j + bbbits - 1; /* logb(rv) */
578    if (i < Emin) /* denormal */
579      j += P - Emin;
580    else
581      j = P + 1 - bbbits;
582#else /*Avoid_Underflow*/
583#ifdef Sudden_Underflow
584#ifdef IBM
585    j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
586#else
587    j = P + 1 - bbbits;
588#endif
589#else /*Sudden_Underflow*/
590    j = bbe;
591    i = j + bbbits - 1; /* logb(rv) */
592    if (i < Emin) /* denormal */
593      j += P - Emin;
594    else
595      j = P + 1 - bbbits;
596#endif /*Sudden_Underflow*/
597#endif /*Avoid_Underflow*/
598    bb2 += j;
599    bd2 += j;
600#ifdef Avoid_Underflow
601    bd2 += scale;
602#endif
603    i = bb2 < bd2 ? bb2 : bd2;
604    if (i > bs2)
605      i = bs2;
606    if (i > 0) {
607      bb2 -= i;
608      bd2 -= i;
609      bs2 -= i;
610      }
611    if (bb5 > 0) {
612      bs = pow5mult(bs, bb5);
613      if (bs == NULL)
614        goto ovfl;
615      bb1 = mult(bs, bb);
616      if (bb1 == NULL)
617        goto ovfl;
618      Bfree(bb);
619      bb = bb1;
620      }
621    if (bb2 > 0) {
622      bb = lshift(bb, bb2);
623      if (bb == NULL)
624        goto ovfl;
625      }
626    if (bd5 > 0) {
627      bd = pow5mult(bd, bd5);
628      if (bd == NULL)
629        goto ovfl;
630      }
631    if (bd2 > 0) {
632      bd = lshift(bd, bd2);
633      if (bd == NULL)
634        goto ovfl;
635      }
636    if (bs2 > 0) {
637      bs = lshift(bs, bs2);
638      if (bs == NULL)
639        goto ovfl;
640      }
641    delta = diff(bb, bd);
642    if (delta == NULL)
643      goto ovfl;
644    dsign = delta->sign;
645    delta->sign = 0;
646    i = cmp(delta, bs);
647#ifdef Honor_FLT_ROUNDS
648    if (rounding != 1) {
649      if (i < 0) {
650        /* Error is less than an ulp */
651        if (!delta->x[0] && delta->wds <= 1) {
652          /* exact */
653#ifdef SET_INEXACT
654          inexact = 0;
655#endif
656          break;
657          }
658        if (rounding) {
659          if (dsign) {
660            adj = 1.;
661            goto apply_adj;
662            }
663          }
664        else if (!dsign) {
665          adj = -1.;
666          if (!word1(rv)
667           && !(word0(rv) & Frac_mask)) {
668            y = word0(rv) & Exp_mask;
669#ifdef Avoid_Underflow
670            if (!scale || y > 2*P*Exp_msk1)
671#else
672            if (y)
673#endif
674              {
675              delta = lshift(delta,Log2P);
676              if (cmp(delta, bs) <= 0)
677              adj = -0.5;
678              }
679            }
680 apply_adj:
681#ifdef Avoid_Underflow
682          if (scale && (y = word0(rv) & Exp_mask)
683            <= 2*P*Exp_msk1)
684            word0(adj) += (2*P+1)*Exp_msk1 - y;
685#else
686#ifdef Sudden_Underflow
687          if ((word0(rv) & Exp_mask) <=
688              P*Exp_msk1) {
689            word0(rv) += P*Exp_msk1;
690            dval(rv) += adj*ulp(dval(rv));
691            word0(rv) -= P*Exp_msk1;
692            }
693          else
694#endif /*Sudden_Underflow*/
695#endif /*Avoid_Underflow*/
696          dval(rv) += adj*ulp(dval(rv));
697          }
698        break;
699        }
700      adj = ratio(delta, bs);
701      if (adj < 1.)
702        adj = 1.;
703      if (adj <= 0x7ffffffe) {
704        /* adj = rounding ? ceil(adj) : floor(adj); */
705        y = adj;
706        if (y != adj) {
707          if (!((rounding>>1) ^ dsign))
708            y++;
709          adj = y;
710          }
711        }
712#ifdef Avoid_Underflow
713      if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
714        word0(adj) += (2*P+1)*Exp_msk1 - y;
715#else
716#ifdef Sudden_Underflow
717      if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
718        word0(rv) += P*Exp_msk1;
719        adj *= ulp(dval(rv));
720        if (dsign)
721          dval(rv) += adj;
722        else
723          dval(rv) -= adj;
724        word0(rv) -= P*Exp_msk1;
725        goto cont;
726        }
727#endif /*Sudden_Underflow*/
728#endif /*Avoid_Underflow*/
729      adj *= ulp(dval(rv));
730      if (dsign)
731        dval(rv) += adj;
732      else
733        dval(rv) -= adj;
734      goto cont;
735      }
736#endif /*Honor_FLT_ROUNDS*/
737
738    if (i < 0) {
739      /* Error is less than half an ulp -- check for
740       * special case of mantissa a power of two.
741       */
742      if (dsign || word1(rv) || word0(rv) & Bndry_mask
743#ifdef IEEE_Arith
744#ifdef Avoid_Underflow
745       || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
746#else
747       || (word0(rv) & Exp_mask) <= Exp_msk1
748#endif
749#endif
750        ) {
751#ifdef SET_INEXACT
752        if (!delta->x[0] && delta->wds <= 1)
753          inexact = 0;
754#endif
755        break;
756        }
757      if (!delta->x[0] && delta->wds <= 1) {
758        /* exact result */
759#ifdef SET_INEXACT
760        inexact = 0;
761#endif
762        break;
763        }
764      delta = lshift(delta,Log2P);
765      if (cmp(delta, bs) > 0)
766        goto drop_down;
767      break;
768      }
769    if (i == 0) {
770      /* exactly half-way between */
771      if (dsign) {
772        if ((word0(rv) & Bndry_mask1) == Bndry_mask1
773         &&  word1(rv) == (
774#ifdef Avoid_Underflow
775      (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
776    ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
777#endif
778               0xffffffff)) {
779          /*boundary case -- increment exponent*/
780          word0(rv) = (word0(rv) & Exp_mask)
781            + Exp_msk1
782#ifdef IBM
783            | Exp_msk1 >> 4
784#endif
785            ;
786          word1(rv) = 0;
787#ifdef Avoid_Underflow
788          dsign = 0;
789#endif
790          break;
791          }
792        }
793      else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
794 drop_down:
795        /* boundary case -- decrement exponent */
796#ifdef Sudden_Underflow /*{{*/
797        L = word0(rv) & Exp_mask;
798#ifdef IBM
799        if (L <  Exp_msk1)
800#else
801#ifdef Avoid_Underflow
802        if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
803#else
804        if (L <= Exp_msk1)
805#endif /*Avoid_Underflow*/
806#endif /*IBM*/
807          goto undfl;
808        L -= Exp_msk1;
809#else /*Sudden_Underflow}{*/
810#ifdef Avoid_Underflow
811        if (scale) {
812          L = word0(rv) & Exp_mask;
813          if (L <= (2*P+1)*Exp_msk1) {
814            if (L > (P+2)*Exp_msk1)
815              /* round even ==> */
816              /* accept rv */
817              break;
818            /* rv = smallest denormal */
819            goto undfl;
820            }
821          }
822#endif /*Avoid_Underflow*/
823        L = (word0(rv) & Exp_mask) - Exp_msk1;
824#endif /*Sudden_Underflow}*/
825        word0(rv) = (UINT32)(L | Bndry_mask1);
826        word1(rv) = 0xffffffffU;
827#ifdef IBM
828        goto cont;
829#else
830        break;
831#endif
832        }
833#ifndef ROUND_BIASED
834      if (!(word1(rv) & LSB))
835        break;
836#endif
837      if (dsign)
838        dval(rv) += ulp(dval(rv));
839#ifndef ROUND_BIASED
840      else {
841        dval(rv) -= ulp(dval(rv));
842#ifndef Sudden_Underflow
843        if (!dval(rv))
844          goto undfl;
845#endif
846        }
847#ifdef Avoid_Underflow
848      dsign = 1 - dsign;
849#endif
850#endif
851      break;
852      }
853    if ((aadj = ratio(delta, bs)) <= 2.) {
854      if (dsign)
855        aadj = aadj1 = 1.;
856      else if (word1(rv) || word0(rv) & Bndry_mask) {
857#ifndef Sudden_Underflow
858        if (word1(rv) == Tiny1 && !word0(rv))
859          goto undfl;
860#endif
861        aadj = 1.;
862        aadj1 = -1.;
863        }
864      else {
865        /* special case -- power of FLT_RADIX to be */
866        /* rounded down... */
867
868        if (aadj < 2./FLT_RADIX)
869          aadj = 1./FLT_RADIX;
870        else
871          aadj *= 0.5;
872        aadj1 = -aadj;
873        }
874      }
875    else {
876      aadj *= 0.5;
877      aadj1 = dsign ? aadj : -aadj;
878#ifdef Check_FLT_ROUNDS
879      switch(Rounding) {
880        case 2: /* towards +infinity */
881          aadj1 -= 0.5;
882          break;
883        case 0: /* towards 0 */
884        case 3: /* towards -infinity */
885          aadj1 += 0.5;
886        }
887#else
888      if (Flt_Rounds == 0)
889        aadj1 += 0.5;
890#endif /*Check_FLT_ROUNDS*/
891      }
892    y = word0(rv) & Exp_mask;
893
894    /* Check for overflow */
895
896    if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
897      dval(rv0) = dval(rv);
898      word0(rv) -= P*Exp_msk1;
899      adj = aadj1 * ulp(dval(rv));
900      dval(rv) += adj;
901      if ((word0(rv) & Exp_mask) >=
902          Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
903        if (word0(rv0) == Big0 && word1(rv0) == Big1)
904          goto ovfl;
905        word0(rv) = Big0;
906        word1(rv) = Big1;
907        goto cont;
908        }
909      else
910        word0(rv) += P*Exp_msk1;
911      }
912    else {
913#ifdef Avoid_Underflow
914      if (scale && y <= 2*P*Exp_msk1) {
915        if (aadj <= 0x7fffffff) {
916          if ((z = (uint32_t)aadj) == 0)
917            z = 1;
918          aadj = (double)z;
919          aadj1 = dsign ? aadj : -aadj;
920          }
921        word0(aadj1) += (UINT32)((2*P+1)*Exp_msk1 - y);
922        }
923      adj = aadj1 * ulp(dval(rv));
924      dval(rv) += adj;
925#else
926#ifdef Sudden_Underflow
927      if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
928        dval(rv0) = dval(rv);
929        word0(rv) += P*Exp_msk1;
930        adj = aadj1 * ulp(dval(rv));
931        dval(rv) += adj;
932#ifdef IBM
933        if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
934#else
935        if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
936#endif
937          {
938          if (word0(rv0) == Tiny0
939           && word1(rv0) == Tiny1)
940            goto undfl;
941          word0(rv) = Tiny0;
942          word1(rv) = Tiny1;
943          goto cont;
944          }
945        else
946          word0(rv) -= P*Exp_msk1;
947        }
948      else {
949        adj = aadj1 * ulp(dval(rv));
950        dval(rv) += adj;
951        }
952#else /*Sudden_Underflow*/
953      /* Compute adj so that the IEEE rounding rules will
954       * correctly round rv + adj in some half-way cases.
955       * If rv * ulp(rv) is denormalized (i.e.,
956       * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
957       * trouble from bits lost to denormalization;
958       * example: 1.2e-307 .
959       */
960      if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
961        aadj1 = (double)(int)(aadj + 0.5);
962        if (!dsign)
963          aadj1 = -aadj1;
964        }
965      adj = aadj1 * ulp(dval(rv));
966      dval(rv) += adj;
967#endif /*Sudden_Underflow*/
968#endif /*Avoid_Underflow*/
969      }
970    z = word0(rv) & Exp_mask;
971#ifndef SET_INEXACT
972#ifdef Avoid_Underflow
973    if (!scale)
974#endif
975    if (y == z) {
976      /* Can we stop now? */
977      L = (Long)aadj;
978      aadj -= L;
979      /* The tolerances below are conservative. */
980      if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
981        if (aadj < .4999999 || aadj > .5000001)
982          break;
983        }
984      else if (aadj < .4999999/FLT_RADIX)
985        break;
986      }
987#endif
988 cont:
989    Bfree(bb);
990    Bfree(bd);
991    Bfree(bs);
992    Bfree(delta);
993    }
994#ifdef SET_INEXACT
995  if (inexact) {
996    if (!oldinexact) {
997      word0(rv0) = Exp_1 + (70 << Exp_shift);
998      word1(rv0) = 0;
999      dval(rv0) += 1.;
1000      }
1001    }
1002  else if (!oldinexact)
1003    clear_inexact();
1004#endif
1005#ifdef Avoid_Underflow
1006  if (scale) {
1007    word0(rv0) = Exp_1 - 2*P*Exp_msk1;
1008    word1(rv0) = 0;
1009    dval(rv) *= dval(rv0);
1010#ifndef NO_ERRNO
1011    /* try to avoid the bug of testing an 8087 register value */
1012    if (word0(rv) == 0 && word1(rv) == 0)
1013      errno = ERANGE;
1014#endif
1015    }
1016#endif /* Avoid_Underflow */
1017#ifdef SET_INEXACT
1018  if (inexact && !(word0(rv) & Exp_mask)) {
1019    /* set underflow bit */
1020    dval(rv0) = 1e-300;
1021    dval(rv0) *= dval(rv0);
1022    }
1023#endif
1024 retfree:
1025  Bfree(bb);
1026  Bfree(bd);
1027  Bfree(bs);
1028  Bfree(bd0);
1029  Bfree(delta);
1030 ret:
1031  if (se)
1032    *se = __UNCONST(s);
1033  return sign ? -dval(rv) : dval(rv);
1034}
1035
1036