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