diff options
-rw-r--r-- | NEWS | 3 | ||||
-rw-r--r-- | doc/mpc.texi | 2 | ||||
-rw-r--r-- | src/pow.c | 32 | ||||
-rw-r--r-- | tests/pow.dat | 17 | ||||
-rw-r--r-- | tests/pow_ui.dat | 2 |
5 files changed, 51 insertions, 5 deletions
@@ -1,6 +1,9 @@ Changes in version 1.0: - First release as a GNU package - License change: LGPLv3+ for code, GFDLv1.3+ for documentation + - 0^0, which returned (NaN,NaN) previously, now returns (1,0), + with the sign of the imaginary part of the result being the opposite of + the sign of the imaginary part of the exponent - Bug fixes: - div and norm now return a value indicating the effective rounding direction, as the other functions. diff --git a/doc/mpc.texi b/doc/mpc.texi index 26e4a36..dfa74ec 100644 --- a/doc/mpc.texi +++ b/doc/mpc.texi @@ -947,6 +947,8 @@ to @var{rnd}. For @code{mpc_pow_d}, @code{mpc_pow_ld}, @code{mpc_pow_si}, @code{mpc_pow_ui}, @code{mpc_pow_z} and @code{mpc_pow_fr}, the imaginary part of @var{op2} is considered as +0. +When both @var{op1} and @var{op2} are zero, the result has real part 1, +and imaginary part 0, with sign being the opposite of that of @var{op2}. @end deftypefun @deftypefun int mpc_exp (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd}) @@ -489,10 +489,36 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) if (y_real && mpfr_zero_p (mpc_realref(y))) /* case y zero */ { - if (x_real && mpfr_zero_p (mpc_realref(x))) /* 0^0 = NaN +i*NaN */ + if (x_real && mpfr_zero_p (mpc_realref(x))) { - mpfr_set_nan (mpc_realref(z)); - mpfr_set_nan (mpc_imagref(z)); + /* we define 0^0 to be (1, -sign(Im(y))*0) since the real part is + coherent with MPFR where 0^0 gives 1, and the sign of the + imaginary part is the most frequent one obtained for random + tiny x and y, for example with the following Sage program: + D = dict() + for sa in [-1,1]: + for sb in [-1,1]: + for sc in [-1,1]: + for sd in [-1,1]: + D[(sa,sb,sc,sd)] = [0,0] + def tiny(): + return (random() - 0.5) * 1e-3 + for n in range(10^4): + a = tiny() + b = tiny() + c = tiny() + d = tiny() + x = a+I*b + y = c+I*d + z = x^y + t = sign(a),sign(b),sign(c),sign(d) + if sign(z.imag()) == 1: + D[t] = [D[t][0],D[t][1]+1] + else: + D[t] = [D[t][0]+1,D[t][1]] */ + mpfr_set_ui (mpc_realref(z), 1, GMP_RNDZ); + mpfr_set_zero (mpc_imagref(z), + mpfr_signbit (mpc_imagref(y)) ? 1 : -1); return 0; } else /* x^0 = 1 +/- i*0 even for x=NaN see algorithms.tex for the diff --git a/tests/pow.dat b/tests/pow.dat index 10b6ee3..6331496 100644 --- a/tests/pow.dat +++ b/tests/pow.dat @@ -133,7 +133,22 @@ 0 0 53 0 53 0 53 -inf 53 -1 53 -inf 53 +1 N N 0 0 53 0 53 0 53 -inf 53 -1 53 -inf 53 -1 N N -0 0 53 nan 53 nan 53 +0 53 +0 53 +0 53 +0 N N +0 0 53 1 53 -0 53 +0 53 +0 53 +0 53 +0 N N +0 0 53 1 53 -0 53 +0 53 +0 53 -0 53 +0 N N +0 0 53 1 53 -0 53 +0 53 -0 53 +0 53 +0 N N +0 0 53 1 53 -0 53 +0 53 -0 53 -0 53 +0 N N +0 0 53 1 53 -0 53 -0 53 +0 53 +0 53 +0 N N +0 0 53 1 53 -0 53 -0 53 +0 53 -0 53 +0 N N +0 0 53 1 53 -0 53 -0 53 -0 53 +0 53 +0 N N +0 0 53 1 53 -0 53 -0 53 -0 53 -0 53 +0 N N +0 0 53 1 53 +0 53 +0 53 +0 53 +0 53 -0 N N +0 0 53 1 53 +0 53 +0 53 +0 53 -0 53 -0 N N +0 0 53 1 53 +0 53 +0 53 -0 53 +0 53 -0 N N +0 0 53 1 53 +0 53 +0 53 -0 53 -0 53 -0 N N +0 0 53 1 53 +0 53 -0 53 +0 53 +0 53 -0 N N +0 0 53 1 53 +0 53 -0 53 +0 53 -0 53 -0 N N +0 0 53 1 53 +0 53 -0 53 -0 53 +0 53 -0 N N +0 0 53 1 53 +0 53 -0 53 -0 53 -0 53 -0 N N 0 0 53 nan 53 nan 53 +0 53 +0 53 +0 53 +1 N N 0 0 53 nan 53 nan 53 +0 53 +0 53 +0 53 -1 N N 0 0 53 0 53 0 53 +0 53 +0 53 +1 53 +0 N N diff --git a/tests/pow_ui.dat b/tests/pow_ui.dat index 26fc41b..9ecde18 100644 --- a/tests/pow_ui.dat +++ b/tests/pow_ui.dat @@ -34,7 +34,7 @@ 0 0 53 +inf 53 nan 53 -inf 53 +1 +1 N N 0 0 53 +inf 53 nan 53 -inf 53 -1 +1 N N -0 0 53 nan 53 nan 53 +0 53 +0 +0 N N +0 0 53 +1 53 -0 53 +0 53 +0 +0 N N 0 0 53 0 53 0 53 +0 53 +0 +1 N N 0 0 53 +1 53 +0 53 +0 53 +1 +0 N N |