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