| 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 |
|
|---|
