| ... | ... | @@ -16,7 +16,7 @@ |
| 16 | 16 | #define shunget(f) ((f)->shlim>=0 ? (void)(f)->rpos-- : (void)0) |
| 17 | 17 | |
| 18 | 18 | #define sh_fromstring(f, s) \ |
| 19 | | ((f)->buf = (f)->rpos = (void *)(s), (f)->rend = (void*)-1) |
| 19 | ((f)->buf = (f)->rpos = (void *)(s), (f)->rend = (void*)-1) |
| 20 | 20 | |
| 21 | 21 | #define LD_B1B_DIG 4 |
| 22 | 22 | #define LD_B1B_MAX 10384593, 717069655, 257060992, 658440191 |
| ... | ... | @@ -50,147 +50,147 @@ |
| 50 | 50 | |
| 51 | 51 | #if __BYTE_ORDER == __LITTLE_ENDIAN |
| 52 | 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; |
| 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 | 64 | }; |
| 65 | 65 | #elif __BYTE_ORDER == __BIG_ENDIAN |
| 66 | 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; |
| 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 | 78 | }; |
| 79 | 79 | #error Unsupported endian |
| 80 | 80 | #endif |
| 81 | 81 | |
| 82 | 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; |
| 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 | 109 | }; |
| 110 | 110 | |
| 111 | 111 | static void __shlim(struct MuslFILE *f, off_t lim) |
| 112 | 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; |
| 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 | 120 | } |
| 121 | 121 | |
| 122 | 122 | static int __toread(struct MuslFILE *f) |
| 123 | 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; |
| 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 | 133 | } |
| 134 | 134 | |
| 135 | 135 | static int __uflow(struct MuslFILE *f) |
| 136 | 136 | { |
| 137 | | unsigned char c; |
| 138 | | if (!__toread(f) && f->read(f, &c, 1)==1) return c; |
| 139 | | return EOF; |
| 137 | unsigned char c; |
| 138 | if (!__toread(f) && f->read(f, &c, 1)==1) return c; |
| 139 | return EOF; |
| 140 | 140 | } |
| 141 | 141 | |
| 142 | 142 | static int __shgetc(struct MuslFILE *f) |
| 143 | 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; |
| 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 | 160 | } |
| 161 | 161 | |
| 162 | 162 | static long long scanexp(struct MuslFILE *f, int pok) |
| 163 | 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; |
| 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 | 186 | } |
| 187 | 187 | |
| 188 | 188 | static float128_t copysignf128(float128_t x, float128_t y) |
| 189 | 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; |
| 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 | 194 | } |
| 195 | 195 | |
| 196 | 196 | static void mul_eq_f128_float(float128_t *x, float op_float) { |
| ... | ... | @@ -214,101 +214,101 @@ static float128_t dbl_to_f128(double x) { |
| 214 | 214 | |
| 215 | 215 | static float128_t fmodf128(float128_t x, float128_t y) |
| 216 | 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; |
| 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 | 221 | |
| 222 | 222 | float128_t zero; |
| 223 | 223 | ui32_to_f128M(0, &zero); |
| 224 | 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); |
| 225 | if (f128M_eq(&y, &zero) || f128M_isSignalingNaN(&y) || ex == 0x7fff) { |
| 226 | //return (x*y)/(x*y); |
| 227 | 227 | float128_t x_times_y; |
| 228 | 228 | f128M_mul(&x, &y, &x_times_y); |
| 229 | 229 | float128_t result; |
| 230 | 230 | f128M_div(&x_times_y, &x_times_y, &result); |
| 231 | 231 | return result; |
| 232 | 232 | } |
| 233 | | ux.i.se = ex; |
| 234 | | uy.i.se = ey; |
| 235 | | //if (ux.f <= uy.f) { |
| 233 | ux.i.se = ex; |
| 234 | uy.i.se = ey; |
| 235 | //if (ux.f <= uy.f) { |
| 236 | 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; |
| 237 | //if (ux.f == uy.f) { |
| 238 | if (f128M_eq(&ux.f, &uy.f)) { |
| 239 | //return 0*x; |
| 240 | 240 | float128_t result; |
| 241 | 241 | f128M_mul(&zero, &x, &result); |
| 242 | 242 | return result; |
| 243 | 243 | } |
| 244 | | return x; |
| 245 | | } |
| 244 | return x; |
| 245 | } |
| 246 | 246 | |
| 247 | | /* normalize x and y */ |
| 248 | | if (!ex) { |
| 249 | | //ux.f *= 0x1p120f; |
| 247 | /* normalize x and y */ |
| 248 | if (!ex) { |
| 249 | //ux.f *= 0x1p120f; |
| 250 | 250 | mul_eq_f128_float(&ux.f, 0x1p120f); |
| 251 | 251 | |
| 252 | | ex = ux.i.se - 120; |
| 253 | | } |
| 254 | | if (!ey) { |
| 255 | | //uy.f *= 0x1p120f; |
| 252 | ex = ux.i.se - 120; |
| 253 | } |
| 254 | if (!ey) { |
| 255 | //uy.f *= 0x1p120f; |
| 256 | 256 | mul_eq_f128_float(&uy.f, 0x1p120f); |
| 257 | 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; |
| 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 | 275 | float128_t result; |
| 276 | 276 | f128M_mul(&zero, &x, &result); |
| 277 | 277 | return result; |
| 278 | 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; |
| 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 | 293 | float128_t result; |
| 294 | 294 | f128M_mul(&zero, &x, &result); |
| 295 | 295 | return result; |
| 296 | 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; |
| 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 | 312 | } |
| 313 | 313 | |
| 314 | 314 | static float128_t int_mul_f128_cast_u32(int sign, uint32_t x0) { |
| ... | ... | @@ -348,7 +348,7 @@ static void mul_eq_f128_int(float128_t *y, int sign) { |
| 348 | 348 | } |
| 349 | 349 | |
| 350 | 350 | static float128_t make_f128(uint64_t hi, uint64_t lo) { |
| 351 | | union ldshape ux; |
| 351 | union ldshape ux; |
| 352 | 352 | ux.i2.hi = hi; |
| 353 | 353 | ux.i2.lo = lo; |
| 354 | 354 | return ux.f; |
| ... | ... | @@ -374,21 +374,21 @@ static void add_eq_f128_dbl(float128_t *a, double b) { |
| 374 | 374 | |
| 375 | 375 | static float128_t scalbnf128(float128_t x, int n) |
| 376 | 376 | { |
| 377 | | union ldshape u; |
| 377 | union ldshape u; |
| 378 | 378 | |
| 379 | | if (n > 16383) { |
| 380 | | //x *= 0x1p16383q; |
| 379 | if (n > 16383) { |
| 380 | //x *= 0x1p16383q; |
| 381 | 381 | mul_eq_f128_f128(&x, make_f128(0x7ffe000000000000, 0x0000000000000000)); |
| 382 | | n -= 16383; |
| 383 | | if (n > 16383) { |
| 384 | | //x *= 0x1p16383q; |
| 382 | n -= 16383; |
| 383 | if (n > 16383) { |
| 384 | //x *= 0x1p16383q; |
| 385 | 385 | mul_eq_f128_f128(&x, make_f128(0x7ffe000000000000, 0x0000000000000000)); |
| 386 | | n -= 16383; |
| 387 | | if (n > 16383) |
| 388 | | n = 16383; |
| 389 | | } |
| 390 | | } else if (n < -16382) { |
| 391 | | //x *= 0x1p-16382q * 0x1p113q; |
| 386 | n -= 16383; |
| 387 | if (n > 16383) |
| 388 | n = 16383; |
| 389 | } |
| 390 | } else if (n < -16382) { |
| 391 | //x *= 0x1p-16382q * 0x1p113q; |
| 392 | 392 | { |
| 393 | 393 | float128_t mul_result; |
| 394 | 394 | float128_t a = make_f128(0x0001000000000000, 0x0000000000000000); |
| ... | ... | @@ -396,9 +396,9 @@ static float128_t scalbnf128(float128_t x, int n) |
| 396 | 396 | f128M_mul(&a, &b, &mul_result); |
| 397 | 397 | mul_eq_f128_f128(&x, mul_result); |
| 398 | 398 | } |
| 399 | | n += 16382 - 113; |
| 400 | | if (n < -16382) { |
| 401 | | //x *= 0x1p-16382q * 0x1p113q; |
| 399 | n += 16382 - 113; |
| 400 | if (n < -16382) { |
| 401 | //x *= 0x1p-16382q * 0x1p113q; |
| 402 | 402 | { |
| 403 | 403 | float128_t mul_result; |
| 404 | 404 | float128_t a = make_f128(0x0001000000000000, 0x0000000000000000); |
| ... | ... | @@ -406,113 +406,113 @@ static float128_t scalbnf128(float128_t x, int n) |
| 406 | 406 | f128M_mul(&a, &b, &mul_result); |
| 407 | 407 | mul_eq_f128_f128(&x, mul_result); |
| 408 | 408 | } |
| 409 | | n += 16382 - 113; |
| 410 | | if (n < -16382) |
| 411 | | n = -16382; |
| 412 | | } |
| 413 | | } |
| 414 | | //u.f = 1.0; |
| 409 | n += 16382 - 113; |
| 410 | if (n < -16382) |
| 411 | n = -16382; |
| 412 | } |
| 413 | } |
| 414 | //u.f = 1.0; |
| 415 | 415 | ui32_to_f128M(1, &u.f); |
| 416 | | u.i.se = 0x3fff + n; |
| 416 | u.i.se = 0x3fff + n; |
| 417 | 417 | mul_eq_f128_f128(&x, u.f); |
| 418 | 418 | return x; |
| 419 | 419 | } |
| 420 | 420 | |
| 421 | 421 | static float128_t fabsf128(float128_t x) |
| 422 | 422 | { |
| 423 | | union ldshape u = {x}; |
| 423 | union ldshape u = {x}; |
| 424 | 424 | |
| 425 | | u.i.se &= 0x7fff; |
| 426 | | return u.f; |
| 425 | u.i.se &= 0x7fff; |
| 426 | return u.f; |
| 427 | 427 | } |
| 428 | 428 | |
| 429 | 429 | static float128_t decfloat(struct MuslFILE *f, int c, int bits, int emin, int sign, int pok) |
| 430 | 430 | { |
| 431 | | uint32_t x[KMAX]; |
| 432 | | static const uint32_t th[] = { LD_B1B_MAX }; |
| 433 | | int i, j, k, a, z; |
| 434 | | long long lrp=0, dc=0; |
| 435 | | long long e10=0; |
| 436 | | int lnz = 0; |
| 437 | | int gotdig = 0, gotrad = 0; |
| 438 | | int rp; |
| 439 | | int e2; |
| 440 | | int emax = -emin-bits+3; |
| 441 | | int denormal = 0; |
| 442 | | float128_t y; |
| 431 | uint32_t x[KMAX]; |
| 432 | static const uint32_t th[] = { LD_B1B_MAX }; |
| 433 | int i, j, k, a, z; |
| 434 | long long lrp=0, dc=0; |
| 435 | long long e10=0; |
| 436 | int lnz = 0; |
| 437 | int gotdig = 0, gotrad = 0; |
| 438 | int rp; |
| 439 | int e2; |
| 440 | int emax = -emin-bits+3; |
| 441 | int denormal = 0; |
| 442 | float128_t y; |
| 443 | 443 | float128_t zero; |
| 444 | 444 | ui32_to_f128M(0, &zero); |
| 445 | | float128_t frac=zero; |
| 446 | | float128_t bias=zero; |
| 447 | | static const int p10s[] = { 10, 100, 1000, 10000, |
| 448 | | 100000, 1000000, 10000000, 100000000 }; |
| 449 | | |
| 450 | | j=0; |
| 451 | | k=0; |
| 452 | | |
| 453 | | /* Don't let leading zeros consume buffer space */ |
| 454 | | for (; c=='0'; c = shgetc(f)) gotdig=1; |
| 455 | | if (c=='.') { |
| 456 | | gotrad = 1; |
| 457 | | for (c = shgetc(f); c=='0'; c = shgetc(f)) gotdig=1, lrp--; |
| 458 | | } |
| 459 | | |
| 460 | | x[0] = 0; |
| 461 | | for (; c-'0'<10U || c=='.'; c = shgetc(f)) { |
| 462 | | if (c == '.') { |
| 463 | | if (gotrad) break; |
| 464 | | gotrad = 1; |
| 465 | | lrp = dc; |
| 466 | | } else if (k < KMAX-3) { |
| 467 | | dc++; |
| 468 | | if (c!='0') lnz = dc; |
| 469 | | if (j) x[k] = x[k]*10 + c-'0'; |
| 470 | | else x[k] = c-'0'; |
| 471 | | if (++j==9) { |
| 472 | | k++; |
| 473 | | j=0; |
| 474 | | } |
| 475 | | gotdig=1; |
| 476 | | } else { |
| 477 | | dc++; |
| 478 | | if (c!='0') { |
| 479 | | lnz = (KMAX-4)*9; |
| 480 | | x[KMAX-4] |= 1; |
| 481 | | } |
| 482 | | } |
| 483 | | } |
| 484 | | if (!gotrad) lrp=dc; |
| 485 | | |
| 486 | | if (gotdig && (c|32)=='e') { |
| 487 | | e10 = scanexp(f, pok); |
| 488 | | if (e10 == LLONG_MIN) { |
| 489 | | if (pok) { |
| 490 | | shunget(f); |
| 491 | | } else { |
| 492 | | shlim(f, 0); |
| 493 | | return zero; |
| 494 | | } |
| 495 | | e10 = 0; |
| 496 | | } |
| 497 | | lrp += e10; |
| 498 | | } else if (c>=0) { |
| 499 | | shunget(f); |
| 500 | | } |
| 501 | | if (!gotdig) { |
| 502 | | errno = EINVAL; |
| 503 | | shlim(f, 0); |
| 504 | | return zero; |
| 505 | | } |
| 506 | | |
| 507 | | /* Handle zero specially to avoid nasty special cases later */ |
| 508 | | if (!x[0]) { |
| 445 | float128_t frac=zero; |
| 446 | float128_t bias=zero; |
| 447 | static const int p10s[] = { 10, 100, 1000, 10000, |
| 448 | 100000, 1000000, 10000000, 100000000 }; |
| 449 | |
| 450 | j=0; |
| 451 | k=0; |
| 452 | |
| 453 | /* Don't let leading zeros consume buffer space */ |
| 454 | for (; c=='0'; c = shgetc(f)) gotdig=1; |
| 455 | if (c=='.') { |
| 456 | gotrad = 1; |
| 457 | for (c = shgetc(f); c=='0'; c = shgetc(f)) gotdig=1, lrp--; |
| 458 | } |
| 459 | |
| 460 | x[0] = 0; |
| 461 | for (; c-'0'<10U || c=='.'; c = shgetc(f)) { |
| 462 | if (c == '.') { |
| 463 | if (gotrad) break; |
| 464 | gotrad = 1; |
| 465 | lrp = dc; |
| 466 | } else if (k < KMAX-3) { |
| 467 | dc++; |
| 468 | if (c!='0') lnz = dc; |
| 469 | if (j) x[k] = x[k]*10 + c-'0'; |
| 470 | else x[k] = c-'0'; |
| 471 | if (++j==9) { |
| 472 | k++; |
| 473 | j=0; |
| 474 | } |
| 475 | gotdig=1; |
| 476 | } else { |
| 477 | dc++; |
| 478 | if (c!='0') { |
| 479 | lnz = (KMAX-4)*9; |
| 480 | x[KMAX-4] |= 1; |
| 481 | } |
| 482 | } |
| 483 | } |
| 484 | if (!gotrad) lrp=dc; |
| 485 | |
| 486 | if (gotdig && (c|32)=='e') { |
| 487 | e10 = scanexp(f, pok); |
| 488 | if (e10 == LLONG_MIN) { |
| 489 | if (pok) { |
| 490 | shunget(f); |
| 491 | } else { |
| 492 | shlim(f, 0); |
| 493 | return zero; |
| 494 | } |
| 495 | e10 = 0; |
| 496 | } |
| 497 | lrp += e10; |
| 498 | } else if (c>=0) { |
| 499 | shunget(f); |
| 500 | } |
| 501 | if (!gotdig) { |
| 502 | errno = EINVAL; |
| 503 | shlim(f, 0); |
| 504 | return zero; |
| 505 | } |
| 506 | |
| 507 | /* Handle zero specially to avoid nasty special cases later */ |
| 508 | if (!x[0]) { |
| 509 | 509 | //return sign * 0.0; |
| 510 | 510 | return dbl_to_f128(sign * 0.0); |
| 511 | 511 | } |
| 512 | 512 | |
| 513 | | /* Optimize small integers (w/no exponent) and over/under-flow */ |
| 514 | | if (lrp==dc && dc<10 && (bits>30 || x[0]>>bits==0)) { |
| 515 | | //return sign * (float128_t)x[0]; |
| 513 | /* Optimize small integers (w/no exponent) and over/under-flow */ |
| 514 | if (lrp==dc && dc<10 && (bits>30 || x[0]>>bits==0)) { |
| 515 | //return sign * (float128_t)x[0]; |
| 516 | 516 | float128_t sign_f128; |
| 517 | 517 | i32_to_f128M(sign, &sign_f128); |
| 518 | 518 | float128_t x0_f128; |
| ... | ... | @@ -521,131 +521,131 @@ static float128_t decfloat(struct MuslFILE *f, int c, int bits, int emin, int si |
| 521 | 521 | f128M_mul(&sign_f128, &x0_f128, &result); |
| 522 | 522 | return result; |
| 523 | 523 | } |
| 524 | | if (lrp > -emin/2) { |
| 525 | | errno = ERANGE; |
| 526 | | //return sign * LDBL_MAX * LDBL_MAX; |
| 527 | | return zero; |
| 528 | | } |
| 529 | | if (lrp < emin-2*LDBL_MANT_DIG) { |
| 530 | | errno = ERANGE; |
| 531 | | //return sign * LDBL_MIN * LDBL_MIN; |
| 532 | | return zero; |
| 533 | | } |
| 534 | | |
| 535 | | /* Align incomplete final B1B digit */ |
| 536 | | if (j) { |
| 537 | | for (; j<9; j++) x[k]*=10; |
| 538 | | k++; |
| 539 | | j=0; |
| 540 | | } |
| 541 | | |
| 542 | | a = 0; |
| 543 | | z = k; |
| 544 | | e2 = 0; |
| 545 | | rp = lrp; |
| 546 | | |
| 547 | | /* Optimize small to mid-size integers (even in exp. notation) */ |
| 548 | | if (lnz<9 && lnz<=rp && rp < 18) { |
| 549 | | if (rp == 9) { |
| 524 | if (lrp > -emin/2) { |
| 525 | errno = ERANGE; |
| 526 | //return sign * LDBL_MAX * LDBL_MAX; |
| 527 | return zero; |
| 528 | } |
| 529 | if (lrp < emin-2*LDBL_MANT_DIG) { |
| 530 | errno = ERANGE; |
| 531 | //return sign * LDBL_MIN * LDBL_MIN; |
| 532 | return zero; |
| 533 | } |
| 534 | |
| 535 | /* Align incomplete final B1B digit */ |
| 536 | if (j) { |
| 537 | for (; j<9; j++) x[k]*=10; |
| 538 | k++; |
| 539 | j=0; |
| 540 | } |
| 541 | |
| 542 | a = 0; |
| 543 | z = k; |
| 544 | e2 = 0; |
| 545 | rp = lrp; |
| 546 | |
| 547 | /* Optimize small to mid-size integers (even in exp. notation) */ |
| 548 | if (lnz<9 && lnz<=rp && rp < 18) { |
| 549 | if (rp == 9) { |
| 550 | 550 | //return sign * (float128_t)(x[0]); |
| 551 | 551 | return int_mul_f128_cast_u32(sign, x[0]); |
| 552 | 552 | } |
| 553 | | if (rp < 9) { |
| 553 | if (rp < 9) { |
| 554 | 554 | //return sign * (float128_t)(x[0]) / p10s[8-rp]; |
| 555 | 555 | return triple_divide(sign, x[0], p10s[8-rp]); |
| 556 | 556 | } |
| 557 | | int bitlim = bits-3*(int)(rp-9); |
| 558 | | if (bitlim>30 || x[0]>>bitlim==0) |
| 559 | | //return sign * (float128_t)(x[0]) * p10s[rp-10]; |
| 560 | | return triple_multiply(sign, x[0], p10s[rp-10]); |
| 561 | | } |
| 562 | | |
| 563 | | /* Drop trailing zeros */ |
| 564 | | for (; !x[z-1]; z--); |
| 565 | | |
| 566 | | /* Align radix point to B1B digit boundary */ |
| 567 | | if (rp % 9) { |
| 568 | | int rpm9 = rp>=0 ? rp%9 : rp%9+9; |
| 569 | | int p10 = p10s[8-rpm9]; |
| 570 | | uint32_t carry = 0; |
| 571 | | for (k=a; k!=z; k++) { |
| 572 | | uint32_t tmp = x[k] % p10; |
| 573 | | x[k] = x[k]/p10 + carry; |
| 574 | | carry = 1000000000/p10 * tmp; |
| 575 | | if (k==a && !x[k]) { |
| 576 | | a = (a+1 & MASK); |
| 577 | | rp -= 9; |
| 578 | | } |
| 579 | | } |
| 580 | | if (carry) x[z++] = carry; |
| 581 | | rp += 9-rpm9; |
| 582 | | } |
| 583 | | |
| 584 | | /* Upscale until desired number of bits are left of radix point */ |
| 585 | | while (rp < 9*LD_B1B_DIG || (rp == 9*LD_B1B_DIG && x[a]<th[0])) { |
| 586 | | uint32_t carry = 0; |
| 587 | | e2 -= 29; |
| 588 | | for (k=(z-1 & MASK); ; k=(k-1 & MASK)) { |
| 589 | | uint64_t tmp = ((uint64_t)x[k] << 29) + carry; |
| 590 | | if (tmp > 1000000000) { |
| 591 | | carry = tmp / 1000000000; |
| 592 | | x[k] = tmp % 1000000000; |
| 593 | | } else { |
| 594 | | carry = 0; |
| 595 | | x[k] = tmp; |
| 596 | | } |
| 597 | | if (k==(z-1 & MASK) && k!=a && !x[k]) z = k; |
| 598 | | if (k==a) break; |
| 599 | | } |
| 600 | | if (carry) { |
| 601 | | rp += 9; |
| 602 | | a = (a-1 & MASK); |
| 603 | | if (a == z) { |
| 604 | | z = (z-1 & MASK); |
| 605 | | x[z-1 & MASK] |= x[z]; |
| 606 | | } |
| 607 | | x[a] = carry; |
| 608 | | } |
| 609 | | } |
| 610 | | |
| 611 | | /* Downscale until exactly number of bits are left of radix point */ |
| 612 | | for (;;) { |
| 613 | | uint32_t carry = 0; |
| 614 | | int sh = 1; |
| 615 | | for (i=0; i<LD_B1B_DIG; i++) { |
| 616 | | k = (a+i & MASK); |
| 617 | | if (k == z || x[k] < th[i]) { |
| 618 | | i=LD_B1B_DIG; |
| 619 | | break; |
| 620 | | } |
| 621 | | if (x[a+i & MASK] > th[i]) break; |
| 622 | | } |
| 623 | | if (i==LD_B1B_DIG && rp==9*LD_B1B_DIG) break; |
| 624 | | /* FIXME: find a way to compute optimal sh */ |
| 625 | | if (rp > 9+9*LD_B1B_DIG) sh = 9; |
| 626 | | e2 += sh; |
| 627 | | for (k=a; k!=z; k=(k+1 & MASK)) { |
| 628 | | uint32_t tmp = x[k] & (1<<sh)-1; |
| 629 | | x[k] = (x[k]>>sh) + carry; |
| 630 | | carry = (1000000000>>sh) * tmp; |
| 631 | | if (k==a && !x[k]) { |
| 632 | | a = (a+1 & MASK); |
| 633 | | i--; |
| 634 | | rp -= 9; |
| 635 | | } |
| 636 | | } |
| 637 | | if (carry) { |
| 638 | | if ((z+1 & MASK) != a) { |
| 639 | | x[z] = carry; |
| 640 | | z = (z+1 & MASK); |
| 641 | | } else x[z-1 & MASK] |= 1; |
| 642 | | } |
| 643 | | } |
| 644 | | |
| 645 | | /* Assemble desired bits into floating point variable */ |
| 646 | | for (y=zero,i=0; i<LD_B1B_DIG; i++) { |
| 647 | | if ((a+i & MASK)==z) x[(z=(z+1 & MASK))-1] = 0; |
| 648 | | //y = 1000000000.0L * y + x[a+i & MASK]; |
| 557 | int bitlim = bits-3*(int)(rp-9); |
| 558 | if (bitlim>30 || x[0]>>bitlim==0) |
| 559 | //return sign * (float128_t)(x[0]) * p10s[rp-10]; |
| 560 | return triple_multiply(sign, x[0], p10s[rp-10]); |
| 561 | } |
| 562 | |
| 563 | /* Drop trailing zeros */ |
| 564 | for (; !x[z-1]; z--); |
| 565 | |
| 566 | /* Align radix point to B1B digit boundary */ |
| 567 | if (rp % 9) { |
| 568 | int rpm9 = rp>=0 ? rp%9 : rp%9+9; |
| 569 | int p10 = p10s[8-rpm9]; |
| 570 | uint32_t carry = 0; |
| 571 | for (k=a; k!=z; k++) { |
| 572 | uint32_t tmp = x[k] % p10; |
| 573 | x[k] = x[k]/p10 + carry; |
| 574 | carry = 1000000000/p10 * tmp; |
| 575 | if (k==a && !x[k]) { |
| 576 | a = (a+1 & MASK); |
| 577 | rp -= 9; |
| 578 | } |
| 579 | } |
| 580 | if (carry) x[z++] = carry; |
| 581 | rp += 9-rpm9; |
| 582 | } |
| 583 | |
| 584 | /* Upscale until desired number of bits are left of radix point */ |
| 585 | while (rp < 9*LD_B1B_DIG || (rp == 9*LD_B1B_DIG && x[a]<th[0])) { |
| 586 | uint32_t carry = 0; |
| 587 | e2 -= 29; |
| 588 | for (k=(z-1 & MASK); ; k=(k-1 & MASK)) { |
| 589 | uint64_t tmp = ((uint64_t)x[k] << 29) + carry; |
| 590 | if (tmp > 1000000000) { |
| 591 | carry = tmp / 1000000000; |
| 592 | x[k] = tmp % 1000000000; |
| 593 | } else { |
| 594 | carry = 0; |
| 595 | x[k] = tmp; |
| 596 | } |
| 597 | if (k==(z-1 & MASK) && k!=a && !x[k]) z = k; |
| 598 | if (k==a) break; |
| 599 | } |
| 600 | if (carry) { |
| 601 | rp += 9; |
| 602 | a = (a-1 & MASK); |
| 603 | if (a == z) { |
| 604 | z = (z-1 & MASK); |
| 605 | x[z-1 & MASK] |= x[z]; |
| 606 | } |
| 607 | x[a] = carry; |
| 608 | } |
| 609 | } |
| 610 | |
| 611 | /* Downscale until exactly number of bits are left of radix point */ |
| 612 | for (;;) { |
| 613 | uint32_t carry = 0; |
| 614 | int sh = 1; |
| 615 | for (i=0; i<LD_B1B_DIG; i++) { |
| 616 | k = (a+i & MASK); |
| 617 | if (k == z || x[k] < th[i]) { |
| 618 | i=LD_B1B_DIG; |
| 619 | break; |
| 620 | } |
| 621 | if (x[a+i & MASK] > th[i]) break; |
| 622 | } |
| 623 | if (i==LD_B1B_DIG && rp==9*LD_B1B_DIG) break; |
| 624 | /* FIXME: find a way to compute optimal sh */ |
| 625 | if (rp > 9+9*LD_B1B_DIG) sh = 9; |
| 626 | e2 += sh; |
| 627 | for (k=a; k!=z; k=(k+1 & MASK)) { |
| 628 | uint32_t tmp = x[k] & (1<<sh)-1; |
| 629 | x[k] = (x[k]>>sh) + carry; |
| 630 | carry = (1000000000>>sh) * tmp; |
| 631 | if (k==a && !x[k]) { |
| 632 | a = (a+1 & MASK); |
| 633 | i--; |
| 634 | rp -= 9; |
| 635 | } |
| 636 | } |
| 637 | if (carry) { |
| 638 | if ((z+1 & MASK) != a) { |
| 639 | x[z] = carry; |
| 640 | z = (z+1 & MASK); |
| 641 | } else x[z-1 & MASK] |= 1; |
| 642 | } |
| 643 | } |
| 644 | |
| 645 | /* Assemble desired bits into floating point variable */ |
| 646 | for (y=zero,i=0; i<LD_B1B_DIG; i++) { |
| 647 | if ((a+i & MASK)==z) x[(z=(z+1 & MASK))-1] = 0; |
| 648 | //y = 1000000000.0L * y + x[a+i & MASK]; |
| 649 | 649 | float128_t const_f128; |
| 650 | 650 | ui64_to_f128M(1000000000, &const_f128); |
| 651 | 651 | float128_t mul_y; |
| ... | ... | @@ -653,56 +653,56 @@ static float128_t decfloat(struct MuslFILE *f, int c, int bits, int emin, int si |
| 653 | 653 | float128_t x_f128; |
| 654 | 654 | ui32_to_f128M(x[a+i & MASK], &x_f128); |
| 655 | 655 | f128M_add(&mul_y, &x_f128, &y); |
| 656 | | } |
| 656 | } |
| 657 | 657 | |
| 658 | | //y *= sign; |
| 658 | //y *= sign; |
| 659 | 659 | mul_eq_f128_int(&y, sign); |
| 660 | 660 | |
| 661 | | /* Limit precision for denormal results */ |
| 662 | | if (bits > LDBL_MANT_DIG+e2-emin) { |
| 663 | | bits = LDBL_MANT_DIG+e2-emin; |
| 664 | | if (bits<0) bits=0; |
| 665 | | denormal = 1; |
| 666 | | } |
| 667 | | |
| 668 | | /* Calculate bias term to force rounding, move out lower bits */ |
| 669 | | if (bits < LDBL_MANT_DIG) { |
| 670 | | bias = copysignf128(dbl_to_f128(scalbn(1, 2*LDBL_MANT_DIG-bits-1)), y); |
| 671 | | frac = fmodf128(y, dbl_to_f128(scalbn(1, LDBL_MANT_DIG-bits))); |
| 672 | | //y -= frac; |
| 661 | /* Limit precision for denormal results */ |
| 662 | if (bits > LDBL_MANT_DIG+e2-emin) { |
| 663 | bits = LDBL_MANT_DIG+e2-emin; |
| 664 | if (bits<0) bits=0; |
| 665 | denormal = 1; |
| 666 | } |
| 667 | |
| 668 | /* Calculate bias term to force rounding, move out lower bits */ |
| 669 | if (bits < LDBL_MANT_DIG) { |
| 670 | bias = copysignf128(dbl_to_f128(scalbn(1, 2*LDBL_MANT_DIG-bits-1)), y); |
| 671 | frac = fmodf128(y, dbl_to_f128(scalbn(1, LDBL_MANT_DIG-bits))); |
| 672 | //y -= frac; |
| 673 | 673 | { |
| 674 | 674 | float128_t new_value; |
| 675 | 675 | f128M_sub(&y, &frac, &new_value); |
| 676 | 676 | y = new_value; |
| 677 | 677 | } |
| 678 | | //y += bias; |
| 678 | //y += bias; |
| 679 | 679 | { |
| 680 | 680 | float128_t new_value; |
| 681 | 681 | f128M_add(&y, &frac, &new_value); |
| 682 | 682 | y = new_value; |
| 683 | 683 | } |
| 684 | | } |
| 684 | } |
| 685 | 685 | |
| 686 | | /* Process tail of decimal input so it can affect rounding */ |
| 687 | | if ((a+i & MASK) != z) { |
| 688 | | uint32_t t = x[a+i & MASK]; |
| 689 | | if (t < 500000000 && (t || (a+i+1 & MASK) != z)) { |
| 690 | | //frac += 0.25*sign; |
| 686 | /* Process tail of decimal input so it can affect rounding */ |
| 687 | if ((a+i & MASK) != z) { |
| 688 | uint32_t t = x[a+i & MASK]; |
| 689 | if (t < 500000000 && (t || (a+i+1 & MASK) != z)) { |
| 690 | //frac += 0.25*sign; |
| 691 | 691 | add_eq_f128_dbl(&frac, 0.25*sign); |
| 692 | 692 | } else if (t > 500000000) { |
| 693 | | //frac += 0.75*sign; |
| 693 | //frac += 0.75*sign; |
| 694 | 694 | add_eq_f128_dbl(&frac, 0.75*sign); |
| 695 | 695 | } else if (t == 500000000) { |
| 696 | | if ((a+i+1 & MASK) == z) { |
| 697 | | //frac += 0.5*sign; |
| 696 | if ((a+i+1 & MASK) == z) { |
| 697 | //frac += 0.5*sign; |
| 698 | 698 | add_eq_f128_dbl(&frac, 0.5*sign); |
| 699 | 699 | } else { |
| 700 | | //frac += 0.75*sign; |
| 700 | //frac += 0.75*sign; |
| 701 | 701 | add_eq_f128_dbl(&frac, 0.75*sign); |
| 702 | 702 | } |
| 703 | | } |
| 704 | | //if (LDBL_MANT_DIG-bits >= 2 && !fmodf128(frac, 1)) |
| 705 | | if (LDBL_MANT_DIG-bits >= 2) { |
| 703 | } |
| 704 | //if (LDBL_MANT_DIG-bits >= 2 && !fmodf128(frac, 1)) |
| 705 | if (LDBL_MANT_DIG-bits >= 2) { |
| 706 | 706 | float128_t one; |
| 707 | 707 | ui32_to_f128M(1, &one); |
| 708 | 708 | float128_t mod_result = fmodf128(frac, one); |
| ... | ... | @@ -711,29 +711,29 @@ static float128_t decfloat(struct MuslFILE *f, int c, int bits, int emin, int si |
| 711 | 711 | add_eq_f128_dbl(&frac, 1.0); |
| 712 | 712 | } |
| 713 | 713 | } |
| 714 | | } |
| 714 | } |
| 715 | 715 | |
| 716 | | //y += frac; |
| 716 | //y += frac; |
| 717 | 717 | { |
| 718 | 718 | float128_t new_value; |
| 719 | 719 | f128M_add(&y, &frac, &new_value); |
| 720 | 720 | y = new_value; |
| 721 | 721 | } |
| 722 | | //y -= bias; |
| 722 | //y -= bias; |
| 723 | 723 | { |
| 724 | 724 | float128_t new_value; |
| 725 | 725 | f128M_sub(&y, &bias, &new_value); |
| 726 | 726 | y = new_value; |
| 727 | 727 | } |
| 728 | 728 | |
| 729 | | if ((e2+LDBL_MANT_DIG & INT_MAX) > emax-5) { |
| 730 | | //if (fabsf128(y) >= 0x1p113) |
| 729 | if ((e2+LDBL_MANT_DIG & INT_MAX) > emax-5) { |
| 730 | //if (fabsf128(y) >= 0x1p113) |
| 731 | 731 | float128_t abs_y = fabsf128(y); |
| 732 | 732 | float128_t mant_f128 = make_f128(0x4070000000000000, 0x0000000000000000); |
| 733 | | if (!f128M_lt(&abs_y, &mant_f128)) { |
| 734 | | if (denormal && bits==LDBL_MANT_DIG+e2-emin) |
| 735 | | denormal = 0; |
| 736 | | //y *= 0.5; |
| 733 | if (!f128M_lt(&abs_y, &mant_f128)) { |
| 734 | if (denormal && bits==LDBL_MANT_DIG+e2-emin) |
| 735 | denormal = 0; |
| 736 | //y *= 0.5; |
| 737 | 737 | { |
| 738 | 738 | float128_t point_5 = dbl_to_f128(0.5); |
| 739 | 739 | float128_t new_value; |
| ... | ... | @@ -741,13 +741,13 @@ static float128_t decfloat(struct MuslFILE *f, int c, int bits, int emin, int si |
| 741 | 741 | y = new_value; |
| 742 | 742 | } |
| 743 | 743 | |
| 744 | | e2++; |
| 745 | | } |
| 746 | | if (e2+LDBL_MANT_DIG>emax || (denormal && !f128M_eq(&frac, &zero))) |
| 747 | | errno = ERANGE; |
| 748 | | } |
| 744 | e2++; |
| 745 | } |
| 746 | if (e2+LDBL_MANT_DIG>emax || (denormal && !f128M_eq(&frac, &zero))) |
| 747 | errno = ERANGE; |
| 748 | } |
| 749 | 749 | |
| 750 | | return scalbnf128(y, e2); |
| 750 | return scalbnf128(y, e2); |
| 751 | 751 | } |
| 752 | 752 | |
| 753 | 753 | static float128_t hexfloat(struct MuslFILE *f, int bits, int emin, int sign, int pok) |
| ... | ... | @@ -760,42 +760,42 @@ static float128_t hexfloat(struct MuslFILE *f, int bits, int emin, int sign, int |
| 760 | 760 | ui32_to_f128M(16, &sixteen); |
| 761 | 761 | float128_t point_5 = dbl_to_f128(0.5); |
| 762 | 762 | |
| 763 | | uint32_t x = 0; |
| 764 | | float128_t y = zero; |
| 765 | | float128_t scale = one; |
| 766 | | float128_t bias = zero; |
| 767 | | int gottail = 0, gotrad = 0, gotdig = 0; |
| 768 | | long long rp = 0; |
| 769 | | long long dc = 0; |
| 770 | | long long e2 = 0; |
| 771 | | int d; |
| 772 | | int c; |
| 773 | | |
| 774 | | c = shgetc(f); |
| 775 | | |
| 776 | | /* Skip leading zeros */ |
| 777 | | for (; c=='0'; c = shgetc(f)) gotdig = 1; |
| 778 | | |
| 779 | | if (c=='.') { |
| 780 | | gotrad = 1; |
| 781 | | c = shgetc(f); |
| 782 | | /* Count zeros after the radix point before significand */ |
| 783 | | for (rp=0; c=='0'; c = shgetc(f), rp--) gotdig = 1; |
| 784 | | } |
| 785 | | |
| 786 | | for (; c-'0'<10U || (c|32)-'a'<6U || c=='.'; c = shgetc(f)) { |
| 787 | | if (c=='.') { |
| 788 | | if (gotrad) break; |
| 789 | | rp = dc; |
| 790 | | gotrad = 1; |
| 791 | | } else { |
| 792 | | gotdig = 1; |
| 793 | | if (c > '9') d = (c|32)+10-'a'; |
| 794 | | else d = c-'0'; |
| 795 | | if (dc<8) { |
| 796 | | x = x*16 + d; |
| 797 | | } else if (dc < LDBL_MANT_DIG/4+1) { |
| 798 | | //y += d*(scale/=16); |
| 763 | uint32_t x = 0; |
| 764 | float128_t y = zero; |
| 765 | float128_t scale = one; |
| 766 | float128_t bias = zero; |
| 767 | int gottail = 0, gotrad = 0, gotdig = 0; |
| 768 | long long rp = 0; |
| 769 | long long dc = 0; |
| 770 | long long e2 = 0; |
| 771 | int d; |
| 772 | int c; |
| 773 | |
| 774 | c = shgetc(f); |
| 775 | |
| 776 | /* Skip leading zeros */ |
| 777 | for (; c=='0'; c = shgetc(f)) gotdig = 1; |
| 778 | |
| 779 | if (c=='.') { |
| 780 | gotrad = 1; |
| 781 | c = shgetc(f); |
| 782 | /* Count zeros after the radix point before significand */ |
| 783 | for (rp=0; c=='0'; c = shgetc(f), rp--) gotdig = 1; |
| 784 | } |
| 785 | |
| 786 | for (; c-'0'<10U || (c|32)-'a'<6U || c=='.'; c = shgetc(f)) { |
| 787 | if (c=='.') { |
| 788 | if (gotrad) break; |
| 789 | rp = dc; |
| 790 | gotrad = 1; |
| 791 | } else { |
| 792 | gotdig = 1; |
| 793 | if (c > '9') d = (c|32)+10-'a'; |
| 794 | else d = c-'0'; |
| 795 | if (dc<8) { |
| 796 | x = x*16 + d; |
| 797 | } else if (dc < LDBL_MANT_DIG/4+1) { |
| 798 | //y += d*(scale/=16); |
| 799 | 799 | { |
| 800 | 800 | float128_t divided; |
| 801 | 801 | f128M_div(&scale, &sixteen, ÷d); |
| ... | ... | @@ -808,8 +808,8 @@ static float128_t hexfloat(struct MuslFILE *f, int bits, int emin, int sign, int |
| 808 | 808 | f128M_add(&y, &add_op, &new_y); |
| 809 | 809 | y = new_y; |
| 810 | 810 | } |
| 811 | | } else if (d && !gottail) { |
| 812 | | //y += 0.5*scale; |
| 811 | } else if (d && !gottail) { |
| 812 | //y += 0.5*scale; |
| 813 | 813 | { |
| 814 | 814 | float128_t add_op; |
| 815 | 815 | f128M_mul(&point_5, &scale, &add_op); |
| ... | ... | @@ -817,60 +817,60 @@ static float128_t hexfloat(struct MuslFILE *f, int bits, int emin, int sign, int |
| 817 | 817 | f128M_add(&y, &add_op, &new_y); |
| 818 | 818 | y = new_y; |
| 819 | 819 | } |
| 820 | | gottail = 1; |
| 821 | | } |
| 822 | | dc++; |
| 823 | | } |
| 824 | | } |
| 825 | | if (!gotdig) { |
| 826 | | shunget(f); |
| 827 | | if (pok) { |
| 828 | | shunget(f); |
| 829 | | if (gotrad) shunget(f); |
| 830 | | } else { |
| 831 | | shlim(f, 0); |
| 832 | | } |
| 833 | | //return sign * 0.0; |
| 820 | gottail = 1; |
| 821 | } |
| 822 | dc++; |
| 823 | } |
| 824 | } |
| 825 | if (!gotdig) { |
| 826 | shunget(f); |
| 827 | if (pok) { |
| 828 | shunget(f); |
| 829 | if (gotrad) shunget(f); |
| 830 | } else { |
| 831 | shlim(f, 0); |
| 832 | } |
| 833 | //return sign * 0.0; |
| 834 | 834 | return dbl_to_f128(sign * 0.0); |
| 835 | | } |
| 836 | | if (!gotrad) rp = dc; |
| 837 | | while (dc<8) x *= 16, dc++; |
| 838 | | if ((c|32)=='p') { |
| 839 | | e2 = scanexp(f, pok); |
| 840 | | if (e2 == LLONG_MIN) { |
| 841 | | if (pok) { |
| 842 | | shunget(f); |
| 843 | | } else { |
| 844 | | shlim(f, 0); |
| 845 | | return zero; |
| 846 | | } |
| 847 | | e2 = 0; |
| 848 | | } |
| 849 | | } else { |
| 850 | | shunget(f); |
| 851 | | } |
| 852 | | e2 += 4*rp - 32; |
| 853 | | |
| 854 | | if (!x) { |
| 835 | } |
| 836 | if (!gotrad) rp = dc; |
| 837 | while (dc<8) x *= 16, dc++; |
| 838 | if ((c|32)=='p') { |
| 839 | e2 = scanexp(f, pok); |
| 840 | if (e2 == LLONG_MIN) { |
| 841 | if (pok) { |
| 842 | shunget(f); |
| 843 | } else { |
| 844 | shlim(f, 0); |
| 845 | return zero; |
| 846 | } |
| 847 | e2 = 0; |
| 848 | } |
| 849 | } else { |
| 850 | shunget(f); |
| 851 | } |
| 852 | e2 += 4*rp - 32; |
| 853 | |
| 854 | if (!x) { |
| 855 | 855 | //return sign * 0.0; |
| 856 | 856 | return dbl_to_f128(sign * 0.0); |
| 857 | 857 | } |
| 858 | | if (e2 > -emin) { |
| 859 | | errno = ERANGE; |
| 860 | | //return sign * LDBL_MAX * LDBL_MAX; |
| 861 | | return zero; |
| 862 | | } |
| 863 | | if (e2 < emin-2*LDBL_MANT_DIG) { |
| 864 | | errno = ERANGE; |
| 865 | | //return sign * LDBL_MIN * LDBL_MIN; |
| 866 | | return zero; |
| 867 | | } |
| 868 | | |
| 869 | | while (x < 0x80000000) { |
| 870 | | //if (y>=0.5) |
| 871 | | if (!f128M_lt(&y, &point_5)) { |
| 872 | | x += x + 1; |
| 873 | | //y += y - 1; |
| 858 | if (e2 > -emin) { |
| 859 | errno = ERANGE; |
| 860 | //return sign * LDBL_MAX * LDBL_MAX; |
| 861 | return zero; |
| 862 | } |
| 863 | if (e2 < emin-2*LDBL_MANT_DIG) { |
| 864 | errno = ERANGE; |
| 865 | //return sign * LDBL_MIN * LDBL_MIN; |
| 866 | return zero; |
| 867 | } |
| 868 | |
| 869 | while (x < 0x80000000) { |
| 870 | //if (y>=0.5) |
| 871 | if (!f128M_lt(&y, &point_5)) { |
| 872 | x += x + 1; |
| 873 | //y += y - 1; |
| 874 | 874 | { |
| 875 | 875 | float128_t minus_one; |
| 876 | 876 | f128M_sub(&y, &one, &minus_one); |
| ... | ... | @@ -878,33 +878,33 @@ static float128_t hexfloat(struct MuslFILE *f, int bits, int emin, int sign, int |
| 878 | 878 | f128M_add(&y, &minus_one, &new_y); |
| 879 | 879 | y = new_y; |
| 880 | 880 | } |
| 881 | | } else { |
| 882 | | x += x; |
| 883 | | //y += y; |
| 881 | } else { |
| 882 | x += x; |
| 883 | //y += y; |
| 884 | 884 | { |
| 885 | 885 | float128_t new_y; |
| 886 | 886 | f128M_add(&y, &y, &new_y); |
| 887 | 887 | y = new_y; |
| 888 | 888 | } |
| 889 | | } |
| 890 | | e2--; |
| 891 | | } |
| 889 | } |
| 890 | e2--; |
| 891 | } |
| 892 | 892 | |
| 893 | | if (bits > 32+e2-emin) { |
| 894 | | bits = 32+e2-emin; |
| 895 | | if (bits<0) bits=0; |
| 896 | | } |
| 893 | if (bits > 32+e2-emin) { |
| 894 | bits = 32+e2-emin; |
| 895 | if (bits<0) bits=0; |
| 896 | } |
| 897 | 897 | |
| 898 | | if (bits < LDBL_MANT_DIG) { |
| 898 | if (bits < LDBL_MANT_DIG) { |
| 899 | 899 | float128_t sign_f128; |
| 900 | 900 | i32_to_f128M(sign, &sign_f128); |
| 901 | | bias = copysignf128(dbl_to_f128(scalbn(1, 32+LDBL_MANT_DIG-bits-1)), sign_f128); |
| 901 | bias = copysignf128(dbl_to_f128(scalbn(1, 32+LDBL_MANT_DIG-bits-1)), sign_f128); |
| 902 | 902 | } |
| 903 | 903 | |
| 904 | | //if (bits<32 && y && !(x&1)) x++, y=0; |
| 905 | | if (bits<32 && !f128M_eq(&y, &zero) && !(x&1)) x++, y=zero; |
| 904 | //if (bits<32 && y && !(x&1)) x++, y=0; |
| 905 | if (bits<32 && !f128M_eq(&y, &zero) && !(x&1)) x++, y=zero; |
| 906 | 906 | |
| 907 | | //y = bias + sign*(float128_t)x + sign*y; |
| 907 | //y = bias + sign*(float128_t)x + sign*y; |
| 908 | 908 | { |
| 909 | 909 | float128_t x_f128; |
| 910 | 910 | ui32_to_f128M(x, &x_f128); |
| ... | ... | @@ -920,25 +920,25 @@ static float128_t hexfloat(struct MuslFILE *f, int bits, int emin, int sign, int |
| 920 | 920 | f128M_add(&bias_op, &sign_mul_y, &new_y); |
| 921 | 921 | y = new_y; |
| 922 | 922 | } |
| 923 | | //y -= bias; |
| 923 | //y -= bias; |
| 924 | 924 | { |
| 925 | 925 | float128_t new_y; |
| 926 | 926 | f128M_sub(&y, &bias, &new_y); |
| 927 | 927 | y = new_y; |
| 928 | 928 | } |
| 929 | 929 | |
| 930 | | if (f128M_eq(&y, &zero)) errno = ERANGE; |
| 930 | if (f128M_eq(&y, &zero)) errno = ERANGE; |
| 931 | 931 | |
| 932 | | return scalbnf128(y, e2); |
| 932 | return scalbnf128(y, e2); |
| 933 | 933 | } |
| 934 | 934 | |
| 935 | 935 | static int isspace(int c) |
| 936 | 936 | { |
| 937 | | return c == ' ' || (unsigned)c-'\t' < 5; |
| 937 | return c == ' ' || (unsigned)c-'\t' < 5; |
| 938 | 938 | } |
| 939 | 939 | |
| 940 | 940 | static inline float128_t makeInf128() { |
| 941 | | union ldshape ux; |
| 941 | union ldshape ux; |
| 942 | 942 | ux.i2.hi = 0x7fff000000000000UL; |
| 943 | 943 | ux.i2.lo = 0x0UL; |
| 944 | 944 | return ux.f; |
| ... | ... | @@ -946,7 +946,7 @@ static inline float128_t makeInf128() { |
| 946 | 946 | |
| 947 | 947 | static inline float128_t makeNaN128() { |
| 948 | 948 | uint64_t rand = 0UL; |
| 949 | | union ldshape ux; |
| 949 | union ldshape ux; |
| 950 | 950 | ux.i2.hi = 0x7fff000000000000UL | (rand & 0xffffffffffffUL); |
| 951 | 951 | ux.i2.lo = 0x0UL; |
| 952 | 952 | return ux.f; |
| ... | ... | @@ -954,86 +954,86 @@ static inline float128_t makeNaN128() { |
| 954 | 954 | |
| 955 | 955 | float128_t __floatscan(struct MuslFILE *f, int prec, int pok) |
| 956 | 956 | { |
| 957 | | int sign = 1; |
| 958 | | size_t i; |
| 959 | | int bits = LDBL_MANT_DIG; |
| 960 | | int emin = LDBL_MIN_EXP-bits; |
| 961 | | int c; |
| 962 | | |
| 963 | | while (isspace((c=shgetc(f)))); |
| 964 | | |
| 965 | | if (c=='+' || c=='-') { |
| 966 | | sign -= 2*(c=='-'); |
| 967 | | c = shgetc(f); |
| 968 | | } |
| 969 | | |
| 970 | | for (i=0; i<8 && (c|32)=="infinity"[i]; i++) |
| 971 | | if (i<7) c = shgetc(f); |
| 972 | | if (i==3 || i==8 || (i>3 && pok)) { |
| 973 | | if (i!=8) { |
| 974 | | shunget(f); |
| 975 | | if (pok) for (; i>3; i--) shunget(f); |
| 976 | | } |
| 977 | | //return sign * INFINITY; |
| 957 | int sign = 1; |
| 958 | size_t i; |
| 959 | int bits = LDBL_MANT_DIG; |
| 960 | int emin = LDBL_MIN_EXP-bits; |
| 961 | int c; |
| 962 | |
| 963 | while (isspace((c=shgetc(f)))); |
| 964 | |
| 965 | if (c=='+' || c=='-') { |
| 966 | sign -= 2*(c=='-'); |
| 967 | c = shgetc(f); |
| 968 | } |
| 969 | |
| 970 | for (i=0; i<8 && (c|32)=="infinity"[i]; i++) |
| 971 | if (i<7) c = shgetc(f); |
| 972 | if (i==3 || i==8 || (i>3 && pok)) { |
| 973 | if (i!=8) { |
| 974 | shunget(f); |
| 975 | if (pok) for (; i>3; i--) shunget(f); |
| 976 | } |
| 977 | //return sign * INFINITY; |
| 978 | 978 | float128_t sign_f128; |
| 979 | 979 | i32_to_f128M(sign, &sign_f128); |
| 980 | 980 | float128_t infinity_f128 = makeInf128(); |
| 981 | 981 | float128_t result; |
| 982 | 982 | f128M_mul(&sign_f128, &infinity_f128, &result); |
| 983 | 983 | return result; |
| 984 | | } |
| 985 | | if (!i) for (i=0; i<3 && (c|32)=="nan"[i]; i++) |
| 986 | | if (i<2) c = shgetc(f); |
| 987 | | if (i==3) { |
| 988 | | if (shgetc(f) != '(') { |
| 989 | | shunget(f); |
| 990 | | return makeNaN128(); |
| 991 | | } |
| 992 | | for (i=1; ; i++) { |
| 993 | | c = shgetc(f); |
| 994 | | if (c-'0'<10U || c-'A'<26U || c-'a'<26U || c=='_') |
| 995 | | continue; |
| 996 | | if (c==')') return makeNaN128(); |
| 997 | | shunget(f); |
| 998 | | if (!pok) { |
| 999 | | errno = EINVAL; |
| 1000 | | shlim(f, 0); |
| 984 | } |
| 985 | if (!i) for (i=0; i<3 && (c|32)=="nan"[i]; i++) |
| 986 | if (i<2) c = shgetc(f); |
| 987 | if (i==3) { |
| 988 | if (shgetc(f) != '(') { |
| 989 | shunget(f); |
| 990 | return makeNaN128(); |
| 991 | } |
| 992 | for (i=1; ; i++) { |
| 993 | c = shgetc(f); |
| 994 | if (c-'0'<10U || c-'A'<26U || c-'a'<26U || c=='_') |
| 995 | continue; |
| 996 | if (c==')') return makeNaN128(); |
| 997 | shunget(f); |
| 998 | if (!pok) { |
| 999 | errno = EINVAL; |
| 1000 | shlim(f, 0); |
| 1001 | 1001 | float128_t zero; |
| 1002 | 1002 | ui32_to_f128M(0, &zero); |
| 1003 | | return zero; |
| 1004 | | } |
| 1005 | | while (i--) shunget(f); |
| 1006 | | return makeNaN128(); |
| 1007 | | } |
| 1008 | | return makeNaN128(); |
| 1009 | | } |
| 1010 | | |
| 1011 | | if (i) { |
| 1012 | | shunget(f); |
| 1013 | | errno = EINVAL; |
| 1014 | | shlim(f, 0); |
| 1003 | return zero; |
| 1004 | } |
| 1005 | while (i--) shunget(f); |
| 1006 | return makeNaN128(); |
| 1007 | } |
| 1008 | return makeNaN128(); |
| 1009 | } |
| 1010 | |
| 1011 | if (i) { |
| 1012 | shunget(f); |
| 1013 | errno = EINVAL; |
| 1014 | shlim(f, 0); |
| 1015 | 1015 | float128_t zero; |
| 1016 | 1016 | ui32_to_f128M(0, &zero); |
| 1017 | | return zero; |
| 1018 | | } |
| 1019 | | |
| 1020 | | if (c=='0') { |
| 1021 | | c = shgetc(f); |
| 1022 | | if ((c|32) == 'x') |
| 1023 | | return hexfloat(f, bits, emin, sign, pok); |
| 1024 | | shunget(f); |
| 1025 | | c = '0'; |
| 1026 | | } |
| 1027 | | |
| 1028 | | return decfloat(f, c, bits, emin, sign, pok); |
| 1017 | return zero; |
| 1018 | } |
| 1019 | |
| 1020 | if (c=='0') { |
| 1021 | c = shgetc(f); |
| 1022 | if ((c|32) == 'x') |
| 1023 | return hexfloat(f, bits, emin, sign, pok); |
| 1024 | shunget(f); |
| 1025 | c = '0'; |
| 1026 | } |
| 1027 | |
| 1028 | return decfloat(f, c, bits, emin, sign, pok); |
| 1029 | 1029 | } |
| 1030 | 1030 | |
| 1031 | 1031 | float128_t parse_f128(const char *s, char **p) { |
| 1032 | | struct MuslFILE f; |
| 1033 | | sh_fromstring(&f, s); |
| 1034 | | shlim(&f, 0); |
| 1035 | | float128_t y = __floatscan(&f, 2, 1); |
| 1036 | | off_t cnt = shcnt(&f); |
| 1037 | | if (p) *p = cnt ? (char *)s + cnt : (char *)s; |
| 1038 | | return y; |
| 1032 | struct MuslFILE f; |
| 1033 | sh_fromstring(&f, s); |
| 1034 | shlim(&f, 0); |
| 1035 | float128_t y = __floatscan(&f, 2, 1); |
| 1036 | off_t cnt = shcnt(&f); |
| 1037 | if (p) *p = cnt ? (char *)s + cnt : (char *)s; |
| 1038 | return y; |
| 1039 | 1039 | } |
Andrew Kelley <andrew@ziglang.org> 2019-03-23 13:46:50-04:00