improve speed of asin for 1+i*y with tiny y

git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@1402 211d60ee-9f03-0410-a15a-8952a2c7a4e4

3 files changed, 74 insertions, 1 deletions

diff --git a/doc/algorithms.tex b/doc/algorithms.tex
index 669d651..d9acd15 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -17,6 +17,7 @@
 \newcommand {\Ulp}{{\operatorname {ulp}}}
 \DeclareMathOperator{\Exp}{\operatorname {Exp}}
 \newcommand {\atantwo}{\operatorname {atan2}}
+\newcommand {\asin}{\operatorname {asin}}
 \newcommand{\error}{\operatorname {error}}
 \newcommand{\relerror}{\operatorname {relerror}}
 \newcommand{\Norm}{\operatorname {N}}
@@ -1142,6 +1143,11 @@ k_R=\left\{
 \right.
 \end{equation*}
 
+\subsection {\texttt {mpc\_asin}}
+
+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}$.
+
 \subsection {\texttt {mpc\_pow}}
 
 The main issue for the power function is to be able to recognize when the
diff --git a/src/asin.c b/src/asin.c
index e796461..a92421a 100644
--- a/src/asin.c
+++ b/src/asin.c
@@ -1,6 +1,6 @@
 /* mpc_asin -- arcsine of a complex number.
 
-Copyright (C) 2009, 2010, 2011, 2012 INRIA
+Copyright (C) 2009, 2010, 2011, 2012, 2013 INRIA
 
 This file is part of GNU MPC.
 
@@ -20,6 +20,67 @@ along with this program. If not, see http://www.gnu.org/licenses/ .
 
 #include "mpc-impl.h"
 
+/* Special case op = 1 + i*y for tiny y (see algorithms.tex).
+   Return 0 if special formula fails, otherwise put in z1 the approximate
+   value which needs to be converted to rop.
+   z1 is a temporary variable with enough precision.
+ */
+static int
+mpc_asin_special (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, mpc_ptr z1)
+{
+  mpfr_exp_t ey = mpfr_get_exp (mpc_imagref (op));
+  mpfr_t abs_y;
+  mpfr_prec_t p;
+  int inex;
+
+  /* |Re(asin(1+i*y)) - pi/2| <= y^(1/2) */
+  if (ey >= 0 || ((-ey) / 2 < mpfr_get_prec (mpc_realref (z1))))
+    return 0;
+  if (mpfr_get_prec (mpc_imagref (z1)) > (-ey) + 2) /* p > -ey + 2 */
+    return 0;
+
+  mpfr_const_pi (mpc_realref (z1), MPFR_RNDN);
+  mpfr_div_2exp (mpc_realref (z1), mpc_realref (z1), 1, MPFR_RNDN); /* exact */
+  p = mpfr_get_prec (mpc_realref (z1));
+  /* if z1 has precision p, the error on z1 is 1/2*ulp(z1) = 2^(-p) so far,
+     and since ey <= -2p, then y^(1/2) <= 1/2*ulp(z1) too, thus the total
+     error is bounded by ulp(z1) */
+  if (!mpfr_can_round (mpc_realref(z1), p, MPFR_RNDN, MPFR_RNDZ,
+                       mpfr_get_prec (mpc_realref(rop)) +
+                       (MPC_RND_RE(rnd) == MPFR_RNDN)))
+    return 0;
+  
+  /* |Im(asin(1+i*y)) - y^(1/2)| <= (1/12) * y^(3/2) for y >= 0 (err >= 0)
+     |Im(asin(1-i*y)) + y^(1/2)| <= (1/12) * y^(3/2) for y >= 0 (err <= 0) */
+  abs_y[0] = mpc_imagref (op)[0];
+  if (mpfr_signbit (mpc_imagref (op)))
+    MPFR_CHANGE_SIGN (abs_y);
+  inex = mpfr_sqrt (mpc_imagref (z1), abs_y, MPFR_RNDN); /* error <= 1/2 ulp */
+  if (mpfr_signbit (mpc_imagref (op)))
+    MPFR_CHANGE_SIGN (mpc_imagref (z1));
+  /* If z1 has precision p, the error on z1 is 1/2*ulp(z1) = 2^(-p) so far,
+     and (1/12) * y^(3/2) <= (1/8) * y * y^(1/2) <= 2^(ey-3)*2^p*ulp(y^(1/2))
+     thus for p+ey-3 <= -1 we have (1/12) * y^(3/2) <= (1/2) * ulp(y^(1/2)),
+     and the total error is bounded by ulp(z1).
+     Note: if y^(1/2) is exactly representable, and ends with many zeroes,
+     then mpfr_can_round below will fail; however in that case we know that
+     Im(asin(1+i*y)) is away from +/-y^(1/2) wrt zero. */
+  if (inex == 0) /* enlarge im(z1) so that the inexact flag is correct */
+    {
+      if (mpfr_signbit (mpc_imagref (op)))
+        mpfr_nextbelow (mpc_imagref (z1));
+      else
+        mpfr_nextabove (mpc_imagref (z1));
+      return 1;
+    }
+  p = mpfr_get_prec (mpc_imagref (z1));
+  if (!mpfr_can_round (mpc_imagref(z1), p - 1, MPFR_RNDA, MPFR_RNDZ,
+                      mpfr_get_prec (mpc_imagref(rop)) +
+                      (MPC_RND_IM(rnd) == MPFR_RNDN)))
+    return 0;
+  return 1;
+}
+
 int
 mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
 {
@@ -153,6 +214,10 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
     mpfr_set_prec (mpc_realref(z1), p);
     mpfr_set_prec (mpc_imagref(z1), p);
 
+    /* try special code for 1+i*y with tiny y */
+    if (loop == 1 && mpc_asin_special (rop, op, rnd, z1))
+      break;
+
     /* z1 <- z^2 */
     mpc_sqr (z1, op, MPC_RNDNN);
     /* err(x) <= 1/2 ulp(x), err(y) <= 1/2 ulp(y) */
diff --git a/tests/asin.dat b/tests/asin.dat
index 9793fcc..e93e8c5 100644
--- a/tests/asin.dat
+++ b/tests/asin.dat
@@ -124,3 +124,5 @@
 + - 53 0x189BF9EC7FCD5Bp-54 53 0x1206ECFA94614Bp-50 53 17 53 42 N N
 - + 2 1.5 2 6 2 96 2 0x1p-8 N N
 - - 8 0xC9p-7 8 0x15p-2 2 96 2 0x1p-8 N N
+- - 53 0x3243f6a8885a3p-49 53 0x1p-500 53 1 53 0x1p-1000 N N
+- + 53 0x3243f6a8885a3p-49 53 -0x1p-500 53 1 53 -0x1p-1000 N N