improved asin (try series for |z| < 1, and special case with tiny Im(z))

1 files changed, 180 insertions, 52 deletions

diff --git a/src/asin.c b/src/asin.c
index e69e2b0..9b09634 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, 2013, 2014 INRIA
+Copyright (C) 2009-2020 INRIA
 
 This file is part of GNU MPC.
 
@@ -80,62 +80,182 @@ 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. */
+/* Put in s an approximation of asin(z) using:
+   asin z = z + 1/2*z^3/3 + (1*3)/(2*4)*z^5/5 + ...
+   Assume |Re(z)|, |Im(z)| < 1/2.
+   Return non-zero if we can get the correct result by rounding s:
+   mpc_set (rop, s, ...) */
 static int
-mpc_asin_special2 (mpc_ptr rop, mpc_srcptr op, mpc_ptr z1)
+mpc_asin_series (mpc_srcptr rop, mpc_ptr s, mpc_srcptr z, mpc_rnd_t rnd)
 {
-  /* 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)
+  mpc_t w, t;
+  unsigned long k, err, kx, ky;
+  mpfr_prec_t p;
+  mpfr_exp_t ex, ey, e;
+
+  /* assume z = (x,y) with |x|,|y| < 2^(-e) with e >= 1, see the error
+     analysis in algorithms.tex */
+  ex = mpfr_get_exp (mpc_realref (z));
+  MPC_ASSERT(ex <= -1);
+  ey = mpfr_get_exp (mpc_imagref (z));
+  MPC_ASSERT(ey <= -1);
+  e = (ex >= ey) ? ex : ey;
+  e = -e;
+  /* now e >= 1 */
+
+  p = mpfr_get_prec (mpc_realref (s)); /* working precision */
+  MPC_ASSERT(mpfr_get_prec (mpc_imagref (s)) == p);
+
+  mpc_init2 (w, p);
+  mpc_init2 (t, p);
+  mpc_set (s, z, MPC_RNDNN);
+  mpc_sqr (w, z, MPC_RNDNN);
+  mpc_set (t, z, MPC_RNDNN);
+  for (k = 1; ;k++)
+    {
+      mpfr_exp_t exp_x, exp_y;
+      mpc_mul (t, t, w, MPC_RNDNN);
+      mpc_mul_ui (t, t, (2 * k - 1) * (2 * k - 1), MPC_RNDNN);
+      mpc_div_ui (t, t, (2 * k) * (2 * k + 1), MPC_RNDNN);
+      exp_x = mpfr_get_exp (mpc_realref (s));
+      exp_y = mpfr_get_exp (mpc_imagref (s));
+      if (mpfr_get_exp (mpc_realref (t)) < exp_x - p &&
+          mpfr_get_exp (mpc_imagref (t)) < exp_y - p)
+        /* Re(t) < 1/2 ulp(Re(s)) and Im(t) < 1/2 ulp(Im(s)),
+           thus adding t to s will not change s */
+        break;
+      mpc_add (s, s, t, MPC_RNDNN);
+    }
+  mpc_clear (w);
+  mpc_clear (t);
+  /* check (2k-1)^2 is exactly representable */
+  MPC_ASSERT(2 * k - 1 <= ULONG_MAX / (2 * k - 1));
+  /* maximal absolute error on Re(s),Im(s) is:
+     (5k-3)k/2*2^(-1-p) for e=1
+     5k/2*2^(-e-p) for e >= 2 */
+  if (e == 1)
+    {
+      MPC_ASSERT(5 * k - 3 <= ULONG_MAX / k);
+      kx = (5 * k - 3) * k;
+    }
+  else
+    kx = 5 * k;
+  kx = (kx + 1) / 2; /* takes into account the 1/2 factor in both cases */
+  /* now (5k-3)k/2 <= kx for e=1, and 5k/2 <= kx for e >= 2, thus
+     the maximal absolute error on Re(s),Im(s) is bounded by kx*2^(-e-p) */
+
+  e = -e;
+  ky = kx;
+
+  /* for the real part, convert the maximal absolute error kx*2^(e-p) into
+     relative error */
+  ex = mpfr_get_exp (mpc_realref (s));
+  /* ulp(Re(s)) = 2^(ex+1-p) */
+  if (ex+1 > e) /* divide kx by 2^(ex+1-e) */
+    while (ex+1 > e)
+      {
+        kx = (kx + 1) / 2;
+        ex --;
+      }
+  else /* multiply kx by 2^(e-(ex+1)) */
+    kx <<= e - (ex+1);
+  /* now the rounding error is bounded by kx*ulp(Re(s)), add the
+     mathematical error which is bounded by ulp(Re(s)): the first neglected
+     term is less than 1/2*ulp(Re(s)), and each term decreases by at least
+     a factor 2, since |z^2| <= 1/2. */
+  kx ++;
+  for (err = 0; kx > 2; err ++, kx = (kx + 1) / 2);
+  /* can we round Re(s) with error less than 2^(EXP(Re(s))-err) ? */
+  if (!mpfr_can_round (mpc_realref (s), p - err, MPFR_RNDN, MPFR_RNDZ,
+                       mpfr_get_prec (mpc_realref (rop)) +
+                       (MPC_RND_RE(rnd) == MPFR_RNDN)))
     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);
