From e39ce820f272a35c183b3921372d826c1251ee51 Mon Sep 17 00:00:00 2001 From: Paul Zimmermann Date: Fri, 14 Feb 2020 17:29:08 +0100 Subject: improve asin for tiny inputs (not only with Re(op)=1) --- doc/algorithms.tex | 43 +++++++++++++++++++++++++++++++++++-- src/asin.c | 62 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 103 insertions(+), 2 deletions(-) diff --git a/doc/algorithms.tex b/doc/algorithms.tex index 8005417..df858ad 100644 --- a/doc/algorithms.tex +++ b/doc/algorithms.tex @@ -1186,9 +1186,48 @@ k_R=\left\{ \subsection {\texttt {mpc\_asin}} -For $0 \leq y \leq 1$, we have $|\Re \asin (1 \pm iy) - \pi/2| +\paragraph{Tiny imaginary part.} +First, in case the real part of the argument is exactly $1$, +for $0 \leq y \leq 1$, we have $|\Re \asin (1 \pm iy) - \pi/2| \leq y^{1/2}$, and $|\Im \asin (1 \pm iy) \mp y^{1/2}| \leq (1/12) y^{3/2}$. +In the general case, +formula 4.4.40 from \cite{AbSt73} gives for $|z| < 1$: +\begin{equation} \label{eq_asin} +\asin z = z + \frac{1}{2} \frac{z^3}{3} + \frac{1 \cdot 3}{2 \cdot 4} + \frac{z^5}{5} + \frac{1 \cdot 3 \cdot 5}{2 \cdot 4 \cdot 6} \frac{z^7}{7} + + \cdots +\end{equation} +thus we have $\asin z = z + z^3 t$ where +\[ t = \frac{1}{2} \frac{1}{3} + \frac{1 \cdot 3}{2 \cdot 4} + \frac{z^2}{5} + \frac{1 \cdot 3 \cdot 5}{2 \cdot 4 \cdot 6} \frac{z^4}{7} + + \cdots, \] +thus +\[ |t| \leq \frac{1}{2} \left( \frac{1}{3} + \frac{|z|^2}{5} + + \frac{|z|^7}{7} + \cdots \right) \] +% var('k'); 1/2*sum((1/2)^(2*k)/(2*k+3),k,0,infinity) +We assume $|z| \leq 1/2$. +Since $\sum_{k \geq 0} (1/2)^{2k}/(2k+3) += 2 \log 3 - 2 \leq 0.198$, thus $|t| < 0.198$, +we deduce for the real part: +\[ |\Re\asin z - \Re z| < 0.198 |z|^3. \] + +For the imaginary part, +we first prove by induction that $|\Im z^k| \leq k |z|^{k-1} |\Im z|$ +for $k \geq 1$. This is trivial for $k=1$. +Now assume it is true for $k$: +\[ |\Im z^{k+1}| \leq |\Im z^k| |\Re z| + |\Im z| |\Re z^k| + \leq k |z|^{k-1} |\Im z| |z| + |\Im z| |z|^k + \leq (k+1) |z|^k |\Im z|. \] +Applying this to each term of Eq.~(\ref{eq_asin}) we get: +\[ |\Im\asin z - \Im z| \leq \frac{1}{2} |z^2| |\Im z| + + \frac{1 \cdot 3}{2 \cdot 4} |z^4| |\Im z| + + \frac{1 \cdot 3 \cdot 5}{2 \cdot 4 \cdot 6} |z^6| |\Im z| + \cdots \] +Still for $|z| \leq 1/2$, we get: +\[ |\Im\asin z - \Im z| \leq \frac{1}{2} |z^2| |\Im z| + + \frac{1}{2} |z^4| |\Im z| + \frac{1}{2} |z^6| |\Im z| + \cdots + \leq \frac{2}{3} |z^2| |\Im z|. \] + \subsection {\texttt {mpc\_pow}} The main issue for the power function is to be able to recognize when the @@ -2079,7 +2118,7 @@ bits~$0$, have the following shape: If $y_k' > \sqrt 2 \cdot 2^{p-1}$, they encode exactly $(y_k')^2$. If $y_k' < \sqrt 2 \cdot 2^{p-1}$, the $2 p - 1$ leading bits encode $(y_k')^2$, and they are followed by a digit~$0$, so the $2 p$ leading -bits encode $2 (y_k')^2$. +bits encode $2 (y_k')^2$. Since $n_1 = n_2$, in particular their leading $2 p$ bits coincide; so either $(y_1')^2 = (y_2')^2$, $(y_1')^2 = 2 (y_2')^2$ or diff --git a/src/asin.c b/src/asin.c index 4cf2353..e69e2b0 100644 --- a/src/asin.c +++ b/src/asin.c @@ -80,6 +80,64 @@ mpc_asin_special (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, mpc_ptr z1) return 1; } +/* Special case for |z| <= 1/2 and tiny imaginary part (see algorithms.tex). + Return 0 if the special case fails (this cannot be confused with an exact + ternary value since this cannot happen). + Otherwise puts an approximation in z1, so that rounding it in rop will get + the correct result. */ +static int +mpc_asin_special2 (mpc_ptr rop, mpc_srcptr op, mpc_ptr z1) +{ + /* for |z| <= 1/2, we know that: + |Re(asin z) - Re(z)| < 0.198 |z|^3 + |Im(asin z) - Im(z)| < 0.667 |z|^2 |Im(z)| + thus if 0.667 |z|^2 < 1/2 ulp(1) we can round to Re(z) and Im(z) */ + mpfr_exp_t ex = mpfr_get_exp (mpc_realref (op)); + mpfr_exp_t ey = mpfr_get_exp (mpc_imagref (op)); + mpfr_prec_t px = mpfr_get_prec (mpc_realref (rop)); + mpfr_prec_t py = mpfr_get_prec (mpc_imagref (rop)); + int sign_re, sign_im, sign_re_t, sign_im_t; + mpc_t t; + + /* |x| < 2^ex and |y| < 2^ey thus |z^2| < 2^(2*max(ex,ey)+1) */ + ex = (ex > ey) ? 2 * ex + 1 : 2 * ey + 1; + /* 0.667 |z|^2 < 1/2 ulp(1) if 2^ex <= 2^-p thus ex <= -p */ + if (ex + px > 0 || ex + py > 0) + return 0; /* special case fails */ + + /* if the precision of the input is larger than that of the output, + we fail since we should deal with a double-rounding, which is not + implemented below */ + if (mpfr_get_prec (mpc_realref (op)) > px || + mpfr_get_prec (mpc_imagref (op)) > py) + return 0; + + mpc_init2 (t, 32); + mpc_pow_ui (t, op, 3, MPC_RNDNN); + sign_re_t = MPFR_SIGN (mpc_realref (t)); + sign_im_t = MPFR_SIGN (mpc_imagref (t)); + mpc_clear (t); + + /* approximate the real part: since asin z = z + z^3/6 + ..., the real part + has to be rounded according to the sign of Re(t) */ + mpfr_set (mpc_realref (z1), mpc_realref (op), MPFR_RNDN); + sign_re = MPFR_SIGN (mpc_realref (op)); + if (sign_re * sign_re_t > 0) + MPFR_ADD_ONE_ULP (mpc_realref (z1)); + else + MPFR_SUB_ONE_ULP (mpc_realref (z1)); + + /* approximate the imaginary part */ + mpfr_set (mpc_imagref (z1), mpc_imagref (op), MPFR_RNDN); + sign_im = MPFR_SIGN (mpc_imagref (op)); + if (sign_im * sign_im_t > 0) + MPFR_ADD_ONE_ULP (mpc_imagref (z1)); + else + MPFR_SUB_ONE_ULP (mpc_imagref (z1)); + + return 1; +} + int mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { @@ -223,6 +281,10 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpc_asin_special (rop, op, rnd, z1)) break; + /* try special code for 1+i*y with tiny y */ + if (loop == 1 && mpc_asin_special2 (rop, op, z1)) + break; + /* z1 <- z^2 */ mpc_sqr (z1, op, MPC_RNDNN); /* err(x) <= 1/2 ulp(x), err(y) <= 1/2 ulp(y) */ -- cgit v1.2.1