source: trunk/libs/newlib/src/newlib/libc/stdlib/strtod.c @ 567

Last change on this file since 567 was 444, checked in by satin@…, 6 years ago

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