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 mpmul(struct t *t) { return mpd2(t, mpfr_mul); }
 354 int mpdiv(struct t *t) { return mpd2(t, mpfr_div); }
 355
 356 int mpacos(struct t *t) { return mpd1(t, mpfr_acos); }
 357 int mpacosf(struct t *t) { return mpf1(t, mpfr_acos); }
 358 int mpacosl(struct t *t) { return mpl1(t, mpfr_acos); }
 359 int mpacosh(struct t *t) { return mpd1(t, mpfr_acosh); }
 360 int mpacoshf(struct t *t) { return mpf1(t, mpfr_acosh); }
 361 int mpacoshl(struct t *t) { return mpl1(t, mpfr_acosh); }
 362 int mpasin(struct t *t) { return mpd1(t, mpfr_asin); }
 363 int mpasinf(struct t *t) { return mpf1(t, mpfr_asin); }
 364 int mpasinl(struct t *t) { return mpl1(t, mpfr_asin); }
 365 int mpasinh(struct t *t) { return mpd1(t, mpfr_asinh); }
 366 int mpasinhf(struct t *t) { return mpf1(t, mpfr_asinh); }
 367 int mpasinhl(struct t *t) { return mpl1(t, mpfr_asinh); }
 368 int mpatan(struct t *t) { return mpd1(t, mpfr_atan); }
 369 int mpatanf(struct t *t) { return mpf1(t, mpfr_atan); }
 370 int mpatanl(struct t *t) { return mpl1(t, mpfr_atan); }
 371 int mpatan2(struct t *t) { return mpd2(t, mpfr_atan2); }
 372 int mpatan2f(struct t *t) { return mpf2(t, mpfr_atan2); }
 373 int mpatan2l(struct t *t) { return mpl2(t, mpfr_atan2); }
 374 int mpatanh(struct t *t) { return mpd1(t, mpfr_atanh); }
 375 int mpatanhf(struct t *t) { return mpf1(t, mpfr_atanh); }
 376 int mpatanhl(struct t *t) { return mpl1(t, mpfr_atanh); }
 377 int mpcbrt(struct t *t) { return mpd1(t, mpfr_cbrt); }
 378 int mpcbrtf(struct t *t) { return mpf1(t, mpfr_cbrt); }
 379 int mpcbrtl(struct t *t) { return mpl1(t, mpfr_cbrt); }
 380 int mpceil(struct t *t) { return mpd1(t, wrap_ceil); }
 381 int mpceilf(struct t *t) { return mpf1(t, wrap_ceil); }
 382 int mpceill(struct t *t) { return mpl1(t, wrap_ceil); }
 383 int mpcopysign(struct t *t) { return mpd2(t, mpfr_copysign); }
 384 int mpcopysignf(struct t *t) { return mpf2(t, mpfr_copysign); }
 385 int mpcopysignl(struct t *t) { return mpl2(t, mpfr_copysign); }
 386 int mpcos(struct t *t) { return mpd1(t, mpfr_cos); }
 387 int mpcosf(struct t *t) { return mpf1(t, mpfr_cos); }
 388 int mpcosl(struct t *t) { return mpl1(t, mpfr_cos); }
 389 int mpcosh(struct t *t) { return mpd1(t, mpfr_cosh); }
 390 int mpcoshf(struct t *t) { return mpf1(t, mpfr_cosh); }
 391 int mpcoshl(struct t *t) { return mpl1(t, mpfr_cosh); }
 392 int mperf(struct t *t) { return mpd1(t, mpfr_erf); }
 393 int mperff(struct t *t) { return mpf1(t, mpfr_erf); }
 394 int mperfl(struct t *t) { return mpl1(t, mpfr_erf); }
 395 int mperfc(struct t *t) { return mpd1(t, mpfr_erfc); }
 396 int mperfcf(struct t *t) { return mpf1(t, mpfr_erfc); }
 397 int mperfcl(struct t *t) { return mpl1(t, mpfr_erfc); }
 398 int mpexp(struct t *t) { return mpd1(t, mpfr_exp); }
 399 int mpexpf(struct t *t) { return mpf1(t, mpfr_exp); }
 400 int mpexpl(struct t *t) { return mpl1(t, mpfr_exp); }
 401 int mpexp2(struct t *t) { return mpd1(t, mpfr_exp2); }
 402 int mpexp2f(struct t *t) { return mpf1(t, mpfr_exp2); }
 403 int mpexp2l(struct t *t) { return mpl1(t, mpfr_exp2); }
 404 int mpexpm1(struct t *t) { return mpd1(t, mpfr_expm1); }
 405 int mpexpm1f(struct t *t) { return mpf1(t, mpfr_expm1); }
 406 int mpexpm1l(struct t *t) { return mpl1(t, mpfr_expm1); }
 407 int mpfabs(struct t *t) { return mpd1(t, mpfr_abs); }
 408 int mpfabsf(struct t *t) { return mpf1(t, mpfr_abs); }
 409 int mpfabsl(struct t *t) { return mpl1(t, mpfr_abs); }
 410 int mpfdim(struct t *t) { return mpd2(t, mpfr_dim); }
 411 int mpfdimf(struct t *t) { return mpf2(t, mpfr_dim); }
 412 int mpfdiml(struct t *t) { return mpl2(t, mpfr_dim); }
 413 int mpfloor(struct t *t) { return mpd1(t, wrap_floor); }
 414 int mpfloorf(struct t *t) { return mpf1(t, wrap_floor); }
 415 int mpfloorl(struct t *t) { return mpl1(t, wrap_floor); }
 416 int mpfmax(struct t *t) { return mpd2(t, mpfr_max); }
 417 int mpfmaxf(struct t *t) { return mpf2(t, mpfr_max); }
 418 int mpfmaxl(struct t *t) { return mpl2(t, mpfr_max); }
 419 int mpfmin(struct t *t) { return mpd2(t, mpfr_min); }
 420 int mpfminf(struct t *t) { return mpf2(t, mpfr_min); }
 421 int mpfminl(struct t *t) { return mpl2(t, mpfr_min); }
 422 int mpfmod(struct t *t) { return mpd2(t, mpfr_fmod); }
 423 int mpfmodf(struct t *t) { return mpf2(t, mpfr_fmod); }
 424 int mpfmodl(struct t *t) { return mpl2(t, mpfr_fmod); }
 425 int mphypot(struct t *t) { return mpd2(t, mpfr_hypot); }
 426 int mphypotf(struct t *t) { return mpf2(t, mpfr_hypot); }
 427 int mphypotl(struct t *t) { return mpl2(t, mpfr_hypot); }
 428 int mplgamma(struct t *t) { return mpd1(t, wrap_lgamma) || (t->i = mplgamma_sign, 0); }
 429 int mplgammaf(struct t *t) { return mpf1(t, wrap_lgamma) || (t->i = mplgamma_sign, 0); }
 430 int mplgammal(struct t *t) { return mpl1(t, wrap_lgamma) || (t->i = mplgamma_sign, 0); }
 431 int mplog(struct t *t) { return mpd1(t, mpfr_log); }
 432 int mplogf(struct t *t) { return mpf1(t, mpfr_log); }
 433 int mplogl(struct t *t) { return mpl1(t, mpfr_log); }
 434 int mplog10(struct t *t) { return mpd1(t, mpfr_log10); }
 435 int mplog10f(struct t *t) { return mpf1(t, mpfr_log10); }
 436 int mplog10l(struct t *t) { return mpl1(t, mpfr_log10); }
 437 int mplog1p(struct t *t) { return mpd1(t, mpfr_log1p); }
 438 int mplog1pf(struct t *t) { return mpf1(t, mpfr_log1p); }
 439 int mplog1pl(struct t *t) { return mpl1(t, mpfr_log1p); }
 440 int mplog2(struct t *t) { return mpd1(t, mpfr_log2); }
 441 int mplog2f(struct t *t) { return mpf1(t, mpfr_log2); }
 442 int mplog2l(struct t *t) { return mpl1(t, mpfr_log2); }
 443 int mplogb(struct t *t)
 444 {
 445         MPFR_DECL_INIT(mx, 53);
 446
 447         t->dy = 0;
 448         t->e = 0;
 449         if (t->x == 0) {
 450                 t->y = -INFINITY;
 451                 t->e |= DIVBYZERO;
 452                 return 0;
 453         }
 454         if (isinf(t->x)) {
 455                 t->y = INFINITY;
 456                 return 0;
 457         }
 458         if (isnan(t->x)) {
 459                 t->y = t->x;
 460                 return 0;
 461         }
 462         mpfr_set_d(mx, t->x, MPFR_RNDN);
 463         t->y = mpfr_get_exp(mx) - 1;
 464         return 0;
 465 }
 466 int mplogbf(struct t *t)
 467 {
 468         MPFR_DECL_INIT(mx, 24);
 469
 470         t->dy = 0;
 471         t->e = 0;
 472         if (t->x == 0) {
 473                 t->y = -INFINITY;
 474                 t->e |= DIVBYZERO;
 475                 return 0;
 476         }
 477         if (isinf(t->x)) {
 478                 t->y = INFINITY;
 479                 return 0;
 480         }
 481         if (isnan(t->x)) {
 482                 t->y = t->x;
 483                 return 0;
 484         }
 485         mpfr_set_flt(mx, t->x, MPFR_RNDN);
 486         t->y = mpfr_get_exp(mx) - 1;
 487         return 0;
 488 }
 489 int mplogbl(struct t *t)
 490 {
 491         MPFR_DECL_INIT(mx, 64);
 492
 493         t->dy = 0;
 494         t->e = 0;
 495         if (t->x == 0) {
 496                 t->y = -INFINITY;
 497                 t->e |= DIVBYZERO;
 498                 return 0;
 499         }
 500         if (isinf(t->x)) {
 501                 t->y = INFINITY;
 502                 return 0;
 503         }
 504         if (isnan(t->x)) {
 505                 t->y = t->x;
 506                 return 0;
 507         }
 508         mpfr_set_ld(mx, t->x, MPFR_RNDN);
 509         t->y = mpfr_get_exp(mx) - 1;
 510         return 0;
 511 }
 512 int mpnearbyint(struct t *t) { return mpd1(t, wrap_nearbyint) || (t->e&=~INEXACT, 0); }
 513 int mpnearbyintf(struct t *t) { return mpf1(t, wrap_nearbyint) || (t->e&=~INEXACT, 0); }
 514 int mpnearbyintl(struct t *t) { return mpl1(t, wrap_nearbyint) || (t->e&=~INEXACT, 0); }
 515 // TODO: hard to implement with mpfr
 516 int mpnextafter(struct t *t)
 517 {
 518         feclearexcept(FE_ALL_EXCEPT);
 519         t->y = nextafter(t->x, t->x2);
 520         t->e = getexcept();
 521         t->dy = 0;
 522         return 0;
 523 }
 524 int mpnextafterf(struct t *t)
 525 {
 526         feclearexcept(FE_ALL_EXCEPT);
 527         t->y = nextafterf(t->x, t->x2);
 528         t->e = getexcept();
 529         t->dy = 0;
 530         return 0;
 531 }
 532 int mpnextafterl(struct t *t)
 533 {
 534         feclearexcept(FE_ALL_EXCEPT);
 535         t->y = nextafterl(t->x, t->x2);
 536         t->e = getexcept();
 537         t->dy = 0;
 538         return 0;
 539 }
 540 int mpnexttoward(struct t *t)
 541 {
 542         feclearexcept(FE_ALL_EXCEPT);
 543         t->y = nexttoward(t->x, t->x2);
 544         t->e = getexcept();
 545         t->dy = 0;
 546         return 0;
 547 }
 548 int mpnexttowardf(struct t *t)
 549 {
 550         feclearexcept(FE_ALL_EXCEPT);
 551         t->y = nexttowardf(t->x, t->x2);
 552         t->e = getexcept();
 553         t->dy = 0;
 554         return 0;
 555 }
 556 int mpnexttowardl(struct t *t) { return mpnextafterl(t); }
 557 int mppow(struct t *t) { return mpd2(t, mpfr_pow); }
 558 int mppowf(struct t *t) { return mpf2(t, mpfr_pow); }
 559 int mppowl(struct t *t) { return mpl2(t, mpfr_pow); }
 560 int mpremainder(struct t *t) { return mpd2(t, mpfr_remainder); }
 561 int mpremainderf(struct t *t) { return mpf2(t, mpfr_remainder); }
 562 int mpremainderl(struct t *t) { return mpl2(t, mpfr_remainder); }
 563 int mprint(struct t *t) { return mpd1(t, mpfr_rint); }
 564 int mprintf(struct t *t) { return mpf1(t, mpfr_rint); }
 565 int mprintl(struct t *t) { return mpl1(t, mpfr_rint); }
 566 int mpround(struct t *t) { return mpd1(t, wrap_round); }
 567 int mproundf(struct t *t) { return mpf1(t, wrap_round); }
 568 int mproundl(struct t *t) { return mpl1(t, wrap_round); }
 569 int mpsin(struct t *t) { return mpd1(t, mpfr_sin); }
 570 int mpsinf(struct t *t) { return mpf1(t, mpfr_sin); }
 571 int mpsinl(struct t *t) { return mpl1(t, mpfr_sin); }
 572 int mpsinh(struct t *t) { return mpd1(t, mpfr_sinh); }
 573 int mpsinhf(struct t *t) { return mpf1(t, mpfr_sinh); }
 574 int mpsinhl(struct t *t) { return mpl1(t, mpfr_sinh); }
 575 int mpsqrt(struct t *t) { return mpd1(t, mpfr_sqrt); }
 576 int mpsqrtf(struct t *t) { return mpf1(t, mpfr_sqrt); }
 577 int mpsqrtl(struct t *t) { return mpl1(t, mpfr_sqrt); }
 578 int mptan(struct t *t) { return mpd1(t, mpfr_tan); }
 579 int mptanf(struct t *t) { return mpf1(t, mpfr_tan); }
 580 int mptanl(struct t *t) { return mpl1(t, mpfr_tan); }
 581 int mptanh(struct t *t) { return mpd1(t, mpfr_tanh); }
 582 int mptanhf(struct t *t) { return mpf1(t, mpfr_tanh); }
 583 int mptanhl(struct t *t) { return mpl1(t, mpfr_tanh); }
 584 // TODO: tgamma(2) raises wrong flags
 585 int mptgamma(struct t *t) { return mpd1(t, mpfr_gamma); }
 586 int mptgammaf(struct t *t) { return mpf1(t, mpfr_gamma); }
 587 int mptgammal(struct t *t) { return mpl1(t, mpfr_gamma); }
 588 int mptrunc(struct t *t) { return mpd1(t, wrap_trunc); }
 589 int mptruncf(struct t *t) { return mpf1(t, wrap_trunc); }
 590 int mptruncl(struct t *t) { return mpl1(t, wrap_trunc); }
 591 int mpj0(struct t *t) { return mpd1(t, mpfr_j0); }
 592 int mpj1(struct t *t) { return mpd1(t, mpfr_j1); }
 593 int mpy0(struct t *t) { return mpd1(t, mpfr_y0); }
 594 int mpy1(struct t *t) { return mpd1(t, mpfr_y1); }
 595 // TODO: non standard functions
 596 int mpscalb(struct t *t)
 597 {
 598         setupfenv(t->r);
 599         t->y = scalb(t->x, t->x2);
 600         t->e = getexcept();
 601         t->dy = 0; // wrong
 602         return 0;
 603 }
 604 int mpscalbf(struct t *t)
 605 {
 606         setupfenv(t->r);
 607         t->y = scalbf(t->x, t->x2);
 608         t->e = getexcept();
 609         t->dy = 0; // wrong
 610         return 0;
 611 }
 612 int mpj0f(struct t *t) { return mpf1(t, mpfr_j0); }
 613 int mpj0l(struct t *t) { return mpl1(t, mpfr_j0); }
 614 int mpj1f(struct t *t) { return mpf1(t, mpfr_j1); }
 615 int mpj1l(struct t *t) { return mpl1(t, mpfr_j1); }
 616 int mpy0f(struct t *t) { return mpf1(t, mpfr_y0); }
 617 int mpy0l(struct t *t) { return mpl1(t, mpfr_y0); }
 618 int mpy1f(struct t *t) { return mpf1(t, mpfr_y1); }
 619 int mpy1l(struct t *t) { return mpl1(t, mpfr_y1); }
 620 int mpexp10(struct t *t) { return mpd1(t, wrap_pow10); }
 621 int mpexp10f(struct t *t) { return mpf1(t, wrap_pow10); }
 622 int mpexp10l(struct t *t) { return mpl1(t, wrap_pow10); }
 623 int mppow10(struct t *t) { return mpd1(t, wrap_pow10); }
 624 int mppow10f(struct t *t) { return mpf1(t, wrap_pow10); }
 625 int mppow10l(struct t *t) { return mpl1(t, wrap_pow10); }
 626
 627 int mpfrexp(struct t *t)
 628 {
 629         mpfr_exp_t e;
 630         int k;
 631         MPFR_DECL_INIT(mx, 53);
 632
 633         t->dy = 0;
 634         t->y = 0;
 635         mpfr_clear_flags();
 636         mpfr_set_d(mx, t->x, MPFR_RNDN);
 637         k = mpfr_frexp(&e, mx, mx, t->r);
 638         t->y = mpfr_get_d(mx, MPFR_RNDN);
 639         t->i = e;
 640         t->e = eflags(isnan(t->x));
 641         return 0;
 642 }
 643
 644 int mpfrexpf(struct t *t)
 645 {
 646         mpfr_exp_t e;
 647         int k;
 648         MPFR_DECL_INIT(mx, 24);
 649
 650         t->dy = 0;
 651         t->y = 0;
 652         mpfr_clear_flags();
 653         mpfr_set_flt(mx, t->x, MPFR_RNDN);
 654         k = mpfr_frexp(&e, mx, mx, t->r);
 655         t->y = mpfr_get_flt(mx, MPFR_RNDN);
 656         t->i = e;
 657         t->e = eflags(isnan(t->x));
 658         return 0;
 659 }
 660
 661 int mpfrexpl(struct t *t)
 662 {
 663         mpfr_exp_t e;
 664         int k;
 665         MPFR_DECL_INIT(mx, 64);
 666
 667         t->dy = 0;
 668         t->y = 0;
 669         mpfr_clear_flags();
 670         mpfr_set_ld(mx, t->x, MPFR_RNDN);
 671         k = mpfr_frexp(&e, mx, mx, t->r);
 672         t->y = mpfr_get_ld(mx, MPFR_RNDN);
 673         t->i = e;
 674         t->e = eflags(isnan(t->x));
 675         return 0;
 676 }
 677
 678 int mpldexp(struct t *t)
 679 {
 680         int k;
 681         MPFR_DECL_INIT(mx, 53);
 682
 683         t->dy = 0;
 684         t->y = 0;
 685         mpfr_clear_flags();
 686         mpfr_set_d(mx, t->x, MPFR_RNDN);
 687         k = mpfr_mul_2si(mx, mx, t->i, t->r);
 688         adjust(mx, mx, k, t->r, DBL);
 689         t->y = mpfr_get_d(mx, MPFR_RNDN);
 690         t->e = eflags(isnan(t->x));
 691         return 0;
 692 }
 693
 694 int mpldexpf(struct t *t)
 695 {
 696         int k;
 697         MPFR_DECL_INIT(mx, 24);
 698
 699         t->dy = 0;
 700         t->y = 0;
 701         mpfr_clear_flags();
 702         mpfr_set_flt(mx, t->x, MPFR_RNDN);
 703         k = mpfr_mul_2si(mx, mx, t->i, t->r);
 704         adjust(mx, mx, k, t->r, FLT);
 705         t->y = mpfr_get_flt(mx, MPFR_RNDN);
 706         t->e = eflags(isnan(t->x));
 707         return 0;
 708 }
 709
 710 int mpldexpl(struct t *t)
 711 {
 712         int k;
 713         MPFR_DECL_INIT(mx, 64);
 714
 715         t->dy = 0;
 716         t->y = 0;
 717         mpfr_clear_flags();
 718         mpfr_set_ld(mx, t->x, MPFR_RNDN);
 719         k = mpfr_mul_2si(mx, mx, t->i, t->r);
 720         adjust(mx, mx, k, t->r, LDBL);
 721         t->y = mpfr_get_ld(mx, MPFR_RNDN);
 722         t->e = eflags(isnan(t->x));
 723         return 0;
 724 }
 725
 726 int mpscalbn(struct t *t) { return mpldexp(t); }
 727 int mpscalbnf(struct t *t) { return mpldexpf(t); }
 728 int mpscalbnl(struct t *t) { return mpldexpl(t); }
 729 int mpscalbln(struct t *t) { return mpldexp(t); }
 730 int mpscalblnf(struct t *t) { return mpldexpf(t); }
 731 int mpscalblnl(struct t *t) { return mpldexpl(t); }
 732
 733 int mplgamma_r(struct t *t) { return mplgamma(t); }
 734 int mplgammaf_r(struct t *t) { return mplgammaf(t); }
 735 int mplgammal_r(struct t *t) { return mplgammal(t); }
 736
 737 int mpilogb(struct t *t)
 738 {
 739         MPFR_DECL_INIT(mx, 53);
 740
 741         mpfr_set_d(mx, t->x, MPFR_RNDN);
 742         t->i = mpfr_get_exp(mx) - 1;
 743         t->e = 0;
 744         if (isinf(t->x) || isnan(t->x) || t->x == 0)
 745                 t->e = INVALID;
 746         return 0;
 747 }
 748 int mpilogbf(struct t *t)
 749 {
 750         MPFR_DECL_INIT(mx, 24);
 751
 752         mpfr_set_flt(mx, t->x, MPFR_RNDN);
 753         t->i = mpfr_get_exp(mx) - 1;
 754         t->e = 0;
 755         if (isinf(t->x) || isnan(t->x) || t->x == 0)
 756                 t->e = INVALID;
 757         return 0;
 758 }
 759 int mpilogbl(struct t *t)
 760 {
 761         MPFR_DECL_INIT(mx, 64);
 762
 763         mpfr_set_ld(mx, t->x, MPFR_RNDN);
 764         t->i = mpfr_get_exp(mx) - 1;
 765         t->e = 0;
 766         if (isinf(t->x) || isnan(t->x) || t->x == 0)
 767                 t->e = INVALID;
 768         return 0;
 769 }
 770
 771 // TODO: ll* is hard to do with mpfr
 772 #define mp_f_i(n) \
 773 int mp##n(struct t *t) \
 774 { \
 775         setupfenv(t->r); \
 776         t->i = n(t->x); \
 777         t->e = getexcept(); \
 778         if (t->e & INVALID) \
 779                 t->i = 0; \
 780         return 0; \
 781 }
 782
 783 mp_f_i(llrint)
 784 mp_f_i(llrintf)
 785 mp_f_i(llrintl)
 786 mp_f_i(lrint)
 787 mp_f_i(lrintf)
 788 mp_f_i(lrintl)
 789 mp_f_i(llround)
 790 mp_f_i(llroundf)
 791 mp_f_i(llroundl)
 792 mp_f_i(lround)
 793 mp_f_i(lroundf)
 794 mp_f_i(lroundl)
 795
 796 int mpmodf(struct t *t)
 797 {
 798         int e, r;
 799
 800         r = mpd1(t, wrap_trunc);
 801         if (r)
 802                 return r;
 803         t->y2 = t->y;
 804         t->dy2 = t->dy;
 805         e = t->e & ~INEXACT;
 806         r = mpd1(t, mpfr_frac);
 807         t->e |= e;
 808         return r;
 809 }
 810
 811 int mpmodff(struct t *t)
 812 {
 813         int e, r;
 814
 815         r = mpf1(t, wrap_trunc);
 816         if (r)
 817                 return r;
 818         t->y2 = t->y;
 819         t->dy2 = t->dy;
 820         e = t->e & ~INEXACT;
 821         r = mpf1(t, mpfr_frac);
 822         t->e |= e;
 823         return r;
 824 }
 825
 826 int mpmodfl(struct t *t)
 827 {
 828         int e, r;
 829
 830         r = mpl1(t, wrap_trunc);
 831         if (r)
 832                 return r;
 833         t->y2 = t->y;
 834         t->dy2 = t->dy;
 835         e = t->e & ~INEXACT;
 836         r = mpl1(t, mpfr_frac);
 837         t->e |= e;
 838         return r;
 839 }
 840
 841 int mpsincos(struct t *t)
 842 {
 843         int e, r;
 844
 845         r = mpd1(t, mpfr_cos);
 846         if (r)
 847                 return r;
 848         t->y2 = t->y;
 849         t->dy2 = t->dy;
 850         e = t->e;
 851         r = mpd1(t, mpfr_sin);
 852         t->e |= e;
 853         return r;
 854 }
 855
 856 int mpsincosf(struct t *t)
 857 {
 858         int e, r;
 859
 860         r = mpf1(t, mpfr_cos);
 861         if (r)
 862                 return r;
 863         t->y2 = t->y;
 864         t->dy2 = t->dy;
 865         e = t->e;
 866         r = mpf1(t, mpfr_sin);
 867         t->e |= e;
 868         return r;
 869 }
 870
 871 int mpsincosl(struct t *t)
 872 {
 873         int e, r;
 874
 875         r = mpl1(t, mpfr_cos);
 876         if (r)
 877                 return r;
 878         t->y2 = t->y;
 879         t->dy2 = t->dy;
 880         e = t->e;
 881         r = mpl1(t, mpfr_sin);
 882         t->e |= e;
 883         return r;
 884 }
 885
 886 int mpremquo(struct t *t) { return mpd2(t, wrap_remquo) || (t->i = mpremquo_q, 0); }
 887 int mpremquof(struct t *t) { return mpf2(t, wrap_remquo) || (t->i = mpremquo_q, 0); }
 888 int mpremquol(struct t *t) { return mpl2(t, wrap_remquo) || (t->i = mpremquo_q, 0); }
 889
 890 int mpfma(struct t *t)
 891 {
 892         int tn;
 893         int r = rmap(t->r);
 894         MPFR_DECL_INIT(mx, 53);
 895         MPFR_DECL_INIT(mx2, 53);
 896         MPFR_DECL_INIT(mx3, 53);
 897         MPFR_DECL_INIT(my, 128);
 898
 899         mpfr_clear_flags();
 900         mpfr_set_d(mx, t->x, MPFR_RNDN);
 901         mpfr_set_d(mx2, t->x2, MPFR_RNDN);
 902         mpfr_set_d(mx3, t->x3, MPFR_RNDN);
 903         tn = mpfr_fma(my, mx, mx2, mx3, r);
 904         gend(t, my, tn, r);
 905         return 0;
 906 }
 907
 908 int mpfmaf(struct t *t)
 909 {
 910         int tn;
 911         int r = rmap(t->r);
 912         MPFR_DECL_INIT(mx, 24);
 913         MPFR_DECL_INIT(mx2, 24);
 914         MPFR_DECL_INIT(mx3, 24);
 915         MPFR_DECL_INIT(my, 128);
 916
 917         mpfr_clear_flags();
 918         mpfr_set_flt(mx, t->x, MPFR_RNDN);
 919         mpfr_set_flt(mx2, t->x2, MPFR_RNDN);
 920         mpfr_set_flt(mx3, t->x3, MPFR_RNDN);
 921         tn = mpfr_fma(my, mx, mx2, mx3, r);
 922         genf(t, my, tn, r);
 923         return 0;
 924 }
 925
 926 int mpfmal(struct t *t)
 927 {
 928 #if LDBL_MANT_DIG == 53
 929         return mpfma(t);
 930 #elif LDBL_MANT_DIG == 64
 931         int tn;
 932         int r = rmap(t->r);
 933         MPFR_DECL_INIT(mx, 64);
 934         MPFR_DECL_INIT(mx2, 64);
 935         MPFR_DECL_INIT(mx3, 64);
 936         MPFR_DECL_INIT(my, 128);
 937
 938         mpfr_clear_flags();
 939         mpfr_set_ld(mx, t->x, MPFR_RNDN);
 940         mpfr_set_ld(mx2, t->x2, MPFR_RNDN);
 941         mpfr_set_ld(mx3, t->x3, MPFR_RNDN);
 942         tn = mpfr_fma(my, mx, mx2, mx3, r);
 943         genl(t, my, tn, r);
 944         return 0;
 945 #else
 946         return -1;
 947 #endif
 948 }
 949
 950 int mpjn(struct t *t) { mpbessel_n = t->i; return mpd1(t, wrap_jn); }
 951 int mpjnf(struct t *t) { mpbessel_n = t->i; return mpf1(t, wrap_jn); }
 952 int mpjnl(struct t *t) { mpbessel_n = t->i; return mpl1(t, wrap_jn); }
 953 int mpyn(struct t *t) { mpbessel_n = t->i; return mpd1(t, wrap_yn); }
 954 int mpynf(struct t *t) { mpbessel_n = t->i; return mpf1(t, wrap_yn); }
 955 int mpynl(struct t *t) { mpbessel_n = t->i; return mpl1(t, wrap_yn); }
 956