]>
Commit | Line | Data |
---|---|---|
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 | |
63 | static 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 | |
87 | double\r | |
88 | strtod(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 |