diff options
-rw-r--r-- | doc/algorithms.tex | 29 | ||||
-rw-r--r-- | src/pow.c | 25 | ||||
-rw-r--r-- | tests/Makefile.am | 2 | ||||
-rw-r--r-- | tests/pow.dat | 12 |
4 files changed, 66 insertions, 2 deletions
diff --git a/doc/algorithms.tex b/doc/algorithms.tex index 4dc9095..a1d09e2 100644 --- a/doc/algorithms.tex +++ b/doc/algorithms.tex @@ -1175,6 +1175,35 @@ A relevant reference is \cite{BrDiJeLeMeMuReStTo09}, especially Section 4.5 which discusses complex floating-point numbers, and gives error bounds for multiplication, division and square root. +\paragraph{Sign of zeroes.} +When the output value has a zero real or imaginary part, its sign should be +decided. When the inputs also have a zero real or imaginary part, we +consider all possible limits, and if all those limits give the same sign, +we take this as the sign of the zero part. +Otherwise, we round to $+0$, except for rounding toward $-\infty$, where we +round to $-0$. + +Example~: consider $x = 1 + 0 i$ and $y = -0 + 0 i$, we consider that $x$ +is the limit of $1 + \epsilon i$ for $\epsilon > 0$ tending to zero. +Similarly, $y = -\delta + \gamma i$ with $\delta, \gamma > 0$. +Now $\log x \approx \epsilon^2/2 + \epsilon i$, thus +$y \log x \approx (-\epsilon^2 \delta/2 - \epsilon \gamma) ++ i (\epsilon^2 \gamma/2 - epsilon \delta)$. +Thus if we neglect terms of order $3$ or more, +whatever the relative growth of $\epsilon, \delta, \gamma$, the real and +imaginary parts of $y \log x$ are negative, thus we decide $x^y$ is rounded +to $1 - 0 i$. + +For $x = 1 - 0 i$ and $y = 0 + 0 i$, we write +$x = 1 - \epsilon i$ and $y = \delta + \gamma i$, +which gives $\log x \approx \epsilon^2/2 - \epsilon i$ +and $y \log x \approx \epsilon \gamma - i \epsilon \delta + O(\epsilon^2)$, +thus we also round $x^y$ to $1 - 0 i$. + +% (1 -0)^(-0 -0): +% x = 1 - epsilon i, y = -delta-gamma*i +% y log(x) = -epsilon gamma + i epsilon delta + O(epsilon^2) + \bibliographystyle{acm} \bibliography{algorithms} @@ -473,8 +473,33 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) if (y_real && (mpfr_integer_p (MPC_RE(y)) || mpfr_cmp_ui (MPC_RE(x), 0) >= 0)) { + int s1, s2; + /* if x = x + eps*i and y = y + delta*i, then for x > 0: + log(x) ~ log(x) + eps/x*i thus + y*log(x) ~ y*log(x) + [y/x*eps + log(x)*delta]*i + For x < 0 and y integer: + log(x) ~ log|x| + [Pi + eps/x]*i thus + y*log(x) ~ y*log|x| + [y*Pi + y/x*eps + log|x|*delta]*i + thus if y is even we get the same term as for x > 0 for the + imaginary part, otherwise its opposite. */ + s1 = mpfr_signbit (MPC_RE(x)) ? -1 : 1; + s1 *= mpfr_signbit (MPC_RE(y)) ? -1 : 1; /* sign of y/x */ + s1 *= MPFR_SIGN(MPC_IM(x)); /* sign of y/x*eps */ + s2 = mpfr_get_exp (MPC_RE(x)) >= 1 ? 1 : -1; /* sign of log|x| */ + s2 *= MPFR_SIGN(MPC_IM(y)); /* sign of log(x)*delta */ + if (s1 != s2) + { + if (MPC_RND_IM(rnd) == GMP_RNDD) + s1 = -1; + else + s1 = 1; /* take +0 as arbitrary sign */ + } + if (MPFR_SIGN(MPC_RE(x)) < 0 && is_odd (MPC_RE(y), 0)) + s1 = -s1; ret = mpfr_pow (MPC_RE(z), MPC_RE(x), MPC_RE(y), MPC_RND_RE(rnd)); ret = MPC_INEX(ret, mpfr_set_ui (MPC_IM(z), 0, MPC_RND_IM(rnd))); + if (s1 == -1) + mpfr_neg (MPC_IM(z), MPC_IM(z), MPC_RND_IM(rnd)); goto end; } diff --git a/tests/Makefile.am b/tests/Makefile.am index b01ca1f..d9e8394 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -8,7 +8,7 @@ check_PROGRAMS = tabs tadd tadd_fr tadd_ui targ tconj tcos tcosh tdiv \ tdiv_2exp tdiv_fr tdiv_ui texp tfr_div tfr_sub tget_version timag tio_str \ tlog tmul tmul_2exp tmul_fr tmul_i tmul_si tmul_ui tneg tnorm tpow tprec \ tproj treal treimref tset tsin tsinh tsqr tsqrt tstrtoc tsub tsub_fr tsub_ui \ -ttan ttanh tui_div tui_ui_sub +ttan ttanh tui_div tui_ui_sub check_LTLIBRARIES=libmpc-tests.la libmpc_tests_la_SOURCES=mpc-tests.h random.c tgeneric.c read_data.c \ diff --git a/tests/pow.dat b/tests/pow.dat index da4f19e..d58865a 100644 --- a/tests/pow.dat +++ b/tests/pow.dat @@ -161,8 +161,13 @@ 0 0 53 +inf 53 nan 53 +0 53 +0 53 -inf 53 -1 N N 0 0 53 +1 53 +0 53 +0 53 +1 53 +0 53 +0 N N -# (1-0*I)^(+0+0*I) -> 1-0*I +# zero cases from algorithms.tex +# (1+0*i)^(-0+0*i) -> 1-0*i +0 0 53 +1 53 -0 53 +1 53 +0 53 -0 53 +0 N N +# (1-0*i)^(+0+0*i) -> 1-0*i 0 0 53 +1 53 -0 53 +1 53 -0 53 +0 53 +0 N N +# (1-0*i)^(-0-0*i) -> 1+0*i +0 0 53 +1 53 +0 53 +1 53 -0 53 -0 53 -0 N N # exact cases # (-4)^(1/4) = 1+i @@ -176,6 +181,11 @@ 0 0 2 -0x1p-18 2 0x1p-18 2 -0x1p14 2 0 3 -0x5p-2 3 0 N N # e=2 n=5 0 0 2 -0x1p-13 2 0x1p-13 2 -0x1p10 2 0 3 -0x5p-2 3 0 N N +# (+2 +0)^(-3 -0) -> (-1/8 -0) +# x = 2 + epsilon*i, y = -3 - delta*i +# log(x) = log(2) + epsilon/2*i + O(epsilon^2) +# y*log(x) = [-3*log(2) + o(1)] + [-3*epsilon/2-delta*log(2)]*i +0 0 2 0x1p-3 2 -0 2 2 2 +0 2 -3 2 -0 N N 0 0 2 +0 2 -2 2 +0 2 0x1p-1 2 -1 2 -0 N N 0 - 2 +0 3 -5 2 +0 53 0xCCCCCCCCCCCCDp-54 2 -1 2 -0 N N - 0 5 -25 2 +0 2 +0 53 0xCCCCCCCCCCCCDp-54 2 -2 2 -0 N N |