Modified abs to be more efficient when real and imaginary part have very

different orders of magnitude; detect when the absolute value is
essentially one of the parts, and detect the cases where a Taylor
expansion of order 1 of the square root is sufficient.
Added test code for the absolute value.
Changed version and modification date after release 0.4.4.
(AE)


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

5 files changed, 313 insertions, 12 deletions

diff --git a/abs.c b/abs.c
index 74e91db..beb1c78 100644
--- a/abs.c
+++ b/abs.c
@@ -1,6 +1,6 @@
 /* mpc_abs -- Absolute value of a complex number.
 
-Copyright (C) 2002, 2004 Andreas Enge, Paul Zimmermann
+Copyright (C) 2002, 2004, 2005 Andreas Enge, Paul Zimmermann
 
 This file is part of the MPC Library.
 
@@ -19,13 +19,14 @@ along with the MPC Library; see the file COPYING.LIB.  If not, write to
 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */
 
+#include <stdio.h>
 #include "gmp.h"
 #include "mpfr.h"
 #include "mpc.h"
 #include "mpc-impl.h"
 
 int
-mpc_abs (mpfr_ptr a, mpc_srcptr b, mp_rnd_t rnd)
+mpc_abs_basic (mpfr_ptr a, mpc_srcptr b, mp_rnd_t rnd)
 {
   mpfr_t u, v;
   mp_prec_t prec;
@@ -44,14 +45,131 @@ mpc_abs (mpfr_ptr a, mpc_srcptr b, mp_rnd_t rnd)
       mpfr_set_prec (v, prec);
 
       /* first compute norm(b)^2 */
-      inexact = mpfr_mul (u, MPC_RE(b), MPC_RE(b), GMP_RNDD); /* err<=1ulp */
-      inexact |= mpfr_mul (v, MPC_IM(b), MPC_IM(b), GMP_RNDD); /* err<=1ulp*/
+      inexact = mpfr_mul (u, MPC_RE(b), MPC_RE(b), GMP_RNDU); /* err<=1ulp */
+      inexact |= mpfr_mul (v, MPC_IM(b), MPC_IM(b), GMP_RNDU); /* err<=1ulp*/
 
-      inexact |= mpfr_add (u, u, v, GMP_RNDD);            /* err <= 3 ulps */
-      inexact |= mpfr_sqrt (u, u, GMP_RNDD); /* err <= 4 ulps */
+      inexact |= mpfr_add (u, u, v, GMP_RNDU);            /* err <= 3 ulps */
+      inexact |= mpfr_sqrt (u, u, GMP_RNDU); /* err <= 4 ulps */
     }
   while (inexact != 0 &&
-         mpfr_can_round (u, prec - 2, GMP_RNDD, rnd, MPFR_PREC(a)) == 0);
+         mpfr_can_round (u, prec - 2, GMP_RNDU, rnd, MPFR_PREC(a)) == 0);
+
+  inexact |= mpfr_set (a, u, rnd);
+
+  mpfr_clear (u);
+  mpfr_clear (v);
+
+  return inexact;
+}
+
+
+int
+mpc_abs (mpfr_ptr a, mpc_srcptr b, mp_rnd_t rnd)
+   /* a more sophisticated version treating special cases separately */
+{
+  mpfr_t u, v;
+  mp_prec_t prec;
+  int inexact;
+
+  if (MPFR_IS_ZERO (MPC_RE (b)))
+     return mpfr_abs (a, MPC_IM (b), rnd);
+  else if (MPFR_IS_ZERO (MPC_IM (b)))
+     return mpfr_abs (a, MPC_RE (b), rnd);
+  
+  prec = MPFR_PREC(a);
+  if (prec >= MPFR_PREC (MPC_RE (b))
+        && 2 * (MPFR_EXP (MPC_RE (b)) - MPFR_EXP (MPC_IM (b)) - 1) >
+           (signed long) prec)
+     /* try to detect the case that the imaginary part is negligible;
+        this is the case, independently of the rounding mode, whenever
+        setting the result to the real part is exact, and the imaginary
+        part contributes less than 1/2 ulp                              */
+  {
+     mpfr_abs (a, MPC_RE (b), GMP_RNDN); /* exact */
+     if (rnd == GMP_RNDU)
+     {
+        mpfr_add_one_ulp (a, GMP_RNDN);
+        return 1;
+     }
+     else
+        return -1;
+  }
+  else if (prec >= MPFR_PREC (MPC_IM (b))
+        && 2 * (MPFR_EXP (MPC_IM (b)) - MPFR_EXP (MPC_RE (b)) - 1) >
+           (signed long) prec)
+  {
+     mpfr_abs (a, MPC_IM (b), GMP_RNDN); /* exact */
+     if (rnd == GMP_RNDU)
+     {
+        mpfr_add_one_ulp (a, GMP_RNDN);
+        return 1;
+     }
+     else
+        return -1;
+  }
+  
+  mpfr_init2 (u, 2*prec);
+  mpfr_init2 (v, 2*prec);
+  
+  do
+    {
+      prec += mpc_ceil_log2 (prec) + 3;
+
+      mpfr_set_prec (u, prec);
+      mpfr_set_prec (v, prec);
+
+      if (4*(MPFR_EXP(MPC_RE(b)) - MPFR_EXP(MPC_IM(b))) > (signed long) prec)
+         /* |Re (b)| >> |Im (b)|, use the Taylor series sqrt (a^2 + b^2)
+            = |a| + b^2 / (2 |a|) - e with 0 <= e <= 1 ulp of the target;
+            so we round towards infinity                                 */
+      {
+         inexact = 1;
+         mpfr_mul (u, MPC_IM (b), MPC_IM (b), GMP_RNDU); /* err <= 1 ulp */
+         if (MPFR_SIGN (MPC_RE (b)) > 0)
+         {
+            mpfr_div (u, u, MPC_RE (b), GMP_RNDU); /* err <= 3 ulp */
+            mpfr_div_2ui (u, u, 1, GMP_RNDU);
+            mpfr_add (u, u, MPC_RE (b), GMP_RNDU); /* err <= 4 ulp */
+         }
+         else
+         {
+            mpfr_div (u, u, MPC_RE (b), GMP_RNDD);
+            mpfr_div_2ui (u, u, 1, GMP_RNDD);
+            mpfr_add (u, u, MPC_RE (b), GMP_RNDD);
+            mpfr_neg (u, u, GMP_RNDU);
+         }
+         /* err <= 5 ulps (one for Taylor approximation) */
+      }
+      else if (4*(MPFR_EXP(MPC_IM(b)) - MPFR_EXP(MPC_RE(b))) >
+               (signed long) prec)
+      {
+         inexact = 1;
+         mpfr_mul (u, MPC_RE (b), MPC_RE (b), GMP_RNDU);
+         if (MPFR_SIGN (MPC_IM (b)) > 0)
+         {
+            mpfr_div (u, u, MPC_IM (b), GMP_RNDU);
+            mpfr_div_2ui (u, u, 1, GMP_RNDU);
+            mpfr_add (u, u, MPC_IM (b), GMP_RNDU);
+         }
+         else
+         {
+            mpfr_div (u, u, MPC_IM (b), GMP_RNDD);
+            mpfr_div_2ui (u, u, 1, GMP_RNDD);
+            mpfr_add (u, u, MPC_IM (b), GMP_RNDD);
+            mpfr_neg (u, u, GMP_RNDU);
+         }
+      }
+      else
+      {
+         /* first compute norm(b)^2 */
+         inexact = mpfr_mul (u, MPC_RE(b), MPC_RE(b), GMP_RNDU); /* err<=1ulp */
+         inexact |= mpfr_mul (v, MPC_IM(b), MPC_IM(b), GMP_RNDU); /* err<=1ulp*/
+         inexact |= mpfr_add (u, u, v, GMP_RNDU);            /* err <= 3 ulps */
+         inexact |= mpfr_sqrt (u, u, GMP_RNDU); /* err <= 4 ulps */
+      }
+   }
+  while (inexact != 0 &&
+         mpfr_can_round (u, prec - 3, GMP_RNDU, rnd, MPFR_PREC(a)) == 0);
 
   inexact |= mpfr_set (a, u, rnd);
 
diff --git a/makefile b/makefile
index ffd4bdd..f2f0bf7 100644
--- a/makefile
+++ b/makefile
@@ -34,7 +34,7 @@ MPFR=$(GMP)
 
 ######################## do not edit below this line ##########################
 
-VERSION=0.4.4
+VERSION=0.4.5
 
 .SUFFIXES: .c .o
 
@@ -53,7 +53,9 @@ libmpc.a: $(OBJECTS)
 	$(AR) cru libmpc.a $(OBJECTS)
 	$(RANLIB) libmpc.a
 
-check: test tmul tsqr tdiv texp
+check: test tmul tsqr tdiv texp tabs
+	@echo Testing mpc_abs
+	./tabs
 	@echo Testing all functions
 	$(RM) -f mpc_test
 	./test
@@ -83,6 +85,9 @@ tdiv: tdiv.c libmpc.a
 texp: texp.c libmpc.a
 	@echo "Building texp" && ((test -e $(GMP)/lib/libgmp.a && test -e $(MPFR)/lib/libmpfr.a && $(CC) $(CFLAGS) $(INCLUDES) texp.c -o texp ./libmpc.a $(MPFR)/lib/libmpfr.a $(GMP)/lib/libgmp.a) || $(CC) $(CFLAGS) $(INCLUDES) -L. -L$(MPFR)/lib -L$(GMP)/lib texp.c -o texp -lmpc -lmpfr -lgmp)
 
+tabs: tabs.c libmpc.a
+	@echo "Building tabs" && ((test -e $(GMP)/lib/libgmp.a && test -e $(MPFR)/lib/libmpfr.a && $(CC) $(CFLAGS) $(INCLUDES) tabs.c -o tabs ./libmpc.a $(MPFR)/lib/libmpfr.a $(GMP)/lib/libgmp.a) || $(CC) $(CFLAGS) $(INCLUDES) -L. -L$(MPFR)/lib -L$(GMP)/lib tabs.c -o tabs -lmpc -lmpfr -lgmp)
+
 clean:
 	$(RM) *.o *~ libmpc.a test tmul tsqr tdiv texp mpc-$(VERSION).tar.gz mpc.aux mpc.cp mpc.cps mpc.dvi mpc.fn mpc.fns mpc.ky mpc.log mpc.pg mpc.ps mpc.toc mpc.tp mpc.vr mpc.vrs
 
diff --git a/mpc-impl.h b/mpc-impl.h
index bbffaf7..b2809b0 100644
--- a/mpc-impl.h
+++ b/mpc-impl.h
@@ -1,6 +1,6 @@
 /* mpc-impl.h -- Internal include file for mpc.
 
-Copyright (C) 2002, 2004 Andreas Enge, Paul Zimmermann
+Copyright (C) 2002, 2004, 2005 Andreas Enge, Paul Zimmermann
 
 This file is part of the MPC Library.
 
@@ -132,6 +132,7 @@ extern "C" {
   int  mpc_mul_naive     __MPC_PROTO ((mpc_ptr, mpc_srcptr, mpc_srcptr, mpc_rnd_t));
   int  mpc_mul_karatsuba __MPC_PROTO ((mpc_ptr, mpc_srcptr, mpc_srcptr, mpc_rnd_t));
   unsigned long  mpc_ceil_log2 __MPC_PROTO ((unsigned long));
+  int  mpc_abs_basic __MPC_PROTO ((mpfr_ptr, mpc_srcptr, mpc_rnd_t));
 
   /* forgotten in mpfr.h from GMP 4.1 */
 #ifdef NEED_CMP_ABS_PROTO
diff --git a/mpc.texi b/mpc.texi
index 377c0b4..b3d068a 100644
--- a/mpc.texi
+++ b/mpc.texi
@@ -9,8 +9,8 @@
 @end iftex
 @comment %**end of header
 
-@set VERSION 0.4.4
-@set DATE {November 2004}
+@set VERSION 0.4.5
+@set DATE {January 2005}
 
 @ifinfo
 @format
diff --git a/tabs.c b/tabs.c
new file mode 100644
index 0000000..f4642ba
--- /dev/null
+++ b/tabs.c
@@ -0,0 +1,177 @@
+/* tsqr -- test file for mpc_sqr.
+
+Copyright (C) 2005 Andreas Enge
+
+This file is part of the MPC Library.
+
+The MPC Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The MPC Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the MPC Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "mpfr.h"
+#include "mpc.h"
+#include "mpc-impl.h"
+
+#define SIGN(x) (((x) == 0 ? 0 : ((x) > 0 ? 1 : -1)))
+
+void cmpabs (mpc_srcptr, mpfr_rnd_t);
+void testabs (mpc_ptr);
+
+
+void cmpabs (mpc_srcptr x, mpfr_rnd_t rnd)
+   /* computes the absolute value of x with the two functions defined in   */
+   /* abs.c with rounding mode rnd and compares the results and return     */
+   /* values. We only consider cases where all precisions are identical.   */
+   /* We also compute the result with four times the precision and check   */
+   /* whether the rounding is correct. Error reports in this part of the   */
+   /* algorithm might still be wrong, though, since there are two          */
+   /* consecutive roundings.                                               */
+{
+   mpfr_t res1, res2, res3;
+   int    inexact1, inexact2;
+
+   mpfr_init2 (res1, MPC_MAX_PREC (x));
+   mpfr_init2 (res2, MPC_MAX_PREC (x));
+   mpfr_init2 (res3, 4 * MPC_MAX_PREC (x));
+
+   inexact1 = mpc_abs (res1, x, rnd);
+   inexact2 = mpc_abs_basic (res2, x, rnd);
+
+   if (mpfr_cmp (res1, res2))
+   {
+      fprintf (stderr, "abs and abs_basic differ for rnd=%s \nx=",
+               mpfr_print_rnd_mode(rnd));
+      mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+      fprintf (stderr, "\nabs       gives ");
+      mpfr_out_str (stderr, 2, 0, res1, GMP_RNDN);
+      fprintf (stderr, "\nabs_basic gives ");
+      mpfr_out_str (stderr, 2, 0, res2, GMP_RNDN);
+      fprintf (stderr, "\n");
+      exit (1);
+   }
+   if (SIGN (inexact1) != SIGN (inexact2))
+   {
+     fprintf (stderr, "The return values of abs and abs_basic differ for rnd=%s \nx=  ",
+              mpfr_print_rnd_mode(rnd));
+     mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+     fprintf (stderr, "\n|x| = ");
+     mpfr_out_str (stderr, 2, 0, res1, GMP_RNDN);
+     fprintf (stderr, "\nabs       gives %i", inexact1);
+     fprintf (stderr, "\nabs_basic gives %i", inexact2);
+     fprintf (stderr, "\n");
+     exit (1);
+   }
+   if (   (rnd == GMP_RNDU && (inexact1 < 0 || inexact2 < 0))
+       || ((rnd == GMP_RNDD || rnd == GMP_RNDZ) 
+            && (inexact1 > 0 || inexact2 > 0)))
+      {
+        fprintf (stderr, "Incorrect return value in abs or abs_basic for rnd=%s \nx=  ",
+                 mpfr_print_rnd_mode(rnd));
+        mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+        fprintf (stderr, "\n|x| = ");
+        mpfr_out_str (stderr, 2, 0, res1, GMP_RNDN);
+        fprintf (stderr, "\nabs       gives %i", inexact1);
+        fprintf (stderr, "\nabs_basic gives %i", inexact2);
+        fprintf (stderr, "\n");
+        exit (1);
+     }
+      
+
+   mpc_abs (res3, x, rnd);
+   mpfr_set (res2, res3, rnd);
+   if (mpfr_cmp (res1, res2))
+   {
+      fprintf (stderr, "rounding in abs might be incorrect for rnd=%s \nx=",
+               mpfr_print_rnd_mode(rnd));
+      mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+      fprintf (stderr, "\nabs                     gives ");
+      mpfr_out_str (stderr, 2, 0, res1, GMP_RNDN);
+      fprintf (stderr, "\nabs quadruple precision gives ");
+      mpfr_out_str (stderr, 2, 0, res3, GMP_RNDN);
+      fprintf (stderr, "\nand is rounded to                 ");
+      mpfr_out_str (stderr, 2, 0, res2, GMP_RNDN);
+      fprintf (stderr, "\n");
+      exit (1);
+  }
+
+   mpfr_clear (res1);
+   mpfr_clear (res2);
+   mpfr_clear (res3);
+}
+
+
+void testabs (mpc_ptr x)
+   /* executes the test for all rounding modes and all possible signs for */
+   /* the real and the imaginary part                                     */
+      
+{
+   mpfr_rnd_t rnd;
+   int i, j;
+   
+   for (i = 0; i < 2; i++)
+   {
+      mpfr_neg (MPC_RE (x), MPC_RE (x), GMP_RNDN);
+      for (j = 0; j < 2; j++)
+      {
+         mpfr_neg (MPC_IM (x), MPC_IM (x), GMP_RNDN);
+         for (rnd = 0; rnd < 4; rnd++)
+            cmpabs (x, rnd);
+      }
+   }
+}
+
+
+int main()
+{
+  mpc_t  x;
+  mpfr_t eps;
+  mp_prec_t prec;
+  int i;
+
+  mpc_init2 (x, 1000);
+  mpfr_init2 (eps, 10);
+
+  /* The most interesting cases arise when real and imaginary part differ */
+  /* much, since only then abs and abs_basic actually differ.             */
+  /* Check 1, 1+epsilon, (1+epsilon) + epsilon*i and the same values with */
+  /* real and imaginary part swapped for epsilon small powers of 2.       */
+  mpfr_set_ui (eps, 1, GMP_RNDN);
+  for (i = 0; i < 12; i++)
+  {
+     mpfr_div_2ui (eps, eps, 200, GMP_RNDN);
+     for (prec = 2; prec < 1000; prec += 44)
+     {
+        mpc_set_prec (x, prec);
+        mpc_set_ui_ui (x, 1, 0, MPC_RNDNN);
+        testabs (x);
+        mpfr_set (MPC_IM (x), eps, GMP_RNDN);
+        testabs (x);
+        mpfr_add_ui (MPC_RE (x), eps, 1, GMP_RNDN);
+        testabs (x);
+        mpc_set_ui_ui (x, 0, 1, MPC_RNDNN);
+        testabs (x);
+        mpfr_set (MPC_RE (x), eps, GMP_RNDN);
+        testabs (x);
+        mpfr_add_ui (MPC_IM (x), eps, 1, GMP_RNDN);
+        testabs (x);
+     }
+  }
+  
+  mpc_clear (x);
+
+  return 0;
+}