Details | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
6099 | serge | 1 | /**************************************************************** |
2 | |||
3 | The author of this software is David M. Gay. |
||
4 | |||
5 | Copyright (C) 1998-2001 by Lucent Technologies |
||
6 | All Rights Reserved |
||
7 | |||
8 | Permission to use, copy, modify, and distribute this software and |
||
9 | its documentation for any purpose and without fee is hereby |
||
10 | granted, provided that the above copyright notice appear in all |
||
11 | copies and that both that the copyright notice and this |
||
12 | permission notice and warranty disclaimer appear in supporting |
||
13 | documentation, and that the name of Lucent or any of its entities |
||
14 | not be used in advertising or publicity pertaining to |
||
15 | distribution of the software without specific, written prior |
||
16 | permission. |
||
17 | |||
18 | LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, |
||
19 | INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. |
||
20 | IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY |
||
21 | SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES |
||
22 | WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER |
||
23 | IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, |
||
24 | ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF |
||
25 | THIS SOFTWARE. |
||
26 | |||
27 | ****************************************************************/ |
||
28 | |||
29 | /* Please send bug reports to David M. Gay (dmg at acm dot org, |
||
30 | * with " at " changed at "@" and " dot " changed to "."). */ |
||
31 | |||
32 | #include <_ansi.h> |
||
33 | #include |
||
34 | #include |
||
35 | #include |
||
36 | #include "mprec.h" |
||
37 | #include "gdtoa.h" |
||
38 | #include "gd_qnan.h" |
||
39 | |||
40 | #include "locale.h" |
||
41 | |||
42 | #if defined (_HAVE_LONG_DOUBLE) && !defined (_LDBL_EQ_DBL) |
||
43 | |||
44 | #define USE_LOCALE |
||
45 | |||
46 | static const int |
||
47 | fivesbits[] = { 0, 3, 5, 7, 10, 12, 14, 17, 19, 21, |
||
48 | 24, 26, 28, 31, 33, 35, 38, 40, 42, 45, |
||
49 | 47, 49, 52 |
||
50 | #ifdef VAX |
||
51 | , 54, 56 |
||
52 | #endif |
||
53 | }; |
||
54 | |||
55 | static _Bigint * |
||
56 | #ifdef KR_headers |
||
57 | sum(p, a, b) struct _reent *p; _Bigint *a; _Bigint *b; |
||
58 | #else |
||
59 | sum(struct _reent *p, _Bigint *a, _Bigint *b) |
||
60 | #endif |
||
61 | { |
||
62 | _Bigint *c; |
||
63 | __ULong carry, *xc, *xa, *xb, *xe, y; |
||
64 | #ifdef Pack_32 |
||
65 | __ULong z; |
||
66 | #endif |
||
67 | |||
68 | if (a->_wds < b->_wds) { |
||
69 | c = b; b = a; a = c; |
||
70 | } |
||
71 | c = Balloc(p, a->_k); |
||
72 | c->_wds = a->_wds; |
||
73 | carry = 0; |
||
74 | xa = a->_x; |
||
75 | xb = b->_x; |
||
76 | xc = c->_x; |
||
77 | xe = xc + b->_wds; |
||
78 | #ifdef Pack_32 |
||
79 | do { |
||
80 | y = (*xa & 0xffff) + (*xb & 0xffff) + carry; |
||
81 | carry = (y & 0x10000) >> 16; |
||
82 | z = (*xa++ >> 16) + (*xb++ >> 16) + carry; |
||
83 | carry = (z & 0x10000) >> 16; |
||
84 | Storeinc(xc, z, y); |
||
85 | } |
||
86 | while(xc < xe); |
||
87 | xe += a->_wds - b->_wds; |
||
88 | while(xc < xe) { |
||
89 | y = (*xa & 0xffff) + carry; |
||
90 | carry = (y & 0x10000) >> 16; |
||
91 | z = (*xa++ >> 16) + carry; |
||
92 | carry = (z & 0x10000) >> 16; |
||
93 | Storeinc(xc, z, y); |
||
94 | } |
||
95 | #else |
||
96 | do { |
||
97 | y = *xa++ + *xb++ + carry; |
||
98 | carry = (y & 0x10000) >> 16; |
||
99 | *xc++ = y & 0xffff; |
||
100 | } |
||
101 | while(xc < xe); |
||
102 | xe += a->_wds - b->_wds; |
||
103 | while(xc < xe) { |
||
104 | y = *xa++ + carry; |
||
105 | carry = (y & 0x10000) >> 16; |
||
106 | *xc++ = y & 0xffff; |
||
107 | } |
||
108 | #endif |
||
109 | if (carry) { |
||
110 | if (c->_wds == c->_maxwds) { |
||
111 | b = Balloc(p, c->_k + 1); |
||
112 | Bcopy(b, c); |
||
113 | Bfree(p, c); |
||
114 | c = b; |
||
115 | } |
||
116 | c->_x[c->_wds++] = 1; |
||
117 | } |
||
118 | return c; |
||
119 | } |
||
120 | |||
121 | static void |
||
122 | #ifdef KR_headers |
||
123 | rshift(b, k) _Bigint *b; int k; |
||
124 | #else |
||
125 | rshift(_Bigint *b, int k) |
||
126 | #endif |
||
127 | { |
||
128 | __ULong *x, *x1, *xe, y; |
||
129 | int n; |
||
130 | |||
131 | x = x1 = b->_x; |
||
132 | n = k >> kshift; |
||
133 | if (n < b->_wds) { |
||
134 | xe = x + b->_wds; |
||
135 | x += n; |
||
136 | if (k &= kmask) { |
||
137 | n = ULbits - k; |
||
138 | y = *x++ >> k; |
||
139 | while(x < xe) { |
||
140 | *x1++ = (y | (*x << n)) & ALL_ON; |
||
141 | y = *x++ >> k; |
||
142 | } |
||
143 | if ((*x1 = y) !=0) |
||
144 | x1++; |
||
145 | } |
||
146 | else |
||
147 | while(x < xe) |
||
148 | *x1++ = *x++; |
||
149 | } |
||
150 | if ((b->_wds = x1 - b->_x) == 0) |
||
151 | b->_x[0] = 0; |
||
152 | } |
||
153 | |||
154 | static int |
||
155 | #ifdef KR_headers |
||
156 | trailz(b) _Bigint *b; |
||
157 | #else |
||
158 | trailz(_Bigint *b) |
||
159 | #endif |
||
160 | { |
||
161 | __ULong L, *x, *xe; |
||
162 | int n = 0; |
||
163 | |||
164 | x = b->_x; |
||
165 | xe = x + b->_wds; |
||
166 | for(n = 0; x < xe && !*x; x++) |
||
167 | n += ULbits; |
||
168 | if (x < xe) { |
||
169 | L = *x; |
||
170 | n += lo0bits(&L); |
||
171 | } |
||
172 | return n; |
||
173 | } |
||
174 | |||
175 | _Bigint * |
||
176 | #ifdef KR_headers |
||
177 | increment(p, b) struct _reent *p; _Bigint *b; |
||
178 | #else |
||
179 | increment(struct _reent *p, _Bigint *b) |
||
180 | #endif |
||
181 | { |
||
182 | __ULong *x, *xe; |
||
183 | _Bigint *b1; |
||
184 | #ifdef Pack_16 |
||
185 | __ULong carry = 1, y; |
||
186 | #endif |
||
187 | |||
188 | x = b->_x; |
||
189 | xe = x + b->_wds; |
||
190 | #ifdef Pack_32 |
||
191 | do { |
||
192 | if (*x < (__ULong)0xffffffffL) { |
||
193 | ++*x; |
||
194 | return b; |
||
195 | } |
||
196 | *x++ = 0; |
||
197 | } while(x < xe); |
||
198 | #else |
||
199 | do { |
||
200 | y = *x + carry; |
||
201 | carry = y >> 16; |
||
202 | *x++ = y & 0xffff; |
||
203 | if (!carry) |
||
204 | return b; |
||
205 | } while(x < xe); |
||
206 | if (carry) |
||
207 | #endif |
||
208 | { |
||
209 | if (b->_wds >= b->_maxwds) { |
||
210 | b1 = Balloc(p,b->_k+1); |
||
211 | Bcopy(b1,b); |
||
212 | Bfree(p,b); |
||
213 | b = b1; |
||
214 | } |
||
215 | b->_x[b->_wds++] = 1; |
||
216 | } |
||
217 | return b; |
||
218 | } |
||
219 | |||
220 | int |
||
221 | #ifdef KR_headers |
||
222 | decrement(b) _Bigint *b; |
||
223 | #else |
||
224 | decrement(_Bigint *b) |
||
225 | #endif |
||
226 | { |
||
227 | __ULong *x, *xe; |
||
228 | #ifdef Pack_16 |
||
229 | __ULong borrow = 1, y; |
||
230 | #endif |
||
231 | |||
232 | x = b->_x; |
||
233 | xe = x + b->_wds; |
||
234 | #ifdef Pack_32 |
||
235 | do { |
||
236 | if (*x) { |
||
237 | --*x; |
||
238 | break; |
||
239 | } |
||
240 | *x++ = 0xffffffffL; |
||
241 | } |
||
242 | while(x < xe); |
||
243 | #else |
||
244 | do { |
||
245 | y = *x - borrow; |
||
246 | borrow = (y & 0x10000) >> 16; |
||
247 | *x++ = y & 0xffff; |
||
248 | } while(borrow && x < xe); |
||
249 | #endif |
||
250 | return STRTOG_Inexlo; |
||
251 | } |
||
252 | |||
253 | static int |
||
254 | #ifdef KR_headers |
||
255 | all_on(b, n) _Bigint *b; int n; |
||
256 | #else |
||
257 | all_on(_Bigint *b, int n) |
||
258 | #endif |
||
259 | { |
||
260 | __ULong *x, *xe; |
||
261 | |||
262 | x = b->_x; |
||
263 | xe = x + (n >> kshift); |
||
264 | while(x < xe) |
||
265 | if ((*x++ & ALL_ON) != ALL_ON) |
||
266 | return 0; |
||
267 | if (n &= kmask) |
||
268 | return ((*x | (ALL_ON << n)) & ALL_ON) == ALL_ON; |
||
269 | return 1; |
||
270 | } |
||
271 | |||
272 | _Bigint * |
||
273 | #ifdef KR_headers |
||
274 | set_ones(p, b, n) struct _reent *p; _Bigint *b; int n; |
||
275 | #else |
||
276 | set_ones(struct _reent *p, _Bigint *b, int n) |
||
277 | #endif |
||
278 | { |
||
279 | int k; |
||
280 | __ULong *x, *xe; |
||
281 | |||
282 | k = (n + ((1 << kshift) - 1)) >> kshift; |
||
283 | if (b->_k < k) { |
||
284 | Bfree(p,b); |
||
285 | b = Balloc(p,k); |
||
286 | } |
||
287 | k = n >> kshift; |
||
288 | if (n &= kmask) |
||
289 | k++; |
||
290 | b->_wds = k; |
||
291 | x = b->_x; |
||
292 | xe = x + k; |
||
293 | while(x < xe) |
||
294 | *x++ = ALL_ON; |
||
295 | if (n) |
||
296 | x[-1] >>= ULbits - n; |
||
297 | return b; |
||
298 | } |
||
299 | |||
300 | static int |
||
301 | rvOK |
||
302 | #ifdef KR_headers |
||
303 | (p, d, fpi, exp, bits, exact, rd, irv) |
||
304 | struct _reent *p; double d; FPI *fpi; Long *exp; __ULong *bits; int exact, rd, *irv; |
||
305 | #else |
||
306 | (struct _reent *p, double d, FPI *fpi, Long *exp, __ULong *bits, int exact, int rd, int *irv) |
||
307 | #endif |
||
308 | { |
||
309 | _Bigint *b; |
||
310 | __ULong carry, inex, lostbits; |
||
311 | int bdif, e, j, k, k1, nb, rv; |
||
312 | |||
313 | carry = rv = 0; |
||
314 | b = d2b(p, d, &e, &bdif); |
||
315 | bdif -= nb = fpi->nbits; |
||
316 | e += bdif; |
||
317 | if (bdif <= 0) { |
||
318 | if (exact) |
||
319 | goto trunc; |
||
320 | goto ret; |
||
321 | } |
||
322 | if (P == nb) { |
||
323 | if ( |
||
324 | #ifndef IMPRECISE_INEXACT |
||
325 | exact && |
||
326 | #endif |
||
327 | fpi->rounding == |
||
328 | #ifdef RND_PRODQUOT |
||
329 | FPI_Round_near |
||
330 | #else |
||
331 | Flt_Rounds |
||
332 | #endif |
||
333 | ) goto trunc; |
||
334 | goto ret; |
||
335 | } |
||
336 | switch(rd) { |
||
337 | case 1: |
||
338 | goto trunc; |
||
339 | case 2: |
||
340 | break; |
||
341 | default: /* round near */ |
||
342 | k = bdif - 1; |
||
343 | if (k < 0) |
||
344 | goto trunc; |
||
345 | if (!k) { |
||
346 | if (!exact) |
||
347 | goto ret; |
||
348 | if (b->_x[0] & 2) |
||
349 | break; |
||
350 | goto trunc; |
||
351 | } |
||
352 | if (b->_x[k>>kshift] & ((__ULong)1 << (k & kmask))) |
||
353 | break; |
||
354 | goto trunc; |
||
355 | } |
||
356 | /* "break" cases: round up 1 bit, then truncate; bdif > 0 */ |
||
357 | carry = 1; |
||
358 | trunc: |
||
359 | inex = lostbits = 0; |
||
360 | if (bdif > 0) { |
||
361 | if ( (lostbits = any_on(b, bdif)) !=0) |
||
362 | inex = STRTOG_Inexlo; |
||
363 | rshift(b, bdif); |
||
364 | if (carry) { |
||
365 | inex = STRTOG_Inexhi; |
||
366 | b = increment(p, b); |
||
367 | if ( (j = nb & kmask) !=0) |
||
368 | j = ULbits - j; |
||
369 | if (hi0bits(b->_x[b->_wds - 1]) != j) { |
||
370 | if (!lostbits) |
||
371 | lostbits = b->_x[0] & 1; |
||
372 | rshift(b, 1); |
||
373 | e++; |
||
374 | } |
||
375 | } |
||
376 | } |
||
377 | else if (bdif < 0) |
||
378 | b = lshift(p, b, -bdif); |
||
379 | if (e < fpi->emin) { |
||
380 | k = fpi->emin - e; |
||
381 | e = fpi->emin; |
||
382 | if (k > nb || fpi->sudden_underflow) { |
||
383 | b->_wds = inex = 0; |
||
384 | *irv = STRTOG_Underflow | STRTOG_Inexlo; |
||
385 | } |
||
386 | else { |
||
387 | k1 = k - 1; |
||
388 | if (k1 > 0 && !lostbits) |
||
389 | lostbits = any_on(b, k1); |
||
390 | if (!lostbits && !exact) |
||
391 | goto ret; |
||
392 | lostbits |= |
||
393 | carry = b->_x[k1>>kshift] & (1 << (k1 & kmask)); |
||
394 | rshift(b, k); |
||
395 | *irv = STRTOG_Denormal; |
||
396 | if (carry) { |
||
397 | b = increment(p, b); |
||
398 | inex = STRTOG_Inexhi | STRTOG_Underflow; |
||
399 | } |
||
400 | else if (lostbits) |
||
401 | inex = STRTOG_Inexlo | STRTOG_Underflow; |
||
402 | } |
||
403 | } |
||
404 | else if (e > fpi->emax) { |
||
405 | e = fpi->emax + 1; |
||
406 | *irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; |
||
407 | #ifndef NO_ERRNO |
||
408 | errno = ERANGE; |
||
409 | #endif |
||
410 | b->_wds = inex = 0; |
||
411 | } |
||
412 | *exp = e; |
||
413 | copybits(bits, nb, b); |
||
414 | *irv |= inex; |
||
415 | rv = 1; |
||
416 | ret: |
||
417 | Bfree(p,b); |
||
418 | return rv; |
||
419 | } |
||
420 | |||
421 | static int |
||
422 | #ifdef KR_headers |
||
423 | mantbits(d) double d; |
||
424 | #else |
||
425 | mantbits(U d) |
||
426 | #endif |
||
427 | { |
||
428 | __ULong L; |
||
429 | #ifdef VAX |
||
430 | L = word1(d) << 16 | word1(d) >> 16; |
||
431 | if (L) |
||
432 | #else |
||
433 | if ( (L = word1(d)) !=0) |
||
434 | #endif |
||
435 | return P - lo0bits(&L); |
||
436 | #ifdef VAX |
||
437 | L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11; |
||
438 | #else |
||
439 | L = word0(d) | Exp_msk1; |
||
440 | #endif |
||
441 | return P - 32 - lo0bits(&L); |
||
442 | } |
||
443 | |||
444 | int |
||
445 | _strtodg_r |
||
446 | #ifdef KR_headers |
||
447 | (p, s00, se, fpi, exp, bits) |
||
448 | struct _reent *p; const char *s00; char **se; FPI *fpi; Long *exp; __ULong *bits; |
||
449 | #else |
||
450 | (struct _reent *p, const char *s00, char **se, FPI *fpi, Long *exp, __ULong *bits) |
||
451 | #endif |
||
452 | { |
||
453 | int abe, abits, asub; |
||
454 | int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm; |
||
455 | int dsign, e, e1, e2, emin, esign, finished, i, inex, irv; |
||
456 | int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign; |
||
457 | int sudden_underflow; |
||
458 | const char *s, *s0, *s1; |
||
459 | //double adj, adj0, rv, tol; |
||
460 | double adj0, tol; |
||
461 | U adj, rv; |
||
462 | Long L; |
||
463 | __ULong y, z; |
||
464 | _Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0; |
||
465 | |||
466 | irv = STRTOG_Zero; |
||
467 | denorm = sign = nz0 = nz = 0; |
||
468 | dval(rv) = 0.; |
||
469 | rvb = 0; |
||
470 | nbits = fpi->nbits; |
||
471 | for(s = s00;;s++) switch(*s) { |
||
472 | case '-': |
||
473 | sign = 1; |
||
474 | /* no break */ |
||
475 | case '+': |
||
476 | if (*++s) |
||
477 | goto break2; |
||
478 | /* no break */ |
||
479 | case 0: |
||
480 | sign = 0; |
||
481 | irv = STRTOG_NoNumber; |
||
482 | s = s00; |
||
483 | goto ret; |
||
484 | case '\t': |
||
485 | case '\n': |
||
486 | case '\v': |
||
487 | case '\f': |
||
488 | case '\r': |
||
489 | case ' ': |
||
490 | continue; |
||
491 | default: |
||
492 | goto break2; |
||
493 | } |
||
494 | break2: |
||
495 | if (*s == '0') { |
||
496 | #ifndef NO_HEX_FP |
||
497 | switch(s[1]) { |
||
498 | case 'x': |
||
499 | case 'X': |
||
500 | irv = gethex(p, &s, fpi, exp, &rvb, sign); |
||
501 | if (irv == STRTOG_NoNumber) { |
||
502 | s = s00; |
||
503 | sign = 0; |
||
504 | } |
||
505 | goto ret; |
||
506 | } |
||
507 | #endif |
||
508 | nz0 = 1; |
||
509 | while(*++s == '0') ; |
||
510 | if (!*s) |
||
511 | goto ret; |
||
512 | } |
||
513 | sudden_underflow = fpi->sudden_underflow; |
||
514 | s0 = s; |
||
515 | y = z = 0; |
||
516 | for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) |
||
517 | if (nd < 9) |
||
518 | y = 10*y + c - '0'; |
||
519 | else if (nd < 16) |
||
520 | z = 10*z + c - '0'; |
||
521 | nd0 = nd; |
||
522 | #ifdef USE_LOCALE |
||
523 | if (strncmp (s, _localeconv_r (p)->decimal_point, |
||
524 | strlen (_localeconv_r (p)->decimal_point)) == 0) |
||
525 | #else |
||
526 | if (c == '.') |
||
527 | #endif |
||
528 | { |
||
529 | decpt = 1; |
||
530 | #ifdef USE_LOCALE |
||
531 | c = *(s += strlen (_localeconv_r (p)->decimal_point)); |
||
532 | #else |
||
533 | c = *++s; |
||
534 | #endif |
||
535 | if (!nd) { |
||
536 | for(; c == '0'; c = *++s) |
||
537 | nz++; |
||
538 | if (c > '0' && c <= '9') { |
||
539 | s0 = s; |
||
540 | nf += nz; |
||
541 | nz = 0; |
||
542 | goto have_dig; |
||
543 | } |
||
544 | goto dig_done; |
||
545 | } |
||
546 | for(; c >= '0' && c <= '9'; c = *++s) { |
||
547 | have_dig: |
||
548 | nz++; |
||
549 | if (c -= '0') { |
||
550 | nf += nz; |
||
551 | for(i = 1; i < nz; i++) |
||
552 | if (nd++ < 9) |
||
553 | y *= 10; |
||
554 | else if (nd <= DBL_DIG + 1) |
||
555 | z *= 10; |
||
556 | if (nd++ < 9) |
||
557 | y = 10*y + c; |
||
558 | else if (nd <= DBL_DIG + 1) |
||
559 | z = 10*z + c; |
||
560 | nz = 0; |
||
561 | } |
||
562 | } |
||
563 | } |
||
564 | dig_done: |
||
565 | e = 0; |
||
566 | if (c == 'e' || c == 'E') { |
||
567 | if (!nd && !nz && !nz0) { |
||
568 | irv = STRTOG_NoNumber; |
||
569 | s = s00; |
||
570 | goto ret; |
||
571 | } |
||
572 | s00 = s; |
||
573 | esign = 0; |
||
574 | switch(c = *++s) { |
||
575 | case '-': |
||
576 | esign = 1; |
||
577 | case '+': |
||
578 | c = *++s; |
||
579 | } |
||
580 | if (c >= '0' && c <= '9') { |
||
581 | while(c == '0') |
||
582 | c = *++s; |
||
583 | if (c > '0' && c <= '9') { |
||
584 | L = c - '0'; |
||
585 | s1 = s; |
||
586 | while((c = *++s) >= '0' && c <= '9') |
||
587 | L = 10*L + c - '0'; |
||
588 | if (s - s1 > 8 || L > 19999) |
||
589 | /* Avoid confusion from exponents |
||
590 | * so large that e might overflow. |
||
591 | */ |
||
592 | e = 19999; /* safe for 16 bit ints */ |
||
593 | else |
||
594 | e = (int)L; |
||
595 | if (esign) |
||
596 | e = -e; |
||
597 | } |
||
598 | else |
||
599 | e = 0; |
||
600 | } |
||
601 | else |
||
602 | s = s00; |
||
603 | } |
||
604 | if (!nd) { |
||
605 | if (!nz && !nz0) { |
||
606 | #ifdef INFNAN_CHECK |
||
607 | /* Check for Nan and Infinity */ |
||
608 | if (!decpt) |
||
609 | switch(c) { |
||
610 | case 'i': |
||
611 | case 'I': |
||
612 | if (match(&s,"nf")) { |
||
613 | --s; |
||
614 | if (!match(&s,"inity")) |
||
615 | ++s; |
||
616 | irv = STRTOG_Infinite; |
||
617 | goto infnanexp; |
||
618 | } |
||
619 | break; |
||
620 | case 'n': |
||
621 | case 'N': |
||
622 | if (match(&s, "an")) { |
||
623 | irv = STRTOG_NaN; |
||
624 | *exp = fpi->emax + 1; |
||
625 | #ifndef No_Hex_NaN |
||
626 | if (*s == '(') /*)*/ |
||
627 | irv = hexnan(&s, fpi, bits); |
||
628 | #endif |
||
629 | goto infnanexp; |
||
630 | } |
||
631 | } |
||
632 | #endif /* INFNAN_CHECK */ |
||
633 | irv = STRTOG_NoNumber; |
||
634 | s = s00; |
||
635 | } |
||
636 | goto ret; |
||
637 | } |
||
638 | |||
639 | irv = STRTOG_Normal; |
||
640 | e1 = e -= nf; |
||
641 | rd = 0; |
||
642 | switch(fpi->rounding & 3) { |
||
643 | case FPI_Round_up: |
||
644 | rd = 2 - sign; |
||
645 | break; |
||
646 | case FPI_Round_zero: |
||
647 | rd = 1; |
||
648 | break; |
||
649 | case FPI_Round_down: |
||
650 | rd = 1 + sign; |
||
651 | } |
||
652 | |||
653 | /* Now we have nd0 digits, starting at s0, followed by a |
||
654 | * decimal point, followed by nd-nd0 digits. The number we're |
||
655 | * after is the integer represented by those digits times |
||
656 | * 10**e */ |
||
657 | |||
658 | if (!nd0) |
||
659 | nd0 = nd; |
||
660 | k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; |
||
661 | dval(rv) = y; |
||
662 | if (k > 9) |
||
663 | dval(rv) = tens[k - 9] * dval(rv) + z; |
||
664 | bd0 = 0; |
||
665 | if (nbits <= P && nd <= DBL_DIG) { |
||
666 | if (!e) { |
||
667 | if (rvOK(p, dval(rv), fpi, exp, bits, 1, rd, &irv)) |
||
668 | goto ret; |
||
669 | } |
||
670 | else if (e > 0) { |
||
671 | if (e <= Ten_pmax) { |
||
672 | #ifdef VAX |
||
673 | goto vax_ovfl_check; |
||
674 | #else |
||
675 | i = fivesbits[e] + mantbits(rv) <= P; |
||
676 | /* rv = */ rounded_product(dval(rv), tens[e]); |
||
677 | if (rvOK(p, dval(rv), fpi, exp, bits, i, rd, &irv)) |
||
678 | goto ret; |
||
679 | e1 -= e; |
||
680 | goto rv_notOK; |
||
681 | #endif |
||
682 | } |
||
683 | i = DBL_DIG - nd; |
||
684 | if (e <= Ten_pmax + i) { |
||
685 | /* A fancier test would sometimes let us do |
||
686 | * this for larger i values. |
||
687 | */ |
||
688 | e2 = e - i; |
||
689 | e1 -= i; |
||
690 | dval(rv) *= tens[i]; |
||
691 | #ifdef VAX |
||
692 | /* VAX exponent range is so narrow we must |
||
693 | * worry about overflow here... |
||
694 | */ |
||
695 | vax_ovfl_check: |
||
696 | dval(adj) = dval(rv); |
||
697 | word0(adj) -= P*Exp_msk1; |
||
698 | /* adj = */ rounded_product(dval(adj), tens[e2]); |
||
699 | if ((word0(adj) & Exp_mask) |
||
700 | > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) |
||
701 | goto rv_notOK; |
||
702 | word0(adj) += P*Exp_msk1; |
||
703 | dval(rv) = dval(adj); |
||
704 | #else |
||
705 | /* rv = */ rounded_product(dval(rv), tens[e2]); |
||
706 | #endif |
||
707 | if (rvOK(p, dval(rv), fpi, exp, bits, 0, rd, &irv)) |
||
708 | goto ret; |
||
709 | e1 -= e2; |
||
710 | } |
||
711 | } |
||
712 | #ifndef Inaccurate_Divide |
||
713 | else if (e >= -Ten_pmax) { |
||
714 | /* rv = */ rounded_quotient(dval(rv), tens[-e]); |
||
715 | if (rvOK(p, dval(rv), fpi, exp, bits, 0, rd, &irv)) |
||
716 | goto ret; |
||
717 | e1 -= e; |
||
718 | } |
||
719 | #endif |
||
720 | } |
||
721 | rv_notOK: |
||
722 | e1 += nd - k; |
||
723 | |||
724 | /* Get starting approximation = rv * 10**e1 */ |
||
725 | |||
726 | e2 = 0; |
||
727 | if (e1 > 0) { |
||
728 | if ( (i = e1 & 15) !=0) |
||
729 | dval(rv) *= tens[i]; |
||
730 | if (e1 &= ~15) { |
||
731 | e1 >>= 4; |
||
732 | while(e1 >= (1 << n_bigtens-1)) { |
||
733 | e2 += ((word0(rv) & Exp_mask) |
||
734 | >> Exp_shift1) - Bias; |
||
735 | word0(rv) &= ~Exp_mask; |
||
736 | word0(rv) |= Bias << Exp_shift1; |
||
737 | dval(rv) *= bigtens[n_bigtens-1]; |
||
738 | e1 -= 1 << n_bigtens-1; |
||
739 | } |
||
740 | e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias; |
||
741 | word0(rv) &= ~Exp_mask; |
||
742 | word0(rv) |= Bias << Exp_shift1; |
||
743 | for(j = 0; e1 > 0; j++, e1 >>= 1) |
||
744 | if (e1 & 1) |
||
745 | dval(rv) *= bigtens[j]; |
||
746 | } |
||
747 | } |
||
748 | else if (e1 < 0) { |
||
749 | e1 = -e1; |
||
750 | if ( (i = e1 & 15) !=0) |
||
751 | dval(rv) /= tens[i]; |
||
752 | if (e1 &= ~15) { |
||
753 | e1 >>= 4; |
||
754 | while(e1 >= (1 << n_bigtens-1)) { |
||
755 | e2 += ((word0(rv) & Exp_mask) |
||
756 | >> Exp_shift1) - Bias; |
||
757 | word0(rv) &= ~Exp_mask; |
||
758 | word0(rv) |= Bias << Exp_shift1; |
||
759 | dval(rv) *= tinytens[n_bigtens-1]; |
||
760 | e1 -= 1 << n_bigtens-1; |
||
761 | } |
||
762 | e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias; |
||
763 | word0(rv) &= ~Exp_mask; |
||
764 | word0(rv) |= Bias << Exp_shift1; |
||
765 | for(j = 0; e1 > 0; j++, e1 >>= 1) |
||
766 | if (e1 & 1) |
||
767 | dval(rv) *= tinytens[j]; |
||
768 | } |
||
769 | } |
||
770 | #ifdef IBM |
||
771 | /* e2 is a correction to the (base 2) exponent of the return |
||
772 | * value, reflecting adjustments above to avoid overflow in the |
||
773 | * native arithmetic. For native IBM (base 16) arithmetic, we |
||
774 | * must multiply e2 by 4 to change from base 16 to 2. |
||
775 | */ |
||
776 | e2 <<= 2; |
||
777 | #endif |
||
778 | rvb = d2b(p, dval(rv), &rve, &rvbits); /* rv = rvb * 2^rve */ |
||
779 | rve += e2; |
||
780 | if ((j = rvbits - nbits) > 0) { |
||
781 | rshift(rvb, j); |
||
782 | rvbits = nbits; |
||
783 | rve += j; |
||
784 | } |
||
785 | bb0 = 0; /* trailing zero bits in rvb */ |
||
786 | e2 = rve + rvbits - nbits; |
||
787 | if (e2 > fpi->emax + 1) |
||
788 | goto huge; |
||
789 | rve1 = rve + rvbits - nbits; |
||
790 | if (e2 < (emin = fpi->emin)) { |
||
791 | denorm = 1; |
||
792 | j = rve - emin; |
||
793 | if (j > 0) { |
||
794 | rvb = lshift(p, rvb, j); |
||
795 | rvbits += j; |
||
796 | } |
||
797 | else if (j < 0) { |
||
798 | rvbits += j; |
||
799 | if (rvbits <= 0) { |
||
800 | if (rvbits < -1) { |
||
801 | ufl: |
||
802 | rvb->_wds = 0; |
||
803 | rvb->_x[0] = 0; |
||
804 | *exp = emin; |
||
805 | irv = STRTOG_Underflow | STRTOG_Inexlo; |
||
806 | goto ret; |
||
807 | } |
||
808 | rvb->_x[0] = rvb->_wds = rvbits = 1; |
||
809 | } |
||
810 | else |
||
811 | rshift(rvb, -j); |
||
812 | } |
||
813 | rve = rve1 = emin; |
||
814 | if (sudden_underflow && e2 + 1 < emin) |
||
815 | goto ufl; |
||
816 | } |
||
817 | |||
818 | /* Now the hard part -- adjusting rv to the correct value.*/ |
||
819 | |||
820 | /* Put digits into bd: true value = bd * 10^e */ |
||
821 | |||
822 | bd0 = s2b(p, s0, nd0, nd, y); |
||
823 | |||
824 | for(;;) { |
||
825 | bd = Balloc(p,bd0->_k); |
||
826 | Bcopy(bd, bd0); |
||
827 | bb = Balloc(p,rvb->_k); |
||
828 | Bcopy(bb, rvb); |
||
829 | bbbits = rvbits - bb0; |
||
830 | bbe = rve + bb0; |
||
831 | bs = i2b(p, 1); |
||
832 | |||
833 | if (e >= 0) { |
||
834 | bb2 = bb5 = 0; |
||
835 | bd2 = bd5 = e; |
||
836 | } |
||
837 | else { |
||
838 | bb2 = bb5 = -e; |
||
839 | bd2 = bd5 = 0; |
||
840 | } |
||
841 | if (bbe >= 0) |
||
842 | bb2 += bbe; |
||
843 | else |
||
844 | bd2 -= bbe; |
||
845 | bs2 = bb2; |
||
846 | j = nbits + 1 - bbbits; |
||
847 | i = bbe + bbbits - nbits; |
||
848 | if (i < emin) /* denormal */ |
||
849 | j += i - emin; |
||
850 | bb2 += j; |
||
851 | bd2 += j; |
||
852 | i = bb2 < bd2 ? bb2 : bd2; |
||
853 | if (i > bs2) |
||
854 | i = bs2; |
||
855 | if (i > 0) { |
||
856 | bb2 -= i; |
||
857 | bd2 -= i; |
||
858 | bs2 -= i; |
||
859 | } |
||
860 | if (bb5 > 0) { |
||
861 | bs = pow5mult(p, bs, bb5); |
||
862 | bb1 = mult(p, bs, bb); |
||
863 | Bfree(p,bb); |
||
864 | bb = bb1; |
||
865 | } |
||
866 | bb2 -= bb0; |
||
867 | if (bb2 > 0) |
||
868 | bb = lshift(p, bb, bb2); |
||
869 | else if (bb2 < 0) |
||
870 | rshift(bb, -bb2); |
||
871 | if (bd5 > 0) |
||
872 | bd = pow5mult(p, bd, bd5); |
||
873 | if (bd2 > 0) |
||
874 | bd = lshift(p, bd, bd2); |
||
875 | if (bs2 > 0) |
||
876 | bs = lshift(p, bs, bs2); |
||
877 | asub = 1; |
||
878 | inex = STRTOG_Inexhi; |
||
879 | delta = diff(p, bb, bd); |
||
880 | if (delta->_wds <= 1 && !delta->_x[0]) |
||
881 | break; |
||
882 | dsign = delta->_sign; |
||
883 | delta->_sign = finished = 0; |
||
884 | L = 0; |
||
885 | i = cmp(delta, bs); |
||
886 | if (rd && i <= 0) { |
||
887 | irv = STRTOG_Normal; |
||
888 | if ( (finished = dsign ^ (rd&1)) !=0) { |
||
889 | if (dsign != 0) { |
||
890 | irv |= STRTOG_Inexhi; |
||
891 | goto adj1; |
||
892 | } |
||
893 | irv |= STRTOG_Inexlo; |
||
894 | if (rve1 == emin) |
||
895 | goto adj1; |
||
896 | for(i = 0, j = nbits; j >= ULbits; |
||
897 | i++, j -= ULbits) { |
||
898 | if (rvb->_x[i] & ALL_ON) |
||
899 | goto adj1; |
||
900 | } |
||
901 | if (j > 1 && lo0bits(rvb->_x + i) < j - 1) |
||
902 | goto adj1; |
||
903 | rve = rve1 - 1; |
||
904 | rvb = set_ones(p, rvb, rvbits = nbits); |
||
905 | break; |
||
906 | } |
||
907 | irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi; |
||
908 | break; |
||
909 | } |
||
910 | if (i < 0) { |
||
911 | /* Error is less than half an ulp -- check for |
||
912 | * special case of mantissa a power of two. |
||
913 | */ |
||
914 | irv = dsign |
||
915 | ? STRTOG_Normal | STRTOG_Inexlo |
||
916 | : STRTOG_Normal | STRTOG_Inexhi; |
||
917 | if (dsign || bbbits > 1 || denorm || rve1 == emin) |
||
918 | break; |
||
919 | delta = lshift(p, delta,1); |
||
920 | if (cmp(delta, bs) > 0) { |
||
921 | irv = STRTOG_Normal | STRTOG_Inexlo; |
||
922 | goto drop_down; |
||
923 | } |
||
924 | break; |
||
925 | } |
||
926 | if (i == 0) { |
||
927 | /* exactly half-way between */ |
||
928 | if (dsign) { |
||
929 | if (denorm && all_on(rvb, rvbits)) { |
||
930 | /*boundary case -- increment exponent*/ |
||
931 | rvb->_wds = 1; |
||
932 | rvb->_x[0] = 1; |
||
933 | rve = emin + nbits - (rvbits = 1); |
||
934 | irv = STRTOG_Normal | STRTOG_Inexhi; |
||
935 | denorm = 0; |
||
936 | break; |
||
937 | } |
||
938 | irv = STRTOG_Normal | STRTOG_Inexlo; |
||
939 | } |
||
940 | else if (bbbits == 1) { |
||
941 | irv = STRTOG_Normal; |
||
942 | drop_down: |
||
943 | /* boundary case -- decrement exponent */ |
||
944 | if (rve1 == emin) { |
||
945 | irv = STRTOG_Normal | STRTOG_Inexhi; |
||
946 | if (rvb->_wds == 1 && rvb->_x[0] == 1) |
||
947 | sudden_underflow = 1; |
||
948 | break; |
||
949 | } |
||
950 | rve -= nbits; |
||
951 | rvb = set_ones(p, rvb, rvbits = nbits); |
||
952 | break; |
||
953 | } |
||
954 | else |
||
955 | irv = STRTOG_Normal | STRTOG_Inexhi; |
||
956 | if (bbbits < nbits && !denorm || !(rvb->_x[0] & 1)) |
||
957 | break; |
||
958 | if (dsign) { |
||
959 | rvb = increment(p, rvb); |
||
960 | j = kmask & (ULbits - (rvbits & kmask)); |
||
961 | if (hi0bits(rvb->_x[rvb->_wds - 1]) != j) |
||
962 | rvbits++; |
||
963 | irv = STRTOG_Normal | STRTOG_Inexhi; |
||
964 | } |
||
965 | else { |
||
966 | if (bbbits == 1) |
||
967 | goto undfl; |
||
968 | decrement(rvb); |
||
969 | irv = STRTOG_Normal | STRTOG_Inexlo; |
||
970 | } |
||
971 | break; |
||
972 | } |
||
973 | if ((dval(adj) = ratio(delta, bs)) <= 2.) { |
||
974 | adj1: |
||
975 | inex = STRTOG_Inexlo; |
||
976 | if (dsign) { |
||
977 | asub = 0; |
||
978 | inex = STRTOG_Inexhi; |
||
979 | } |
||
980 | else if (denorm && bbbits <= 1) { |
||
981 | undfl: |
||
982 | rvb->_wds = 0; |
||
983 | rve = emin; |
||
984 | irv = STRTOG_Underflow | STRTOG_Inexlo; |
||
985 | break; |
||
986 | } |
||
987 | adj0 = dval(adj) = 1.; |
||
988 | } |
||
989 | else { |
||
990 | adj0 = dval(adj) *= 0.5; |
||
991 | if (dsign) { |
||
992 | asub = 0; |
||
993 | inex = STRTOG_Inexlo; |
||
994 | } |
||
995 | if (dval(adj) < 2147483647.) { |
||
996 | L = adj0; |
||
997 | adj0 -= L; |
||
998 | switch(rd) { |
||
999 | case 0: |
||
1000 | if (adj0 >= .5) |
||
1001 | goto inc_L; |
||
1002 | break; |
||
1003 | case 1: |
||
1004 | if (asub && adj0 > 0.) |
||
1005 | goto inc_L; |
||
1006 | break; |
||
1007 | case 2: |
||
1008 | if (!asub && adj0 > 0.) { |
||
1009 | inc_L: |
||
1010 | L++; |
||
1011 | inex = STRTOG_Inexact - inex; |
||
1012 | } |
||
1013 | } |
||
1014 | dval(adj) = L; |
||
1015 | } |
||
1016 | } |
||
1017 | y = rve + rvbits; |
||
1018 | |||
1019 | /* adj *= ulp(dval(rv)); */ |
||
1020 | /* if (asub) rv -= adj; else rv += adj; */ |
||
1021 | |||
1022 | if (!denorm && rvbits < nbits) { |
||
1023 | rvb = lshift(p, rvb, j = nbits - rvbits); |
||
1024 | rve -= j; |
||
1025 | rvbits = nbits; |
||
1026 | } |
||
1027 | ab = d2b(p, dval(adj), &abe, &abits); |
||
1028 | if (abe < 0) |
||
1029 | rshift(ab, -abe); |
||
1030 | else if (abe > 0) |
||
1031 | ab = lshift(p, ab, abe); |
||
1032 | rvb0 = rvb; |
||
1033 | if (asub) { |
||
1034 | /* rv -= adj; */ |
||
1035 | j = hi0bits(rvb->_x[rvb->_wds-1]); |
||
1036 | rvb = diff(p, rvb, ab); |
||
1037 | k = rvb0->_wds - 1; |
||
1038 | if (denorm) |
||
1039 | /* do nothing */; |
||
1040 | else if (rvb->_wds <= k |
||
1041 | || hi0bits( rvb->_x[k]) > |
||
1042 | hi0bits(rvb0->_x[k])) { |
||
1043 | /* unlikely; can only have lost 1 high bit */ |
||
1044 | if (rve1 == emin) { |
||
1045 | --rvbits; |
||
1046 | denorm = 1; |
||
1047 | } |
||
1048 | else { |
||
1049 | rvb = lshift(p, rvb, 1); |
||
1050 | --rve; |
||
1051 | --rve1; |
||
1052 | L = finished = 0; |
||
1053 | } |
||
1054 | } |
||
1055 | } |
||
1056 | else { |
||
1057 | rvb = sum(p, rvb, ab); |
||
1058 | k = rvb->_wds - 1; |
||
1059 | if (k >= rvb0->_wds |
||
1060 | || hi0bits(rvb->_x[k]) < hi0bits(rvb0->_x[k])) { |
||
1061 | if (denorm) { |
||
1062 | if (++rvbits == nbits) |
||
1063 | denorm = 0; |
||
1064 | } |
||
1065 | else { |
||
1066 | rshift(rvb, 1); |
||
1067 | rve++; |
||
1068 | rve1++; |
||
1069 | L = 0; |
||
1070 | } |
||
1071 | } |
||
1072 | } |
||
1073 | Bfree(p,ab); |
||
1074 | Bfree(p,rvb0); |
||
1075 | if (finished) |
||
1076 | break; |
||
1077 | |||
1078 | z = rve + rvbits; |
||
1079 | if (y == z && L) { |
||
1080 | /* Can we stop now? */ |
||
1081 | tol = dval(adj) * 5e-16; /* > max rel error */ |
||
1082 | dval(adj) = adj0 - .5; |
||
1083 | if (dval(adj) < -tol) { |
||
1084 | if (adj0 > tol) { |
||
1085 | irv |= inex; |
||
1086 | break; |
||
1087 | } |
||
1088 | } |
||
1089 | else if (dval(adj) > tol && adj0 < 1. - tol) { |
||
1090 | irv |= inex; |
||
1091 | break; |
||
1092 | } |
||
1093 | } |
||
1094 | bb0 = denorm ? 0 : trailz(rvb); |
||
1095 | Bfree(p,bb); |
||
1096 | Bfree(p,bd); |
||
1097 | Bfree(p,bs); |
||
1098 | Bfree(p,delta); |
||
1099 | } |
||
1100 | if (!denorm && (j = nbits - rvbits)) { |
||
1101 | if (j > 0) |
||
1102 | rvb = lshift(p, rvb, j); |
||
1103 | else |
||
1104 | rshift(rvb, -j); |
||
1105 | rve -= j; |
||
1106 | } |
||
1107 | *exp = rve; |
||
1108 | Bfree(p,bb); |
||
1109 | Bfree(p,bd); |
||
1110 | Bfree(p,bs); |
||
1111 | Bfree(p,bd0); |
||
1112 | Bfree(p,delta); |
||
1113 | if (rve > fpi->emax) { |
||
1114 | huge: |
||
1115 | rvb->_wds = 0; |
||
1116 | irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; |
||
1117 | #ifndef NO_ERRNO |
||
1118 | errno = ERANGE; |
||
1119 | #endif |
||
1120 | infnanexp: |
||
1121 | *exp = fpi->emax + 1; |
||
1122 | } |
||
1123 | ret: |
||
1124 | if (denorm) { |
||
1125 | if (sudden_underflow) { |
||
1126 | rvb->_wds = 0; |
||
1127 | irv = STRTOG_Underflow | STRTOG_Inexlo; |
||
1128 | } |
||
1129 | else { |
||
1130 | irv = (irv & ~STRTOG_Retmask) | |
||
1131 | (rvb->_wds > 0 ? STRTOG_Denormal : STRTOG_Zero); |
||
1132 | if (irv & STRTOG_Inexact) |
||
1133 | irv |= STRTOG_Underflow; |
||
1134 | } |
||
1135 | } |
||
1136 | if (se) |
||
1137 | *se = (char *)s; |
||
1138 | if (sign) |
||
1139 | irv |= STRTOG_Neg; |
||
1140 | if (rvb) { |
||
1141 | copybits(bits, nbits, rvb); |
||
1142 | Bfree(p,rvb); |
||
1143 | } |
||
1144 | return irv; |
||
1145 | } |
||
1146 | |||
1147 | #endif /* _HAVE_LONG_DOUBLE && !_LDBL_EQ_DBL */>>>=>>>>=>=>>>>=>=>>>>>>=>>>=><=>><>><>><>><>>><>><>><>><>=>=>=>=>=>>=>=>=>=>>=>>>=>=>>>=>><>><>><>>>><>>=>>>><>><>>>>>>>>>>><>>>>>>>> |