+  /* same for the imaginary part */
+  ey = mpfr_get_exp (mpc_imagref (s));
+  /* we take for e the exponent of Im(z), which amounts to divide the error by
+     2^delta where delta is the exponent difference between Re(z) and Im(z)
+     (see algorithms.tex) */
+  e = mpfr_get_exp (mpc_imagref (z));
+  /* ulp(Im(s)) = 2^(ey+1-p) */
+  if (ey+1 > e) /* divide ky by 2^(ey+1-e) */
+    while (ey+1 > e)
+      {
+        ky = (ky + 1) / 2;
+        ey --;
+      }
+  else /* multiply ky by 2^(e-(ey+1)) */
+    ky <<= e - (ey+1);
+  /* now the rounding error is bounded by ky*ulp(Im(s)), add the
+     mathematical error which is bounded by ulp(Im(s)): the first neglected
+     term is less than 1/2*ulp(Im(s)), and each term decreases by at least
+     a factor 2, since |z^2| <= 1/2. */
+  ky ++;
+  for (err = 0; ky > 2; err ++, ky = (ky + 1) / 2);
+  /* can we round Im(s) with error less than 2^(EXP(Im(s))-err) ? */
+  return mpfr_can_round (mpc_imagref (s), p - err, MPFR_RNDN, MPFR_RNDZ,
+                         mpfr_get_prec (mpc_imagref (rop)) +
+                         (MPC_RND_IM(rnd) == MPFR_RNDN));
+}
 
