]>
git.proxmox.com Git - mirror_edk2.git/blob - StdLib/LibC/gdtoa/strtod.c
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.
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.
12 *****************************************************************
14 The author of this software is David M. Gay.
16 Copyright (C) 1998-2001 by Lucent Technologies
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
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
39 Please send bug reports to David M. Gay (dmg at acm dot org,
40 with " at " changed at "@" and " dot " changed to ".").
42 *****************************************************************
44 NetBSD: strtod.c,v 1.4.14.1 2008/04/08 21:10:55 jdc Exp
46 #include <LibConfig.h>
59 #define Avoid_Underflow
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. */
63 static CONST
double tinytens
[] = { 1e-16, 1e-32, 1e-64, 1e-128,
64 9007199254740992.e
-256
69 #ifdef Honor_FLT_ROUNDS
70 #define Rounding rounding
71 #undef Check_FLT_ROUNDS
72 #define Check_FLT_ROUNDS
74 #define Rounding Flt_Rounds
77 //#ifndef __HAVE_LONG_DOUBLE
78 //__strong_alias(_strtold, strtod)
79 //__weak_alias(strtold, _strtold)
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) */
88 strtod(CONST
char *s00
, char **se
)
90 #ifdef Avoid_Underflow
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
;
99 Bigint
*bb
= NULL
, *bb1
, *bd0
;
100 Bigint
*bd
= NULL
, *bs
= NULL
, *delta
= NULL
; /* pacify gcc */
102 int inexact
, oldinexact
;
104 #ifdef Honor_FLT_ROUNDS
108 sign
= nz0
= nz
= decpt
= 0;
136 static FPI fpi
= { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI
};
143 #if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD)
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;
152 switch((i
= gethex(&s
, &fpi
, &expt
, &bb
, sign
)) & STRTOG_Retmask
) {
153 case STRTOG_NoNumber
:
161 copybits(bits
, fpi
.nbits
, bb
);
164 ULtod((/* LINTED */(U
*)&rv
)->L
, bits
, expt
, i
);
177 for(nd
= nf
= 0; (c
= *s
) >= '0' && c
<= '9'; nd
++, s
++)
184 if (c
== *localeconv()->decimal_point
)
192 for(; c
== '0'; c
= *++s
)
194 if (c
> '0' && c
<= '9') {
202 for(; c
>= '0' && c
<= '9'; c
= *++s
) {
207 for(i
= 1; i
< nz
; i
++)
210 else if (nd
<= DBL_DIG
+ 1)
214 else if (nd
<= DBL_DIG
+ 1)
222 if (c
== 'e' || c
== 'E') {
223 if (!nd
&& !nz
&& !nz0
) {
235 if (c
>= '0' && c
<= '9') {
238 if (c
> '0' && c
<= '9') {
241 while((c
= *++s
) >= '0' && c
<= '9')
243 if (s
- s1
> 8 || L
> 19999)
244 /* Avoid confusion from exponents
245 * so large that e might overflow.
247 e
= 19999; /* safe for 16 bit ints */
262 /* Check for Nan and Infinity */
265 static FPI fpinan
= /* only 52 explicit bits */
266 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI
};
272 if (match(&s
,"nf")) {
274 if (!match(&s
,"inity"))
276 word0(rv
) = 0x7ff00000;
283 if (match(&s
, "an")) {
286 && hexnan(&s
, &fpinan
, bits
)
288 word0(rv
) = (UINT32
)(0x7ff00000U
| bits
[1]);
289 word1(rv
) = (UINT32
)bits
[0];
293 word0(rv
) = NAN_WORD0
;
294 word1(rv
) = NAN_WORD1
;
301 #endif /* INFNAN_CHECK */
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
317 k
= nd
< DBL_DIG
+ 1 ? nd
: DBL_DIG
+ 1;
318 dval(rv
) = (double)y
;
322 oldinexact
= get_inexact();
324 dval(rv
) = tens
[k
- 9] * dval(rv
) + z
;
329 #ifndef Honor_FLT_ROUNDS
341 #ifdef Honor_FLT_ROUNDS
342 /* round correctly FLT_ROUNDS = 2 or 3 */
348 /* rv = */ rounded_product(dval(rv
), tens
[e
]);
353 if (e
<= Ten_pmax
+ i
) {
354 /* A fancier test would sometimes let us do
355 * this for larger i values.
357 #ifdef Honor_FLT_ROUNDS
358 /* round correctly FLT_ROUNDS = 2 or 3 */
367 /* VAX exponent range is so narrow we must
368 * worry about overflow here...
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
))
376 word0(rv
) += P
*Exp_msk1
;
378 /* rv = */ rounded_product(dval(rv
), tens
[e
]);
383 #ifndef Inaccurate_Divide
384 else if (e
>= -Ten_pmax
) {
385 #ifdef Honor_FLT_ROUNDS
386 /* round correctly FLT_ROUNDS = 2 or 3 */
392 /* rv = */ rounded_quotient(dval(rv
), tens
[-e
]);
403 oldinexact
= get_inexact();
405 #ifdef Avoid_Underflow
408 #ifdef Honor_FLT_ROUNDS
409 if ((rounding
= Flt_Rounds
) >= 2) {
411 rounding
= rounding
== 2 ? 0 : 2;
417 #endif /*IEEE_Arith*/
419 /* Get starting approximation = rv * 10**e1 */
422 if ( (i
= e1
& 15) !=0)
425 if (e1
> DBL_MAX_10_EXP
) {
430 /* Can't trust HUGE_VAL */
432 #ifdef Honor_FLT_ROUNDS
434 case 0: /* toward 0 */
435 case 3: /* toward -infinity */
440 word0(rv
) = Exp_mask
;
443 #else /*Honor_FLT_ROUNDS*/
444 word0(rv
) = Exp_mask
;
446 #endif /*Honor_FLT_ROUNDS*/
448 /* set overflow bit */
450 dval(rv0
) *= dval(rv0
);
455 #endif /*IEEE_Arith*/
460 e1
= (unsigned int)e1
>> 4;
461 for(j
= 0; e1
> 1; j
++, e1
= (unsigned int)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
))
470 if (z
> Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
)) {
471 /* set to largest number */
472 /* (Can't trust DBL_MAX) */
477 word0(rv
) += P
*Exp_msk1
;
482 if ( (i
= e1
& 15) !=0)
485 if (e1
>= 1 << n_bigtens
)
487 #ifdef Avoid_Underflow
490 for(j
= 0; e1
> 0; j
++, e1
= (unsigned int)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 */
499 word0(rv
) = (P
+2)*Exp_msk1
;
501 word0(rv
) &= 0xffffffff << (j
-32);
504 word1(rv
) &= 0xffffffff << j
;
507 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
509 dval(rv
) *= tinytens
[j
];
510 /* The last multiplication could underflow. */
511 dval(rv0
) = dval(rv
);
512 dval(rv
) *= tinytens
[j
];
514 dval(rv
) = 2.*dval(rv0
);
515 dval(rv
) *= tinytens
[j
];
527 #ifndef Avoid_Underflow
530 /* The refinement below will clean
531 * this approximation up.
538 /* Now the hard part -- adjusting rv to the correct value.*/
540 /* Put digits into bd: true value = bd * 10^e */
542 bd0
= s2b(s0
, nd0
, nd
, y
);
551 bb
= d2b(dval(rv
), &bbe
, &bbbits
); /* rv = bb * 2^bbe */
571 #ifdef Honor_FLT_ROUNDS
575 #ifdef Avoid_Underflow
577 i
= j
+ bbbits
- 1; /* logb(rv) */
578 if (i
< Emin
) /* denormal */
582 #else /*Avoid_Underflow*/
583 #ifdef Sudden_Underflow
585 j
= 1 + 4*P
- 3 - bbbits
+ ((bbe
+ bbbits
- 1) & 3);
589 #else /*Sudden_Underflow*/
591 i
= j
+ bbbits
- 1; /* logb(rv) */
592 if (i
< Emin
) /* denormal */
596 #endif /*Sudden_Underflow*/
597 #endif /*Avoid_Underflow*/
600 #ifdef Avoid_Underflow
603 i
= bb2
< bd2
? bb2
: bd2
;
612 bs
= pow5mult(bs
, bb5
);
622 bb
= lshift(bb
, bb2
);
627 bd
= pow5mult(bd
, bd5
);
632 bd
= lshift(bd
, bd2
);
637 bs
= lshift(bs
, bs2
);
641 delta
= diff(bb
, bd
);
647 #ifdef Honor_FLT_ROUNDS
650 /* Error is less than an ulp */
651 if (!delta
->x
[0] && delta
->wds
<= 1) {
667 && !(word0(rv
) & Frac_mask
)) {
668 y
= word0(rv
) & Exp_mask
;
669 #ifdef Avoid_Underflow
670 if (!scale
|| y
> 2*P
*Exp_msk1
)
675 delta
= lshift(delta
,Log2P
);
676 if (cmp(delta
, bs
) <= 0)
681 #ifdef Avoid_Underflow
682 if (scale
&& (y
= word0(rv
) & Exp_mask
)
684 word0(adj
) += (2*P
+1)*Exp_msk1
- y
;
686 #ifdef Sudden_Underflow
687 if ((word0(rv
) & Exp_mask
) <=
689 word0(rv
) += P
*Exp_msk1
;
690 dval(rv
) += adj
*ulp(dval(rv
));
691 word0(rv
) -= P
*Exp_msk1
;
694 #endif /*Sudden_Underflow*/
695 #endif /*Avoid_Underflow*/
696 dval(rv
) += adj
*ulp(dval(rv
));
700 adj
= ratio(delta
, bs
);
703 if (adj
<= 0x7ffffffe) {
704 /* adj = rounding ? ceil(adj) : floor(adj); */
707 if (!((rounding
>>1) ^ dsign
))
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
;
716 #ifdef Sudden_Underflow
717 if ((word0(rv
) & Exp_mask
) <= P
*Exp_msk1
) {
718 word0(rv
) += P
*Exp_msk1
;
719 adj
*= ulp(dval(rv
));
724 word0(rv
) -= P
*Exp_msk1
;
727 #endif /*Sudden_Underflow*/
728 #endif /*Avoid_Underflow*/
729 adj
*= ulp(dval(rv
));
736 #endif /*Honor_FLT_ROUNDS*/
739 /* Error is less than half an ulp -- check for
740 * special case of mantissa a power of two.
742 if (dsign
|| word1(rv
) || word0(rv
) & Bndry_mask
744 #ifdef Avoid_Underflow
745 || (word0(rv
) & Exp_mask
) <= (2*P
+1)*Exp_msk1
747 || (word0(rv
) & Exp_mask
) <= Exp_msk1
752 if (!delta
->x
[0] && delta
->wds
<= 1)
757 if (!delta
->x
[0] && delta
->wds
<= 1) {
764 delta
= lshift(delta
,Log2P
);
765 if (cmp(delta
, bs
) > 0)
770 /* exactly half-way between */
772 if ((word0(rv
) & Bndry_mask1
) == Bndry_mask1
774 #ifdef Avoid_Underflow
775 (scale
&& (y
= word0(rv
) & Exp_mask
) <= 2*P
*Exp_msk1
)
776 ? (0xffffffff & (0xffffffff << (2*P
+1-(y
>>Exp_shift
)))) :
779 /*boundary case -- increment exponent*/
780 word0(rv
) = (word0(rv
) & Exp_mask
)
787 #ifdef Avoid_Underflow
793 else if (!(word0(rv
) & Bndry_mask
) && !word1(rv
)) {
795 /* boundary case -- decrement exponent */
796 #ifdef Sudden_Underflow /*{{*/
797 L
= word0(rv
) & Exp_mask
;
801 #ifdef Avoid_Underflow
802 if (L
<= (scale
? (2*P
+1)*Exp_msk1
: Exp_msk1
))
805 #endif /*Avoid_Underflow*/
809 #else /*Sudden_Underflow}{*/
810 #ifdef Avoid_Underflow
812 L
= word0(rv
) & Exp_mask
;
813 if (L
<= (2*P
+1)*Exp_msk1
) {
814 if (L
> (P
+2)*Exp_msk1
)
818 /* rv = smallest denormal */
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
;
834 if (!(word1(rv
) & LSB
))
838 dval(rv
) += ulp(dval(rv
));
841 dval(rv
) -= ulp(dval(rv
));
842 #ifndef Sudden_Underflow
847 #ifdef Avoid_Underflow
853 if ((aadj
= ratio(delta
, bs
)) <= 2.) {
856 else if (word1(rv
) || word0(rv
) & Bndry_mask
) {
857 #ifndef Sudden_Underflow
858 if (word1(rv
) == Tiny1
&& !word0(rv
))
865 /* special case -- power of FLT_RADIX to be */
866 /* rounded down... */
868 if (aadj
< 2./FLT_RADIX
)
877 aadj1
= dsign
? aadj
: -aadj
;
878 #ifdef Check_FLT_ROUNDS
880 case 2: /* towards +infinity */
883 case 0: /* towards 0 */
884 case 3: /* towards -infinity */
890 #endif /*Check_FLT_ROUNDS*/
892 y
= word0(rv
) & Exp_mask
;
894 /* Check for overflow */
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
));
901 if ((word0(rv
) & Exp_mask
) >=
902 Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
)) {
903 if (word0(rv0
) == Big0
&& word1(rv0
) == Big1
)
910 word0(rv
) += P
*Exp_msk1
;
913 #ifdef Avoid_Underflow
914 if (scale
&& y
<= 2*P
*Exp_msk1
) {
915 if (aadj
<= 0x7fffffff) {
916 if ((z
= (uint32_t)aadj
) == 0)
919 aadj1
= dsign
? aadj
: -aadj
;
921 word0(aadj1
) += (UINT32
)((2*P
+1)*Exp_msk1
- y
);
923 adj
= aadj1
* ulp(dval(rv
));
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
));
933 if ((word0(rv
) & Exp_mask
) < P
*Exp_msk1
)
935 if ((word0(rv
) & Exp_mask
) <= P
*Exp_msk1
)
938 if (word0(rv0
) == Tiny0
939 && word1(rv0
) == Tiny1
)
946 word0(rv
) -= P
*Exp_msk1
;
949 adj
= aadj1
* ulp(dval(rv
));
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 .
960 if (y
<= (P
-1)*Exp_msk1
&& aadj
> 1.) {
961 aadj1
= (double)(int)(aadj
+ 0.5);
965 adj
= aadj1
* ulp(dval(rv
));
967 #endif /*Sudden_Underflow*/
968 #endif /*Avoid_Underflow*/
970 z
= word0(rv
) & Exp_mask
;
972 #ifdef Avoid_Underflow
976 /* Can we stop now? */
979 /* The tolerances below are conservative. */
980 if (dsign
|| word1(rv
) || word0(rv
) & Bndry_mask
) {
981 if (aadj
< .4999999 || aadj
> .5000001)
984 else if (aadj
< .4999999/FLT_RADIX
)
997 word0(rv0
) = Exp_1
+ (70 << Exp_shift
);
1002 else if (!oldinexact
)
1005 #ifdef Avoid_Underflow
1007 word0(rv0
) = Exp_1
- 2*P
*Exp_msk1
;
1009 dval(rv
) *= dval(rv0
);
1011 /* try to avoid the bug of testing an 8087 register value */
1012 if (word0(rv
) == 0 && word1(rv
) == 0)
1016 #endif /* Avoid_Underflow */
1018 if (inexact
&& !(word0(rv
) & Exp_mask
)) {
1019 /* set underflow bit */
1021 dval(rv0
) *= dval(rv0
);
1033 return sign
? -dval(rv
) : dval(rv
);