From f04c4ccf618dfb01fe6d878f602006a643f13071 Mon Sep 17 00:00:00 2001 From: Andreas Enge Date: Mon, 28 Nov 2022 12:48:56 +0100 Subject: [PATCH] Fix bug in mpc_pow_fr. * tests/pow_fr.dat: Correct test and add more tests. * src/pow.c (mpc_pow): Correct sign of zero part in result for c*(1+I) or c*(1-I) raised to an even positive integer. * doc/algorithms.tex: Add comment concerning the sign. * NEWS: Add entry. --- NEWS | 3 + doc/algorithms.tex | 9 +++ src/pow.c | 146 +++++++++++++++++++++++++-------------------- tests/pow_fr.dat | 7 ++- 4 files changed, 99 insertions(+), 66 deletions(-) diff --git a/NEWS b/NEWS index ba37c9d..6926338 100644 --- a/NEWS +++ b/NEWS @@ -6,6 +6,9 @@ Changes in version 1.3.0: - New experimental function: mpc_eta_fund - Bug fixes: - mpc_asin for asin(z) with small |Re(z)| and tiny |Im(z)| + - mpc_pow_fr: sign of zero part of result when the base has up to sign + the same real and imaginary part, and the exponent is an even + positive integer - mpc_fma Changes in version 1.2.1, released in October 2020: diff --git a/doc/algorithms.tex b/doc/algorithms.tex index aaff643..4cdbb09 100644 --- a/doc/algorithms.tex +++ b/doc/algorithms.tex @@ -2126,6 +2126,15 @@ the determined cases: where $\sigma_2$ (resp $\rho_1$, $\rho_2$) is the sign of $x_2$ (resp. $y_1$, $y_2$) and with the convention $0^0=+1$. +FIXME: This misses the cases $(x_1 \pm x_1 i)^{y_1 + 0 i}$ for $y_1 \geq 2$ +even and $x_1 \neq 0$; the case (d)(ii) above applies when $x_2 > 0$, but +misses $x_2 < 0$. In any case, the sign of the imaginary part +(for $y_1$ divisible by~$4$) +or the real part +(for $y_1$ divisible by~$2$, but not by~$4$) +cannot be determined. + + \subsection {\texttt {mpc\_pow\_ui}} \label {ssec:mpcpowui} diff --git a/src/pow.c b/src/pow.c index 0fc6932..2bab8b8 100644 --- a/src/pow.c +++ b/src/pow.c @@ -1,6 +1,6 @@ /* mpc_pow -- Raise a complex number to the power of another complex number. -Copyright (C) 2009, 2010, 2011, 2012, 2014, 2015, 2016, 2018, 2020 INRIA +Copyright (C) 2009, 2010, 2011, 2012, 2014, 2015, 2016, 2018, 2020, 2022 INRIA This file is part of GNU MPC. @@ -481,7 +481,8 @@ is_odd (mpfr_srcptr y, mpfr_exp_t k) int mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) { - int ret = -2, loop, x_real, x_imag, y_real, z_real = 0, z_imag = 0; + int ret = -2, loop, x_real, x_imag, y_real, z_real = 0, z_imag = 0, + ramified = 0; mpc_t t, u; mpfr_prec_t p, pr, pi, maxprec; int saved_underflow, saved_overflow; @@ -640,6 +641,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) if (mpfr_cmpabs (mpc_realref(x), mpc_imagref(x)) == 0 && y_real && mpfr_integer_p (mpc_realref(y)) && is_odd (mpc_realref(y), 0) == 0) { + ramified = 1; if (is_odd (mpc_realref(y), -1)) /* y/2 is odd */ z_imag = 1; else @@ -764,78 +766,94 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) if (z_real) { - /* When the result is real (see algorithm.tex for details), + /* When the result is real (see algorithm.tex for details) and + x=x1 * (1 \pm i), y a positive integer divisible by 4, then + Im(x^y) = 0i with a sign that cannot be determined (and is thus + chosen as _1). Otherwise, Im(x^y) = + sign(imag(y))*0i, if |x| > 1 + sign(imag(x))*sign(real(y))*0i, if |x| = 1 - sign(imag(y))*0i, if |x| < 1 */ - mpfr_t n; - int inex, cx1; - int sign_zi, sign_rex, sign_imx; - /* cx1 < 0 if |x| < 1 - cx1 = 0 if |x| = 1 - cx1 > 0 if |x| > 1 - */ - - sign_rex = mpfr_signbit (mpc_realref (x)); - sign_imx = mpfr_signbit (mpc_imagref (x)); - mpfr_init (n); - inex = mpc_norm (n, x, MPFR_RNDN); - cx1 = mpfr_cmp_ui (n, 1); - if (cx1 == 0 && inex != 0) - cx1 = -inex; - - sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0) - || (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y))) - || (cx1 > 0 && mpfr_signbit (mpc_imagref (y))); - - /* copy RE(y) to n since if z==y we will destroy Re(y) below */ - mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y))); - mpfr_set (n, mpc_realref (y), MPFR_RNDN); - ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd)); - if (y_real && (x_real || x_imag)) - { - /* FIXME: with y_real we assume Im(y) is really 0, which is the case - for example when y comes from pow_fr, but in case Im(y) is +0 or - -0, we might get different results */ - mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd)); - fix_sign (z, sign_rex, sign_imx, n); - ret = MPC_INEX(ret, 0); /* imaginary part is exact */ - } - else - { - inex_im = mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd)); - ret = MPC_INEX (ret, inex_im); - /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */ - if (MPC_RND_IM (rnd) == MPFR_RNDD || sign_zi) - mpc_conj (z, z, MPC_RNDNN); - } + if (ramified) + ret = MPC_INEX ( + mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd)), + mpfr_set_ui (mpc_imagref (z), 0, MPFR_RNDN)); + else { + mpfr_t n; + int inex, cx1; + int sign_zi, sign_rex, sign_imx; + /* cx1 < 0 if |x| < 1 + cx1 = 0 if |x| = 1 + cx1 > 0 if |x| > 1 + */ + + sign_rex = mpfr_signbit (mpc_realref (x)); + sign_imx = mpfr_signbit (mpc_imagref (x)); + mpfr_init (n); + inex = mpc_norm (n, x, MPFR_RNDN); + cx1 = mpfr_cmp_ui (n, 1); + if (cx1 == 0 && inex != 0) + cx1 = -inex; + + sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0) + || (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y))) + || (cx1 > 0 && mpfr_signbit (mpc_imagref (y))); + + /* copy RE(y) to n since if z==y we will destroy Re(y) below */ + mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y))); + mpfr_set (n, mpc_realref (y), MPFR_RNDN); + ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd)); + if (y_real && (x_real || x_imag)) + { + /* FIXME: with y_real we assume Im(y) is really 0, which is the case + for example when y comes from pow_fr, but in case Im(y) is +0 or + -0, we might get different results */ + mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd)); + fix_sign (z, sign_rex, sign_imx, n); + ret = MPC_INEX(ret, 0); /* imaginary part is exact */ + } + else + { + inex_im = mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd)); + ret = MPC_INEX (ret, inex_im); + /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */ + if (MPC_RND_IM (rnd) == MPFR_RNDD || sign_zi) + mpc_conj (z, z, MPC_RNDNN); + } - mpfr_clear (n); + mpfr_clear (n); + } } else if (z_imag) { - ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd)); - /* if z is imaginary and y real, then x cannot be real */ - if (y_real && x_imag) - { - int sign_rex = mpfr_signbit (mpc_realref (x)); - - /* If z overlaps with y we set Re(z) before checking Re(y) below, - but in that case y=0, which was dealt with above. */ - mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd)); - /* Note: fix_sign only does something when y is an integer, - then necessarily y = 1 or 3 (mod 4), and in that case the - sign of Im(x) is irrelevant. */ - fix_sign (z, sign_rex, 0, mpc_realref (y)); - ret = MPC_INEX(0, ret); - } + if (ramified) + ret = MPC_INEX ( + mpfr_set_ui (mpc_realref (z), 0, MPFR_RNDN), + mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_RE(rnd))); else - { - inex_re = mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd)); - ret = MPC_INEX(inex_re, ret); - } + { + ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd)); + /* if z is imaginary and y real, then x cannot be real */ + if (y_real && x_imag) + { + int sign_rex = mpfr_signbit (mpc_realref (x)); + + /* If z overlaps with y we set Re(z) before checking Re(y) below, + but in that case y=0, which was dealt with above. */ + mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd)); + /* Note: fix_sign only does something when y is an integer, + then necessarily y = 1 or 3 (mod 4), and in that case the + sign of Im(x) is irrelevant. */ + fix_sign (z, sign_rex, 0, mpc_realref (y)); + ret = MPC_INEX(0, ret); + } + else + { + inex_re = mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd)); + ret = MPC_INEX(inex_re, ret); + } + } } else ret = mpc_set (z, u, rnd); diff --git a/tests/pow_fr.dat b/tests/pow_fr.dat index 66162c9..632c39a 100644 --- a/tests/pow_fr.dat +++ b/tests/pow_fr.dat @@ -1,6 +1,6 @@ # Data file for mpc_pow_fr. # -# Copyright (C) 2011 INRIA +# Copyright (C) 2011, 2022 INRIA # # This file is part of GNU MPC. # @@ -74,4 +74,7 @@ 0 - 8 -0 8 -0x89p-10 2 +0 2 0x3p-2 3 7 N N # issue revealed by random tests (with GMP_CHECK_RANDOMIZE=1669437260) -- 0 2 -0x1p-28 2 +0 2 0x1.8p-8 2 0x1.8p-8 2 4 N N +- 0 2 -0x3p-29 2 +0 2 0x1.8p-8 2 0x1.8p-8 2 4 N N +- 0 2 -0x3p-29 2 +0 2 0x1.8p-8 2 -0x1.8p-8 2 4 N N +0 - 2 +0 2 0x1p-14 2 0x1.8p-8 2 0x1.8p-8 2 2 N N +0 + 2 +0 2 -0x1p-14 2 0x1.8p-8 2 -0x1.8p-8 2 2 N N -- 2.33.0