9 case RN: return MPFR_RNDN;
10 case RZ: return MPFR_RNDZ;
11 case RD: return MPFR_RNDD;
12 case RU: return MPFR_RNDU;
17 enum {FLT, DBL, LDBL};
18 static const int emin[] = {
23 static const int emax[] = {
31 mpfr_out_str(stdout, 10, 0, x, MPFR_RNDN);
36 round x into y considering x is already rounded (t = up or down)
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
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)
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);
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);
68 return mpfr_set(y, x, r);
71 static int adjust(mpfr_t mr, mpfr_t my, int t, int r, int type)
74 //printf("adj %d\n", t);
76 t = adjust_round(mr, my, t, r);
77 //printf("rnd %d\n", t);
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);
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);
99 //static int eflags(mpfr_t mr, mpfr_t my, int t)
100 static int eflags(int naninput)
104 if (mpfr_inexflag_p())
106 // if (mpfr_underflow_p() && (t || mpfr_cmp(mr, my) != 0))
107 if (mpfr_underflow_p() && i)
109 if (mpfr_overflow_p())
113 if (!naninput && (mpfr_nanflag_p() || mpfr_erangeflag_p()))
118 static void genf(struct t *p, mpfr_t my, int t, int r)
120 MPFR_DECL_INIT(mr, 24);
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));
127 if (!isfinite(p->y)) {
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
141 static int mpf1(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, mpfr_rnd_t))
145 MPFR_DECL_INIT(mx, 24);
146 MPFR_DECL_INIT(my, 128);
149 mpfr_set_flt(mx, p->x, MPFR_RNDN);
156 static int mpf2(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t))
160 MPFR_DECL_INIT(mx, 24);
161 MPFR_DECL_INIT(mx2, 24);
162 MPFR_DECL_INIT(my, 128);
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);
172 static void gend(struct t *p, mpfr_t my, int t, int r)
174 MPFR_DECL_INIT(mr, 53);
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));
181 if (!isfinite(p->y)) {
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
195 static int mpd1(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, mpfr_rnd_t))
199 MPFR_DECL_INIT(mx, 53);
200 MPFR_DECL_INIT(my, 128);
203 mpfr_set_d(mx, p->x, MPFR_RNDN);
210 static int mpd2(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t))
214 MPFR_DECL_INIT(mx, 53);
215 MPFR_DECL_INIT(mx2, 53);
216 MPFR_DECL_INIT(my, 128);
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);
226 #if LDBL_MANT_DIG == 64
227 static void genl(struct t *p, mpfr_t my, int t, int r)
229 MPFR_DECL_INIT(mr, 64);
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));
236 if (!isfinite(p->y)) {
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
251 static int mpl1(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, mpfr_rnd_t))
253 #if LDBL_MANT_DIG == 53
255 #elif LDBL_MANT_DIG == 64
258 MPFR_DECL_INIT(mx, 64);
259 MPFR_DECL_INIT(my, 128);
262 mpfr_set_ld(mx, p->x, MPFR_RNDN);
272 static int mpl2(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, const mpfr_t, mpfr_rnd_t))
274 #if LDBL_MANT_DIG == 53
276 #elif LDBL_MANT_DIG == 64
279 MPFR_DECL_INIT(mx, 64);
280 MPFR_DECL_INIT(mx2, 64);
281 MPFR_DECL_INIT(my, 128);
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);
295 static int mplgamma_sign;
296 static int wrap_lgamma(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
298 return mpfr_lgamma(my, &mplgamma_sign, mx, r);
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)
303 return mpfr_remquo(my, &mpremquo_q, mx, mx2, r);
305 static int mpbessel_n;
306 static int wrap_jn(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
308 return mpfr_jn(my, mpbessel_n, mx, r);
310 static int wrap_yn(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
312 return mpfr_yn(my, mpbessel_n, mx, r);
314 static int wrap_ceil(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
316 return mpfr_ceil(my, mx);
318 static int wrap_floor(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
320 return mpfr_floor(my, mx);
322 static int wrap_round(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
324 return mpfr_round(my, mx);
326 static int wrap_trunc(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
328 return mpfr_trunc(my, mx);
330 static int wrap_nearbyint(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
332 int i = mpfr_rint(my, mx, r);
333 mpfr_clear_inexflag();
336 static int wrap_pow10(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
338 return mpfr_ui_pow(my, 10, mx, r);
342 static int wrap_sinpi(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
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);
350 int mpsinpi(struct t *t) { return mpd1(t, wrap_sinpi); }
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); }
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)
451 MPFR_DECL_INIT(mx, 53);
468 mpfr_set_d(mx, t->x, MPFR_RNDN);
469 t->y = mpfr_get_exp(mx) - 1;
472 int mplogbf(struct t *t)
474 MPFR_DECL_INIT(mx, 24);
491 mpfr_set_flt(mx, t->x, MPFR_RNDN);
492 t->y = mpfr_get_exp(mx) - 1;
495 int mplogbl(struct t *t)
497 MPFR_DECL_INIT(mx, 64);
514 mpfr_set_ld(mx, t->x, MPFR_RNDN);
515 t->y = mpfr_get_exp(mx) - 1;
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)
524 feclearexcept(FE_ALL_EXCEPT);
525 t->y = nextafter(t->x, t->x2);
530 int mpnextafterf(struct t *t)
532 feclearexcept(FE_ALL_EXCEPT);
533 t->y = nextafterf(t->x, t->x2);
538 int mpnextafterl(struct t *t)
540 feclearexcept(FE_ALL_EXCEPT);
541 t->y = nextafterl(t->x, t->x2);
546 int mpnexttoward(struct t *t)
548 feclearexcept(FE_ALL_EXCEPT);
549 t->y = nexttoward(t->x, t->x2);
554 int mpnexttowardf(struct t *t)
556 feclearexcept(FE_ALL_EXCEPT);
557 t->y = nexttowardf(t->x, t->x2);
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)
605 t->y = scalb(t->x, t->x2);
610 int mpscalbf(struct t *t)
613 t->y = scalbf(t->x, t->x2);
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); }
633 int mpfrexp(struct t *t)
637 MPFR_DECL_INIT(mx, 53);
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);
646 t->e = eflags(isnan(t->x));
650 int mpfrexpf(struct t *t)
654 MPFR_DECL_INIT(mx, 24);
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);
663 t->e = eflags(isnan(t->x));
667 int mpfrexpl(struct t *t)
671 MPFR_DECL_INIT(mx, 64);
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);
680 t->e = eflags(isnan(t->x));
684 int mpldexp(struct t *t)
687 MPFR_DECL_INIT(mx, 53);
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));
700 int mpldexpf(struct t *t)
703 MPFR_DECL_INIT(mx, 24);
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));
716 int mpldexpl(struct t *t)
719 MPFR_DECL_INIT(mx, 64);
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));
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); }
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); }
743 int mpilogb(struct t *t)
745 MPFR_DECL_INIT(mx, 53);
747 mpfr_set_d(mx, t->x, MPFR_RNDN);
748 t->i = mpfr_get_exp(mx) - 1;
750 if (isinf(t->x) || isnan(t->x) || t->x == 0)
754 int mpilogbf(struct t *t)
756 MPFR_DECL_INIT(mx, 24);
758 mpfr_set_flt(mx, t->x, MPFR_RNDN);
759 t->i = mpfr_get_exp(mx) - 1;
761 if (isinf(t->x) || isnan(t->x) || t->x == 0)
765 int mpilogbl(struct t *t)
767 MPFR_DECL_INIT(mx, 64);
769 mpfr_set_ld(mx, t->x, MPFR_RNDN);
770 t->i = mpfr_get_exp(mx) - 1;
772 if (isinf(t->x) || isnan(t->x) || t->x == 0)
777 // TODO: ll* is hard to do with mpfr
779 int mp##n(struct t *t) \
783 t->e = getexcept(); \
784 if (t->e & INVALID) \
802 int mpmodf(struct t *t)
806 r = mpd1(t, wrap_trunc);
812 r = mpd1(t, mpfr_frac);
817 int mpmodff(struct t *t)
821 r = mpf1(t, wrap_trunc);
827 r = mpf1(t, mpfr_frac);
832 int mpmodfl(struct t *t)
836 r = mpl1(t, wrap_trunc);
842 r = mpl1(t, mpfr_frac);
847 int mpsincos(struct t *t)
851 r = mpd1(t, mpfr_cos);
857 r = mpd1(t, mpfr_sin);
862 int mpsincosf(struct t *t)
866 r = mpf1(t, mpfr_cos);
872 r = mpf1(t, mpfr_sin);
877 int mpsincosl(struct t *t)
881 r = mpl1(t, mpfr_cos);
887 r = mpl1(t, mpfr_sin);
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); }
896 int mpfma(struct t *t)
900 MPFR_DECL_INIT(mx, 53);
901 MPFR_DECL_INIT(mx2, 53);
902 MPFR_DECL_INIT(mx3, 53);
903 MPFR_DECL_INIT(my, 128);
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);
914 int mpfmaf(struct t *t)
918 MPFR_DECL_INIT(mx, 24);
919 MPFR_DECL_INIT(mx2, 24);
920 MPFR_DECL_INIT(mx3, 24);
921 MPFR_DECL_INIT(my, 128);
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);
932 int mpfmal(struct t *t)
934 #if LDBL_MANT_DIG == 53
936 #elif LDBL_MANT_DIG == 64
939 MPFR_DECL_INIT(mx, 64);
940 MPFR_DECL_INIT(mx2, 64);
941 MPFR_DECL_INIT(mx3, 64);
942 MPFR_DECL_INIT(my, 128);
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);
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); }