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