improve asin for tiny inputs (not only with Re(op)=1)

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) */