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
); }
