1 // ---------------------------------------------------------------------- 2 // Copyright 2005-2014 Rich Felker, et al. 3 // Copyright 2014 The Android Open Source Project. 4 // 5 // Permission is hereby granted, free of charge, to any person obtaining 6 // a copy of this software and associated documentation files (the 7 // "Software"), to deal in the Software without restriction, including 8 // without limitation the rights to use, copy, modify, merge, publish, 9 // distribute, sublicense, and/or sell copies of the Software, and to 10 // permit persons to whom the Software is furnished to do so, subject to 11 // the following conditions: 12 // 13 // The above copyright notice and this permission notice shall be 14 // included in all copies or substantial portions of the Software. 15 // 16 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 17 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 18 // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 19 // IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY 20 // CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, 21 // TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE 22 // SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 23 // ---------------------------------------------------------------------- 24 25 #include <stdint.h> 26 #include <stdio.h> 27 #include <math.h> 28 #include <float.h> 29 #include <limits.h> 30 #include <errno.h> 31 #include <ctype.h> 32 33 #include "shgetc.h" 34 #include "floatscan.h" 35 36 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 37 38 #define LD_B1B_DIG 2 39 #define LD_B1B_MAX 9007199, 254740991 40 #define KMAX 128 41 42 #else /* LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 */ 43 44 #define LD_B1B_DIG 3 45 #define LD_B1B_MAX 18, 446744073, 709551615 46 #define KMAX 2048 47 48 #endif 49 50 #define MASK (KMAX-1) 51 52 #define CONCAT2(x,y) x ## y 53 #define CONCAT(x,y) CONCAT2(x,y) 54 55 // ANDROID: The original Musl sources use a FILE* handle, but redefine the 56 // type to match our custom fake_file_t from shgetc.h 57 #undef FILE 58 #define FILE struct fake_file_t 59 60 static long long scanexp(FILE *f, int pok) 61 { 62 int c; 63 int x; 64 long long y; 65 int neg = 0; 66 67 c = shgetc(f); 68 if (c=='+' || c=='-') { 69 neg = (c=='-'); 70 c = shgetc(f); 71 if (c-'0'>=10U && pok) shunget(f); 72 } 73 if (c-'0'>=10U) { 74 shunget(f); 75 return LLONG_MIN; 76 } 77 for (x=0; c-'0'<10U && x<INT_MAX/10; c = shgetc(f)) 78 x = 10*x + c-'0'; 79 for (y=x; c-'0'<10U && y<LLONG_MAX/100; c = shgetc(f)) 80 y = 10*y + c-'0'; 81 for (; c-'0'<10U; c = shgetc(f)); 82 shunget(f); 83 return neg ? -y : y; 84 } 85 86 87 static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int pok) 88 { 89 uint32_t x[KMAX]; 90 static const uint32_t th[] = { LD_B1B_MAX }; 91 int i, j, k, a, z; 92 long long lrp=0, dc=0; 93 long long e10=0; 94 int lnz = 0; 95 int gotdig = 0, gotrad = 0; 96 int rp; 97 int e2; 98 int emax = -emin-bits+3; 99 int denormal = 0; 100 long double y; 101 long double frac=0; 102 long double bias=0; 103 static const int p10s[] = { 10, 100, 1000, 10000, 104 100000, 1000000, 10000000, 100000000 }; 105 106 j=0; 107 k=0; 108 109 /* Don't let leading zeros consume buffer space */ 110 for (; c=='0'; c = shgetc(f)) gotdig=1; 111 if (c=='.') { 112 gotrad = 1; 113 for (c = shgetc(f); c=='0'; c = shgetc(f)) gotdig=1, lrp--; 114 } 115 116 x[0] = 0; 117 for (; c-'0'<10U || c=='.'; c = shgetc(f)) { 118 if (c == '.') { 119 if (gotrad) break; 120 gotrad = 1; 121 lrp = dc; 122 } else if (k < KMAX-3) { 123 dc++; 124 if (c!='0') lnz = dc; 125 if (j) x[k] = x[k]*10 + c-'0'; 126 else x[k] = c-'0'; 127 if (++j==9) { 128 k++; 129 j=0; 130 } 131 gotdig=1; 132 } else { 133 dc++; 134 if (c!='0') x[KMAX-4] |= 1; 135 } 136 } 137 if (!gotrad) lrp=dc; 138 139 if (gotdig && (c|32)=='e') { 140 e10 = scanexp(f, pok); 141 if (e10 == LLONG_MIN) { 142 if (pok) { 143 shunget(f); 144 } else { 145 shlim(f, 0); 146 return 0; 147 } 148 e10 = 0; 149 } 150 lrp += e10; 151 } else if (c>=0) { 152 shunget(f); 153 } 154 if (!gotdig) { 155 errno = EINVAL; 156 shlim(f, 0); 157 return 0; 158 } 159 160 /* Handle zero specially to avoid nasty special cases later */ 161 if (!x[0]) return sign * 0.0; 162 163 /* Optimize small integers (w/no exponent) and over/under-flow */ 164 if (lrp==dc && dc<10 && (bits>30 || x[0]>>bits==0)) 165 return sign * (long double)x[0]; 166 if (lrp > -emin/2) { 167 errno = ERANGE; 168 return sign * LDBL_MAX * LDBL_MAX; 169 } 170 if (lrp < emin-2*LDBL_MANT_DIG) { 171 errno = ERANGE; 172 return sign * LDBL_MIN * LDBL_MIN; 173 } 174 175 /* Align incomplete final B1B digit */ 176 if (j) { 177 for (; j<9; j++) x[k]*=10; 178 k++; 179 j=0; 180 } 181 182 a = 0; 183 z = k; 184 e2 = 0; 185 rp = lrp; 186 187 /* Optimize small to mid-size integers (even in exp. notation) */ 188 if (lnz<9 && lnz<=rp && rp < 18) { 189 if (rp == 9) return sign * (long double)x[0]; 190 if (rp < 9) return sign * (long double)x[0] / p10s[8-rp]; 191 int bitlim = bits-3*(int)(rp-9); 192 if (bitlim>30 || x[0]>>bitlim==0) 193 return sign * (long double)x[0] * p10s[rp-10]; 194 } 195 196 /* Align radix point to B1B digit boundary */ 197 if (rp % 9) { 198 int rpm9 = rp>=0 ? rp%9 : rp%9+9; 199 int p10 = p10s[8-rpm9]; 200 uint32_t carry = 0; 201 for (k=a; k!=z; k++) { 202 uint32_t tmp = x[k] % p10; 203 x[k] = x[k]/p10 + carry; 204 carry = 1000000000/p10 * tmp; 205 if (k==a && !x[k]) { 206 a = (a+1 & MASK); 207 rp -= 9; 208 } 209 } 210 if (carry) x[z++] = carry; 211 rp += 9-rpm9; 212 } 213 214 /* Upscale until desired number of bits are left of radix point */ 215 while (rp < 9*LD_B1B_DIG || (rp == 9*LD_B1B_DIG && x[a]<th[0])) { 216 uint32_t carry = 0; 217 e2 -= 29; 218 for (k=(z-1 & MASK); ; k=(k-1 & MASK)) { 219 uint64_t tmp = ((uint64_t)x[k] << 29) + carry; 220 if (tmp > 1000000000) { 221 carry = tmp / 1000000000; 222 x[k] = tmp % 1000000000; 223 } else { 224 carry = 0; 225 x[k] = tmp; 226 } 227 if (k==(z-1 & MASK) && k!=a && !x[k]) z = k; 228 if (k==a) break; 229 } 230 if (carry) { 231 rp += 9; 232 a = (a-1 & MASK); 233 if (a == z) { 234 z = (z-1 & MASK); 235 x[z-1 & MASK] |= x[z]; 236 } 237 x[a] = carry; 238 } 239 } 240 241 /* Downscale until exactly number of bits are left of radix point */ 242 for (;;) { 243 uint32_t carry = 0; 244 int sh = 1; 245 for (i=0; i<LD_B1B_DIG; i++) { 246 k = (a+i & MASK); 247 if (k == z || x[k] < th[i]) { 248 i=LD_B1B_DIG; 249 break; 250 } 251 if (x[a+i & MASK] > th[i]) break; 252 } 253 if (i==LD_B1B_DIG && rp==9*LD_B1B_DIG) break; 254 /* FIXME: find a way to compute optimal sh */ 255 if (rp > 9+9*LD_B1B_DIG) sh = 9; 256 e2 += sh; 257 for (k=a; k!=z; k=(k+1 & MASK)) { 258 uint32_t tmp = x[k] & (1<<sh)-1; 259 x[k] = (x[k]>>sh) + carry; 260 carry = (1000000000>>sh) * tmp; 261 if (k==a && !x[k]) { 262 a = (a+1 & MASK); 263 i--; 264 rp -= 9; 265 } 266 } 267 if (carry) { 268 if ((z+1 & MASK) != a) { 269 x[z] = carry; 270 z = (z+1 & MASK); 271 } else x[z-1 & MASK] |= 1; 272 } 273 } 274 275 /* Assemble desired bits into floating point variable */ 276 for (y=i=0; i<LD_B1B_DIG; i++) { 277 if ((a+i & MASK)==z) x[(z=(z+1 & MASK))-1] = 0; 278 y = 1000000000.0L * y + x[a+i & MASK]; 279 } 280 281 y *= sign; 282 283 /* Limit precision for denormal results */ 284 if (bits > LDBL_MANT_DIG+e2-emin) { 285 bits = LDBL_MANT_DIG+e2-emin; 286 if (bits<0) bits=0; 287 denormal = 1; 288 } 289 290 /* Calculate bias term to force rounding, move out lower bits */ 291 if (bits < LDBL_MANT_DIG) { 292 bias = copysignl(scalbn(1, 2*LDBL_MANT_DIG-bits-1), y); 293 frac = fmodl(y, scalbn(1, LDBL_MANT_DIG-bits)); 294 y -= frac; 295 y += bias; 296 } 297 298 /* Process tail of decimal input so it can affect rounding */ 299 if ((a+i & MASK) != z) { 300 uint32_t t = x[a+i & MASK]; 301 if (t < 500000000 && (t || (a+i+1 & MASK) != z)) 302 frac += 0.25*sign; 303 else if (t > 500000000) 304 frac += 0.75*sign; 305 else if (t == 500000000) { 306 if ((a+i+1 & MASK) == z) 307 frac += 0.5*sign; 308 else 309 frac += 0.75*sign; 310 } 311 if (LDBL_MANT_DIG-bits >= 2 && !fmodl(frac, 1)) 312 frac++; 313 } 314 315 y += frac; 316 y -= bias; 317 318 if ((e2+LDBL_MANT_DIG & INT_MAX) > emax-5) { 319 if (fabs(y) >= CONCAT(0x1p, LDBL_MANT_DIG)) { 320 if (denormal && bits==LDBL_MANT_DIG+e2-emin) 321 denormal = 0; 322 y *= 0.5; 323 e2++; 324 } 325 if (e2+LDBL_MANT_DIG>emax || (denormal && frac)) 326 errno = ERANGE; 327 } 328 329 return scalbnl(y, e2); 330 } 331 332 static long double hexfloat(FILE *f, int bits, int emin, int sign, int pok) 333 { 334 uint32_t x = 0; 335 long double y = 0; 336 long double scale = 1; 337 long double bias = 0; 338 int gottail = 0, gotrad = 0, gotdig = 0; 339 long long rp = 0; 340 long long dc = 0; 341 long long e2 = 0; 342 int d; 343 int c; 344 345 c = shgetc(f); 346 347 /* Skip leading zeros */ 348 for (; c=='0'; c = shgetc(f)) gotdig = 1; 349 350 if (c=='.') { 351 gotrad = 1; 352 c = shgetc(f); 353 /* Count zeros after the radix point before significand */ 354 for (rp=0; c=='0'; c = shgetc(f), rp--) gotdig = 1; 355 } 356 357 for (; c-'0'<10U || (c|32)-'a'<6U || c=='.'; c = shgetc(f)) { 358 if (c=='.') { 359 if (gotrad) break; 360 rp = dc; 361 gotrad = 1; 362 } else { 363 gotdig = 1; 364 if (c > '9') d = (c|32)+10-'a'; 365 else d = c-'0'; 366 if (dc<8) { 367 x = x*16 + d; 368 } else if (dc < LDBL_MANT_DIG/4+1) { 369 y += d*(scale/=16); 370 } else if (d && !gottail) { 371 y += 0.5*scale; 372 gottail = 1; 373 } 374 dc++; 375 } 376 } 377 if (!gotdig) { 378 shunget(f); 379 if (pok) { 380 shunget(f); 381 if (gotrad) shunget(f); 382 } else { 383 shlim(f, 0); 384 } 385 return sign * 0.0; 386 } 387 if (!gotrad) rp = dc; 388 while (dc<8) x *= 16, dc++; 389 if ((c|32)=='p') { 390 e2 = scanexp(f, pok); 391 if (e2 == LLONG_MIN) { 392 if (pok) { 393 shunget(f); 394 } else { 395 shlim(f, 0); 396 return 0; 397 } 398 e2 = 0; 399 } 400 } else { 401 shunget(f); 402 } 403 e2 += 4*rp - 32; 404 405 if (!x) return sign * 0.0; 406 if (e2 > -emin) { 407 errno = ERANGE; 408 return sign * LDBL_MAX * LDBL_MAX; 409 } 410 if (e2 < emin-2*LDBL_MANT_DIG) { 411 errno = ERANGE; 412 return sign * LDBL_MIN * LDBL_MIN; 413 } 414 415 while (x < 0x80000000) { 416 if (y>=0.5) { 417 x += x + 1; 418 y += y - 1; 419 } else { 420 x += x; 421 y += y; 422 } 423 e2--; 424 } 425 426 if (bits > 32+e2-emin) { 427 bits = 32+e2-emin; 428 if (bits<0) bits=0; 429 } 430 431 if (bits < LDBL_MANT_DIG) 432 bias = copysignl(scalbn(1, 32+LDBL_MANT_DIG-bits-1), sign); 433 434 if (bits<32 && y && !(x&1)) x++, y=0; 435 436 y = bias + sign*(long double)x + sign*y; 437 y -= bias; 438 439 if (!y) errno = ERANGE; 440 441 return scalbnl(y, e2); 442 } 443 444 long double __floatscan(FILE *f, int prec, int pok) 445 { 446 int sign = 1; 447 size_t i; 448 int bits; 449 int emin; 450 int c; 451 452 switch (prec) { 453 case 0: 454 bits = FLT_MANT_DIG; 455 emin = FLT_MIN_EXP-bits; 456 break; 457 case 1: 458 bits = DBL_MANT_DIG; 459 emin = DBL_MIN_EXP-bits; 460 break; 461 case 2: 462 bits = LDBL_MANT_DIG; 463 emin = LDBL_MIN_EXP-bits; 464 break; 465 default: 466 return 0; 467 } 468 469 while (isspace((c=shgetc(f)))); 470 471 if (c=='+' || c=='-') { 472 sign -= 2*(c=='-'); 473 c = shgetc(f); 474 } 475 476 for (i=0; i<8 && (c|32)=="infinity"[i]; i++) 477 if (i<7) c = shgetc(f); 478 if (i==3 || i==8 || (i>3 && pok)) { 479 if (i!=8) { 480 shunget(f); 481 if (pok) for (; i>3; i--) shunget(f); 482 } 483 return sign * INFINITY; 484 } 485 if (!i) for (i=0; i<3 && (c|32)=="nan"[i]; i++) 486 if (i<2) c = shgetc(f); 487 if (i==3) { 488 if (shgetc(f) != '(') { 489 shunget(f); 490 return NAN; 491 } 492 for (i=1; ; i++) { 493 c = shgetc(f); 494 if (c-'0'<10U || c-'A'<26U || c-'a'<26U || c=='_') 495 continue; 496 if (c==')') return NAN; 497 shunget(f); 498 if (!pok) { 499 errno = EINVAL; 500 shlim(f, 0); 501 return 0; 502 } 503 while (i--) shunget(f); 504 return NAN; 505 } 506 return NAN; 507 } 508 509 if (i) { 510 shunget(f); 511 errno = EINVAL; 512 shlim(f, 0); 513 return 0; 514 } 515 516 if (c=='0') { 517 c = shgetc(f); 518 if ((c|32) == 'x') 519 return hexfloat(f, bits, emin, sign, pok); 520 shunget(f); 521 c = '0'; 522 } 523 524 return decfloat(f, c, bits, emin, sign, pok); 525 } 526
