OSDN Git Service

2005-03-05 Andreas Tobler <a.tobler@schweiz.ch>
[pf3gnuchains/gcc-fork.git] / libjava / java / lang / strtod.c
1 /*
2 FUNCTION
3         <<strtod>>, <<strtodf>>---string to double or float
4
5 INDEX
6         strtod
7 INDEX
8         _strtod_r
9 INDEX
10         strtodf
11
12 ANSI_SYNOPSIS
13         #include <stdlib.h>
14         double strtod(const char *<[str]>, char **<[tail]>);
15         float strtodf(const char *<[str]>, char **<[tail]>);
16
17         double _strtod_r(void *<[reent]>, 
18                          const char *<[str]>, char **<[tail]>);
19
20 TRAD_SYNOPSIS
21         #include <stdlib.h>
22         double strtod(<[str]>,<[tail]>)
23         char *<[str]>;
24         char **<[tail]>;
25
26         float strtodf(<[str]>,<[tail]>)
27         char *<[str]>;
28         char **<[tail]>;
29
30         double _strtod_r(<[reent]>,<[str]>,<[tail]>)
31         char *<[reent]>;
32         char *<[str]>;
33         char **<[tail]>;
34
35 DESCRIPTION
36         The function <<strtod>> parses the character string <[str]>,
37         producing a substring which can be converted to a double
38         value.  The substring converted is the longest initial
39         subsequence of <[str]>, beginning with the first
40         non-whitespace character, that has the format:
41         .[+|-]<[digits]>[.][<[digits]>][(e|E)[+|-]<[digits]>] 
42         The substring contains no characters if <[str]> is empty, consists
43         entirely of whitespace, or if the first non-whitespace
44         character is something other than <<+>>, <<->>, <<.>>, or a
45         digit. If the substring is empty, no conversion is done, and
46         the value of <[str]> is stored in <<*<[tail]>>>.  Otherwise,
47         the substring is converted, and a pointer to the final string
48         (which will contain at least the terminating null character of
49         <[str]>) is stored in <<*<[tail]>>>.  If you want no
50         assignment to <<*<[tail]>>>, pass a null pointer as <[tail]>.
51         <<strtodf>> is identical to <<strtod>> except for its return type.
52
53         This implementation returns the nearest machine number to the
54         input decimal string.  Ties are broken by using the IEEE
55         round-even rule.
56
57         The alternate function <<_strtod_r>> is a reentrant version.
58         The extra argument <[reent]> is a pointer to a reentrancy structure.
59
60 RETURNS
61         <<strtod>> returns the converted substring value, if any.  If
62         no conversion could be performed, 0 is returned.  If the
63         correct value is out of the range of representable values,
64         plus or minus <<HUGE_VAL>> is returned, and <<ERANGE>> is
65         stored in errno. If the correct value would cause underflow, 0
66         is returned and <<ERANGE>> is stored in errno.
67
68 Supporting OS subroutines required: <<close>>, <<fstat>>, <<isatty>>,
69 <<lseek>>, <<read>>, <<sbrk>>, <<write>>.
70 */
71
72 /****************************************************************
73  *
74  * The author of this software is David M. Gay.
75  *
76  * Copyright (c) 1991 by AT&T.
77  *
78  * Permission to use, copy, modify, and distribute this software for any
79  * purpose without fee is hereby granted, provided that this entire notice
80  * is included in all copies of any software which is or includes a copy
81  * or modification of this software and in all copies of the supporting
82  * documentation for such software.
83  *
84  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
85  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
86  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
87  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
88  *
89  ***************************************************************/
90
91 /* Please send bug reports to
92         David M. Gay
93         AT&T Bell Laboratories, Room 2C-463
94         600 Mountain Avenue
95         Murray Hill, NJ 07974-2070
96         U.S.A.
97         dmg@research.att.com or research!dmg
98  */
99
100 #include <string.h>
101 #include <float.h>
102 #include <errno.h>
103 #include "mprec.h"
104
105 double
106 _DEFUN (_strtod_r, (ptr, s00, se),
107         struct _Jv_reent *ptr _AND
108         _CONST char *s00 _AND
109         char **se)
110 {
111   int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, e1, esign, i, j,
112     k, nd, nd0, nf, nz, nz0, sign;
113   int digits = 0;  /* Number of digits found in fraction part. */
114   long e;
115   _CONST char *s, *s0, *s1;
116   double aadj, aadj1, adj;
117   long L;
118   unsigned long y, z;
119   union double_union rv, rv0;
120
121   _Jv_Bigint *bb = NULL, *bb1, *bd = NULL, *bd0, *bs = NULL, *delta = NULL;
122   sign = nz0 = nz = 0;
123   rv.d = 0.;
124   for (s = s00;; s++)
125     switch (*s)
126       {
127       case '-':
128         sign = 1;
129         /* no break */
130       case '+':
131         if (*++s)
132           goto break2;
133         /* no break */
134       case 0:
135         s = s00;
136         goto ret;
137       case '\t':
138       case '\n':
139       case '\v':
140       case '\f':
141       case '\r':
142       case ' ':
143         continue;
144       default:
145         goto break2;
146       }
147 break2:
148   if (*s == '0')
149     {
150       digits++;
151       nz0 = 1;
152       while (*++s == '0')
153         digits++;
154       if (!*s)
155         goto ret;
156     }
157   s0 = s;
158   y = z = 0;
159   for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
160     {
161       digits++;
162       if (nd < 9)
163         y = 10 * y + c - '0';
164       else if (nd < 16)
165         z = 10 * z + c - '0';
166     }
167   nd0 = nd;
168   if (c == '.')
169     {
170       c = *++s;
171       if (!nd)
172         {
173           for (; c == '0'; c = *++s)
174             {
175               digits++;
176               nz++;
177             }
178           if (c > '0' && c <= '9')
179             {
180               digits++;
181               s0 = s;
182               nf += nz;
183               nz = 0;
184               goto have_dig;
185             }
186           goto dig_done;
187         }
188       for (; c >= '0' && c <= '9'; c = *++s)
189         {
190           digits++;
191         have_dig:
192           nz++;
193           if (c -= '0')
194             {
195               nf += nz;
196               for (i = 1; i < nz; i++)
197                 if (nd++ < 9)
198                   y *= 10;
199                 else if (nd <= DBL_DIG + 1)
200                   z *= 10;
201               if (nd++ < 9)
202                 y = 10 * y + c;
203               else if (nd <= DBL_DIG + 1)
204                 z = 10 * z + c;
205               nz = 0;
206             }
207         }
208     }
209 dig_done:
210   e = 0;
211   if (c == 'e' || c == 'E')
212     {
213       if (!nd && !nz && !nz0)
214         {
215           s = s00;
216           goto ret;
217         }
218       s00 = s;
219       esign = 0;
220       switch (c = *++s)
221         {
222         case '-':
223           esign = 1;
224         case '+':
225           c = *++s;
226         }
227       if (c >= '0' && c <= '9')
228         {
229           while (c == '0')
230             c = *++s;
231           if (c > '0' && c <= '9')
232             {
233               e = c - '0';
234               s1 = s;
235               while ((c = *++s) >= '0' && c <= '9')
236                 e = 10 * e + c - '0';
237               if (s - s1 > 8)
238                 /* Avoid confusion from exponents
239                  * so large that e might overflow.
240                  */
241                 e = 9999999L;
242               if (esign)
243                 e = -e;
244             }
245         }
246       else
247         {
248           /* No exponent after an 'E' : that's an error. */
249           ptr->_errno = EINVAL;
250           e = 0;
251           s = s00;
252           goto ret;
253         }
254     }
255   if (!nd)
256     {
257       if (!nz && !nz0)
258         s = s00;
259       goto ret;
260     }
261   e1 = e -= nf;
262
263   /* Now we have nd0 digits, starting at s0, followed by a
264    * decimal point, followed by nd-nd0 digits.  The number we're
265    * after is the integer represented by those digits times
266    * 10**e */
267
268   if (!nd0)
269     nd0 = nd;
270   k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
271   rv.d = y;
272   if (k > 9)
273     rv.d = tens[k - 9] * rv.d + z;
274   bd0 = 0;
275   if (nd <= DBL_DIG
276 #ifndef RND_PRODQUOT
277       && FLT_ROUNDS == 1
278 #endif
279     )
280     {
281       if (!e)
282         goto ret;
283       if (e > 0)
284         {
285           if (e <= Ten_pmax)
286             {
287 #ifdef VAX
288               goto vax_ovfl_check;
289 #else
290               /* rv.d = */ rounded_product (rv.d, tens[e]);
291               goto ret;
292 #endif
293             }
294           i = DBL_DIG - nd;
295           if (e <= Ten_pmax + i)
296             {
297               /* A fancier test would sometimes let us do
298                                  * this for larger i values.
299                                  */
300               e -= i;
301               rv.d *= tens[i];
302 #ifdef VAX
303               /* VAX exponent range is so narrow we must
304                * worry about overflow here...
305                */
306             vax_ovfl_check:
307               word0 (rv) -= P * Exp_msk1;
308               /* rv.d = */ rounded_product (rv.d, tens[e]);
309               if ((word0 (rv) & Exp_mask)
310                   > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P))
311                 goto ovfl;
312               word0 (rv) += P * Exp_msk1;
313 #else
314               /* rv.d = */ rounded_product (rv.d, tens[e]);
315 #endif
316               goto ret;
317             }
318         }
319 #ifndef Inaccurate_Divide
320       else if (e >= -Ten_pmax)
321         {
322           /* rv.d = */ rounded_quotient (rv.d, tens[-e]);
323           goto ret;
324         }
325 #endif
326     }
327   e1 += nd - k;
328
329   /* Get starting approximation = rv.d * 10**e1 */
330
331   if (e1 > 0)
332     {
333       if ((i = e1 & 15))
334         rv.d *= tens[i];
335
336       if (e1 &= ~15)
337         {
338           if (e1 > DBL_MAX_10_EXP)
339             {
340             ovfl:
341               ptr->_errno = ERANGE;
342
343               /* Force result to IEEE infinity. */
344               word0 (rv) = Exp_mask;
345               word1 (rv) = 0;
346
347               if (bd0)
348                 goto retfree;
349               goto ret;
350             }
351           if (e1 >>= 4)
352             {
353               for (j = 0; e1 > 1; j++, e1 >>= 1)
354                 if (e1 & 1)
355                   rv.d *= bigtens[j];
356               /* The last multiplication could overflow. */
357               word0 (rv) -= P * Exp_msk1;
358               rv.d *= bigtens[j];
359               if ((z = word0 (rv) & Exp_mask)
360                   > Exp_msk1 * (DBL_MAX_EXP + Bias - P))
361                 goto ovfl;
362               if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P))
363                 {
364                   /* set to largest number */
365                   /* (Can't trust DBL_MAX) */
366                   word0 (rv) = Big0;
367 #ifndef _DOUBLE_IS_32BITS
368                   word1 (rv) = Big1;
369 #endif
370                 }
371               else
372                 word0 (rv) += P * Exp_msk1;
373             }
374
375         }
376     }
377   else if (e1 < 0)
378     {
379       e1 = -e1;
380       if ((i = e1 & 15))
381         rv.d /= tens[i];
382       if (e1 &= ~15)
383         {
384           e1 >>= 4;
385           if (e1 >= 1 << n_bigtens)
386             goto undfl;
387           for (j = 0; e1 > 1; j++, e1 >>= 1)
388             if (e1 & 1)
389               rv.d *= tinytens[j];
390           /* The last multiplication could underflow. */
391           rv0.d = rv.d;
392           rv.d *= tinytens[j];
393           if (!rv.d)
394             {
395               rv.d = 2. * rv0.d;
396               rv.d *= tinytens[j];
397               if (!rv.d)
398                 {
399                 undfl:
400                   rv.d = 0.;
401                   ptr->_errno = ERANGE;
402                   if (bd0)
403                     goto retfree;
404                   goto ret;
405                 }
406 #ifndef _DOUBLE_IS_32BITS
407               word0 (rv) = Tiny0;
408               word1 (rv) = Tiny1;
409 #else
410               word0 (rv) = Tiny1;
411 #endif
412               /* The refinement below will clean
413                * this approximation up.
414                */
415             }
416         }
417     }
418
419   /* Now the hard part -- adjusting rv to the correct value.*/
420
421   /* Put digits into bd: true value = bd * 10^e */
422
423   bd0 = s2b (ptr, s0, nd0, nd, y);
424
425   for (;;)
426     {
427       bd = Balloc (ptr, bd0->_k);
428       Bcopy (bd, bd0);
429       bb = d2b (ptr, rv.d, &bbe, &bbbits);      /* rv.d = bb * 2^bbe */
430       bs = i2b (ptr, 1);
431
432       if (e >= 0)
433         {
434           bb2 = bb5 = 0;
435           bd2 = bd5 = e;
436         }
437       else
438         {
439           bb2 = bb5 = -e;
440           bd2 = bd5 = 0;
441         }
442       if (bbe >= 0)
443         bb2 += bbe;
444       else
445         bd2 -= bbe;
446       bs2 = bb2;
447 #ifdef Sudden_Underflow
448 #ifdef IBM
449       j = 1 + 4 * P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
450 #else
451       j = P + 1 - bbbits;
452 #endif
453 #else
454       i = bbe + bbbits - 1;     /* logb(rv.d) */
455       if (i < Emin)             /* denormal */
456         j = bbe + (P - Emin);
457       else
458         j = P + 1 - bbbits;
459 #endif
460       bb2 += j;
461       bd2 += j;
462       i = bb2 < bd2 ? bb2 : bd2;
463       if (i > bs2)
464         i = bs2;
465       if (i > 0)
466         {
467           bb2 -= i;
468           bd2 -= i;
469           bs2 -= i;
470         }
471       if (bb5 > 0)
472         {
473           bs = pow5mult (ptr, bs, bb5);
474           bb1 = mult (ptr, bs, bb);
475           Bfree (ptr, bb);
476           bb = bb1;
477         }
478       if (bb2 > 0)
479         bb = lshift (ptr, bb, bb2);
480       if (bd5 > 0)
481         bd = pow5mult (ptr, bd, bd5);
482       if (bd2 > 0)
483         bd = lshift (ptr, bd, bd2);
484       if (bs2 > 0)
485         bs = lshift (ptr, bs, bs2);
486       delta = diff (ptr, bb, bd);
487       dsign = delta->_sign;
488       delta->_sign = 0;
489       i = cmp (delta, bs);
490       if (i < 0)
491         {
492           /* Error is less than half an ulp -- check for
493            * special case of mantissa a power of two.
494            */
495           if (dsign || word1 (rv) || word0 (rv) & Bndry_mask)
496             break;
497           delta = lshift (ptr, delta, Log2P);
498           if (cmp (delta, bs) > 0)
499             goto drop_down;
500           break;
501         }
502       if (i == 0)
503         {
504           /* exactly half-way between */
505           if (dsign)
506             {
507               if ((word0 (rv) & Bndry_mask1) == Bndry_mask1
508                   && word1 (rv) == 0xffffffff)
509                 {
510                   /*boundary case -- increment exponent*/
511                   word0 (rv) = (word0 (rv) & Exp_mask)
512                     + Exp_msk1
513 #ifdef IBM
514                     | Exp_msk1 >> 4
515 #endif
516                     ;
517 #ifndef _DOUBLE_IS_32BITS
518                   word1 (rv) = 0;
519 #endif
520                   break;
521                 }
522             }
523           else if (!(word0 (rv) & Bndry_mask) && !word1 (rv))
524             {
525             drop_down:
526               /* boundary case -- decrement exponent */
527 #ifdef Sudden_Underflow
528               L = word0 (rv) & Exp_mask;
529 #ifdef IBM
530               if (L < Exp_msk1)
531 #else
532               if (L <= Exp_msk1)
533 #endif
534                 goto undfl;
535               L -= Exp_msk1;
536 #else
537               L = (word0 (rv) & Exp_mask) - Exp_msk1;
538 #endif
539               word0 (rv) = L | Bndry_mask1;
540 #ifndef _DOUBLE_IS_32BITS
541               word1 (rv) = 0xffffffff;
542 #endif
543 #ifdef IBM
544               goto cont;
545 #else
546               break;
547 #endif
548             }
549 #ifndef ROUND_BIASED
550           if (!(word1 (rv) & LSB))
551             break;
552 #endif
553           if (dsign)
554             rv.d += ulp (rv.d);
555 #ifndef ROUND_BIASED
556           else
557             {
558               rv.d -= ulp (rv.d);
559 #ifndef Sudden_Underflow
560               if (!rv.d)
561                 goto undfl;
562 #endif
563             }
564 #endif
565           break;
566         }
567       if ((aadj = ratio (delta, bs)) <= 2.)
568         {
569           if (dsign)
570             aadj = aadj1 = 1.;
571           else if (word1 (rv) || word0 (rv) & Bndry_mask)
572             {
573 #ifndef Sudden_Underflow
574               if (word1 (rv) == Tiny1 && !word0 (rv))
575                 goto undfl;
576 #endif
577               aadj = 1.;
578               aadj1 = -1.;
579             }
580           else
581             {
582               /* special case -- power of FLT_RADIX to be */
583               /* rounded down... */
584
585               if (aadj < 2. / FLT_RADIX)
586                 aadj = 1. / FLT_RADIX;
587               else
588                 aadj *= 0.5;
589               aadj1 = -aadj;
590             }
591         }
592       else
593         {
594           aadj *= 0.5;
595           aadj1 = dsign ? aadj : -aadj;
596 #ifdef Check_FLT_ROUNDS
597           switch (FLT_ROUNDS)
598             {
599             case 2:             /* towards +infinity */
600               aadj1 -= 0.5;
601               break;
602             case 0:             /* towards 0 */
603             case 3:             /* towards -infinity */
604               aadj1 += 0.5;
605             }
606 #else
607           if (FLT_ROUNDS == 0)
608             aadj1 += 0.5;
609 #endif
610         }
611       y = word0 (rv) & Exp_mask;
612
613       /* Check for overflow */
614
615       if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1))
616         {
617           rv0.d = rv.d;
618           word0 (rv) -= P * Exp_msk1;
619           adj = aadj1 * ulp (rv.d);
620           rv.d += adj;
621           if ((word0 (rv) & Exp_mask) >=
622               Exp_msk1 * (DBL_MAX_EXP + Bias - P))
623             {
624               if (word0 (rv0) == Big0 && word1 (rv0) == Big1)
625                 goto ovfl;
626 #ifdef _DOUBLE_IS_32BITS
627               word0 (rv) = Big1;
628 #else
629               word0 (rv) = Big0;
630               word1 (rv) = Big1;
631 #endif
632               goto cont;
633             }
634           else
635             word0 (rv) += P * Exp_msk1;
636         }
637       else
638         {
639 #ifdef Sudden_Underflow
640           if ((word0 (rv) & Exp_mask) <= P * Exp_msk1)
641             {
642               rv0.d = rv.d;
643               word0 (rv) += P * Exp_msk1;
644               adj = aadj1 * ulp (rv.d);
645               rv.d += adj;
646 #ifdef IBM
647               if ((word0 (rv) & Exp_mask) < P * Exp_msk1)
648 #else
649               if ((word0 (rv) & Exp_mask) <= P * Exp_msk1)
650 #endif
651                 {
652                   if (word0 (rv0) == Tiny0
653                       && word1 (rv0) == Tiny1)
654                     goto undfl;
655                   word0 (rv) = Tiny0;
656                   word1 (rv) = Tiny1;
657                   goto cont;
658                 }
659               else
660                 word0 (rv) -= P * Exp_msk1;
661             }
662           else
663             {
664               adj = aadj1 * ulp (rv.d);
665               rv.d += adj;
666             }
667 #else
668           /* Compute adj so that the IEEE rounding rules will
669            * correctly round rv.d + adj in some half-way cases.
670            * If rv.d * ulp(rv.d) is denormalized (i.e.,
671            * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
672            * trouble from bits lost to denormalization;
673            * example: 1.2e-307 .
674            */
675           if (y <= (P - 1) * Exp_msk1 && aadj >= 1.)
676             {
677               aadj1 = (double) (int) (aadj + 0.5);
678               if (!dsign)
679                 aadj1 = -aadj1;
680             }
681           adj = aadj1 * ulp (rv.d);
682           rv.d += adj;
683 #endif
684         }
685       z = word0 (rv) & Exp_mask;
686       if (y == z)
687         {
688           /* Can we stop now? */
689           L = aadj;
690           aadj -= L;
691           /* The tolerances below are conservative. */
692           if (dsign || word1 (rv) || word0 (rv) & Bndry_mask)
693             {
694               if (aadj < .4999999 || aadj > .5000001)
695                 break;
696             }
697           else if (aadj < .4999999 / FLT_RADIX)
698             break;
699         }
700     cont:
701       Bfree (ptr, bb);
702       Bfree (ptr, bd);
703       Bfree (ptr, bs);
704       Bfree (ptr, delta);
705     }
706 retfree:
707   Bfree (ptr, bb);
708   Bfree (ptr, bd);
709   Bfree (ptr, bs);
710   Bfree (ptr, bd0);
711   Bfree (ptr, delta);
712 ret:
713   if (se)
714     *se = (char *) s;
715   if (digits == 0)
716     ptr->_errno = EINVAL;
717   return sign ? -rv.d : rv.d;
718 }
719