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

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