| ... | @@ -0,0 +1,1038 @@ |
| 1 | // Code ported from musl libc 8f12c4e110acb3bbbdc8abfb3a552c3ced718039 |
| 2 | // and then modified to use softfloat and to assume f128 for everything |
| 3 | |
| 4 | #include "parse_f128.h" |
| 5 | #include "softfloat.h" |
| 6 | #include <stddef.h> |
| 7 | #include <sys/types.h> |
| 8 | #include <errno.h> |
| 9 | #include <limits.h> |
| 10 | #include <string.h> |
| 11 | #include <math.h> |
| 12 | |
| 13 | #define shcnt(f) ((f)->shcnt + ((f)->rpos - (f)->buf)) |
| 14 | #define shlim(f, lim) __shlim((f), (lim)) |
| 15 | #define shgetc(f) (((f)->rpos != (f)->shend) ? *(f)->rpos++ : __shgetc(f)) |
| 16 | #define shunget(f) ((f)->shlim>=0 ? (void)(f)->rpos-- : (void)0) |
| 17 | |
| 18 | #define sh_fromstring(f, s) \ |
| 19 | ((f)->buf = (f)->rpos = (void *)(s), (f)->rend = (void*)-1) |
| 20 | |
| 21 | #define LD_B1B_DIG 4 |
| 22 | #define LD_B1B_MAX 10384593, 717069655, 257060992, 658440191 |
| 23 | #define KMAX 2048 |
| 24 | |
| 25 | #define MASK (KMAX-1) |
| 26 | |
| 27 | #define CONCAT2(x,y) x ## y |
| 28 | #define CONCAT(x,y) CONCAT2(x,y) |
| 29 | |
| 30 | #define F_PERM 1 |
| 31 | #define F_NORD 4 |
| 32 | #define F_NOWR 8 |
| 33 | #define F_EOF 16 |
| 34 | #define F_ERR 32 |
| 35 | #define F_SVB 64 |
| 36 | #define F_APP 128 |
| 37 | |
| 38 | #define EOF (-1) |
| 39 | |
| 40 | #define LDBL_MANT_DIG 113 |
| 41 | #define LDBL_MIN_EXP (-16381) |
| 42 | #define LDBL_MAX_EXP 16384 |
| 43 | |
| 44 | #define LDBL_DIG 33 |
| 45 | #define LDBL_MIN_10_EXP (-4931) |
| 46 | #define LDBL_MAX_10_EXP 4932 |
| 47 | |
| 48 | #define DECIMAL_DIG 36 |
| 49 | |
| 50 | |
| 51 | #if __BYTE_ORDER == __LITTLE_ENDIAN |
| 52 | union ldshape { |
| 53 | float128_t f; |
| 54 | struct { |
| 55 | uint64_t lo; |
| 56 | uint32_t mid; |
| 57 | uint16_t top; |
| 58 | uint16_t se; |
| 59 | } i; |
| 60 | struct { |
| 61 | uint64_t lo; |
| 62 | uint64_t hi; |
| 63 | } i2; |
| 64 | }; |
| 65 | #elif __BYTE_ORDER == __BIG_ENDIAN |
| 66 | union ldshape { |
| 67 | float128_t f; |
| 68 | struct { |
| 69 | uint16_t se; |
| 70 | uint16_t top; |
| 71 | uint32_t mid; |
| 72 | uint64_t lo; |
| 73 | } i; |
| 74 | struct { |
| 75 | uint64_t hi; |
| 76 | uint64_t lo; |
| 77 | } i2; |
| 78 | }; |
| 79 | #error Unsupported endian |
| 80 | #endif |
| 81 | |
| 82 | struct MuslFILE { |
| 83 | unsigned flags; |
| 84 | unsigned char *rpos, *rend; |
| 85 | int (*close)(struct MuslFILE *); |
| 86 | unsigned char *wend, *wpos; |
| 87 | unsigned char *mustbezero_1; |
| 88 | unsigned char *wbase; |
| 89 | size_t (*read)(struct MuslFILE *, unsigned char *, size_t); |
| 90 | size_t (*write)(struct MuslFILE *, const unsigned char *, size_t); |
| 91 | off_t (*seek)(struct MuslFILE *, off_t, int); |
| 92 | unsigned char *buf; |
| 93 | size_t buf_size; |
| 94 | struct MuslFILE *prev, *next; |
| 95 | int fd; |
| 96 | int pipe_pid; |
| 97 | long lockcount; |
| 98 | int mode; |
| 99 | volatile int lock; |
| 100 | int lbf; |
| 101 | void *cookie; |
| 102 | off_t off; |
| 103 | char *getln_buf; |
| 104 | void *mustbezero_2; |
| 105 | unsigned char *shend; |
| 106 | off_t shlim, shcnt; |
| 107 | struct MuslFILE *prev_locked, *next_locked; |
| 108 | struct __locale_struct *locale; |
| 109 | }; |
| 110 | |
| 111 | static void __shlim(struct MuslFILE *f, off_t lim) |
| 112 | { |
| 113 | f->shlim = lim; |
| 114 | f->shcnt = f->buf - f->rpos; |
| 115 | /* If lim is nonzero, rend must be a valid pointer. */ |
| 116 | if (lim && f->rend - f->rpos > lim) |
| 117 | f->shend = f->rpos + lim; |
| 118 | else |
| 119 | f->shend = f->rend; |
| 120 | } |
| 121 | |
| 122 | static int __toread(struct MuslFILE *f) |
| 123 | { |
| 124 | f->mode |= f->mode-1; |
| 125 | if (f->wpos != f->wbase) f->write(f, 0, 0); |
| 126 | f->wpos = f->wbase = f->wend = 0; |
| 127 | if (f->flags & F_NORD) { |
| 128 | f->flags |= F_ERR; |
| 129 | return EOF; |
| 130 | } |
| 131 | f->rpos = f->rend = f->buf + f->buf_size; |
| 132 | return (f->flags & F_EOF) ? EOF : 0; |
| 133 | } |
| 134 | |
| 135 | static int __uflow(struct MuslFILE *f) |
| 136 | { |
| 137 | unsigned char c; |
| 138 | if (!__toread(f) && f->read(f, &c, 1)==1) return c; |
| 139 | return EOF; |
| 140 | } |
| 141 | |
| 142 | static int __shgetc(struct MuslFILE *f) |
| 143 | { |
| 144 | int c; |
| 145 | off_t cnt = shcnt(f); |
| 146 | if (f->shlim && cnt >= f->shlim || (c=__uflow(f)) < 0) { |
| 147 | f->shcnt = f->buf - f->rpos + cnt; |
| 148 | f->shend = f->rpos; |
| 149 | f->shlim = -1; |
| 150 | return EOF; |
| 151 | } |
| 152 | cnt++; |
| 153 | if (f->shlim && f->rend - f->rpos > f->shlim - cnt) |
| 154 | f->shend = f->rpos + (f->shlim - cnt); |
| 155 | else |
| 156 | f->shend = f->rend; |
| 157 | f->shcnt = f->buf - f->rpos + cnt; |
| 158 | if (f->rpos[-1] != c) f->rpos[-1] = c; |
| 159 | return c; |
| 160 | } |
| 161 | |
| 162 | static long long scanexp(struct MuslFILE *f, int pok) |
| 163 | { |
| 164 | int c; |
| 165 | int x; |
| 166 | long long y; |
| 167 | int neg = 0; |
| 168 | |
| 169 | c = shgetc(f); |
| 170 | if (c=='+' || c=='-') { |
| 171 | neg = (c=='-'); |
| 172 | c = shgetc(f); |
| 173 | if (c-'0'>=10U && pok) shunget(f); |
| 174 | } |
| 175 | if (c-'0'>=10U) { |
| 176 | shunget(f); |
| 177 | return LLONG_MIN; |
| 178 | } |
| 179 | for (x=0; c-'0'<10U && x<INT_MAX/10; c = shgetc(f)) |
| 180 | x = 10*x + c-'0'; |
| 181 | for (y=x; c-'0'<10U && y<LLONG_MAX/100; c = shgetc(f)) |
| 182 | y = 10*y + c-'0'; |
| 183 | for (; c-'0'<10U; c = shgetc(f)); |
| 184 | shunget(f); |
| 185 | return neg ? -y : y; |
| 186 | } |
| 187 | |
| 188 | static float128_t copysignf128(float128_t x, float128_t y) |
| 189 | { |
| 190 | union ldshape ux = {x}, uy = {y}; |
| 191 | ux.i.se &= 0x7fff; |
| 192 | ux.i.se |= uy.i.se & 0x8000; |
| 193 | return ux.f; |
| 194 | } |
| 195 | |
| 196 | static void mul_eq_f128_float(float128_t *x, float op_float) { |
| 197 | //x *= 0x1p120f; |
| 198 | float32_t op_f32; |
| 199 | memcpy(&op_f32, &op_float, sizeof(float)); |
| 200 | float128_t op_f128; |
| 201 | f32_to_f128M(op_f32, &op_f128); |
| 202 | float128_t new_value; |
| 203 | f128M_mul(x, &op_f128, &new_value); |
| 204 | *x = new_value; |
| 205 | } |
| 206 | |
| 207 | static float128_t dbl_to_f128(double x) { |
| 208 | float64_t x_f64; |
| 209 | memcpy(&x_f64, &x, sizeof(double)); |
| 210 | float128_t result; |
| 211 | f64_to_f128M(x_f64, &result); |
| 212 | return result; |
| 213 | } |
| 214 | |
| 215 | static float128_t fmodf128(float128_t x, float128_t y) |
| 216 | { |
| 217 | union ldshape ux = {x}, uy = {y}; |
| 218 | int ex = ux.i.se & 0x7fff; |
| 219 | int ey = uy.i.se & 0x7fff; |
| 220 | int sx = ux.i.se & 0x8000; |
| 221 | |
| 222 | float128_t zero; |
| 223 | ui32_to_f128M(0, &zero); |
| 224 | // if (y == 0 || isnan(y) || ex == 0x7fff) |
| 225 | if (f128M_eq(&y, &zero) || f128M_isSignalingNaN(&y) || ex == 0x7fff) { |
| 226 | //return (x*y)/(x*y); |
| 227 | float128_t x_times_y; |
| 228 | f128M_mul(&x, &y, &x_times_y); |
| 229 | float128_t result; |
| 230 | f128M_div(&x_times_y, &x_times_y, &result); |
| 231 | return result; |
| 232 | } |
| 233 | ux.i.se = ex; |
| 234 | uy.i.se = ey; |
| 235 | //if (ux.f <= uy.f) { |
| 236 | if (f128M_le(&ux.f, &uy.f)) { |
| 237 | //if (ux.f == uy.f) { |
| 238 | if (f128M_eq(&ux.f, &uy.f)) { |
| 239 | //return 0*x; |
| 240 | float128_t result; |
| 241 | f128M_mul(&zero, &x, &result); |
| 242 | return result; |
| 243 | } |
| 244 | return x; |
| 245 | } |
| 246 | |
| 247 | /* normalize x and y */ |
| 248 | if (!ex) { |
| 249 | //ux.f *= 0x1p120f; |
| 250 | mul_eq_f128_float(&ux.f, 0x1p120f); |
| 251 | |
| 252 | ex = ux.i.se - 120; |
| 253 | } |
| 254 | if (!ey) { |
| 255 | //uy.f *= 0x1p120f; |
| 256 | mul_eq_f128_float(&uy.f, 0x1p120f); |
| 257 | |
| 258 | ey = uy.i.se - 120; |
| 259 | } |
| 260 | |
| 261 | /* x mod y */ |
| 262 | uint64_t hi, lo, xhi, xlo, yhi, ylo; |
| 263 | xhi = (ux.i2.hi & -1ULL>>16) | 1ULL<<48; |
| 264 | yhi = (uy.i2.hi & -1ULL>>16) | 1ULL<<48; |
| 265 | xlo = ux.i2.lo; |
| 266 | ylo = uy.i2.lo; |
| 267 | for (; ex > ey; ex--) { |
| 268 | hi = xhi - yhi; |
| 269 | lo = xlo - ylo; |
| 270 | if (xlo < ylo) |
| 271 | hi -= 1; |
| 272 | if (hi >> 63 == 0) { |
| 273 | if ((hi|lo) == 0) { |
| 274 | //return 0*x; |
| 275 | float128_t result; |
| 276 | f128M_mul(&zero, &x, &result); |
| 277 | return result; |
| 278 | } |
| 279 | xhi = 2*hi + (lo>>63); |
| 280 | xlo = 2*lo; |
| 281 | } else { |
| 282 | xhi = 2*xhi + (xlo>>63); |
| 283 | xlo = 2*xlo; |
| 284 | } |
| 285 | } |
| 286 | hi = xhi - yhi; |
| 287 | lo = xlo - ylo; |
| 288 | if (xlo < ylo) |
| 289 | hi -= 1; |
| 290 | if (hi >> 63 == 0) { |
| 291 | if ((hi|lo) == 0) { |
| 292 | //return 0*x; |
| 293 | float128_t result; |
| 294 | f128M_mul(&zero, &x, &result); |
| 295 | return result; |
| 296 | } |
| 297 | xhi = hi; |
| 298 | xlo = lo; |
| 299 | } |
| 300 | for (; xhi >> 48 == 0; xhi = 2*xhi + (xlo>>63), xlo = 2*xlo, ex--); |
| 301 | ux.i2.hi = xhi; |
| 302 | ux.i2.lo = xlo; |
| 303 | |
| 304 | /* scale result */ |
| 305 | if (ex <= 0) { |
| 306 | ux.i.se = (ex+120)|sx; |
| 307 | //ux.f *= 0x1p-120f; |
| 308 | mul_eq_f128_float(&ux.f, 0x1p-120f); |
| 309 | } else |
| 310 | ux.i.se = ex|sx; |
| 311 | return ux.f; |
| 312 | } |
| 313 | |
| 314 | static float128_t int_mul_f128_cast_u32(int sign, uint32_t x0) { |
| 315 | float128_t x0_f128; |
| 316 | ui32_to_f128M(x0, &x0_f128); |
| 317 | float128_t sign_f128; |
| 318 | i32_to_f128M(sign, &sign_f128); |
| 319 | float128_t result; |
| 320 | f128M_mul(&sign_f128, &x0_f128, &result); |
| 321 | return result; |
| 322 | } |
| 323 | |
| 324 | static float128_t triple_divide(int sign, uint32_t x0, int p10s) { |
| 325 | float128_t part1 = int_mul_f128_cast_u32(sign, x0); |
| 326 | float128_t p10s_f128; |
| 327 | i32_to_f128M(p10s, &p10s_f128); |
| 328 | float128_t result; |
| 329 | f128M_div(&part1, &p10s_f128, &result); |
| 330 | return result; |
| 331 | } |
| 332 | |
| 333 | static float128_t triple_multiply(int sign, uint32_t x0, int p10s) { |
| 334 | float128_t part1 = int_mul_f128_cast_u32(sign, x0); |
| 335 | float128_t p10s_f128; |
| 336 | i32_to_f128M(p10s, &p10s_f128); |
| 337 | float128_t result; |
| 338 | f128M_mul(&part1, &p10s_f128, &result); |
| 339 | return result; |
| 340 | } |
| 341 | |
| 342 | static void mul_eq_f128_int(float128_t *y, int sign) { |
| 343 | float128_t sign_f128; |
| 344 | i32_to_f128M(sign, &sign_f128); |
| 345 | float128_t new_value; |
| 346 | f128M_mul(y, &sign_f128, &new_value); |
| 347 | *y = new_value; |
| 348 | } |
| 349 | |
| 350 | static float128_t literal_f128(__float128 x) { |
| 351 | float128_t result; |
| 352 | memcpy(&result, &x, 16); |
| 353 | return result; |
| 354 | } |
| 355 | |
| 356 | static void mul_eq_f128_f128(float128_t *a, float128_t b) { |
| 357 | float128_t new_value; |
| 358 | f128M_mul(a, &b, &new_value); |
| 359 | *a = new_value; |
| 360 | } |
| 361 | |
| 362 | static void add_eq_f128_dbl(float128_t *a, double b) { |
| 363 | float64_t b_f64; |
| 364 | memcpy(&b_f64, &b, sizeof(double)); |
| 365 | |
| 366 | float128_t b_f128; |
| 367 | f64_to_f128M(b_f64, &b_f128); |
| 368 | |
| 369 | float128_t new_value; |
| 370 | f128M_add(a, &b_f128, &new_value); |
| 371 | *a = new_value; |
| 372 | } |
| 373 | |
| 374 | static float128_t scalbnf128(float128_t x, int n) |
| 375 | { |
| 376 | union ldshape u; |
| 377 | |
| 378 | if (n > 16383) { |
| 379 | //x *= 0x1p16383q; |
| 380 | mul_eq_f128_f128(&x, literal_f128(0x1p16383q)); |
| 381 | n -= 16383; |
| 382 | if (n > 16383) { |
| 383 | //x *= 0x1p16383q; |
| 384 | mul_eq_f128_f128(&x, literal_f128(0x1p16383q)); |
| 385 | n -= 16383; |
| 386 | if (n > 16383) |
| 387 | n = 16383; |
| 388 | } |
| 389 | } else if (n < -16382) { |
| 390 | //x *= 0x1p-16382q * 0x1p113q; |
| 391 | { |
| 392 | float128_t mul_result; |
| 393 | float128_t a = literal_f128(0x1p-16382q); |
| 394 | float128_t b = literal_f128(0x1p113q); |
| 395 | f128M_mul(&a, &b, &mul_result); |
| 396 | mul_eq_f128_f128(&x, mul_result); |
| 397 | } |
| 398 | n += 16382 - 113; |
| 399 | if (n < -16382) { |
| 400 | //x *= 0x1p-16382q * 0x1p113q; |
| 401 | { |
| 402 | float128_t mul_result; |
| 403 | float128_t a = literal_f128(0x1p-16382q); |
| 404 | float128_t b = literal_f128(0x1p113q); |
| 405 | f128M_mul(&a, &b, &mul_result); |
| 406 | mul_eq_f128_f128(&x, mul_result); |
| 407 | } |
| 408 | n += 16382 - 113; |
| 409 | if (n < -16382) |
| 410 | n = -16382; |
| 411 | } |
| 412 | } |
| 413 | //u.f = 1.0; |
| 414 | ui32_to_f128M(1, &u.f); |
| 415 | u.i.se = 0x3fff + n; |
| 416 | mul_eq_f128_f128(&x, u.f); |
| 417 | return x; |
| 418 | } |
| 419 | |
| 420 | static float128_t fabsf128(float128_t x) |
| 421 | { |
| 422 | union ldshape u = {x}; |
| 423 | |
| 424 | u.i.se &= 0x7fff; |
| 425 | return u.f; |
| 426 | } |
| 427 | |
| 428 | static float128_t decfloat(struct MuslFILE *f, int c, int bits, int emin, int sign, int pok) |
| 429 | { |
| 430 | uint32_t x[KMAX]; |
| 431 | static const uint32_t th[] = { LD_B1B_MAX }; |
| 432 | int i, j, k, a, z; |
| 433 | long long lrp=0, dc=0; |
| 434 | long long e10=0; |
| 435 | int lnz = 0; |
| 436 | int gotdig = 0, gotrad = 0; |
| 437 | int rp; |
| 438 | int e2; |
| 439 | int emax = -emin-bits+3; |
| 440 | int denormal = 0; |
| 441 | float128_t y; |
| 442 | float128_t zero; |
| 443 | ui32_to_f128M(0, &zero); |
| 444 | float128_t frac=zero; |
| 445 | float128_t bias=zero; |
| 446 | static const int p10s[] = { 10, 100, 1000, 10000, |
| 447 | 100000, 1000000, 10000000, 100000000 }; |
| 448 | |
| 449 | j=0; |
| 450 | k=0; |
| 451 | |
| 452 | /* Don't let leading zeros consume buffer space */ |
| 453 | for (; c=='0'; c = shgetc(f)) gotdig=1; |
| 454 | if (c=='.') { |
| 455 | gotrad = 1; |
| 456 | for (c = shgetc(f); c=='0'; c = shgetc(f)) gotdig=1, lrp--; |
| 457 | } |
| 458 | |
| 459 | x[0] = 0; |
| 460 | for (; c-'0'<10U || c=='.'; c = shgetc(f)) { |
| 461 | if (c == '.') { |
| 462 | if (gotrad) break; |
| 463 | gotrad = 1; |
| 464 | lrp = dc; |
| 465 | } else if (k < KMAX-3) { |
| 466 | dc++; |
| 467 | if (c!='0') lnz = dc; |
| 468 | if (j) x[k] = x[k]*10 + c-'0'; |
| 469 | else x[k] = c-'0'; |
| 470 | if (++j==9) { |
| 471 | k++; |
| 472 | j=0; |
| 473 | } |
| 474 | gotdig=1; |
| 475 | } else { |
| 476 | dc++; |
| 477 | if (c!='0') { |
| 478 | lnz = (KMAX-4)*9; |
| 479 | x[KMAX-4] |= 1; |
| 480 | } |
| 481 | } |
| 482 | } |
| 483 | if (!gotrad) lrp=dc; |
| 484 | |
| 485 | if (gotdig && (c|32)=='e') { |
| 486 | e10 = scanexp(f, pok); |
| 487 | if (e10 == LLONG_MIN) { |
| 488 | if (pok) { |
| 489 | shunget(f); |
| 490 | } else { |
| 491 | shlim(f, 0); |
| 492 | return zero; |
| 493 | } |
| 494 | e10 = 0; |
| 495 | } |
| 496 | lrp += e10; |
| 497 | } else if (c>=0) { |
| 498 | shunget(f); |
| 499 | } |
| 500 | if (!gotdig) { |
| 501 | errno = EINVAL; |
| 502 | shlim(f, 0); |
| 503 | return zero; |
| 504 | } |
| 505 | |
| 506 | /* Handle zero specially to avoid nasty special cases later */ |
| 507 | if (!x[0]) { |
| 508 | //return sign * 0.0; |
| 509 | return dbl_to_f128(sign * 0.0); |
| 510 | } |
| 511 | |
| 512 | /* Optimize small integers (w/no exponent) and over/under-flow */ |
| 513 | if (lrp==dc && dc<10 && (bits>30 || x[0]>>bits==0)) { |
| 514 | //return sign * (float128_t)x[0]; |
| 515 | float128_t sign_f128; |
| 516 | i32_to_f128M(sign, &sign_f128); |
| 517 | float128_t x0_f128; |
| 518 | ui32_to_f128M(x[0], &x0_f128); |
| 519 | float128_t result; |
| 520 | f128M_mul(&sign_f128, &x0_f128, &result); |
| 521 | return result; |
| 522 | } |
| 523 | if (lrp > -emin/2) { |
| 524 | errno = ERANGE; |
| 525 | //return sign * LDBL_MAX * LDBL_MAX; |
| 526 | return zero; |
| 527 | } |
| 528 | if (lrp < emin-2*LDBL_MANT_DIG) { |
| 529 | errno = ERANGE; |
| 530 | //return sign * LDBL_MIN * LDBL_MIN; |
| 531 | return zero; |
| 532 | } |
| 533 | |
| 534 | /* Align incomplete final B1B digit */ |
| 535 | if (j) { |
| 536 | for (; j<9; j++) x[k]*=10; |
| 537 | k++; |
| 538 | j=0; |
| 539 | } |
| 540 | |
| 541 | a = 0; |
| 542 | z = k; |
| 543 | e2 = 0; |
| 544 | rp = lrp; |
| 545 | |
| 546 | /* Optimize small to mid-size integers (even in exp. notation) */ |
| 547 | if (lnz<9 && lnz<=rp && rp < 18) { |
| 548 | if (rp == 9) { |
| 549 | //return sign * (float128_t)(x[0]); |
| 550 | return int_mul_f128_cast_u32(sign, x[0]); |
| 551 | } |
| 552 | if (rp < 9) { |
| 553 | //return sign * (float128_t)(x[0]) / p10s[8-rp]; |
| 554 | return triple_divide(sign, x[0], p10s[8-rp]); |
| 555 | } |
| 556 | int bitlim = bits-3*(int)(rp-9); |
| 557 | if (bitlim>30 || x[0]>>bitlim==0) |
| 558 | //return sign * (float128_t)(x[0]) * p10s[rp-10]; |
| 559 | return triple_multiply(sign, x[0], p10s[rp-10]); |
| 560 | } |
| 561 | |
| 562 | /* Drop trailing zeros */ |
| 563 | for (; !x[z-1]; z--); |
| 564 | |
| 565 | /* Align radix point to B1B digit boundary */ |
| 566 | if (rp % 9) { |
| 567 | int rpm9 = rp>=0 ? rp%9 : rp%9+9; |
| 568 | int p10 = p10s[8-rpm9]; |
| 569 | uint32_t carry = 0; |
| 570 | for (k=a; k!=z; k++) { |
| 571 | uint32_t tmp = x[k] % p10; |
| 572 | x[k] = x[k]/p10 + carry; |
| 573 | carry = 1000000000/p10 * tmp; |
| 574 | if (k==a && !x[k]) { |
| 575 | a = (a+1 & MASK); |
| 576 | rp -= 9; |
| 577 | } |
| 578 | } |
| 579 | if (carry) x[z++] = carry; |
| 580 | rp += 9-rpm9; |
| 581 | } |
| 582 | |
| 583 | /* Upscale until desired number of bits are left of radix point */ |
| 584 | while (rp < 9*LD_B1B_DIG || (rp == 9*LD_B1B_DIG && x[a]<th[0])) { |
| 585 | uint32_t carry = 0; |
| 586 | e2 -= 29; |
| 587 | for (k=(z-1 & MASK); ; k=(k-1 & MASK)) { |
| 588 | uint64_t tmp = ((uint64_t)x[k] << 29) + carry; |
| 589 | if (tmp > 1000000000) { |
| 590 | carry = tmp / 1000000000; |
| 591 | x[k] = tmp % 1000000000; |
| 592 | } else { |
| 593 | carry = 0; |
| 594 | x[k] = tmp; |
| 595 | } |
| 596 | if (k==(z-1 & MASK) && k!=a && !x[k]) z = k; |
| 597 | if (k==a) break; |
| 598 | } |
| 599 | if (carry) { |
| 600 | rp += 9; |
| 601 | a = (a-1 & MASK); |
| 602 | if (a == z) { |
| 603 | z = (z-1 & MASK); |
| 604 | x[z-1 & MASK] |= x[z]; |
| 605 | } |
| 606 | x[a] = carry; |
| 607 | } |
| 608 | } |
| 609 | |
| 610 | /* Downscale until exactly number of bits are left of radix point */ |
| 611 | for (;;) { |
| 612 | uint32_t carry = 0; |
| 613 | int sh = 1; |
| 614 | for (i=0; i<LD_B1B_DIG; i++) { |
| 615 | k = (a+i & MASK); |
| 616 | if (k == z || x[k] < th[i]) { |
| 617 | i=LD_B1B_DIG; |
| 618 | break; |
| 619 | } |
| 620 | if (x[a+i & MASK] > th[i]) break; |
| 621 | } |
| 622 | if (i==LD_B1B_DIG && rp==9*LD_B1B_DIG) break; |
| 623 | /* FIXME: find a way to compute optimal sh */ |
| 624 | if (rp > 9+9*LD_B1B_DIG) sh = 9; |
| 625 | e2 += sh; |
| 626 | for (k=a; k!=z; k=(k+1 & MASK)) { |
| 627 | uint32_t tmp = x[k] & (1<<sh)-1; |
| 628 | x[k] = (x[k]>>sh) + carry; |
| 629 | carry = (1000000000>>sh) * tmp; |
| 630 | if (k==a && !x[k]) { |
| 631 | a = (a+1 & MASK); |
| 632 | i--; |
| 633 | rp -= 9; |
| 634 | } |
| 635 | } |
| 636 | if (carry) { |
| 637 | if ((z+1 & MASK) != a) { |
| 638 | x[z] = carry; |
| 639 | z = (z+1 & MASK); |
| 640 | } else x[z-1 & MASK] |= 1; |
| 641 | } |
| 642 | } |
| 643 | |
| 644 | /* Assemble desired bits into floating point variable */ |
| 645 | for (y=zero,i=0; i<LD_B1B_DIG; i++) { |
| 646 | if ((a+i & MASK)==z) x[(z=(z+1 & MASK))-1] = 0; |
| 647 | //y = 1000000000.0L * y + x[a+i & MASK]; |
| 648 | float128_t const_f128; |
| 649 | ui64_to_f128M(1000000000, &const_f128); |
| 650 | float128_t mul_y; |
| 651 | f128M_mul(&const_f128, &y, &mul_y); |
| 652 | float128_t x_f128; |
| 653 | ui32_to_f128M(x[a+i & MASK], &x_f128); |
| 654 | f128M_add(&mul_y, &x_f128, &y); |
| 655 | } |
| 656 | |
| 657 | //y *= sign; |
| 658 | mul_eq_f128_int(&y, sign); |
| 659 | |
| 660 | /* Limit precision for denormal results */ |
| 661 | if (bits > LDBL_MANT_DIG+e2-emin) { |
| 662 | bits = LDBL_MANT_DIG+e2-emin; |
| 663 | if (bits<0) bits=0; |
| 664 | denormal = 1; |
| 665 | } |
| 666 | |
| 667 | /* Calculate bias term to force rounding, move out lower bits */ |
| 668 | if (bits < LDBL_MANT_DIG) { |
| 669 | bias = copysignf128(dbl_to_f128(scalbn(1, 2*LDBL_MANT_DIG-bits-1)), y); |
| 670 | frac = fmodf128(y, dbl_to_f128(scalbn(1, LDBL_MANT_DIG-bits))); |
| 671 | //y -= frac; |
| 672 | { |
| 673 | float128_t new_value; |
| 674 | f128M_sub(&y, &frac, &new_value); |
| 675 | y = new_value; |
| 676 | } |
| 677 | //y += bias; |
| 678 | { |
| 679 | float128_t new_value; |
| 680 | f128M_add(&y, &frac, &new_value); |
| 681 | y = new_value; |
| 682 | } |
| 683 | } |
| 684 | |
| 685 | /* Process tail of decimal input so it can affect rounding */ |
| 686 | if ((a+i & MASK) != z) { |
| 687 | uint32_t t = x[a+i & MASK]; |
| 688 | if (t < 500000000 && (t || (a+i+1 & MASK) != z)) { |
| 689 | //frac += 0.25*sign; |
| 690 | add_eq_f128_dbl(&frac, 0.25*sign); |
| 691 | } else if (t > 500000000) { |
| 692 | //frac += 0.75*sign; |
| 693 | add_eq_f128_dbl(&frac, 0.75*sign); |
| 694 | } else if (t == 500000000) { |
| 695 | if ((a+i+1 & MASK) == z) { |
| 696 | //frac += 0.5*sign; |
| 697 | add_eq_f128_dbl(&frac, 0.5*sign); |
| 698 | } else { |
| 699 | //frac += 0.75*sign; |
| 700 | add_eq_f128_dbl(&frac, 0.75*sign); |
| 701 | } |
| 702 | } |
| 703 | //if (LDBL_MANT_DIG-bits >= 2 && !fmodf128(frac, 1)) |
| 704 | if (LDBL_MANT_DIG-bits >= 2) { |
| 705 | float128_t one; |
| 706 | ui32_to_f128M(1, &one); |
| 707 | float128_t mod_result = fmodf128(frac, one); |
| 708 | if (f128M_eq(&mod_result, &zero)) { |
| 709 | //frac++; |
| 710 | add_eq_f128_dbl(&frac, 1.0); |
| 711 | } |
| 712 | } |
| 713 | } |
| 714 | |
| 715 | //y += frac; |
| 716 | { |
| 717 | float128_t new_value; |
| 718 | f128M_add(&y, &frac, &new_value); |
| 719 | y = new_value; |
| 720 | } |
| 721 | //y -= bias; |
| 722 | { |
| 723 | float128_t new_value; |
| 724 | f128M_sub(&y, &bias, &new_value); |
| 725 | y = new_value; |
| 726 | } |
| 727 | |
| 728 | if ((e2+LDBL_MANT_DIG & INT_MAX) > emax-5) { |
| 729 | //if (fabsf128(y) >= 0x1p113) |
| 730 | float128_t abs_y = fabsf128(y); |
| 731 | float128_t mant_f128 = literal_f128(0x1p113q); |
| 732 | if (!f128M_lt(&abs_y, &mant_f128)) { |
| 733 | if (denormal && bits==LDBL_MANT_DIG+e2-emin) |
| 734 | denormal = 0; |
| 735 | //y *= 0.5; |
| 736 | { |
| 737 | float128_t point_5 = dbl_to_f128(0.5); |
| 738 | float128_t new_value; |
| 739 | f128M_mul(&y, &point_5, &new_value); |
| 740 | y = new_value; |
| 741 | } |
| 742 | |
| 743 | e2++; |
| 744 | } |
| 745 | if (e2+LDBL_MANT_DIG>emax || (denormal && !f128M_eq(&frac, &zero))) |
| 746 | errno = ERANGE; |
| 747 | } |
| 748 | |
| 749 | return scalbnf128(y, e2); |
| 750 | } |
| 751 | |
| 752 | static float128_t hexfloat(struct MuslFILE *f, int bits, int emin, int sign, int pok) |
| 753 | { |
| 754 | float128_t zero; |
| 755 | ui32_to_f128M(0, &zero); |
| 756 | float128_t one; |
| 757 | ui32_to_f128M(1, &one); |
| 758 | float128_t sixteen; |
| 759 | ui32_to_f128M(16, &sixteen); |
| 760 | float128_t point_5 = dbl_to_f128(0.5); |
| 761 | |
| 762 | uint32_t x = 0; |
| 763 | float128_t y = zero; |
| 764 | float128_t scale = one; |
| 765 | float128_t bias = zero; |
| 766 | int gottail = 0, gotrad = 0, gotdig = 0; |
| 767 | long long rp = 0; |
| 768 | long long dc = 0; |
| 769 | long long e2 = 0; |
| 770 | int d; |
| 771 | int c; |
| 772 | |
| 773 | c = shgetc(f); |
| 774 | |
| 775 | /* Skip leading zeros */ |
| 776 | for (; c=='0'; c = shgetc(f)) gotdig = 1; |
| 777 | |
| 778 | if (c=='.') { |
| 779 | gotrad = 1; |
| 780 | c = shgetc(f); |
| 781 | /* Count zeros after the radix point before significand */ |
| 782 | for (rp=0; c=='0'; c = shgetc(f), rp--) gotdig = 1; |
| 783 | } |
| 784 | |
| 785 | for (; c-'0'<10U || (c|32)-'a'<6U || c=='.'; c = shgetc(f)) { |
| 786 | if (c=='.') { |
| 787 | if (gotrad) break; |
| 788 | rp = dc; |
| 789 | gotrad = 1; |
| 790 | } else { |
| 791 | gotdig = 1; |
| 792 | if (c > '9') d = (c|32)+10-'a'; |
| 793 | else d = c-'0'; |
| 794 | if (dc<8) { |
| 795 | x = x*16 + d; |
| 796 | } else if (dc < LDBL_MANT_DIG/4+1) { |
| 797 | //y += d*(scale/=16); |
| 798 | { |
| 799 | float128_t divided; |
| 800 | f128M_div(&scale, &sixteen, ÷d); |
| 801 | scale = divided; |
| 802 | float128_t d_f128; |
| 803 | i32_to_f128M(d, &d_f128); |
| 804 | float128_t add_op; |
| 805 | f128M_mul(&d_f128, &scale, &add_op); |
| 806 | float128_t new_y; |
| 807 | f128M_add(&y, &add_op, &new_y); |
| 808 | y = new_y; |
| 809 | } |
| 810 | } else if (d && !gottail) { |
| 811 | //y += 0.5*scale; |
| 812 | { |
| 813 | float128_t add_op; |
| 814 | f128M_mul(&point_5, &scale, &add_op); |
| 815 | float128_t new_y; |
| 816 | f128M_add(&y, &add_op, &new_y); |
| 817 | y = new_y; |
| 818 | } |
| 819 | gottail = 1; |
| 820 | } |
| 821 | dc++; |
| 822 | } |
| 823 | } |
| 824 | if (!gotdig) { |
| 825 | shunget(f); |
| 826 | if (pok) { |
| 827 | shunget(f); |
| 828 | if (gotrad) shunget(f); |
| 829 | } else { |
| 830 | shlim(f, 0); |
| 831 | } |
| 832 | //return sign * 0.0; |
| 833 | return dbl_to_f128(sign * 0.0); |
| 834 | } |
| 835 | if (!gotrad) rp = dc; |
| 836 | while (dc<8) x *= 16, dc++; |
| 837 | if ((c|32)=='p') { |
| 838 | e2 = scanexp(f, pok); |
| 839 | if (e2 == LLONG_MIN) { |
| 840 | if (pok) { |
| 841 | shunget(f); |
| 842 | } else { |
| 843 | shlim(f, 0); |
| 844 | return zero; |
| 845 | } |
| 846 | e2 = 0; |
| 847 | } |
| 848 | } else { |
| 849 | shunget(f); |
| 850 | } |
| 851 | e2 += 4*rp - 32; |
| 852 | |
| 853 | if (!x) { |
| 854 | //return sign * 0.0; |
| 855 | return dbl_to_f128(sign * 0.0); |
| 856 | } |
| 857 | if (e2 > -emin) { |
| 858 | errno = ERANGE; |
| 859 | //return sign * LDBL_MAX * LDBL_MAX; |
| 860 | return zero; |
| 861 | } |
| 862 | if (e2 < emin-2*LDBL_MANT_DIG) { |
| 863 | errno = ERANGE; |
| 864 | //return sign * LDBL_MIN * LDBL_MIN; |
| 865 | return zero; |
| 866 | } |
| 867 | |
| 868 | while (x < 0x80000000) { |
| 869 | //if (y>=0.5) |
| 870 | if (!f128M_lt(&y, &point_5)) { |
| 871 | x += x + 1; |
| 872 | //y += y - 1; |
| 873 | { |
| 874 | float128_t minus_one; |
| 875 | f128M_sub(&y, &one, &minus_one); |
| 876 | float128_t new_y; |
| 877 | f128M_add(&y, &minus_one, &new_y); |
| 878 | y = new_y; |
| 879 | } |
| 880 | } else { |
| 881 | x += x; |
| 882 | //y += y; |
| 883 | { |
| 884 | float128_t new_y; |
| 885 | f128M_add(&y, &y, &new_y); |
| 886 | y = new_y; |
| 887 | } |
| 888 | } |
| 889 | e2--; |
| 890 | } |
| 891 | |
| 892 | if (bits > 32+e2-emin) { |
| 893 | bits = 32+e2-emin; |
| 894 | if (bits<0) bits=0; |
| 895 | } |
| 896 | |
| 897 | if (bits < LDBL_MANT_DIG) { |
| 898 | float128_t sign_f128; |
| 899 | i32_to_f128M(sign, &sign_f128); |
| 900 | bias = copysignf128(dbl_to_f128(scalbn(1, 32+LDBL_MANT_DIG-bits-1)), sign_f128); |
| 901 | } |
| 902 | |
| 903 | //if (bits<32 && y && !(x&1)) x++, y=0; |
| 904 | if (bits<32 && !f128M_eq(&y, &zero) && !(x&1)) x++, y=zero; |
| 905 | |
| 906 | //y = bias + sign*(float128_t)x + sign*y; |
| 907 | { |
| 908 | float128_t x_f128; |
| 909 | ui32_to_f128M(x, &x_f128); |
| 910 | float128_t sign_f128; |
| 911 | i32_to_f128M(sign, &sign_f128); |
| 912 | float128_t sign_mul_x; |
| 913 | f128M_mul(&sign_f128, &x_f128, &sign_mul_x); |
| 914 | float128_t sign_mul_y; |
| 915 | f128M_mul(&sign_f128, &y, &sign_mul_y); |
| 916 | float128_t bias_op; |
| 917 | f128M_add(&bias, &sign_mul_x, &bias_op); |
| 918 | float128_t new_y; |
| 919 | f128M_add(&bias_op, &sign_mul_y, &new_y); |
| 920 | y = new_y; |
| 921 | } |
| 922 | //y -= bias; |
| 923 | { |
| 924 | float128_t new_y; |
| 925 | f128M_sub(&y, &bias, &new_y); |
| 926 | y = new_y; |
| 927 | } |
| 928 | |
| 929 | if (f128M_eq(&y, &zero)) errno = ERANGE; |
| 930 | |
| 931 | return scalbnf128(y, e2); |
| 932 | } |
| 933 | |
| 934 | static int isspace(int c) |
| 935 | { |
| 936 | return c == ' ' || (unsigned)c-'\t' < 5; |
| 937 | } |
| 938 | |
| 939 | static inline float128_t makeInf128() { |
| 940 | union ldshape ux; |
| 941 | ux.i2.hi = 0x7fff000000000000UL; |
| 942 | ux.i2.lo = 0x0UL; |
| 943 | return ux.f; |
| 944 | } |
| 945 | |
| 946 | static inline float128_t makeNaN128() { |
| 947 | uint64_t rand = 0UL; |
| 948 | union ldshape ux; |
| 949 | ux.i2.hi = 0x7fff000000000000UL | (rand & 0xffffffffffffUL); |
| 950 | ux.i2.lo = 0x0UL; |
| 951 | return ux.f; |
| 952 | } |
| 953 | |
| 954 | float128_t __floatscan(struct MuslFILE *f, int prec, int pok) |
| 955 | { |
| 956 | int sign = 1; |
| 957 | size_t i; |
| 958 | int bits = LDBL_MANT_DIG; |
| 959 | int emin = LDBL_MIN_EXP-bits; |
| 960 | int c; |
| 961 | |
| 962 | while (isspace((c=shgetc(f)))); |
| 963 | |
| 964 | if (c=='+' || c=='-') { |
| 965 | sign -= 2*(c=='-'); |
| 966 | c = shgetc(f); |
| 967 | } |
| 968 | |
| 969 | for (i=0; i<8 && (c|32)=="infinity"[i]; i++) |
| 970 | if (i<7) c = shgetc(f); |
| 971 | if (i==3 || i==8 || (i>3 && pok)) { |
| 972 | if (i!=8) { |
| 973 | shunget(f); |
| 974 | if (pok) for (; i>3; i--) shunget(f); |
| 975 | } |
| 976 | //return sign * INFINITY; |
| 977 | float128_t sign_f128; |
| 978 | i32_to_f128M(sign, &sign_f128); |
| 979 | float128_t infinity_f128 = makeInf128(); |
| 980 | float128_t result; |
| 981 | f128M_mul(&sign_f128, &infinity_f128, &result); |
| 982 | return result; |
| 983 | } |
| 984 | if (!i) for (i=0; i<3 && (c|32)=="nan"[i]; i++) |
| 985 | if (i<2) c = shgetc(f); |
| 986 | if (i==3) { |
| 987 | if (shgetc(f) != '(') { |
| 988 | shunget(f); |
| 989 | return makeNaN128(); |
| 990 | } |
| 991 | for (i=1; ; i++) { |
| 992 | c = shgetc(f); |
| 993 | if (c-'0'<10U || c-'A'<26U || c-'a'<26U || c=='_') |
| 994 | continue; |
| 995 | if (c==')') return makeNaN128(); |
| 996 | shunget(f); |
| 997 | if (!pok) { |
| 998 | errno = EINVAL; |
| 999 | shlim(f, 0); |
| 1000 | float128_t zero; |
| 1001 | ui32_to_f128M(0, &zero); |
| 1002 | return zero; |
| 1003 | } |
| 1004 | while (i--) shunget(f); |
| 1005 | return makeNaN128(); |
| 1006 | } |
| 1007 | return makeNaN128(); |
| 1008 | } |
| 1009 | |
| 1010 | if (i) { |
| 1011 | shunget(f); |
| 1012 | errno = EINVAL; |
| 1013 | shlim(f, 0); |
| 1014 | float128_t zero; |
| 1015 | ui32_to_f128M(0, &zero); |
| 1016 | return zero; |
| 1017 | } |
| 1018 | |
| 1019 | if (c=='0') { |
| 1020 | c = shgetc(f); |
| 1021 | if ((c|32) == 'x') |
| 1022 | return hexfloat(f, bits, emin, sign, pok); |
| 1023 | shunget(f); |
| 1024 | c = '0'; |
| 1025 | } |
| 1026 | |
| 1027 | return decfloat(f, c, bits, emin, sign, pok); |
| 1028 | } |
| 1029 | |
| 1030 | float128_t parse_f128(const char *restrict s, char **restrict p) { |
| 1031 | struct MuslFILE f; |
| 1032 | sh_fromstring(&f, s); |
| 1033 | shlim(&f, 0); |
| 1034 | float128_t y = __floatscan(&f, 2, 1); |
| 1035 | off_t cnt = shcnt(&f); |
| 1036 | if (p) *p = cnt ? (char *)s + cnt : (char *)s; |
| 1037 | return y; |
| 1038 | } |
Andrew Kelley <andrew@ziglang.org> 2019-03-22 14:56:03-04:00