atan.c: transform test into assertion, add comments

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

3 files changed, 13 insertions, 2 deletions

diff --git a/NEWS b/NEWS
index 73c69eb..5127113 100644
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,7 @@
 Changes in version 1.0:
   - First release as a GNU package
   - License change: LGPLv3+ for code, GFDLv1.3+ for documentation
+  - 100% of all lines are covered by tests
   - Functions renamed:
     mpc_mul_2exp to mpc_mul_2ui, mpc_div_2exp to mpc_div_2ui
   - 0^0, which returned (NaN,NaN) previously, now returns (1,+0).
diff --git a/TODO b/TODO
index 8858a50..df70e2e 100644
--- a/TODO
+++ b/TODO
@@ -12,6 +12,10 @@ mpc-impl.h.
 From Andreas Enge 05 July 2012:
 Add support for rounding mode MPFR_RNDA.
 
+From Andreas Enge and Paul Zimmermann 6 July 2012:
+Improve speed of Im (atan) for x+i*y with small y, for instance by using
+the Taylor series directly.
+
 Bench:
 - from Andreas Enge 9 June 2009:
   Scripts and web page comparing timings with different systems,
diff --git a/src/atan.c b/src/atan.c
index 8aa8553..c0b01a4 100644
--- a/src/atan.c
+++ b/src/atan.c
@@ -332,6 +332,8 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
         mpfr_sub (y, a, b, GMP_RNDU);
 
         if (mpfr_zero_p (y))
+           /* FIXME: happens when x and y have very different magnitudes;
+              could be handled more efficiently                           */
           ok = 0;
         else
           {
@@ -345,8 +347,12 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
               err = (expo < 0) ? 1 : expo + 2;
 
             mpfr_div_2ui (y, y, 2, GMP_RNDN);
-            if (mpfr_zero_p (y))
-              err = 2; /* underflow */
+            MPC_ASSERT (!mpfr_zero_p (y));
+               /* FIXME: underflow. Since the main term of the Taylor series
+                  in y=0 is 1/(x^2+1) * y, this means that y is very small
+                  and/or x very large; but then the mpfr_zero_p (y) above
+                  should be true. This needs a proof, or better yet,
+                  special code.                                              */
 
             ok = mpfr_can_round (y, p - err, GMP_RNDU, GMP_RNDD,
                                  prec + (MPC_RND_IM (rnd) == GMP_RNDN));