nsz Git - libc-test/blob - src/math/gen/mp.c

   1 #include <stdio.h>
   2
   3 #include <stdint.h>
   4 #include <mpfr.h>
   5 #include "gen.h"
   6
   7 static int rmap(int r)
   8 {
   9         switch (r) {
  10         case RN: return MPFR_RNDN;
  11         case RZ: return MPFR_RNDZ;
  12         case RD: return MPFR_RNDD;
  13         case RU: return MPFR_RNDU;
  14         }
  15         return -1;
  16 }
  17
  18 void debug(mpfr_t x)
  19 {
  20         mpfr_out_str(stdout, 10, 0, x, MPFR_RNDN);
  21         printf("\n");
  22 }
  23
  24 void mpsetup()
  25 {
  26         mpfr_set_emin(-1073);
  27         mpfr_set_emax(1024);
  28 }
  29 void mpsetupf()
  30 {
  31         mpfr_set_emin(-148);
  32         mpfr_set_emax(128);
  33 }
  34 #if LDBL_MANT_DIG == 64
  35 void mpsetupl()
  36 {
  37         mpfr_set_emin(-16444);
  38         mpfr_set_emax(16384);
  39 }
  40 #endif
  41
  42 /*
  43 round x into y considering x is already rounded (t = up or down)
  44
  45 only cases where adjustment is done:
  46         x=...|1...0, t=up    -> x=nextbelow(x)
  47         x=...|1...0, t=down  -> x=nextabove(x)
  48 where | is the rounding point, ... is 0 or 1 bit patterns
  49 */
  50
  51 // TODO: adjust(y, 0, 2, RN); when prec is 24 (0 vs 0x1p-149f), special case x=0
  52 static int adjust_round(mpfr_t y, mpfr_t x, int t, int r)
  53 {
  54         mp_limb_t *p, *q;
  55         unsigned xp, yp;
  56         int t2;
  57
  58         xp = mpfr_get_prec(x);
  59         yp = mpfr_get_prec(y);
  60         if (yp >= xp || r != MPFR_RNDN || t == 0 || !mpfr_number_p(x) || mpfr_zero_p(x)) {
  61                 t2 = mpfr_set(y, x, r);
  62                 return t2 ? t2 : t;
  63         }
  64         p = x->_mpfr_d;
  65         yp++;
  66         q = p + (xp + mp_bits_per_limb - 1)/mp_bits_per_limb - (yp + mp_bits_per_limb - 1)/mp_bits_per_limb;
  67         if ((*p & 1 << -xp%mp_bits_per_limb) || !(*q & 1 << -yp%mp_bits_per_limb)) {
  68                 t2 = mpfr_set(y, x, r);
  69                 return t2 ? t2 : t;
  70         }
  71         if (t > 0)
  72                 mpfr_nextbelow(x);
  73         else
  74                 mpfr_nextabove(x);
  75         return mpfr_set(y, x, r);
  76 }
  77
  78 static int adjust(mpfr_t mr, mpfr_t my, int t, int r)
  79 {
  80 //      double d, dn, dp;
  81 //printf("adj %d\n", t);
  82 //debug(my);
  83         t = adjust_round(mr, my, t, r);
  84 //printf("rnd %d\n", t);
  85 //debug(mr);
  86         t = mpfr_subnormalize(mr, t, r);
  87 //printf("sub %d\n", t);
  88 //debug(mr);
  89 //      d = mpfr_get_d(mr, r);
  90 //      dn = nextafter(d, INFINITY);
  91 //      dp = nextafter(d, -INFINITY);
  92 //printf("c\n %.21e %a\n %.21e %a\n %.21e %a\n",d,d,dn,dn,dp,dp);
  93 //      dn = nextafterf(d, INFINITY);
  94 //      dp = nextafterf(d, -INFINITY);
  95 //printf("cf\n %.21e %a\n %.21e %a\n %.21e %a\n",d,d,dn,dn,dp,dp);
  96         return t;
  97 }
  98
  99 // TODO
 100 //static int eflags(mpfr_t mr, mpfr_t my, int t)
 101 static int eflags(int naninput)
 102 {
 103         int i = 0;
 104
 105         if (mpfr_inexflag_p())
 106                 i |= FE_INEXACT;
 107 //      if (mpfr_underflow_p() && (t || mpfr_cmp(mr, my) != 0))
 108         if (mpfr_underflow_p() && i)
 109                 i |= FE_UNDERFLOW;
 110         if (mpfr_overflow_p())
 111                 i |= FE_OVERFLOW;
 112         if (mpfr_divby0_p())
 113                 i |= FE_DIVBYZERO;
 114         if (!naninput && (mpfr_nanflag_p() || mpfr_erangeflag_p()))
 115                 i |= FE_INVALID;
 116         return i;
 117 }
 118
 119 static void genf(struct t *p, mpfr_t my, int t, int r)
 120 {
 121         MPFR_DECL_INIT(mr, 24);
 122         int i;
 123
 124         t = adjust(mr, my, t, r);
 125         p->y = mpfr_get_flt(mr, r);
 126         p->e = eflags(isnan(p->x) || isnan(p->x2) || isnan(p->x3));
 127         i = eulpf(p->y);
 128         if (!isfinite(p->y)) {
 129                 p->dy = 0;
 130         } else if (i < 0) {
 131                 mpfr_div_2si(mr, mr, i, MPFR_RNDN);
 132                 mpfr_div_2si(my, my, i, MPFR_RNDN);
 133                 mpfr_sub(my, mr, my, MPFR_RNDN);
 134                 p->dy = mpfr_get_flt(my, r);
 135         } else {
 136                 mpfr_sub(my, mr, my, MPFR_RNDN);
 137                 mpfr_div_2si(my, my, i, MPFR_RNDN);
 138                 p->dy = mpfr_get_flt(my, r);
 139         }
 140 }
 141
 142 static int mpf1(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, mpfr_rnd_t))
 143 {
 144         int tn;
 145         int r = rmap(p->r);
 146         MPFR_DECL_INIT(mx, 24);
 147         MPFR_DECL_INIT(my, 128);
 148
 149         mpsetupf();
 150         mpfr_clear_flags();
 151         mpfr_set_flt(mx, p->x, MPFR_RNDN);
 152         tn = fmp(my, mx, r);
 153         p->x2 = 0;
 154         genf(p, my, tn, r);
 155         if ((p->e & INEXACT) && nextafterf(p->y, 0) == 0) {
 156                 mpfr_set_emin(-(1<<20));
 157                 tn = fmp(my, mx, r);
 158                 mpfr_mul_2si(my, my, 149, MPFR_RNDN);
 159                 p->dy = scalbnl(p->y, 149) - mpfr_get_ld(my, r);
 160                 mpfr_set_emin(-148);
 161         }
 162         return 0;
 163 }
 164
 165 static int mpf2(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t))
 166 {
 167         int tn;
 168         int r = rmap(p->r);
 169         MPFR_DECL_INIT(mx, 24);
 170         MPFR_DECL_INIT(mx2, 24);
 171         MPFR_DECL_INIT(my, 128);
 172
 173         mpsetupf();
 174         mpfr_clear_flags();
 175         mpfr_set_flt(mx, p->x, MPFR_RNDN);
 176         mpfr_set_flt(mx2, p->x2, MPFR_RNDN);
 177         tn = fmp(my, mx, mx2, r);
 178         genf(p, my, tn, r);
 179         if ((p->e & INEXACT) && nextafterf(p->y, 0) == 0) {
 180                 mpfr_set_emin(-(1<<20));
 181                 tn = fmp(my, mx, mx2, r);
 182                 mpfr_mul_2si(my, my, 149, MPFR_RNDN);
 183                 p->dy = scalbnl(p->y, 149) - mpfr_get_ld(my, r);
 184                 mpfr_set_emin(-148);
 185         }
 186         return 0;
 187 }
 188
 189 static void gend(struct t *p, mpfr_t my, int t, int r)
 190 {
 191         MPFR_DECL_INIT(mr, 53);
 192         int i;
 193
 194         t = adjust(mr, my, t, r);
 195         p->y = mpfr_get_d(mr, r);
 196         p->e = eflags(isnan(p->x) || isnan(p->x2) || isnan(p->x3));
 197         i = eulp(p->y);
 198         if (!isfinite(p->y)) {
 199                 p->dy = 0;
 200         } else if (i < 0) {
 201                 mpfr_div_2si(mr, mr, i, MPFR_RNDN);
 202                 mpfr_div_2si(my, my, i, MPFR_RNDN);
 203                 mpfr_sub(my, mr, my, MPFR_RNDN);
 204                 p->dy = mpfr_get_flt(my, r);
 205         } else {
 206                 mpfr_sub(my, mr, my, MPFR_RNDN);
 207                 mpfr_div_2si(my, my, i, MPFR_RNDN);
 208                 p->dy = mpfr_get_flt(my, r);
 209         }
 210 }
 211
 212 static int mpd1(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, mpfr_rnd_t))
 213 {
 214         int tn;
 215         int r = rmap(p->r);
 216         MPFR_DECL_INIT(mx, 53);
 217         MPFR_DECL_INIT(my, 128);
 218
 219         mpsetup();
 220         mpfr_clear_flags();
 221         mpfr_set_d(mx, p->x, MPFR_RNDN);
 222         tn = fmp(my, mx, r);
 223 //printf("underflow: %d\n", mpfr_underflow_p());
 224         p->x2 = 0;
 225         gend(p, my, tn, r);
 226 //printf("dy: %a  %.3f\n", p->dy, p->dy);
 227         if ((p->e & INEXACT) && nextafter(p->y, 0) == 0) {
 228                 mpfr_set_emin(-(1<<20));
 229                 tn = fmp(my, mx, r);
 230                 mpfr_mul_2si(my, my, 1074, MPFR_RNDN);
 231 //debug(my);
 232                 p->dy = scalbnl(p->y, 1074) - mpfr_get_ld(my, r);
 233 //printf("dy: %a  %.3f\n", p->dy, p->dy);
 234                 mpfr_set_emin(-1073);
 235         }
 236         return 0;
 237 }
 238
 239 static int mpd2(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t))
 240 {
 241         int tn;
 242         int r = rmap(p->r);
 243         MPFR_DECL_INIT(mx, 53);
 244         MPFR_DECL_INIT(mx2, 53);
 245         MPFR_DECL_INIT(my, 128);
 246
 247         mpsetup();
 248         mpfr_clear_flags();
 249         mpfr_set_d(mx, p->x, MPFR_RNDN);
 250         mpfr_set_d(mx2, p->x2, MPFR_RNDN);
 251         tn = fmp(my, mx, mx2, r);
 252         gend(p, my, tn, r);
 253         if ((p->e & INEXACT) && nextafter(p->y, 0) == 0) {
 254                 mpfr_set_emin(-(1<<20));
 255                 tn = fmp(my, mx, mx2, r);
 256                 mpfr_mul_2si(my, my, 1074, MPFR_RNDN);
 257                 p->dy = scalbnl(p->y, 1074) - mpfr_get_ld(my, r);
 258                 mpfr_set_emin(-1073);
 259         }
 260         return 0;
 261 }
 262
 263 #if LDBL_MANT_DIG == 64
 264 static void genl(struct t *p, mpfr_t my, int t, int r)
 265 {
 266         MPFR_DECL_INIT(mr, 64);
 267         int i;
 268
 269         t = adjust(mr, my, t, r);
 270         p->y = mpfr_get_ld(mr, r);
 271         p->e = eflags(isnan(p->x) || isnan(p->x2) || isnan(p->x3));
 272         i = eulpl(p->y);
 273         if (!isfinite(p->y)) {
 274                 p->dy = 0;
 275         } else if (i < 0) {
 276                 mpfr_div_2si(mr, mr, i, MPFR_RNDN);
 277                 mpfr_div_2si(my, my, i, MPFR_RNDN);
 278                 mpfr_sub(my, mr, my, MPFR_RNDN);
 279                 p->dy = mpfr_get_flt(my, r);
 280         } else {
 281                 mpfr_sub(my, mr, my, MPFR_RNDN);
 282                 mpfr_div_2si(my, my, i, MPFR_RNDN);
 283                 p->dy = mpfr_get_flt(my, r);
 284         }
 285 }
 286 #endif
 287
 288 static int mpl1(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, mpfr_rnd_t))
 289 {
 290 #if LDBL_MANT_DIG == 53
 291         return mpd1(p, fmp);
 292 #elif LDBL_MANT_DIG == 64
 293         int tn;
 294         int r = rmap(p->r);
 295         MPFR_DECL_INIT(mx, 64);
 296         MPFR_DECL_INIT(my, 128);
 297
 298         mpsetupl();
 299         mpfr_clear_flags();
 300         mpfr_set_ld(mx, p->x, MPFR_RNDN);
 301         tn = fmp(my, mx, r);
 302         p->x2 = 0;
 303         genl(p, my, tn, r);
 304         if ((p->e & INEXACT) && nextafterl(p->y, 0) == 0) {
 305                 mpfr_set_emin(-(1<<20));
 306                 tn = fmp(my, mx, r);
 307                 mpfr_mul_2si(my, my, 16445, MPFR_RNDN);
 308                 p->dy = scalbnl(p->y, 16445) - mpfr_get_ld(my, r);
 309                 mpfr_set_emin(-16444);
 310         }
 311         return 0;
 312 #else
 313         return -1;
 314 #endif
 315 }
 316
 317 static int mpl2(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t))
 318 {
 319 #if LDBL_MANT_DIG == 53
 320         return mpd2(p, fmp);
 321 #elif LDBL_MANT_DIG == 64
 322         int tn;
 323         int r = rmap(p->r);
 324         MPFR_DECL_INIT(mx, 64);
 325         MPFR_DECL_INIT(mx2, 64);
 326         MPFR_DECL_INIT(my, 128);
 327
 328         mpsetupl();
 329         mpfr_clear_flags();
 330         mpfr_set_ld(mx, p->x, MPFR_RNDN);
 331         mpfr_set_ld(mx2, p->x2, MPFR_RNDN);
 332         tn = fmp(my, mx, mx2, r);
 333         genl(p, my, tn, r);
 334         if ((p->e & INEXACT) && nextafterl(p->y, 0) == 0) {
 335                 mpfr_set_emin(-(1<<20));
 336                 tn = fmp(my, mx, mx2, r);
 337                 mpfr_mul_2si(my, my, 16445, MPFR_RNDN);
 338                 p->dy = scalbnl(p->y, 16445) - mpfr_get_ld(my, r);
 339                 mpfr_set_emin(-16444);
 340         }
 341         return 0;
 342 #else
 343         return -1;
 344 #endif
 345 }
 346
 347 // TODO
 348 static int mplgamma_sign;
 349 static int wrap_lgamma(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 350 {
 351         return mpfr_lgamma(my, &mplgamma_sign, mx, r);
 352 }
 353 static long mpremquo_q;
 354 static int wrap_remquo(mpfr_t my, const mpfr_t mx, const mpfr_t mx2, mpfr_rnd_t r)
 355 {
 356         return mpfr_remquo(my, &mpremquo_q, mx, mx2, r);
 357 }
 358 static int mpbessel_n;
 359 static int wrap_jn(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 360 {
 361         return mpfr_jn(my, mpbessel_n, mx, r);
 362 }
 363 static int wrap_yn(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 364 {
 365         return mpfr_yn(my, mpbessel_n, mx, r);
 366 }
 367 static int wrap_ceil(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 368 {
 369         return mpfr_ceil(my, mx);
 370 }
 371 static int wrap_floor(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 372 {
 373         return mpfr_floor(my, mx);
 374 }
 375 static int wrap_round(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 376 {
 377         return mpfr_round(my, mx);
 378 }
 379 static int wrap_trunc(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 380 {
 381         return mpfr_trunc(my, mx);
 382 }
 383 static int wrap_nearbyint(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 384 {
 385         int i = mpfr_rint(my, mx, r);
 386         mpfr_clear_inexflag();
 387         return i;
 388 }
 389 static int wrap_pow10(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 390 {
 391         return mpfr_ui_pow(my, 10, mx, r);
 392 }
 393
 394 int mpadd(struct t *t) { return mpd2(t, mpfr_add); }
 395
 396 int mpacos(struct t *t) { return mpd1(t, mpfr_acos); }
 397 int mpacosf(struct t *t) { return mpf1(t, mpfr_acos); }
 398 int mpacosl(struct t *t) { return mpl1(t, mpfr_acos); }
 399 int mpacosh(struct t *t) { return mpd1(t, mpfr_acosh); }
 400 int mpacoshf(struct t *t) { return mpf1(t, mpfr_acosh); }
 401 int mpacoshl(struct t *t) { return mpl1(t, mpfr_acosh); }
 402 int mpasin(struct t *t) { return mpd1(t, mpfr_asin); }
 403 int mpasinf(struct t *t) { return mpf1(t, mpfr_asin); }
 404 int mpasinl(struct t *t) { return mpl1(t, mpfr_asin); }
 405 int mpasinh(struct t *t) { return mpd1(t, mpfr_asinh); }
 406 int mpasinhf(struct t *t) { return mpf1(t, mpfr_asinh); }
 407 int mpasinhl(struct t *t) { return mpl1(t, mpfr_asinh); }
 408 int mpatan(struct t *t) { return mpd1(t, mpfr_atan); }
 409 int mpatanf(struct t *t) { return mpf1(t, mpfr_atan); }
 410 int mpatanl(struct t *t) { return mpl1(t, mpfr_atan); }
 411 int mpatan2(struct t *t) { return mpd2(t, mpfr_atan2); }
 412 int mpatan2f(struct t *t) { return mpf2(t, mpfr_atan2); }
 413 int mpatan2l(struct t *t) { return mpl2(t, mpfr_atan2); }
 414 int mpatanh(struct t *t) { return mpd1(t, mpfr_atanh); }
 415 int mpatanhf(struct t *t) { return mpf1(t, mpfr_atanh); }
 416 int mpatanhl(struct t *t) { return mpl1(t, mpfr_atanh); }
 417 int mpcbrt(struct t *t) { return mpd1(t, mpfr_cbrt); }
 418 int mpcbrtf(struct t *t) { return mpf1(t, mpfr_cbrt); }
 419 int mpcbrtl(struct t *t) { return mpl1(t, mpfr_cbrt); }
 420 int mpceil(struct t *t) { return mpd1(t, wrap_ceil); }
 421 int mpceilf(struct t *t) { return mpf1(t, wrap_ceil); }
 422 int mpceill(struct t *t) { return mpl1(t, wrap_ceil); }
 423 int mpcopysign(struct t *t) { return mpd2(t, mpfr_copysign); }
 424 int mpcopysignf(struct t *t) { return mpf2(t, mpfr_copysign); }
 425 int mpcopysignl(struct t *t) { return mpl2(t, mpfr_copysign); }
 426 int mpcos(struct t *t) { return mpd1(t, mpfr_cos); }
 427 int mpcosf(struct t *t) { return mpf1(t, mpfr_cos); }
 428 int mpcosl(struct t *t) { return mpl1(t, mpfr_cos); }
 429 int mpcosh(struct t *t) { return mpd1(t, mpfr_cosh); }
 430 int mpcoshf(struct t *t) { return mpf1(t, mpfr_cosh); }
 431 int mpcoshl(struct t *t) { return mpl1(t, mpfr_cosh); }
 432 int mperf(struct t *t) { return mpd1(t, mpfr_erf); }
 433 int mperff(struct t *t) { return mpf1(t, mpfr_erf); }
 434 int mperfl(struct t *t) { return mpl1(t, mpfr_erf); }
 435 int mperfc(struct t *t) { return mpd1(t, mpfr_erfc); }
 436 int mperfcf(struct t *t) { return mpf1(t, mpfr_erfc); }
 437 int mperfcl(struct t *t) { return mpl1(t, mpfr_erfc); }
 438 int mpexp(struct t *t) { return mpd1(t, mpfr_exp); }
 439 int mpexpf(struct t *t) { return mpf1(t, mpfr_exp); }
 440 int mpexpl(struct t *t) { return mpl1(t, mpfr_exp); }
 441 int mpexp2(struct t *t) { return mpd1(t, mpfr_exp2); }
 442 int mpexp2f(struct t *t) { return mpf1(t, mpfr_exp2); }
 443 int mpexp2l(struct t *t) { return mpl1(t, mpfr_exp2); }
 444 int mpexpm1(struct t *t) { return mpd1(t, mpfr_expm1); }
 445 int mpexpm1f(struct t *t) { return mpf1(t, mpfr_expm1); }
 446 int mpexpm1l(struct t *t) { return mpl1(t, mpfr_expm1); }
 447 int mpfabs(struct t *t) { return mpd1(t, mpfr_abs); }
 448 int mpfabsf(struct t *t) { return mpf1(t, mpfr_abs); }
 449 int mpfabsl(struct t *t) { return mpl1(t, mpfr_abs); }
 450 int mpfdim(struct t *t) { return mpd2(t, mpfr_dim); }
 451 int mpfdimf(struct t *t) { return mpf2(t, mpfr_dim); }
 452 int mpfdiml(struct t *t) { return mpl2(t, mpfr_dim); }
 453 int mpfloor(struct t *t) { return mpd1(t, wrap_floor); }
 454 int mpfloorf(struct t *t) { return mpf1(t, wrap_floor); }
 455 int mpfloorl(struct t *t) { return mpl1(t, wrap_floor); }
 456 int mpfmax(struct t *t) { return mpd2(t, mpfr_max); }
 457 int mpfmaxf(struct t *t) { return mpf2(t, mpfr_max); }
 458 int mpfmaxl(struct t *t) { return mpl2(t, mpfr_max); }
 459 int mpfmin(struct t *t) { return mpd2(t, mpfr_min); }
 460 int mpfminf(struct t *t) { return mpf2(t, mpfr_min); }
 461 int mpfminl(struct t *t) { return mpl2(t, mpfr_min); }
 462 int mpfmod(struct t *t) { return mpd2(t, mpfr_fmod); }
 463 int mpfmodf(struct t *t) { return mpf2(t, mpfr_fmod); }
 464 int mpfmodl(struct t *t) { return mpl2(t, mpfr_fmod); }
 465 int mphypot(struct t *t) { return mpd2(t, mpfr_hypot); }
 466 int mphypotf(struct t *t) { return mpf2(t, mpfr_hypot); }
 467 int mphypotl(struct t *t) { return mpl2(t, mpfr_hypot); }
 468 int mplgamma(struct t *t) { return mpd1(t, wrap_lgamma) || (t->i = mplgamma_sign, 0); }
 469 int mplgammaf(struct t *t) { return mpf1(t, wrap_lgamma) || (t->i = mplgamma_sign, 0); }
 470 int mplgammal(struct t *t) { return mpl1(t, wrap_lgamma) || (t->i = mplgamma_sign, 0); }
 471 int mplog(struct t *t) { return mpd1(t, mpfr_log); }
 472 int mplogf(struct t *t) { return mpf1(t, mpfr_log); }
 473 int mplogl(struct t *t) { return mpl1(t, mpfr_log); }
 474 int mplog10(struct t *t) { return mpd1(t, mpfr_log10); }
 475 int mplog10f(struct t *t) { return mpf1(t, mpfr_log10); }
 476 int mplog10l(struct t *t) { return mpl1(t, mpfr_log10); }
 477 int mplog1p(struct t *t) { return mpd1(t, mpfr_log1p); }
 478 int mplog1pf(struct t *t) { return mpf1(t, mpfr_log1p); }
 479 int mplog1pl(struct t *t) { return mpl1(t, mpfr_log1p); }
 480 int mplog2(struct t *t) { return mpd1(t, mpfr_log2); }
 481 int mplog2f(struct t *t) { return mpf1(t, mpfr_log2); }
 482 int mplog2l(struct t *t) { return mpl1(t, mpfr_log2); }
 483 int mplogb(struct t *t)
 484 {
 485         MPFR_DECL_INIT(mx, 53);
 486
 487         t->dy = 0;
 488         t->e = 0;
 489         if (t->x == 0) {
 490                 t->y = -INFINITY;
 491                 t->e |= DIVBYZERO;
 492                 return 0;
 493         }
 494         if (isinf(t->x)) {
 495                 t->y = INFINITY;
 496                 return 0;
 497         }
 498         if (isnan(t->x)) {
 499                 t->y = t->x;
 500                 return 0;
 501         }
 502         mpfr_set_d(mx, t->x, MPFR_RNDN);
 503         t->y = mpfr_get_exp(mx) - 1;
 504         return 0;
 505 }
 506 int mplogbf(struct t *t)
 507 {
 508         MPFR_DECL_INIT(mx, 24);
 509
 510         t->dy = 0;
 511         t->e = 0;
 512         if (t->x == 0) {
 513                 t->y = -INFINITY;
 514                 t->e |= DIVBYZERO;
 515                 return 0;
 516         }
 517         if (isinf(t->x)) {
 518                 t->y = INFINITY;
 519                 return 0;
 520         }
 521         if (isnan(t->x)) {
 522                 t->y = t->x;
 523                 return 0;
 524         }
 525         mpfr_set_flt(mx, t->x, MPFR_RNDN);
 526         t->y = mpfr_get_exp(mx) - 1;
 527         return 0;
 528 }
 529 int mplogbl(struct t *t)
 530 {
 531         MPFR_DECL_INIT(mx, 64);
 532
 533         t->dy = 0;
 534         t->e = 0;
 535         if (t->x == 0) {
 536                 t->y = -INFINITY;
 537                 t->e |= DIVBYZERO;
 538                 return 0;
 539         }
 540         if (isinf(t->x)) {
 541                 t->y = INFINITY;
 542                 return 0;
 543         }
 544         if (isnan(t->x)) {
 545                 t->y = t->x;
 546                 return 0;
 547         }
 548         mpfr_set_ld(mx, t->x, MPFR_RNDN);
 549         t->y = mpfr_get_exp(mx) - 1;
 550         return 0;
 551 }
 552 int mpnearbyint(struct t *t) { return mpd1(t, wrap_nearbyint) || (t->e&=~INEXACT, 0); }
 553 int mpnearbyintf(struct t *t) { return mpf1(t, wrap_nearbyint) || (t->e&=~INEXACT, 0); }
 554 int mpnearbyintl(struct t *t) { return mpl1(t, wrap_nearbyint) || (t->e&=~INEXACT, 0); }
 555 // TODO: hard to implement with mpfr
 556 int mpnextafter(struct t *t)
 557 {
 558         feclearexcept(FE_ALL_EXCEPT);
 559         t->y = nextafter(t->x, t->x2);
 560         t->e = getexcept();
 561         t->dy = 0;
 562         return 0;
 563 }
 564 int mpnextafterf(struct t *t)
 565 {
 566         feclearexcept(FE_ALL_EXCEPT);
 567         t->y = nextafterf(t->x, t->x2);
 568         t->e = getexcept();
 569         t->dy = 0;
 570         return 0;
 571 }
 572 int mpnextafterl(struct t *t)
 573 {
 574         feclearexcept(FE_ALL_EXCEPT);
 575         t->y = nextafterl(t->x, t->x2);
 576         t->e = getexcept();
 577         t->dy = 0;
 578         return 0;
 579 }
 580 int mpnexttoward(struct t *t)
 581 {
 582         feclearexcept(FE_ALL_EXCEPT);
 583         t->y = nexttoward(t->x, t->x2);
 584         t->e = getexcept();
 585         t->dy = 0;
 586         return 0;
 587 }
 588 int mpnexttowardf(struct t *t)
 589 {
 590         feclearexcept(FE_ALL_EXCEPT);
 591         t->y = nexttowardf(t->x, t->x2);
 592         t->e = getexcept();
 593         t->dy = 0;
 594         return 0;
 595 }
 596 int mpnexttowardl(struct t *t) { return mpnextafterl(t); }
 597 int mppow(struct t *t) { return mpd2(t, mpfr_pow); }
 598 int mppowf(struct t *t) { return mpf2(t, mpfr_pow); }
 599 int mppowl(struct t *t) { return mpl2(t, mpfr_pow); }
 600 int mpremainder(struct t *t) { return mpd2(t, mpfr_remainder); }
 601 int mpremainderf(struct t *t) { return mpf2(t, mpfr_remainder); }
 602 int mpremainderl(struct t *t) { return mpl2(t, mpfr_remainder); }
 603 int mprint(struct t *t) { return mpd1(t, mpfr_rint); }
 604 int mprintf(struct t *t) { return mpf1(t, mpfr_rint); }
 605 int mprintl(struct t *t) { return mpl1(t, mpfr_rint); }
 606 int mpround(struct t *t) { return mpd1(t, wrap_round); }
 607 int mproundf(struct t *t) { return mpf1(t, wrap_round); }
 608 int mproundl(struct t *t) { return mpl1(t, wrap_round); }
 609 int mpsin(struct t *t) { return mpd1(t, mpfr_sin); }
 610 int mpsinf(struct t *t) { return mpf1(t, mpfr_sin); }
 611 int mpsinl(struct t *t) { return mpl1(t, mpfr_sin); }
 612 int mpsinh(struct t *t) { return mpd1(t, mpfr_sinh); }
 613 int mpsinhf(struct t *t) { return mpf1(t, mpfr_sinh); }
 614 int mpsinhl(struct t *t) { return mpl1(t, mpfr_sinh); }
 615 int mpsqrt(struct t *t) { return mpd1(t, mpfr_sqrt); }
 616 int mpsqrtf(struct t *t) { return mpf1(t, mpfr_sqrt); }
 617 int mpsqrtl(struct t *t) { return mpl1(t, mpfr_sqrt); }
 618 int mptan(struct t *t) { return mpd1(t, mpfr_tan); }
 619 int mptanf(struct t *t) { return mpf1(t, mpfr_tan); }
 620 int mptanl(struct t *t) { return mpl1(t, mpfr_tan); }
 621 int mptanh(struct t *t) { return mpd1(t, mpfr_tanh); }
 622 int mptanhf(struct t *t) { return mpf1(t, mpfr_tanh); }
 623 int mptanhl(struct t *t) { return mpl1(t, mpfr_tanh); }
 624 // TODO: tgamma(2) raises wrong flags
 625 int mptgamma(struct t *t) { return mpd1(t, mpfr_gamma); }
 626 int mptgammaf(struct t *t) { return mpf1(t, mpfr_gamma); }
 627 int mptgammal(struct t *t) { return mpl1(t, mpfr_gamma); }
 628 int mptrunc(struct t *t) { return mpd1(t, wrap_trunc); }
 629 int mptruncf(struct t *t) { return mpf1(t, wrap_trunc); }
 630 int mptruncl(struct t *t) { return mpl1(t, wrap_trunc); }
 631 int mpj0(struct t *t) { return mpd1(t, mpfr_j0); }
 632 int mpj1(struct t *t) { return mpd1(t, mpfr_j1); }
 633 int mpy0(struct t *t) { return mpd1(t, mpfr_y0); }
 634 int mpy1(struct t *t) { return mpd1(t, mpfr_y1); }
 635 // TODO: non standard functions
 636 int mpscalb(struct t *t)
 637 {
 638         setupfenv(t->r);
 639         t->y = scalb(t->x, t->x2);
 640         t->e = getexcept();
 641         t->dy = 0; // wrong
 642         return 0;
 643 }
 644 int mpscalbf(struct t *t)
 645 {
 646         setupfenv(t->r);
 647         t->y = scalbf(t->x, t->x2);
 648         t->e = getexcept();
 649         t->dy = 0; // wrong
 650         return 0;
 651 }
 652 int mpj0f(struct t *t) { return mpf1(t, mpfr_j0); }
 653 int mpj0l(struct t *t) { return mpl1(t, mpfr_j0); }
 654 int mpj1f(struct t *t) { return mpf1(t, mpfr_j1); }
 655 int mpj1l(struct t *t) { return mpl1(t, mpfr_j1); }
 656 int mpy0f(struct t *t) { return mpf1(t, mpfr_y0); }
 657 int mpy0l(struct t *t) { return mpl1(t, mpfr_y0); }
 658 int mpy1f(struct t *t) { return mpf1(t, mpfr_y1); }
 659 int mpy1l(struct t *t) { return mpl1(t, mpfr_y1); }
 660 int mpexp10(struct t *t) { return mpd1(t, wrap_pow10); }
 661 int mpexp10f(struct t *t) { return mpf1(t, wrap_pow10); }
 662 int mpexp10l(struct t *t) { return mpl1(t, wrap_pow10); }
 663 int mppow10(struct t *t) { return mpd1(t, wrap_pow10); }
 664 int mppow10f(struct t *t) { return mpf1(t, wrap_pow10); }
 665 int mppow10l(struct t *t) { return mpl1(t, wrap_pow10); }
 666
 667 int mpfrexp(struct t *t)
 668 {
 669         mpfr_exp_t e;
 670         int k;
 671         MPFR_DECL_INIT(mx, 53);
 672
 673         t->dy = 0;
 674         t->y = 0;
 675         mpsetup();
 676         mpfr_clear_flags();
 677         mpfr_set_d(mx, t->x, MPFR_RNDN);
 678         k = mpfr_frexp(&e, mx, mx, t->r);
 679         mpfr_subnormalize(mx, k, t->r);
 680         t->y = mpfr_get_d(mx, MPFR_RNDN);
 681         t->i = e;
 682         t->e = eflags(isnan(t->x));
 683         return 0;
 684 }
 685
 686 int mpfrexpf(struct t *t)
 687 {
 688         mpfr_exp_t e;
 689         int k;
 690         MPFR_DECL_INIT(mx, 24);
 691
 692         t->dy = 0;
 693         t->y = 0;
 694         mpsetup();
 695         mpfr_clear_flags();
 696         mpfr_set_flt(mx, t->x, MPFR_RNDN);
 697         k = mpfr_frexp(&e, mx, mx, t->r);
 698         mpfr_subnormalize(mx, k, t->r);
 699         t->y = mpfr_get_flt(mx, MPFR_RNDN);
 700         t->i = e;
 701         t->e = eflags(isnan(t->x));
 702         return 0;
 703 }
 704
 705 int mpfrexpl(struct t *t)
 706 {
 707         mpfr_exp_t e;
 708         int k;
 709         MPFR_DECL_INIT(mx, 64);
 710
 711         t->dy = 0;
 712         t->y = 0;
 713         mpsetup();
 714         mpfr_clear_flags();
 715         mpfr_set_ld(mx, t->x, MPFR_RNDN);
 716         k = mpfr_frexp(&e, mx, mx, t->r);
 717         mpfr_subnormalize(mx, k, t->r);
 718         t->y = mpfr_get_ld(mx, MPFR_RNDN);
 719         t->i = e;
 720         t->e = eflags(isnan(t->x));
 721         return 0;
 722 }
 723
 724 int mpldexp(struct t *t)
 725 {
 726         int k;
 727         MPFR_DECL_INIT(mx, 53);
 728
 729         t->dy = 0;
 730         t->y = 0;
 731         mpsetup();
 732         mpfr_clear_flags();
 733         mpfr_set_d(mx, t->x, MPFR_RNDN);
 734         k = mpfr_mul_2si(mx, mx, t->i, t->r);
 735         mpfr_subnormalize(mx, k, t->r);
 736         t->y = mpfr_get_d(mx, MPFR_RNDN);
 737         t->e = eflags(isnan(t->x));
 738         return 0;
 739 }
 740
 741 int mpldexpf(struct t *t)
 742 {
 743         int k;
 744         MPFR_DECL_INIT(mx, 24);
 745
 746         t->dy = 0;
 747         t->y = 0;
 748         mpsetup();
 749         mpfr_clear_flags();
 750         mpfr_set_flt(mx, t->x, MPFR_RNDN);
 751         k = mpfr_mul_2si(mx, mx, t->i, t->r);
 752         mpfr_subnormalize(mx, k, t->r);
 753         t->y = mpfr_get_flt(mx, MPFR_RNDN);
 754         t->e = eflags(isnan(t->x));
 755         return 0;
 756 }
 757
 758 int mpldexpl(struct t *t)
 759 {
 760         int k;
 761         MPFR_DECL_INIT(mx, 64);
 762
 763         t->dy = 0;
 764         t->y = 0;
 765         mpsetup();
 766         mpfr_clear_flags();
 767         mpfr_set_ld(mx, t->x, MPFR_RNDN);
 768         k = mpfr_mul_2si(mx, mx, t->i, t->r);
 769         mpfr_subnormalize(mx, k, t->r);
 770         t->y = mpfr_get_ld(mx, MPFR_RNDN);
 771         t->e = eflags(isnan(t->x));
 772         return 0;
 773 }
 774
 775 int mpscalbn(struct t *t) { return mpldexp(t); }
 776 int mpscalbnf(struct t *t) { return mpldexpf(t); }
 777 int mpscalbnl(struct t *t) { return mpldexpl(t); }
 778 int mpscalbln(struct t *t) { return mpldexp(t); }
 779 int mpscalblnf(struct t *t) { return mpldexpf(t); }
 780 int mpscalblnl(struct t *t) { return mpldexpl(t); }
 781
 782 int mplgamma_r(struct t *t) { return mplgamma(t); }
 783 int mplgammaf_r(struct t *t) { return mplgammaf(t); }
 784 int mplgammal_r(struct t *t) { return mplgammal(t); }
 785
 786 int mpilogb(struct t *t)
 787 {
 788         MPFR_DECL_INIT(mx, 53);
 789
 790         mpfr_set_d(mx, t->x, MPFR_RNDN);
 791         t->i = mpfr_get_exp(mx) - 1;
 792         t->e = 0;
 793         if (isinf(t->x) || isnan(t->x) || t->x == 0)
 794                 t->e = INVALID;
 795         return 0;
 796 }
 797 int mpilogbf(struct t *t)
 798 {
 799         MPFR_DECL_INIT(mx, 24);
 800
 801         mpfr_set_flt(mx, t->x, MPFR_RNDN);
 802         t->i = mpfr_get_exp(mx) - 1;
 803         t->e = 0;
 804         if (isinf(t->x) || isnan(t->x) || t->x == 0)
 805                 t->e = INVALID;
 806         return 0;
 807 }
 808 int mpilogbl(struct t *t)
 809 {
 810         MPFR_DECL_INIT(mx, 64);
 811
 812         mpfr_set_ld(mx, t->x, MPFR_RNDN);
 813         t->i = mpfr_get_exp(mx) - 1;
 814         t->e = 0;
 815         if (isinf(t->x) || isnan(t->x) || t->x == 0)
 816                 t->e = INVALID;
 817         return 0;
 818 }
 819
 820 // TODO: ll* is hard to do with mpfr
 821 #define mp_f_i(n) \
 822 int mp##n(struct t *t) \
 823 { \
 824         setupfenv(t->r); \
 825         t->i = n(t->x); \
 826         t->e = getexcept(); \
 827         return 0; \
 828 }
 829
 830 mp_f_i(llrint)
 831 mp_f_i(llrintf)
 832 mp_f_i(llrintl)
 833 mp_f_i(lrint)
 834 mp_f_i(lrintf)
 835 mp_f_i(lrintl)
 836 mp_f_i(llround)
 837 mp_f_i(llroundf)
 838 mp_f_i(llroundl)
 839 mp_f_i(lround)
 840 mp_f_i(lroundf)
 841 mp_f_i(lroundl)
 842
 843 int mpmodf(struct t *t)
 844 {
 845         int e, r;
 846
 847         r = mpd1(t, wrap_trunc);
 848         if (r)
 849                 return r;
 850         t->y2 = t->y;
 851         t->dy2 = t->dy;
 852         e = t->e & ~INEXACT;
 853         r = mpd1(t, mpfr_frac);
 854         t->e |= e;
 855         return r;
 856 }
 857
 858 int mpmodff(struct t *t)
 859 {
 860         int e, r;
 861
 862         r = mpf1(t, wrap_trunc);
 863         if (r)
 864                 return r;
 865         t->y2 = t->y;
 866         t->dy2 = t->dy;
 867         e = t->e & ~INEXACT;
 868         r = mpf1(t, mpfr_frac);
 869         t->e |= e;
 870         return r;
 871 }
 872
 873 int mpmodfl(struct t *t)
 874 {
 875         int e, r;
 876
 877         r = mpl1(t, wrap_trunc);
 878         if (r)
 879                 return r;
 880         t->y2 = t->y;
 881         t->dy2 = t->dy;
 882         e = t->e & ~INEXACT;
 883         r = mpl1(t, mpfr_frac);
 884         t->e |= e;
 885         return r;
 886 }
 887
 888 int mpsincos(struct t *t)
 889 {
 890         int e, r;
 891
 892         r = mpd1(t, mpfr_cos);
 893         if (r)
 894                 return r;
 895         t->y2 = t->y;
 896         t->dy2 = t->dy;
 897         e = t->e;
 898         r = mpd1(t, mpfr_sin);
 899         t->e |= e;
 900         return r;
 901 }
 902
 903 int mpsincosf(struct t *t)
 904 {
 905         int e, r;
 906
 907         r = mpf1(t, mpfr_cos);
 908         if (r)
 909                 return r;
 910         t->y2 = t->y;
 911         t->dy2 = t->dy;
 912         e = t->e;
 913         r = mpf1(t, mpfr_sin);
 914         t->e |= e;
 915         return r;
 916 }
 917
 918 int mpsincosl(struct t *t)
 919 {
 920         int e, r;
 921
 922         r = mpl1(t, mpfr_cos);
 923         if (r)
 924                 return r;
 925         t->y2 = t->y;
 926         t->dy2 = t->dy;
 927         e = t->e;
 928         r = mpl1(t, mpfr_sin);
 929         t->e |= e;
 930         return r;
 931 }
 932
 933 int mpremquo(struct t *t) { return mpd2(t, wrap_remquo) || (t->i = mpremquo_q, 0); }
 934 int mpremquof(struct t *t) { return mpf2(t, wrap_remquo) || (t->i = mpremquo_q, 0); }
 935 int mpremquol(struct t *t) { return mpl2(t, wrap_remquo) || (t->i = mpremquo_q, 0); }
 936
 937 int mpfma(struct t *t)
 938 {
 939         int tn;
 940         int r = rmap(t->r);
 941         MPFR_DECL_INIT(mx, 53);
 942         MPFR_DECL_INIT(mx2, 53);
 943         MPFR_DECL_INIT(mx3, 53);
 944         MPFR_DECL_INIT(my, 128);
 945
 946         mpsetup();
 947         mpfr_clear_flags();
 948         mpfr_set_d(mx, t->x, MPFR_RNDN);
 949         mpfr_set_d(mx2, t->x2, MPFR_RNDN);
 950         mpfr_set_d(mx3, t->x3, MPFR_RNDN);
 951         tn = mpfr_fma(my, mx, mx2, mx3, r);
 952         gend(t, my, tn, r);
 953         if ((t->e & INEXACT) && nextafter(t->y, 0) == 0) {
 954                 mpfr_set_emin(-(1<<20));
 955                 tn = mpfr_fma(my, mx, mx2, mx3, r);
 956                 mpfr_mul_2si(my, my, 1074, MPFR_RNDN);
 957                 t->dy = scalbnl(t->y, 1074) - mpfr_get_ld(my, r);
 958                 mpfr_set_emin(-1073);
 959         }
 960         return 0;
 961 }
 962
 963 int mpfmaf(struct t *t)
 964 {
 965         int tn;
 966         int r = rmap(t->r);
 967         MPFR_DECL_INIT(mx, 24);
 968         MPFR_DECL_INIT(mx2, 24);
 969         MPFR_DECL_INIT(mx3, 24);
 970         MPFR_DECL_INIT(my, 128);
 971
 972         mpsetupf();
 973         mpfr_clear_flags();
 974         mpfr_set_flt(mx, t->x, MPFR_RNDN);
 975         mpfr_set_flt(mx2, t->x2, MPFR_RNDN);
 976         mpfr_set_flt(mx3, t->x3, MPFR_RNDN);
 977         tn = mpfr_fma(my, mx, mx2, mx3, r);
 978         genf(t, my, tn, r);
 979         if ((t->e & INEXACT) && nextafterf(t->y, 0) == 0) {
 980                 mpfr_set_emin(-(1<<20));
 981                 tn = mpfr_fma(my, mx, mx2, mx3, r);
 982                 mpfr_mul_2si(my, my, 149, MPFR_RNDN);
 983                 t->dy = scalbnl(t->y, 149) - mpfr_get_ld(my, r);
 984                 mpfr_set_emin(-148);
 985         }
 986         return 0;
 987 }
 988
 989 int mpfmal(struct t *t)
 990 {
 991 #if LDBL_MANT_DIG == 53
 992         return mpfma(t);
 993 #elif LDBL_MANT_DIG == 64
 994         int tn;
 995         int r = rmap(t->r);
 996         MPFR_DECL_INIT(mx, 64);
 997         MPFR_DECL_INIT(mx2, 64);
 998         MPFR_DECL_INIT(mx3, 64);
 999         MPFR_DECL_INIT(my, 128);
1000
1001         mpsetupl();
1002         mpfr_clear_flags();
1003         mpfr_set_ld(mx, t->x, MPFR_RNDN);
1004         mpfr_set_ld(mx2, t->x2, MPFR_RNDN);
1005         mpfr_set_ld(mx3, t->x3, MPFR_RNDN);
1006         tn = mpfr_fma(my, mx, mx2, mx3, r);
1007         genl(t, my, tn, r);
1008         if ((t->e & INEXACT) && nextafterl(t->y, 0) == 0) {
1009                 mpfr_set_emin(-(1<<20));
1010                 tn = mpfr_fma(my, mx, mx2, mx3, r);
1011                 mpfr_mul_2si(my, my, 16445, MPFR_RNDN);
1012                 t->dy = scalbnl(t->y, 16445) - mpfr_get_ld(my, r);
1013                 mpfr_set_emin(-16444);
1014         }
1015         return 0;
1016 #else
1017         return -1;
1018 #endif
1019 }
1020
1021 int mpjn(struct t *t) { mpbessel_n = t->i; return mpd1(t, wrap_jn); }
1022 int mpjnf(struct t *t) { mpbessel_n = t->i; return mpf1(t, wrap_jn); }
1023 int mpjnl(struct t *t) { mpbessel_n = t->i; return mpl1(t, wrap_jn); }
1024 int mpyn(struct t *t) { mpbessel_n = t->i; return mpd1(t, wrap_yn); }
1025 int mpynf(struct t *t) { mpbessel_n = t->i; return mpf1(t, wrap_yn); }
1026 int mpynl(struct t *t) { mpbessel_n = t->i; return mpl1(t, wrap_yn); }
1027