diff options
-rw-r--r-- | doc/algorithms.tex | 186 | ||||
-rw-r--r-- | src/pow.c | 9 | ||||
-rw-r--r-- | tests/pow.dat | 142 |
3 files changed, 327 insertions, 10 deletions
diff --git a/doc/algorithms.tex b/doc/algorithms.tex index a1d09e2..14b29ef 100644 --- a/doc/algorithms.tex +++ b/doc/algorithms.tex @@ -1188,7 +1188,7 @@ 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)$. ++ 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 @@ -1204,6 +1204,190 @@ thus we also round $x^y$ to $1 - 0 i$. % x = 1 - epsilon i, y = -delta-gamma*i % y log(x) = -epsilon gamma + i epsilon delta + O(epsilon^2) + + + +The sign of zero parts are chosen so that they are consistent with the formula +$x^y = \exp(y\log x)$. +Let $\phi \in [-\pi, +\pi]$ the argument of $x = |x| e^{i\phi}$, $y_1$ +(resp. $y_2$) the real (resp. imaginary) part of $y = y_1 + y_2 i$. +Then +\[ +x^y=\exp(A(x,y)) (\cos B(x,y)+\sin B(x,y) i) +\] where +\begin {eqnarray*} +A(x,y) & = & y_1\log|x|-y_2\phi,\\ +B(x,y) & = & y_2\log|x|+y_1\phi. +\end {eqnarray*} +As $|x^y| = \exp(A(x,y))$ is positive, the value of $B(x,y)$ determines the +sign of each part of $x^y$. +Special study is needed around $B(x, y)$ values of the form $k \pi/2$ to +determine the sign of the zero part: let $Q_{x_0}$ (resp. $Q_{y_0}$) the +quadrant where $x_0$ (resp. $y_0$) lies, when $B(x_0,y_0)$ lies around zero, +if $B(x,y)$ stays non negative for $x$ and $y$ in an open domain of $Q_{x_0} +\times Q_{y_0}$ containing $(x_0, y_0)$ then $\sin B(x_0, y_0) = \sin(+0) = ++0$, if it stays non positive then $\sin B(x_0, y_0) = \sin (-0) = -0$. +If $B(x,y)$ tends to $+\pi$ (resp. $-\pi$) staying lesser (resp. greater) than +it when $(x,y)$ tends to $(x_0, y_0)$, then $\sin B(x_0,y_0) = +0$. +Conversely, if $B(x,y)$ tend to $+\pi$ (resp. $-\pi$) remaining greater (resp. +lesser) than it, then $\sin B(x_0,y_0) = -0$. +The same applies to the real part when $B(x,y)$ lies around $\pm \pi/2$. + +First, let us solve $B(x,y) = 0$. +\begin {enumerate} +\item +If +\[ +y_2 \log |x| = y_1 \phi = +0, +\] +then +\[ +x^y = +\exp(A(x,y)) +0i, +\] +the solutions are summarized in the following table: +\[ +\begin {array}{l c r c l l l l} +(x_1 -0i)^{y_1 -0i} & = & x_1^{y_1} +0i +& \text {with} & 0 < x_1 < 1 & & y_1 <0 \\ + +(x_1 -0i)^{-0 -0i} & = & +1 +0i +& \text {with} & |x_1| < 1 \\ + +x^{-0 -0i} & = & +1 +0i +& \text {with} & |x| < 1 & x_2 < 0 \\ + +x^{+0 -0i} & = & +1 +0i +& \text {with} & |x| < 1 & x_2 > 0 \\ + +(x_1 +0i)^{+0 -0i} & = & +1 +0i +& \text {with} & |x_1| < 1\\ + +(x_1 +0i)^{y_1 -0i} & = & x_1^{y_1} +0i +& \text {with} & 0 < x_1 < 1 & & y_1 > 0 \\ + +\hline + +(x_1 -0i)^{y_1 +0i} & = & x_1^{y_1} +0i +& \text {with} & x_1 \geq 1 & & y_1 <0 \\ + +(x_1 -0i)^{-0 +0i} & = & +1 +0i +& \text {with} & |x_1| \geq 1 \\ + +x^{-0 +0i} & = & +1 +0i +& \text {with} & |x| \geq 1 & x_2 < 0 \\ + +x^{+0 +0i} & = & +1 +0i +& \text {with} & |x| \geq 1 & x_2 > 0 \\ + +(x_1 +0i)^{+0 +0i} & = & +1 +0i +& \text {with} & |x_1| \geq 1\\ + +(x_1 +0i)^{y_1 +0i} & = & x_1^{y_1} +0i +& \text {with} & x_1 \geq 1 & & y_1 > 0 \\ + +\hline + +(+1 -0i)^{y_1 +y_2i} & = & +1 +0i +& \text {with} & & & y_1 <0 & y_2 > 0\\ + +(\pm 1 -0i)^{-0 +y_2i} & = & \exp (y_2 \phi) +0i +& \text {with} & & & & y_2 > 0\\ + +x^{-0 +y_2i} & = & \exp (y_2 \phi) +0i +& \text {with} & |x| = 1 & x_2 < 0 & & y_2 > 0 \\ + +x^{+0 +y_2i} & = & \exp (y_2 \phi) +0i +& \text {with} & |x| = 1 & x_2 > 0 & & y_2 > 0\\ + +(\pm1 +0i)^{+0 +y_2i} & = & \exp (y_2 \phi) +0i +& \text {with} & & & & y_2 > 0\\ + +(+1 +0i)^{y_1 +y_2i} & = & +1 +0i +& \text {with} & & & y_1 > 0 & y_2 > 0 +\end {array} +\] + +\item +If +\[ +y_2 \log |x| = y_1 \phi = -0, +\] +then +\[ +x^y = +\exp(A(x,y)) -0i, +\] +the solutions are summarized in the following table: +\[ +\begin {array}{l c r c l l l l} +(+1 +0i)^{y_1 +y_2i} & = & +1 -0i +& \text {with} & & & y_1 < 0 & y_2 < 0\\ + +(\pm 1 +0i)^{-0 +y_2i} & = & \exp(y_2 \phi) -0i +& \text {with} & & & & y_2 <0 \\ + +x^{-0 +y_2i} & = & \exp(y_2 \phi) -0i +& \text {with} & |x| = 1 & x_2 > 0 & & y_2 <0 \\ + +x^{+0 +y_2i} & = & \exp(y_2 \phi) -0i +& \text {with} & |x| = 1 & x_2 < 0 & & y_2 <0 \\ + +(\pm 1 -0i)^{+0 +y_2i} & = & \exp(y_2 \phi) -0i +& \text {with} & & & & y_2 < 0\\ + +(+1 -0i)^{y_1 +y_2i} & = & +1 -0i +& \text {with} & & & y_1 > 0 & y_2 < 0\\ + +\hline + +(x_1 +0i)^{y_1 -0i} & = & x^{y_1} -0i +& \text {with} & x_1 \geq 1 & & y_1 <0 \\ + +(x_1 +0i)^{-0 -0i} & = & +1 -0i +& \text {with} & |x_1| \geq 1 \\ + +x^{-0 -0i} & = & +1 -0i +& \text {with} & |x| \geq 1 & x_2 > 0 \\ + +x^{+0 -0i} & = & +1 -0i +& \text {with} & |x| \geq 1 & x_2 < 0 \\ + +(x_1 -0i)^{+0 -0i} & = & +1 -0i +& \text {with} & |x_1| \geq 1\\ + +(x_1 -0i)^{y_1 -0i} & = & x_1^{y_1} -0i +& \text {with} & x_1 \geq 1 & & y_1 > 0 \\ + +\hline + +(x_1 +0i)^{y_1 +0i} & = & x_1^{y_1} -0i +& \text {with} & 0 < x_1 < 1 & & y_1 <0 \\ + +(x_1 +0i)^{-0 +0i} & = & +1 -0i +& \text {with} & |x_1| < 1 \\ + +x^{-0 +0i} & = & +1 -0i +& \text {with} & |x| < 1 & x_2 > 0 \\ + +x^{+0 +0i} & = & +1 -0i +& \text {with} & |x| < 1 & x_2 < 0 \\ + +(x_1 -0i)^{+0 +0i} & = & +1 -0i +& \text {with} & |x_1| < 1 \\ + +(x_1 -0i)^{y_1 +0i} & = & x_1^{y_1} -0i +& \text {with} & 0 < x_1 < 1 & & y_1 > 0 \\ +\end {array} +\] + +\item +Is it possible that +\[ +\frac{\log |x|}{\phi} = -\frac{y_1}{y_2} \neq 0? +\] +In other words, can $\frac{\log |x|}{\arg (x)}$ be a rational number when $x$ +is a dyadic complex? +\end {enumerate} + \bibliographystyle{acm} \bibliography{algorithms} @@ -463,7 +463,10 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) /* Special case 1^y = 1 */ if (mpfr_cmp_ui (MPC_RE(x), 1) == 0) { - ret = mpc_set_ui (z, 1, rnd); + ret = mpc_set_ui (z, +1, rnd); + if (mpfr_signbit (MPC_IM (y)) + && mpfr_signbit (MPC_RE (y)) != mpfr_signbit (MPC_IM (x))) + mpc_conj (z, z, MPC_RNDNN); goto end; } @@ -635,8 +638,12 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) if (z_real) { + int s = mpfr_signbit (MPC_IM (y)) + && mpfr_signbit (MPC_RE (y)) != mpfr_signbit (MPC_IM (x)); ret = mpfr_set (MPC_RE(z), MPC_RE(u), MPC_RND_RE(rnd)); ret = MPC_INEX(ret, mpfr_set_ui (MPC_IM(z), 0, MPC_RND_IM(rnd))); + if (s) + mpc_conj (z, z, MPC_RNDNN); } else if (z_imag) { diff --git a/tests/pow.dat b/tests/pow.dat index 7baaad0..281853a 100644 --- a/tests/pow.dat +++ b/tests/pow.dat @@ -169,6 +169,128 @@ # (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 +# (x-0i)^(y-0i) = +a +0i when 0 < x < 1 and y <= -0 +0 0 53 4 53 +0 53 +0.5 53 -0 53 -2 53 -0 N N +0 0 53 1 53 +0 53 +0.5 53 -0 53 -0 53 -0 N N +# x^(-0-0i) = +1 +0i when |x| < 1 and Im(x) < 0 +0 0 53 1 53 +0 53 -0.5 53 -0.5 53 -0 53 -0 N N +0 0 53 1 53 +0 53 -0 53 -0.5 53 -0 53 -0 N N +0 0 53 1 53 +0 53 +0 53 -0.5 53 -0 53 -0 N N +0 0 53 1 53 +0 53 +0.5 53 -0.5 53 -0 53 -0 N N +# x^(+0-0i) = +1 +0i when |x| < 1 and Im(x) > 0 +0 0 53 1 53 +0 53 +0.5 53 +0.5 53 +0 53 -0 N N +0 0 53 1 53 +0 53 +0 53 +0.5 53 +0 53 -0 N N +0 0 53 1 53 +0 53 -0 53 +0.5 53 +0 53 -0 N N +0 0 53 1 53 +0 53 -0.5 53 +0.5 53 +0 53 -0 N N +# (x+0i)^(y-0i) = +a +0i when 0 < x < 1 and y >= +0 +0 0 53 1 53 +0 53 +0.5 53 +0 53 +0 53 -0 N N +0 0 53 0.25 53 +0 53 +0.5 53 +0 53 +2 53 -0 N N + +# (x-0i)^(y+0i) = +a +0i when x >= 1 and y <= -0 +0 0 53 0.5 53 +0 53 +2 53 -0 53 -1 53 +0 N N +0 0 53 1 53 +0 53 +3 53 -0 53 -0 53 +0 N N +0 0 53 1 53 +0 53 +1 53 -0 53 -1 53 +0 N N +0 0 53 1 53 +0 53 +1 53 -0 53 -0 53 +0 N N +# x^(-0+0i) = +1 +0i when |x| >= 1 and Im(x) < 0 +0 0 53 1 53 +0 53 +0.5 53 -0.5 53 -0 53 +0 N N +0 0 53 1 53 +0 53 +0 53 -0.5 53 -0 53 +0 N N +0 0 53 1 53 +0 53 -0 53 -0.5 53 -0 53 +0 N N +0 0 53 1 53 +0 53 -0.5 53 -0.5 53 -0 53 +0 N N +# x^(+0+0i) = +1 +0i when |x| >= 1 and Im(x) > 0 +0 0 53 1 53 +0 53 +0.5 53 +0.5 53 +0 53 +0 N N +0 0 53 1 53 +0 53 +0 53 -1 53 +0 53 +0 N N +0 0 53 1 53 +0 53 -0 53 +1 53 +0 53 +0 N N +0 0 53 1 53 +0 53 -0.5 53 +0.5 53 +0 53 +0 N N +# (x+0i)^(y+0i) = +a +0i when x >= 1 and y >= +0 +0 0 53 1 53 +0 53 +2 53 +0 53 +0 53 +0 N N +0 0 53 4 53 +0 53 +2 53 +0 53 +2 53 +0 N N +0 0 53 1 53 +0 53 +1 53 +0 53 +0 53 +0 N N +0 0 53 1 53 +0 53 +1 53 +0 53 +2 53 +0 N N + +# (+1-0i)^(y1+y2 i) = +1 +0i when y1 <= -0 and y2 > 0 +0 0 53 1 53 +0 53 +1 53 -0 53 -2 53 +1 N N +0 0 53 1 53 +0 53 +1 53 -0 53 -1 53 +2 N N +# (+/-1-0i)^(+0+y i) = +a +0i when y > 0 +0 0 53 1 53 +0 53 +1 53 -0 53 -0 53 +1 N N +- 0 53 +0x10BBEEE9177E19p-43 53 +0 53 -1 53 -0 53 -0 53 +2 N N +# x^(-0-y i) = +a +0i when |x| = 1, Im(x) < 0 and y > 0 ++ 0 53 +0x1724046EB0933Ap-48 53 +0 53 -1 53 -0 53 -0 53 +1 N N ++ 0 53 +0x1724046EB0933Ap-48 53 +0 53 -0 53 -1 53 -0 53 +2 N N +- 0 53 +0x1BD4567B975381p-46 53 +0 53 +0 53 -1 53 -0 53 +3 N N +0 0 53 1 53 +0 53 +1 53 -0 53 -0 53 +4 N N +# x^(+0+y i) = +a +0i when |x| = 1, Im(x) > 0 and y > 0 +0 0 53 1 53 +0 53 +1 53 +0 53 +0 53 +1 N N ++ 0 53 +0x1620227B598EF9p-57 53 +0 53 +0 53 +1 53 +0 53 +2 N N +- 0 53 +0x1265D4E92B6B9Bp-59 53 +0 53 -0 53 +1 53 +0 53 +3 N N ++ 0 53 +0x1D4102BC3F7D4Cp-71 53 +0 53 -1 53 +0 53 +0 53 +4 N N +# (+/-1+0i)^(0+yi) = +a +0i when y > 0 ++ 0 53 +0x1E989F5D6DFF5Cp-62 53 +0 53 -1 53 +0 53 +0 53 +2 N N +0 0 53 1 53 +0 53 +1 53 +0 53 +0 53 +2 N N +# (+1 +0i)^(y1+y2 i) = +1 +0i when y1 >= +0 and y2 >0 +0 0 53 1 53 +0 53 +1 53 +0 53 +2 53 +2 N N +0 0 53 1 53 +0 53 +1 53 +0 53 +0 53 +2 N N + +# (+1 +0i)^y = +1 -0i when Re(y) <= -0 and Im(y) < 0 +0 0 53 1 53 -0 53 +1 53 +0 53 -1 53 -1 N N +# (+/-1 +0i)^(-0 +yi) = +a -0i when y < 0 +0 0 53 1 53 -0 53 +1 53 +0 53 -0 53 -1 N N ++ 0 53 +0x1724046EB0933Ap-48 53 -0 53 -1 53 +0 53 -0 53 -1 N N +# x^(-0+y i) = +a -0i when |x| = 1, Im(x) >= +0 and y < 0 ++ 0 53 +0x1724046EB0933Ap-48 53 -0 53 +0 53 +1 53 -0 53 -2 N N +- 0 53 +0x1BD4567B975381p-46 53 -0 53 -0 53 +1 53 -0 53 -3 N N +# x^(+0+y i) = +a -0i when |x| = 1, Im(x) <= -0 and y < 0 +- 0 53 +0x1A9BCC46F767DFp-55 53 -0 53 +0 53 -1 53 +0 53 -1 N N ++ 0 53 +0x1620227B598EF9p-57 53 -0 53 -0 53 -1 53 +0 53 -2 N N +# (+/-1 -0i)^(+0+y i) = +a -0i when y < 0 +0 0 53 1 53 -0 53 +1 53 -0 53 +0 53 -1 N N ++ 0 53 +0x1620227B598EF9p-57 53 -0 53 -1 53 -0 53 +0 53 -1 N N +# (+1 -0i)^y = +1 -0i when Re(y) > 0 and Im(y) < 0 +0 0 53 1 53 -0 53 +1 53 -0 53 +2 53 -3 N N + +# (x +0i)^(y -0i) = +a -0i when x >= 1 and y < 0 +0 0 53 0.5 53 -0 53 +2 53 +0 53 -1 53 -0 N N +0 0 53 1 53 -0 53 +1 53 +0 53 -1 53 -0 N N +# x^(-0 -0i) = +1 -0i when |x| >= 1 and Im (x) >= +0 +0 0 53 1 53 -0 53 +1.5 53 +0 53 -0 53 -0 N N +0 0 53 1 53 -0 53 +1 53 +0 53 -0 53 -0 N N +0 0 53 1 53 -0 53 -1 53 +0 53 -0 53 -0 N N +0 0 53 1 53 -0 53 -1.5 53 +0 53 -0 53 -0 N N +0 0 53 1 53 -0 53 +1.5 53 +4 53 -0 53 -0 N N +0 0 53 1 53 -0 53 +1 53 +4 53 -0 53 -0 N N +0 0 53 1 53 -0 53 -1 53 +4 53 -0 53 -0 N N +0 0 53 1 53 -0 53 -1.5 53 +4 53 -0 53 -0 N N +# x^(+0 -0i) = +1 -0i when |x| >= 1 and Im (x) <= -0 +0 0 53 1 53 -0 53 +1.5 53 -0 53 +0 53 -0 N N +0 0 53 1 53 -0 53 +1 53 -0 53 +0 53 -0 N N +0 0 53 1 53 -0 53 -1 53 -0 53 +0 53 -0 N N +0 0 53 1 53 -0 53 -1.5 53 -0 53 +0 53 -0 N N +0 0 53 1 53 -0 53 +1.5 53 -4 53 +0 53 -0 N N +0 0 53 1 53 -0 53 +1 53 -4 53 +0 53 -0 N N +0 0 53 1 53 -0 53 -1 53 -4 53 +0 53 -0 N N +0 0 53 1 53 -0 53 -1.5 53 -4 53 +0 53 -0 N N +# (x -0i)^(y -0i) = x^y -0i when x >= 1 and y > 0 +0 0 53 9 53 -0 53 +3 53 -0 53 +2 53 -0 N N +0 0 53 1 53 -0 53 +1 53 -0 53 +2 53 -0 N N + +# (x +0i)^(y +0i) = x^y -0i when 0 < x < 1 and y < 0 +0 0 53 2 53 -0 53 +0.5 53 +0 53 -1 53 +0 N N +# x^(-0+0i) = +1 -0i when |x| < 1 and Im(x) >= +0 +0 0 53 1 53 -0 53 -0.5 53 +0 53 -0 53 +0 N N +0 0 53 1 53 -0 53 -0.1 53 +0.3 53 -0 53 +0 N N +0 0 53 1 53 -0 53 -0.0 53 +0.3 53 -0 53 +0 N N +0 0 53 1 53 -0 53 +0.0 53 +0.3 53 -0 53 +0 N N +0 0 53 1 53 -0 53 +0.1 53 +0.3 53 -0 53 +0 N N +0 0 53 1 53 -0 53 +0.5 53 +0 53 -0 53 +0 N N +# x^(+0+0i) = +1 -0i when |x| < 1 and Im(x) <= -0 +0 0 53 1 53 -0 53 -0.5 53 -0 53 +0 53 +0 N N +0 0 53 1 53 -0 53 -0.1 53 -0.3 53 +0 53 +0 N N +0 0 53 1 53 -0 53 -0.0 53 -0.3 53 +0 53 +0 N N +0 0 53 1 53 -0 53 +0.0 53 -0.3 53 +0 53 +0 N N +0 0 53 1 53 -0 53 +0.1 53 -0.3 53 +0 53 +0 N N +0 0 53 1 53 -0 53 +0.5 53 -0 53 +0 53 +0 N N +# (x-0i)^(y+0i) = x^y -0i when 0 < x < 1 and y > 0 +0 0 53 +0.25 53 -0 53 +0.5 53 -0 53 +2 53 +0 N N + # exact cases # (-4)^(1/4) = 1+i 0 0 2 1 2 1 2 -4 2 0 2 0x1p-2 2 0 N N @@ -195,14 +317,18 @@ # x = -2 - epsilon*i, y = -3 - delta*i # log(x) = log(2) - [Pi-epsilon/2]*i + O(epsilon^2) # y*log(x) ~ -3*log(2) + [3*Pi-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 -- - 53 -0x85649E3220691p-63 53 -0x14A25D455A9D0Dp-60 3 5 2 3 2 -3 2 +0 N N -+ 0 53 0xABCC77118461Dp-74 2 +0 3 5 3 5 2 -8 2 0 N N -+ 0 53 -0x127DB86014739Dp-93 2 +0 2 -1 2 -0 2 1 4 -9 N N -+ + 24 0xC1F98Dp-21 24 0x12FF89p-2 24 -7 24 +0 24 0xCFFFF3p-21 24 +0 N N +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 +# undefined zero sign in result +- 0 5 -25 2 0 2 +0 53 0xCCCCCCCCCCCCDp-54 2 -2 2 -0 N N + +- - 53 -0x85649E3220691p-63 53 -0x14A25D455A9D0Dp-60 3 5 2 3 2 -3 2 +0 N N ++ 0 53 0xABCC77118461Dp-74 2 +0 3 5 3 5 2 -8 2 +0 N N + ++ 0 53 -0x127DB86014739Dp-93 2 -0 2 -1 2 -0 2 1 4 -9 N N ++ + 24 0xC1F98Dp-21 24 0x12FF89p-2 24 -7 24 +0 24 0xCFFFF3p-21 24 +0 N N # underflow case - - 24 +0 24 +0 24 2 24 0x44CCCDp-20 24 -0x7FFFF200 24 -0x7FFFF200 N N - 0 53 0x14D55AFA6E0BB0p210433620 53 0 53 +0 53 0x44CCCCFFFFFFFp-48 53 0x5F5E100 53 +0 N N |