1#include <stdio.h> 2#include <stdint.h> 3#include <mpfr.h> 4#include "gen.h" 5 6static 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 17enum {FLT, DBL, LDBL}; 18static const int emin[] = { 19[FLT] = -148, 20[DBL] = -1073, 21[LDBL] = -16444 22}; 23static const int emax[] = { 24[FLT] = 128, 25[DBL] = 1024, 26[LDBL] = 16384 27}; 28 29void debug(mpfr_t x) 30{ 31 mpfr_out_str(stdout, 10, 0, x, MPFR_RNDN); 32 printf("\n"); 33} 34 35/* 36round x into y considering x is already rounded (t = up or down) 37 38only cases where adjustment is done: 39 x=...|1...0, t=up -> x=nextbelow(x) 40 x=...|1...0, t=down -> x=nextabove(x) 41where | 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 45static 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 71static 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) 100static 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 118static 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 141static 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 156static 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 172static 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 195static 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 210static 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 227static 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 251static 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 272static 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 295static int mplgamma_sign; 296static 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} 300static long mpremquo_q; 301static 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} 305static int mpbessel_n; 306static 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} 310static 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} 314static int wrap_ceil(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r) 315{ 316 return mpfr_ceil(my, mx); 317} 318static int wrap_floor(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r) 319{ 320 return mpfr_floor(my, mx); 321} 322static int wrap_round(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r) 323{ 324 return mpfr_round(my, mx); 325} 326static int wrap_trunc(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r) 327{ 328 return mpfr_trunc(my, mx); 329} 330static 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} 336static 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 342static 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} 350int mpsinpi(struct t *t) { return mpd1(t, wrap_sinpi); } 351 352int mpadd(struct t *t) { return mpd2(t, mpfr_add); } 353int mpaddf(struct t *t) { return mpf2(t, mpfr_add); } 354int mpaddl(struct t *t) { return mpl2(t, mpfr_add); } 355int mpmul(struct t *t) { return mpd2(t, mpfr_mul); } 356int mpmulf(struct t *t) { return mpf2(t, mpfr_mul); } 357int mpmull(struct t *t) { return mpl2(t, mpfr_mul); } 358int mpdiv(struct t *t) { return mpd2(t, mpfr_div); } 359int mpdivf(struct t *t) { return mpf2(t, mpfr_div); } 360int mpdivl(struct t *t) { return mpl2(t, mpfr_div); } 361 362int mpacos(struct t *t) { return mpd1(t, mpfr_acos); } 363int mpacosf(struct t *t) { return mpf1(t, mpfr_acos); } 364int mpacosl(struct t *t) { return mpl1(t, mpfr_acos); } 365int mpacosh(struct t *t) { return mpd1(t, mpfr_acosh); } 366int mpacoshf(struct t *t) { return mpf1(t, mpfr_acosh); } 367int mpacoshl(struct t *t) { return mpl1(t, mpfr_acosh); } 368int mpasin(struct t *t) { return mpd1(t, mpfr_asin); } 369int mpasinf(struct t *t) { return mpf1(t, mpfr_asin); } 370int mpasinl(struct t *t) { return mpl1(t, mpfr_asin); } 371int mpasinh(struct t *t) { return mpd1(t, mpfr_asinh); } 372int mpasinhf(struct t *t) { return mpf1(t, mpfr_asinh); } 373int mpasinhl(struct t *t) { return mpl1(t, mpfr_asinh); } 374int mpatan(struct t *t) { return mpd1(t, mpfr_atan); } 375int mpatanf(struct t *t) { return mpf1(t, mpfr_atan); } 376int mpatanl(struct t *t) { return mpl1(t, mpfr_atan); } 377int mpatan2(struct t *t) { return mpd2(t, mpfr_atan2); } 378int mpatan2f(struct t *t) { return mpf2(t, mpfr_atan2); } 379int mpatan2l(struct t *t) { return mpl2(t, mpfr_atan2); } 380int mpatanh(struct t *t) { return mpd1(t, mpfr_atanh); } 381int mpatanhf(struct t *t) { return mpf1(t, mpfr_atanh); } 382int mpatanhl(struct t *t) { return mpl1(t, mpfr_atanh); } 383int mpcbrt(struct t *t) { return mpd1(t, mpfr_cbrt); } 384int mpcbrtf(struct t *t) { return mpf1(t, mpfr_cbrt); } 385int mpcbrtl(struct t *t) { return mpl1(t, mpfr_cbrt); } 386int mpceil(struct t *t) { return mpd1(t, wrap_ceil); } 387int mpceilf(struct t *t) { return mpf1(t, wrap_ceil); } 388int mpceill(struct t *t) { return mpl1(t, wrap_ceil); } 389int mpcopysign(struct t *t) { return mpd2(t, mpfr_copysign); } 390int mpcopysignf(struct t *t) { return mpf2(t, mpfr_copysign); } 391int mpcopysignl(struct t *t) { return mpl2(t, mpfr_copysign); } 392int mpcos(struct t *t) { return mpd1(t, mpfr_cos); } 393int mpcosf(struct t *t) { return mpf1(t, mpfr_cos); } 394int mpcosl(struct t *t) { return mpl1(t, mpfr_cos); } 395int mpcosh(struct t *t) { return mpd1(t, mpfr_cosh); } 396int mpcoshf(struct t *t) { return mpf1(t, mpfr_cosh); } 397int mpcoshl(struct t *t) { return mpl1(t, mpfr_cosh); } 398int mperf(struct t *t) { return mpd1(t, mpfr_erf); } 399int mperff(struct t *t) { return mpf1(t, mpfr_erf); } 400int mperfl(struct t *t) { return mpl1(t, mpfr_erf); } 401int mperfc(struct t *t) { return mpd1(t, mpfr_erfc); } 402int mperfcf(struct t *t) { return mpf1(t, mpfr_erfc); } 403int mperfcl(struct t *t) { return mpl1(t, mpfr_erfc); } 404int mpexp(struct t *t) { return mpd1(t, mpfr_exp); } 405int mpexpf(struct t *t) { return mpf1(t, mpfr_exp); } 406int mpexpl(struct t *t) { return mpl1(t, mpfr_exp); } 407int mpexp2(struct t *t) { return mpd1(t, mpfr_exp2); } 408int mpexp2f(struct t *t) { return mpf1(t, mpfr_exp2); } 409int mpexp2l(struct t *t) { return mpl1(t, mpfr_exp2); } 410int mpexpm1(struct t *t) { return mpd1(t, mpfr_expm1); } 411int mpexpm1f(struct t *t) { return mpf1(t, mpfr_expm1); } 412int mpexpm1l(struct t *t) { return mpl1(t, mpfr_expm1); } 413int mpfabs(struct t *t) { return mpd1(t, mpfr_abs); } 414int mpfabsf(struct t *t) { return mpf1(t, mpfr_abs); } 415int mpfabsl(struct t *t) { return mpl1(t, mpfr_abs); } 416int mpfdim(struct t *t) { return mpd2(t, mpfr_dim); } 417int mpfdimf(struct t *t) { return mpf2(t, mpfr_dim); } 418int mpfdiml(struct t *t) { return mpl2(t, mpfr_dim); } 419int mpfloor(struct t *t) { return mpd1(t, wrap_floor); } 420int mpfloorf(struct t *t) { return mpf1(t, wrap_floor); } 421int mpfloorl(struct t *t) { return mpl1(t, wrap_floor); } 422int mpfmax(struct t *t) { return mpd2(t, mpfr_max); } 423int mpfmaxf(struct t *t) { return mpf2(t, mpfr_max); } 424int mpfmaxl(struct t *t) { return mpl2(t, mpfr_max); } 425int mpfmin(struct t *t) { return mpd2(t, mpfr_min); } 426int mpfminf(struct t *t) { return mpf2(t, mpfr_min); } 427int mpfminl(struct t *t) { return mpl2(t, mpfr_min); } 428int mpfmod(struct t *t) { return mpd2(t, mpfr_fmod); } 429int mpfmodf(struct t *t) { return mpf2(t, mpfr_fmod); } 430int mpfmodl(struct t *t) { return mpl2(t, mpfr_fmod); } 431int mphypot(struct t *t) { return mpd2(t, mpfr_hypot); } 432int mphypotf(struct t *t) { return mpf2(t, mpfr_hypot); } 433int mphypotl(struct t *t) { return mpl2(t, mpfr_hypot); } 434int mplgamma(struct t *t) { return mpd1(t, wrap_lgamma) || (t->i = mplgamma_sign, 0); } 435int mplgammaf(struct t *t) { return mpf1(t, wrap_lgamma) || (t->i = mplgamma_sign, 0); } 436int mplgammal(struct t *t) { return mpl1(t, wrap_lgamma) || (t->i = mplgamma_sign, 0); } 437int mplog(struct t *t) { return mpd1(t, mpfr_log); } 438int mplogf(struct t *t) { return mpf1(t, mpfr_log); } 439int mplogl(struct t *t) { return mpl1(t, mpfr_log); } 440int mplog10(struct t *t) { return mpd1(t, mpfr_log10); } 441int mplog10f(struct t *t) { return mpf1(t, mpfr_log10); } 442int mplog10l(struct t *t) { return mpl1(t, mpfr_log10); } 443int mplog1p(struct t *t) { return mpd1(t, mpfr_log1p); } 444int mplog1pf(struct t *t) { return mpf1(t, mpfr_log1p); } 445int mplog1pl(struct t *t) { return mpl1(t, mpfr_log1p); } 446int mplog2(struct t *t) { return mpd1(t, mpfr_log2); } 447int mplog2f(struct t *t) { return mpf1(t, mpfr_log2); } 448int mplog2l(struct t *t) { return mpl1(t, mpfr_log2); } 449int 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} 472int 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} 495int 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} 518int mpnearbyint(struct t *t) { return mpd1(t, wrap_nearbyint) || (t->e&=~INEXACT, 0); } 519int mpnearbyintf(struct t *t) { return mpf1(t, wrap_nearbyint) || (t->e&=~INEXACT, 0); } 520int mpnearbyintl(struct t *t) { return mpl1(t, wrap_nearbyint) || (t->e&=~INEXACT, 0); } 521// TODO: hard to implement with mpfr 522int 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} 530int 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} 538int 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} 546int 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} 554int 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} 562int mpnexttowardl(struct t *t) { return mpnextafterl(t); } 563int mppow(struct t *t) { return mpd2(t, mpfr_pow); } 564int mppowf(struct t *t) { return mpf2(t, mpfr_pow); } 565int mppowl(struct t *t) { return mpl2(t, mpfr_pow); } 566int mpremainder(struct t *t) { return mpd2(t, mpfr_remainder); } 567int mpremainderf(struct t *t) { return mpf2(t, mpfr_remainder); } 568int mpremainderl(struct t *t) { return mpl2(t, mpfr_remainder); } 569int mprint(struct t *t) { return mpd1(t, mpfr_rint); } 570int mprintf(struct t *t) { return mpf1(t, mpfr_rint); } 571int mprintl(struct t *t) { return mpl1(t, mpfr_rint); } 572int mpround(struct t *t) { return mpd1(t, wrap_round); } 573int mproundf(struct t *t) { return mpf1(t, wrap_round); } 574int mproundl(struct t *t) { return mpl1(t, wrap_round); } 575int mpsin(struct t *t) { return mpd1(t, mpfr_sin); } 576int mpsinf(struct t *t) { return mpf1(t, mpfr_sin); } 577int mpsinl(struct t *t) { return mpl1(t, mpfr_sin); } 578int mpsinh(struct t *t) { return mpd1(t, mpfr_sinh); } 579int mpsinhf(struct t *t) { return mpf1(t, mpfr_sinh); } 580int mpsinhl(struct t *t) { return mpl1(t, mpfr_sinh); } 581int mpsqrt(struct t *t) { return mpd1(t, mpfr_sqrt); } 582int mpsqrtf(struct t *t) { return mpf1(t, mpfr_sqrt); } 583int mpsqrtl(struct t *t) { return mpl1(t, mpfr_sqrt); } 584int mptan(struct t *t) { return mpd1(t, mpfr_tan); } 585int mptanf(struct t *t) { return mpf1(t, mpfr_tan); } 586int mptanl(struct t *t) { return mpl1(t, mpfr_tan); } 587int mptanh(struct t *t) { return mpd1(t, mpfr_tanh); } 588int mptanhf(struct t *t) { return mpf1(t, mpfr_tanh); } 589int mptanhl(struct t *t) { return mpl1(t, mpfr_tanh); } 590// TODO: tgamma(2) raises wrong flags 591int mptgamma(struct t *t) { return mpd1(t, mpfr_gamma); } 592int mptgammaf(struct t *t) { return mpf1(t, mpfr_gamma); } 593int mptgammal(struct t *t) { return mpl1(t, mpfr_gamma); } 594int mptrunc(struct t *t) { return mpd1(t, wrap_trunc); } 595int mptruncf(struct t *t) { return mpf1(t, wrap_trunc); } 596int mptruncl(struct t *t) { return mpl1(t, wrap_trunc); } 597int mpj0(struct t *t) { return mpd1(t, mpfr_j0); } 598int mpj1(struct t *t) { return mpd1(t, mpfr_j1); } 599int mpy0(struct t *t) { return mpd1(t, mpfr_y0); } 600int mpy1(struct t *t) { return mpd1(t, mpfr_y1); } 601// TODO: non standard functions 602int 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} 610int 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} 618int mpj0f(struct t *t) { return mpf1(t, mpfr_j0); } 619int mpj0l(struct t *t) { return mpl1(t, mpfr_j0); } 620int mpj1f(struct t *t) { return mpf1(t, mpfr_j1); } 621int mpj1l(struct t *t) { return mpl1(t, mpfr_j1); } 622int mpy0f(struct t *t) { return mpf1(t, mpfr_y0); } 623int mpy0l(struct t *t) { return mpl1(t, mpfr_y0); } 624int mpy1f(struct t *t) { return mpf1(t, mpfr_y1); } 625int mpy1l(struct t *t) { return mpl1(t, mpfr_y1); } 626int mpexp10(struct t *t) { return mpd1(t, wrap_pow10); } 627int mpexp10f(struct t *t) { return mpf1(t, wrap_pow10); } 628int mpexp10l(struct t *t) { return mpl1(t, wrap_pow10); } 629int mppow10(struct t *t) { return mpd1(t, wrap_pow10); } 630int mppow10f(struct t *t) { return mpf1(t, wrap_pow10); } 631int mppow10l(struct t *t) { return mpl1(t, wrap_pow10); } 632 633int 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 650int 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 667int 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 684int 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 700int 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 716int 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 732int mpscalbn(struct t *t) { return mpldexp(t); } 733int mpscalbnf(struct t *t) { return mpldexpf(t); } 734int mpscalbnl(struct t *t) { return mpldexpl(t); } 735int mpscalbln(struct t *t) { return mpldexp(t); } 736int mpscalblnf(struct t *t) { return mpldexpf(t); } 737int mpscalblnl(struct t *t) { return mpldexpl(t); } 738 739int mplgamma_r(struct t *t) { return mplgamma(t); } 740int mplgammaf_r(struct t *t) { return mplgammaf(t); } 741int mplgammal_r(struct t *t) { return mplgammal(t); } 742 743int 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} 754int 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} 765int 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) \ 779int 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 789mp_f_i(llrint) 790mp_f_i(llrintf) 791mp_f_i(llrintl) 792mp_f_i(lrint) 793mp_f_i(lrintf) 794mp_f_i(lrintl) 795mp_f_i(llround) 796mp_f_i(llroundf) 797mp_f_i(llroundl) 798mp_f_i(lround) 799mp_f_i(lroundf) 800mp_f_i(lroundl) 801 802int 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 817int 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 832int 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 847int 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 862int 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 877int 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 892int mpremquo(struct t *t) { return mpd2(t, wrap_remquo) || (t->i = mpremquo_q, 0); } 893int mpremquof(struct t *t) { return mpf2(t, wrap_remquo) || (t->i = mpremquo_q, 0); } 894int mpremquol(struct t *t) { return mpl2(t, wrap_remquo) || (t->i = mpremquo_q, 0); } 895 896int 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 914int 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 932int 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 956int mpjn(struct t *t) { mpbessel_n = t->i; return mpd1(t, wrap_jn); } 957int mpjnf(struct t *t) { mpbessel_n = t->i; return mpf1(t, wrap_jn); } 958int mpjnl(struct t *t) { mpbessel_n = t->i; return mpl1(t, wrap_jn); } 959int mpyn(struct t *t) { mpbessel_n = t->i; return mpd1(t, wrap_yn); } 960int mpynf(struct t *t) { mpbessel_n = t->i; return mpf1(t, wrap_yn); } 961int mpynl(struct t *t) { mpbessel_n = t->i; return mpl1(t, wrap_yn); } 962 963