-  /* 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));
+/* Put in s an approximation of asin(z) for |Re(z)| < 1 and tiny Im(y)
+   (see algorithms.tex). Assume z = x + i*y:
+   |Re(asin(z)) - asin(x)| <= Pi*beta^2/(1-beta)
+   |Im(asin(z)) - y/sqrt(1-x^2)| <= Pi/2*beta^2/(1-beta)
+   where beta = |y|/(1-|x|).
+   We assume |x| >= 1/2 > |y| here, thus beta < 1, and the bounds
+   simplify to 16y^2.
+   Assume Re(s) and Im(s) have the same precision.
+   Return non-zero if we can get the correct result by rounding s:
+   mpc_set (rop, s, ...) */
+static int
+mpc_asin_tiny (mpc_srcptr rop, mpc_ptr s, mpc_srcptr z, mpc_rnd_t rnd)
+{
+  mpfr_exp_t ey, e1, e2, es;
+  mpfr_prec_t p = mpfr_get_prec (mpc_realref (s)), err;
 
-  /* 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));
+  ey = mpfr_get_exp (mpc_imagref (z));
 
-  return 1;
+  MPC_ASSERT(mpfr_get_exp (mpc_realref (z)) >= 0); /* |x| >= 1/2 */
+  MPC_ASSERT(ey < 0);  /* |y| < 1/2 */
+
+  e2 = 2 * ey + 4;
+  /* for |x| >= 0.5, |asin x| >= 0.5, thus no need to compute asin(x)
+     if 16y^2 >= 1/2 ulp(1) for the target variable */
+  if (e2 >= - (mpfr_exp_t) mpfr_get_prec (mpc_realref (rop)))
+    return 0;
+
+  /* real part */
+  mpfr_asin (mpc_realref (s), mpc_realref (z), MPFR_RNDN);
+  /* check that we can round with error < 16y^2 */
+  e1 = mpfr_get_exp (mpc_realref (s)) - p;
+  err = (e1 >= e2) ? 0 : e2 - e1;
+  if (!mpfr_can_round (mpc_realref (s), p - err, MPFR_RNDN, MPFR_RNDZ,
+                       mpfr_get_prec (mpc_realref (rop)) +
+                       (MPC_RND_RE(rnd) == MPFR_RNDN)))
+    return 0;
+
+  /* now compute the approximate imaginary part y/sqrt(1-x^2) */
+  mpfr_sqr (mpc_imagref (s), mpc_realref (z), MPFR_RNDN);
+  /* now Im(s) approximates x^2, with 1/4 <= Im(s) <= 1 and absolute error
+     less than 2^-p */
+  mpfr_ui_sub (mpc_imagref (s), 1, mpc_imagref (s), MPFR_RNDN);
+  /* now Im(s) approximates 1-x^2, with 0 <= Im(s) <= 3/4, assuming p >= 2,
+     and absolute error less than 2^(1-p) */
+  mpfr_sqrt (mpc_imagref (s), mpc_imagref (s), MPFR_RNDN); /* sqrt(1-x^2) */
+  es = mpfr_get_exp (mpc_imagref (s));
+  /* now Im(s) approximates sqrt(1-x^2), with 0 <= Im(s) <= 7/8,
+     assuming p >= 3, and absolute error less than 2^-p + 2^(1-p)/s
+     <= 3*2^-p/s, thus relative error less than 3*2^-p/s^2.
+     Since |s| >= 2^(es-1), the relative error is less than 3*2^(-p+2-2*es) */
+  mpfr_div (mpc_imagref (s), mpc_imagref (z), mpc_imagref (s), MPFR_RNDN);
+  /* now Im(s) approximates y/sqrt(1-x^2), with relative error less than
+     (3*2^(2-2*es)+2)*2^-p. Since es <= 0, 2-2*es >= 2, thus the relative
+     error is less than 2^(4-2*es-p) */
+  err = 4 - 2 * es;
+  return mpfr_can_round (mpc_imagref (s), p - err, MPFR_RNDN, MPFR_RNDZ,
+                         mpfr_get_prec (mpc_imagref (rop)) +
+                         (MPC_RND_IM(rnd) == MPFR_RNDN));
 }
 
 int
@@ -281,8 +401,16 @@ 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))
+    /* try special code for small z */
+    if (mpfr_get_exp (mpc_realref (op)) <= -1 &&
+        mpfr_get_exp (mpc_imagref (op)) <= -1 &&
+        mpc_asin_series (rop, z1, op, rnd))
+      break;
+
+    /* try special code for 1/2 <= |x| < 1 and |y| < 1/2 */
+    if (mpfr_get_exp (mpc_realref (op)) == 0 &&
+        mpfr_get_exp (mpc_imagref (op)) <= -1 &&
+        mpc_asin_tiny (rop, z1, op, rnd))
       break;
 
     /* z1 <- z^2 */