X-Git-Url: http://nsz.repo.hu/git/?a=blobdiff_plain;f=src%2Fmath%2Fgen%2Fmp.c;h=be2ed39d6bc7ca9658e032d2a31d3e8557634354;hb=ae0f0fe09b7fc9d44d072c3fd08372991d852b1d;hp=96cd33a7a18e77ea8d169c9b8e909d6d2ce2016b;hpb=f9d17902a35b6403b7c8354845e9f13f882c1c8e;p=libc-test

diff --git a/src/math/gen/mp.c b/src/math/gen/mp.c
index 96cd33a..be2ed39 100644
--- a/src/math/gen/mp.c
+++ b/src/math/gen/mp.c
@@ -123,7 +123,7 @@ static void genf(struct t *p, mpfr_t my, int t, int r)
 
 	t = adjust(mr, my, t, r);
 	p->y = mpfr_get_flt(mr, r);
-	p->e = eflags(isnan(p->x) || isnan(p->x2));
+	p->e = eflags(isnan(p->x) || isnan(p->x2) || isnan(p->x3));
 	i = eulpf(p->y);
 	if (!isfinite(p->y)) {
 		p->dy = 0;
@@ -193,7 +193,7 @@ static void gend(struct t *p, mpfr_t my, int t, int r)
 
 	t = adjust(mr, my, t, r);
 	p->y = mpfr_get_d(mr, r);
-	p->e = eflags(isnan(p->x) || isnan(p->x2));
+	p->e = eflags(isnan(p->x) || isnan(p->x2) || isnan(p->x3));
 	i = eulp(p->y);
 	if (!isfinite(p->y)) {
 		p->dy = 0;
@@ -268,7 +268,7 @@ static void genl(struct t *p, mpfr_t my, int t, int r)
 
 	t = adjust(mr, my, t, r);
 	p->y = mpfr_get_ld(mr, r);
-	p->e = eflags(isnan(p->x) || isnan(p->x2));
+	p->e = eflags(isnan(p->x) || isnan(p->x2) || isnan(p->x3));
 	i = eulpl(p->y);
 	if (!isfinite(p->y)) {
 		p->dy = 0;
@@ -327,8 +327,8 @@ static int mpl2(struct t *p, int (*fmp)(mpfr_t, const mpfr_t, const mpfr_t, mpfr
 
 	mpsetupl();
 	mpfr_clear_flags();
-	mpfr_set_d(mx, p->x, MPFR_RNDN);
-	mpfr_set_d(mx2, p->x2, MPFR_RNDN);
+	mpfr_set_ld(mx, p->x, MPFR_RNDN);
+	mpfr_set_ld(mx2, p->x2, MPFR_RNDN);
 	tn = fmp(my, mx, mx2, r);
 	genl(p, my, tn, r);
 	if ((p->e & INEXACT) && nextafterl(p->y, 0) == 0) {
@@ -350,6 +350,20 @@ static int wrap_lgamma(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 {
 	return mpfr_lgamma(my, &mplgamma_sign, mx, r);
 }
+static long mpremquo_q;
+static int wrap_remquo(mpfr_t my, const mpfr_t mx, const mpfr_t mx2, mpfr_rnd_t r)
+{
+	return mpfr_remquo(my, &mpremquo_q, mx, mx2, r);
+}
+static int mpbessel_n;
+static int wrap_jn(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
+{
+	return mpfr_jn(my, mpbessel_n, mx, r);
+}
+static int wrap_yn(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
+{
+	return mpfr_yn(my, mpbessel_n, mx, r);
+}
 static int wrap_ceil(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 {
 	return mpfr_ceil(my, mx);
@@ -377,6 +391,8 @@ static int wrap_pow10(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 	return mpfr_ui_pow(my, 10, mx, r);
 }
 
+int mpadd(struct t *t) { return mpd2(t, mpfr_add); }
+
 int mpacos(struct t *t) { return mpd1(t, mpfr_acos); }
 int mpacosf(struct t *t) { return mpf1(t, mpfr_acos); }
 int mpacosl(struct t *t) { return mpl1(t, mpfr_acos); }
@@ -468,39 +484,116 @@ int mplogb(struct t *t)
 {
 	MPFR_DECL_INIT(mx, 53);
 
-	mpfr_set_d(mx, t->x, MPFR_RNDN);
-	t->y = mpfr_get_exp(mx) + 1;
 	t->dy = 0;
 	t->e = 0;
+	if (t->x == 0) {
+		t->y = -INFINITY;
+		t->e |= DIVBYZERO;
+		return 0;
+	}
+	if (isinf(t->x)) {
+		t->y = INFINITY;
+		return 0;
+	}
+	if (isnan(t->x)) {
+		t->y = t->x;
+		return 0;
+	}
+	mpfr_set_d(mx, t->x, MPFR_RNDN);
+	t->y = mpfr_get_exp(mx) - 1;
 	return 0;
 }
 int mplogbf(struct t *t)
 {
 	MPFR_DECL_INIT(mx, 24);
 
-	mpfr_set_flt(mx, t->x, MPFR_RNDN);
-	t->y = mpfr_get_exp(mx) + 1;
 	t->dy = 0;
 	t->e = 0;
+	if (t->x == 0) {
+		t->y = -INFINITY;
+		t->e |= DIVBYZERO;
+		return 0;
+	}
+	if (isinf(t->x)) {
+		t->y = INFINITY;
+		return 0;
+	}
+	if (isnan(t->x)) {
+		t->y = t->x;
+		return 0;
+	}
+	mpfr_set_flt(mx, t->x, MPFR_RNDN);
+	t->y = mpfr_get_exp(mx) - 1;
 	return 0;
 }
 int mplogbl(struct t *t)
 {
 	MPFR_DECL_INIT(mx, 64);
 
-	mpfr_set_ld(mx, t->x, MPFR_RNDN);
-	t->y = mpfr_get_exp(mx) + 1;
 	t->dy = 0;
 	t->e = 0;
+	if (t->x == 0) {
+		t->y = -INFINITY;
+		t->e |= DIVBYZERO;
+		return 0;
+	}
+	if (isinf(t->x)) {
+		t->y = INFINITY;
+		return 0;
+	}
+	if (isnan(t->x)) {
+		t->y = t->x;
+		return 0;
+	}
+	mpfr_set_ld(mx, t->x, MPFR_RNDN);
+	t->y = mpfr_get_exp(mx) - 1;
+	return 0;
+}
+int mpnearbyint(struct t *t) { return mpd1(t, wrap_nearbyint) || (t->e&=~INEXACT, 0); }
+int mpnearbyintf(struct t *t) { return mpf1(t, wrap_nearbyint) || (t->e&=~INEXACT, 0); }
+int mpnearbyintl(struct t *t) { return mpl1(t, wrap_nearbyint) || (t->e&=~INEXACT, 0); }
+// TODO: hard to implement with mpfr
+int mpnextafter(struct t *t)
+{
+	feclearexcept(FE_ALL_EXCEPT);
+	t->y = nextafter(t->x, t->x2);
+	t->e = getexcept();
+	t->dy = 0;
+	return 0;
+}
+int mpnextafterf(struct t *t)
+{
+	feclearexcept(FE_ALL_EXCEPT);
+	t->y = nextafterf(t->x, t->x2);
+	t->e = getexcept();
+	t->dy = 0;
+	return 0;
+}
+int mpnextafterl(struct t *t)
+{
+	feclearexcept(FE_ALL_EXCEPT);
+	t->y = nextafterl(t->x, t->x2);
+	t->e = getexcept();
+	t->dy = 0;
+	return 0;
+}
+int mpnexttoward(struct t *t)
+{
+	feclearexcept(FE_ALL_EXCEPT);
+	t->y = nexttoward(t->x, t->x2);
+	t->e = getexcept();
+	t->dy = 0;
+	return 0;
+}
+int mpnexttowardf(struct t *t)
+{
+	feclearexcept(FE_ALL_EXCEPT);
+	t->y = nexttowardf(t->x, t->x2);
+	t->e = getexcept();
+	t->dy = 0;
 	return 0;
 }
-int mpnearbyint(struct t *t) { return mpd1(t, wrap_nearbyint); }
-int mpnearbyintf(struct t *t) { return mpf1(t, wrap_nearbyint); }
-int mpnearbyintl(struct t *t) { return mpl1(t, wrap_nearbyint); }
-int mpnextafter(struct t *t) { return -1; }
-int mpnextafterf(struct t *t) { return -1; }
-int mpnextafterl(struct t *t) { return -1; }
-int mpnexttowardl(struct t *t) { return -1; }
+int mpnexttowardl(struct t *t) { return mpnextafterl(t); }
 int mppow(struct t *t) { return mpd2(t, mpfr_pow); }
 int mppowf(struct t *t) { return mpf2(t, mpfr_pow); }
 int mppowl(struct t *t) { return mpl2(t, mpfr_pow); }
@@ -528,6 +621,7 @@ int mptanl(struct t *t) { return mpl1(t, mpfr_tan); }
 int mptanh(struct t *t) { return mpd1(t, mpfr_tanh); }
 int mptanhf(struct t *t) { return mpf1(t, mpfr_tanh); }
 int mptanhl(struct t *t) { return mpl1(t, mpfr_tanh); }
+// TODO: tgamma(2) raises wrong flags
 int mptgamma(struct t *t) { return mpd1(t, mpfr_gamma); }
 int mptgammaf(struct t *t) { return mpf1(t, mpfr_gamma); }
 int mptgammal(struct t *t) { return mpl1(t, mpfr_gamma); }
@@ -538,8 +632,23 @@ int mpj0(struct t *t) { return mpd1(t, mpfr_j0); }
 int mpj1(struct t *t) { return mpd1(t, mpfr_j1); }
 int mpy0(struct t *t) { return mpd1(t, mpfr_y0); }
 int mpy1(struct t *t) { return mpd1(t, mpfr_y1); }
-int mpscalb(struct t *t) { return -1; }
-int mpscalbf(struct t *t) { return -1; }
+// TODO: non standard functions
+int mpscalb(struct t *t)
+{
+	setupfenv(t->r);
+	t->y = scalb(t->x, t->x2);
+	t->e = getexcept();
+	t->dy = 0; // wrong
+	return 0;
+}
+int mpscalbf(struct t *t)
+{
+	setupfenv(t->r);
+	t->y = scalbf(t->x, t->x2);
+	t->e = getexcept();
+	t->dy = 0; // wrong
+	return 0;
+}
 int mpj0f(struct t *t) { return mpf1(t, mpfr_j0); }
 int mpj0l(struct t *t) { return mpl1(t, mpfr_j0); }
 int mpj1f(struct t *t) { return mpf1(t, mpfr_j1); }
@@ -673,3 +782,246 @@ int mpscalblnl(struct t *t) { return mpldexpl(t); }
 int mplgamma_r(struct t *t) { return mplgamma(t); }
 int mplgammaf_r(struct t *t) { return mplgammaf(t); }
 int mplgammal_r(struct t *t) { return mplgammal(t); }
+
+int mpilogb(struct t *t)
+{
+	MPFR_DECL_INIT(mx, 53);
+
+	mpfr_set_d(mx, t->x, MPFR_RNDN);
+	t->i = mpfr_get_exp(mx) - 1;
+	t->e = 0;
+	if (isinf(t->x) || isnan(t->x) || t->x == 0)
+		t->e = INVALID;
+	return 0;
+}
+int mpilogbf(struct t *t)
+{
+	MPFR_DECL_INIT(mx, 24);
+
+	mpfr_set_flt(mx, t->x, MPFR_RNDN);
+	t->i = mpfr_get_exp(mx) - 1;
+	t->e = 0;
+	if (isinf(t->x) || isnan(t->x) || t->x == 0)
+		t->e = INVALID;
+	return 0;
+}
+int mpilogbl(struct t *t)
+{
+	MPFR_DECL_INIT(mx, 64);
+
+	mpfr_set_ld(mx, t->x, MPFR_RNDN);
+	t->i = mpfr_get_exp(mx) - 1;
+	t->e = 0;
+	if (isinf(t->x) || isnan(t->x) || t->x == 0)
+		t->e = INVALID;
+	return 0;
+}
+
+// TODO: ll* is hard to do with mpfr
+#define mp_f_i(n) \
+int mp##n(struct t *t) \
+{ \
+	setupfenv(t->r); \
+	t->i = n(t->x); \
+	t->e = getexcept(); \
+	return 0; \
+}
+
+mp_f_i(llrint)
+mp_f_i(llrintf)
+mp_f_i(llrintl)
+mp_f_i(lrint)
+mp_f_i(lrintf)
+mp_f_i(lrintl)
+mp_f_i(llround)
+mp_f_i(llroundf)
+mp_f_i(llroundl)
+mp_f_i(lround)
+mp_f_i(lroundf)
+mp_f_i(lroundl)
+
+int mpmodf(struct t *t)
+{
+	int e, r;
+
+	r = mpd1(t, wrap_trunc);
+	if (r)
+		return r;
+	t->y2 = t->y;
+	t->dy2 = t->dy;
+	e = t->e & ~INEXACT;
+	r = mpd1(t, mpfr_frac);
+	t->e |= e;
+	return r;
+}
+
+int mpmodff(struct t *t)
+{
+	int e, r;
+
+	r = mpf1(t, wrap_trunc);
+	if (r)
+		return r;
+	t->y2 = t->y;
+	t->dy2 = t->dy;
+	e = t->e & ~INEXACT;
+	r = mpf1(t, mpfr_frac);
+	t->e |= e;
+	return r;
+}
+
+int mpmodfl(struct t *t)
+{
+	int e, r;
+
+	r = mpl1(t, wrap_trunc);
+	if (r)
+		return r;
+	t->y2 = t->y;
+	t->dy2 = t->dy;
+	e = t->e & ~INEXACT;
+	r = mpl1(t, mpfr_frac);
+	t->e |= e;
+	return r;
+}
+
+int mpsincos(struct t *t)
+{
+	int e, r;
+
+	r = mpd1(t, mpfr_cos);
+	if (r)
+		return r;
+	t->y2 = t->y;
+	t->dy2 = t->dy;
+	e = t->e;
+	r = mpd1(t, mpfr_sin);
+	t->e |= e;
+	return r;
+}
+
+int mpsincosf(struct t *t)
+{
+	int e, r;
+
+	r = mpf1(t, mpfr_cos);
+	if (r)
+		return r;
+	t->y2 = t->y;
+	t->dy2 = t->dy;
+	e = t->e;
+	r = mpf1(t, mpfr_sin);
+	t->e |= e;
+	return r;
+}
+
+int mpsincosl(struct t *t)
+{
+	int e, r;
+
+	r = mpl1(t, mpfr_cos);
+	if (r)
+		return r;
+	t->y2 = t->y;
+	t->dy2 = t->dy;
+	e = t->e;
+	r = mpl1(t, mpfr_sin);
+	t->e |= e;
+	return r;
+}
+
+int mpremquo(struct t *t) { return mpd2(t, wrap_remquo) || (t->i = mpremquo_q, 0); }
+int mpremquof(struct t *t) { return mpf2(t, wrap_remquo) || (t->i = mpremquo_q, 0); }
+int mpremquol(struct t *t) { return mpl2(t, wrap_remquo) || (t->i = mpremquo_q, 0); }
+
+int mpfma(struct t *t)
+{
+	int tn;
+	int r = rmap(t->r);
+	MPFR_DECL_INIT(mx, 53);
+	MPFR_DECL_INIT(mx2, 53);
+	MPFR_DECL_INIT(mx3, 53);
+	MPFR_DECL_INIT(my, 128);
+
+	mpsetup();
+	mpfr_clear_flags();
+	mpfr_set_d(mx, t->x, MPFR_RNDN);
+	mpfr_set_d(mx2, t->x2, MPFR_RNDN);
+	mpfr_set_d(mx3, t->x3, MPFR_RNDN);
+	tn = mpfr_fma(my, mx, mx2, mx3, r);
+	gend(t, my, tn, r);
+	if ((t->e & INEXACT) && nextafter(t->y, 0) == 0) {
+		mpfr_set_emin(-(1<<20));
+		tn = mpfr_fma(my, mx, mx2, mx3, r);
+		mpfr_mul_2si(my, my, 1074, MPFR_RNDN);
+		t->dy = scalbnl(t->y, 1074) - mpfr_get_ld(my, r);
+		mpfr_set_emin(-1073);
+	}
+	return 0;
+}
+
+int mpfmaf(struct t *t)
+{
+	int tn;
+	int r = rmap(t->r);
+	MPFR_DECL_INIT(mx, 24);
+	MPFR_DECL_INIT(mx2, 24);
+	MPFR_DECL_INIT(mx3, 24);
+	MPFR_DECL_INIT(my, 128);
+
+	mpsetupf();
+	mpfr_clear_flags();
+	mpfr_set_flt(mx, t->x, MPFR_RNDN);
+	mpfr_set_flt(mx2, t->x2, MPFR_RNDN);
+	mpfr_set_flt(mx3, t->x3, MPFR_RNDN);
+	tn = mpfr_fma(my, mx, mx2, mx3, r);
+	genf(t, my, tn, r);
+	if ((t->e & INEXACT) && nextafterf(t->y, 0) == 0) {
+		mpfr_set_emin(-(1<<20));
+		tn = mpfr_fma(my, mx, mx2, mx3, r);
+		mpfr_mul_2si(my, my, 149, MPFR_RNDN);
+		t->dy = scalbnl(t->y, 149) - mpfr_get_ld(my, r);
+		mpfr_set_emin(-148);
+	}
+	return 0;
+}
+
+int mpfmal(struct t *t)
+{
+#if LDBL_MANT_DIG == 53
+	return mpfma(t);
+#elif LDBL_MANT_DIG == 64
+	int tn;
+	int r = rmap(t->r);
+	MPFR_DECL_INIT(mx, 64);
+	MPFR_DECL_INIT(mx2, 64);
+	MPFR_DECL_INIT(mx3, 64);
+	MPFR_DECL_INIT(my, 128);
+
+	mpsetupl();
+	mpfr_clear_flags();
+	mpfr_set_ld(mx, t->x, MPFR_RNDN);
+	mpfr_set_ld(mx2, t->x2, MPFR_RNDN);
+	mpfr_set_ld(mx3, t->x3, MPFR_RNDN);
+	tn = mpfr_fma(my, mx, mx2, mx3, r);
+	genl(t, my, tn, r);
+	if ((t->e & INEXACT) && nextafterl(t->y, 0) == 0) {
+		mpfr_set_emin(-(1<<20));
+		tn = mpfr_fma(my, mx, mx2, mx3, r);
+		mpfr_mul_2si(my, my, 16445, MPFR_RNDN);
+		t->dy = scalbnl(t->y, 16445) - mpfr_get_ld(my, r);
+		mpfr_set_emin(-16444);
+	}
+	return 0;
+#else
+	return -1;
+#endif
+}
+
+int mpjn(struct t *t) { mpbessel_n = t->i; return mpd1(t, wrap_jn); }
+int mpjnf(struct t *t) { mpbessel_n = t->i; return mpf1(t, wrap_jn); }
+int mpjnl(struct t *t) { mpbessel_n = t->i; return mpl1(t, wrap_jn); }
+int mpyn(struct t *t) { mpbessel_n = t->i; return mpd1(t, wrap_yn); }
+int mpynf(struct t *t) { mpbessel_n = t->i; return mpf1(t, wrap_yn); }
+int mpynl(struct t *t) { mpbessel_n = t->i; return mpl1(t, wrap_yn); }
+