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, *bb1, *bd, *bd0, *bs, *delta;
 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