new mpc_check file to check against GNU libc

5 files changed, 468 insertions, 2 deletions

diff --git a/tests/Makefile.am b/tests/Makefile.am
index 2b7d447..ee9fb31 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -37,7 +37,8 @@ check_PROGRAMS = tabs tacos tacosh tadd tadd_fr tadd_si tadd_ui targ	\
   tpow_d tpow_fr tpow_ld tpow_si tpow_ui tpow_z tprec tproj treal	\
   treimref trootofunity                                                 \
   tset tsin tsin_cos tsinh tsqr tsqrt tstrtoc tsub tsub_fr	        \
-  tsub_ui tsum tswap ttan ttanh tui_div tui_ui_sub tget_version exceptions
+  tsub_ui tsum tswap ttan ttanh tui_div tui_ui_sub tget_version exceptions \
+  mpc_check
 
 check_LTLIBRARIES=libmpc-tests.la
 libmpc_tests_la_SOURCES = mpc-tests.h check_data.c clear_parameters.c	\
@@ -65,7 +66,8 @@ DATA_SETS = abs.dat acos.dat acosh.dat add.dat add_fr.dat arg.dat	\
   pow_fr.dat pow_si.dat pow_ui.dat pow_z.dat proj.dat rootofunity.dat   \
   sin.dat sinh.dat                                                      \
   sqr.dat sqrt.dat strtoc.dat sub.dat sub_fr.dat tan.dat tanh.dat
-EXTRA_DIST = data_check.tpl tgeneric.tpl $(DATA_SETS) $(DESCRIPTIONS)
+EXTRA_DIST = data_check.tpl tgeneric.tpl $(DATA_SETS) $(DESCRIPTIONS) \
+             mpc_check_template1.c mpc_check_template2.c mpc_check_template3.c
 
 LOG_COMPILER = $(VALGRIND)
 AM_LOG_FLAGS = $(VALGRIND_OPTS)
diff --git a/tests/mpc_check.c b/tests/mpc_check.c
new file mode 100644
index 0000000..569699e
--- /dev/null
+++ b/tests/mpc_check.c
@@ -0,0 +1,214 @@
+/* mpc_check -- compare mpc functions against the GNU libc implementation
+
+Copyright (C) 2020 INRIA
+
+This file is part of GNU MPC.
+
+GNU MPC 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 3 of the License, or (at your
+option) any later version.
+
+GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ .
+*/
+
+/* the GNU libc provides the following functions (as of 2.31),
+   with 'f' suffix for the float/binary32 version, with no suffix
+   for the double/binary64 version, with 'l' suffix for the long double
+   version, and with 'f128' suffix for the __float128 version:
+
+   cabs     casinh    cexp      csinh
+   cacos    catan     clog      csqrt
+   cacosh   catanh    clog10    ctan
+   carg     ccos      cpow      ctanh
+   casin    ccosh     csin
+*/
+
+#define _GNU_SOURCE /* for clog10 */
+#include <stdlib.h>
+#include <string.h>
+#include <complex.h>
+#include "mpc-tests.h"
+
+gmp_randstate_t state;
+unsigned long seed = 1;
+int verbose = 0;
+
+static unsigned long
+ulp_error (mpfr_t x, mpfr_t y)
+{
+  mpfr_t z;
+  mpfr_prec_t p = mpfr_get_prec (y);
+  unsigned long n;
+
+  if (mpfr_cmp (x, y) == 0)
+    return 0;
+
+  mpfr_init2 (z, p);
+  mpfr_sub (z, x, y, MPFR_RNDN);
+  mpfr_abs (z, z, MPFR_RNDN);
+  /* divide by ulp(y) = 2^(EXP(y) - p) */
+  mpfr_div_2si (z, z, mpfr_get_exp (y) - p, MPFR_RNDN);
+  n = mpfr_get_ui (z, MPFR_RNDZ);
+  mpfr_clear (z);
+  return n;
+}
+
+#define FOO add
+#define CFOO(x,y) (x+y)
+#include "mpc_check_template3.c"
+
+#define FOO sub
+#define CFOO(x,y) (x-y)
+#include "mpc_check_template3.c"
+
+#define FOO mul
+#define CFOO(x,y) (x*y)
+#include "mpc_check_template3.c"
+
+#define FOO div
+#define CFOO(x,y) (x/y)
+#include "mpc_check_template3.c"
+
+#define FOO pow
+#include "mpc_check_template3.c"
+
+#define FOO abs
+#include "mpc_check_template2.c"
+
+#define FOO arg
+#include "mpc_check_template2.c"
+
+#define FOO sqrt
+#include "mpc_check_template1.c"
+
+#define FOO acos
+#include "mpc_check_template1.c"
+
+#define FOO acosh
+#include "mpc_check_template1.c"
+
+#define FOO asin
+#include "mpc_check_template1.c"
+
+#define FOO asinh
+#include "mpc_check_template1.c"
+
+#define FOO atan
+#include "mpc_check_template1.c"
+
+#define FOO atanh
+#include "mpc_check_template1.c"
+
+#define FOO cos
+#include "mpc_check_template1.c"
+
+#define FOO cosh
+#include "mpc_check_template1.c"
+
+#define FOO exp
+#include "mpc_check_template1.c"
+
+#define FOO log
+#include "mpc_check_template1.c"
+
+#define FOO log10
+#include "mpc_check_template1.c"
+
+#define FOO sin
+#include "mpc_check_template1.c"
+
+#define FOO sinh
+#include "mpc_check_template1.c"
+
+#define FOO tan
+#include "mpc_check_template1.c"
+
+#define FOO tanh
+#include "mpc_check_template1.c"
+
+int
+main (int argc, char *argv[])
+{
+  mpfr_prec_t p = 53; /* default precision */
+  unsigned long n = 1000000;
+
+  while (argc >= 2 && argv[1][0] == '-')
+    {
+      if (argc >= 3 && strcmp (argv[1], "-p") == 0)
+        {
+          p = atoi (argv[2]);
+          argc -= 2;
+          argv += 2;
+        }
+      else if (argc >= 3 && strcmp (argv[1], "-seed") == 0)
+        {
+          seed = atoi (argv[2]);
+          argc -= 2;
+          argv += 2;
+        }
+      else if (strcmp (argv[1], "-v") == 0)
+        {
+          verbose ++;
+          argc --;
+          argv ++;
+        }
+      else
+        {
+          fprintf (stderr, "Unknown option %s\n", argv[1]);
+          exit (1);
+        }
+    }
+
+  MPC_ASSERT (p == 53);
+
+  if (p == 53)
+    {
+      mpfr_set_emin (-1073);
+      mpfr_set_emax (1024);
+    }
+
+  gmp_randinit_default (state);
+
+  printf ("Using random seed %lu\n", seed);
+
+  /* (complex,complex) -> complex */
+  test_add (p, n);
+  test_sub (p, n);
+  test_mul (p, n);
+  test_div (p, n);
+  test_pow (p, n);
+
+  /* complex -> real */
+  test_abs (p, n);
+  test_arg (p, n);
+
+  /* complex -> complex */
+  test_sqrt (p, n);
+  test_acos (p, n);
+  test_acosh (p, n);
+  test_asin (p, n);
+  test_asinh (p, n);
+  test_atan (p, n);
+  test_atanh (p, n);
+  test_cos (p, n);
+  test_cosh (p, n);
+  test_exp (p, n);
+  test_log (p, n);
+  test_log10 (p, n);
+  test_sin (p, n);
+  test_sinh (p, n);
+  test_tan (p, n);
+  test_tanh (p, n);
+
+  gmp_randclear (state);
+
+  return 0;
+}
+   
diff --git a/tests/mpc_check_template1.c b/tests/mpc_check_template1.c
new file mode 100644
index 0000000..771054b
--- /dev/null
+++ b/tests/mpc_check_template1.c
@@ -0,0 +1,82 @@
+/* template for complex -> complex function, for example acos
+
+Copyright (C) 2020 INRIA
+
+This file is part of GNU MPC.
+
+GNU MPC 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 3 of the License, or (at your
+option) any later version.
+
+GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ .
+*/
+
+#define CAT1(X,Y) X ## Y
+#define CAT2(X,Y) CAT1(X,Y)
+#define STRINGIFY(x) #x
+#define TOSTRING(x) STRINGIFY(x)
+
+#define FUN CAT2(test_,FOO)
+#define MPC_FOO CAT2(mpc_,FOO)
+#define CFOO CAT2(c,FOO)
+#define BAR TOSTRING(FOO)
+
+static void
+FUN (mpfr_prec_t p, unsigned long n)
+{
+  unsigned long i = 0;
+  mpc_t x, z, t;
+  double complex xx, zz;
+  int inex;
+  unsigned long errors = 0, max_err_re = 0, max_err_im = 0;
+
+  gmp_randseed_ui (state, seed);
+
+  mpc_init2 (x, p);
+  mpc_init2 (z, p);
+  mpc_init2 (t, p);
+  for (i = 0; i < n; i++)
+    {
+      mpc_urandom (x, state);
+      inex = MPC_FOO (z, x, MPC_RNDNN);
+      mpfr_subnormalize (mpc_realref (z), MPC_INEX_RE(inex), MPFR_RNDN);
+      mpfr_subnormalize (mpc_imagref (z), MPC_INEX_IM(inex), MPFR_RNDN);
+      xx = mpc_get_dc (x, MPC_RNDNN);
+      zz = CFOO (xx);
+      mpc_set_dc (t, zz, MPFR_RNDN);
+      if (mpc_cmp (z, t) != 0)
+        {
+          unsigned long err_re = ulp_error (mpc_realref (t), mpc_realref (z));
+          unsigned long err_im = ulp_error (mpc_imagref (t), mpc_imagref (z));
+          if (verbose > 0 && (err_re > max_err_re || err_im > max_err_im))
+            {
+              mpfr_printf ("   mpc_%s and c%s differ by (%lu,%lu) ulp(s)\n      for x=(%Re,%Re)\n",
+                           BAR, BAR, err_re, err_im,
+                           mpc_realref (x), mpc_imagref (x));
+              mpfr_printf ("   mpc_%s gives (%Re,%Re)\n", BAR,
+                           mpc_realref (z), mpc_imagref (z));
+              mpfr_printf ("      c%s gives (%Re,%Re)\n", BAR,
+                           mpc_realref (t), mpc_imagref (t));
+            }
+          errors ++;
+          if (err_re > max_err_re)
+            max_err_re = err_re;
+          if (err_im > max_err_im)
+            max_err_im = err_im;
+        }
+    }
+  mpc_clear (x);
+  mpc_clear (z);
+  mpc_clear (t);
+  printf ("Number of errors for %s: %lu/%lu (max %lu,%lu)\n", BAR,
+          errors, n, max_err_re, max_err_im);
+}
+
+#undef FOO
diff --git a/tests/mpc_check_template2.c b/tests/mpc_check_template2.c
new file mode 100644
index 0000000..fdb26c3
--- /dev/null
+++ b/tests/mpc_check_template2.c
@@ -0,0 +1,77 @@
+/* template for complex -> real function, for example abs
+
+Copyright (C) 2020 INRIA
+
+This file is part of GNU MPC.
+
+GNU MPC 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 3 of the License, or (at your
+option) any later version.
+
+GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ .
+*/
+
+#define CAT1(X,Y) X ## Y
+#define CAT2(X,Y) CAT1(X,Y)
+#define STRINGIFY(x) #x
+#define TOSTRING(x) STRINGIFY(x)
+
+#define FUN CAT2(test_,FOO)
+#define MPC_FOO CAT2(mpc_,FOO)
+#define CFOO CAT2(c,FOO)
+#define BAR TOSTRING(FOO)
+
+static void
+FUN (mpfr_prec_t p, unsigned long n)
+{
+  unsigned long i = 0;
+  mpc_t x;
+  mpfr_t z, t;
+  double complex xx;
+  double zz;
+  int inex;
+  unsigned long errors = 0, max_err = 0;
+
+  gmp_randseed_ui (state, seed);
+
+  mpc_init2 (x, p);
+  mpfr_init2 (z, p);
+  mpfr_init2 (t, p);
+  for (i = 0; i < n; i++)
+    {
+      mpc_urandom (x, state);
+      inex = MPC_FOO (z, x, MPC_RNDNN);
+      mpfr_subnormalize (z, inex, MPFR_RNDN);
+      xx = mpc_get_dc (x, MPC_RNDNN);
+      zz = CFOO (xx);
+      mpfr_set_d (t, zz, MPFR_RNDN);
+      if (mpfr_cmp (z, t) != 0)
+        {
+          unsigned long err = ulp_error (t, z);
+          if (verbose > 0 && err > max_err)
+            {
+              mpfr_printf ("   mpc_%s and c%s differ by %lu ulp(s)\n      for x=(%Re,%Re)\n",
+                           BAR, BAR, err,
+                           mpc_realref (x), mpc_imagref (x));
+              mpfr_printf ("   mpc_%s gives %Re\n", BAR, z);
+              mpfr_printf ("      c%s gives %Re\n", BAR, t);
+            }
+          errors++;
+          if (err > max_err)
+            max_err = err;
+        }
+    }
+  mpc_clear (x);
+  mpfr_clear (z);
+  mpfr_clear (t);
+  printf ("Number of errors for %s: %lu/%lu (max %lu)\n", BAR, errors, n, max_err);
+}
+
+#undef FOO
diff --git a/tests/mpc_check_template3.c b/tests/mpc_check_template3.c
new file mode 100644
index 0000000..506e7ac
--- /dev/null
+++ b/tests/mpc_check_template3.c
@@ -0,0 +1,91 @@
+/* template for (complex,complex) -> complex function, for example add
+
+Copyright (C) 2020 INRIA
+
+This file is part of GNU MPC.
+
+GNU MPC 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 3 of the License, or (at your
+option) any later version.
+
+GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ .
+*/
+
+#define CAT1(X,Y) X ## Y
+#define CAT2(X,Y) CAT1(X,Y)
+#define STRINGIFY(x) #x
+#define TOSTRING(x) STRINGIFY(x)
+
+#define FUN CAT2(test_,FOO)
+#define MPC_FOO CAT2(mpc_,FOO)
+#ifndef CFOO
+#define CFOO CAT2(c,FOO)
+#endif
+#define BAR TOSTRING(FOO)
+
+static void
+FUN (mpfr_prec_t p, unsigned long n)
+{
+  unsigned long i = 0;
+  mpc_t x, y, z, t;
+  double complex xx, yy, zz;
+  int inex;
+  unsigned long errors = 0, max_err_re = 0, max_err_im = 0;
+
+  gmp_randseed_ui (state, seed);
+
+  mpc_init2 (x, p);
+  mpc_init2 (y, p);
+  mpc_init2 (z, p);
+  mpc_init2 (t, p);
+  for (i = 0; i < n; i++)
+    {
+      mpc_urandom (x, state);
+      mpc_urandom (y, state);
+      inex = MPC_FOO (z, x, y, MPC_RNDNN);
+      mpfr_subnormalize (mpc_realref (z), MPC_INEX_RE(inex), MPFR_RNDN);
+      mpfr_subnormalize (mpc_imagref (z), MPC_INEX_IM(inex), MPFR_RNDN);
+      xx = mpc_get_dc (x, MPC_RNDNN);
+      yy = mpc_get_dc (y, MPC_RNDNN);
+      zz = CFOO(xx, yy);
+      mpc_set_dc (t, zz, MPFR_RNDN);
+      if (mpc_cmp (z, t) != 0)
+        {
+          unsigned long err_re = ulp_error (mpc_realref (t), mpc_realref (z));
+          unsigned long err_im = ulp_error (mpc_imagref (t), mpc_imagref (z));
+          if (verbose > 0 && (err_re > max_err_re || err_im > max_err_im))
+            {
+              printf ("   mpc_%s and c%s differ by (%lu,%lu) ulp(s)\n",
+                      BAR, BAR, err_re, err_im);
+              mpfr_printf ("      for x=(%Re,%Re)\n          y=(%Re,%Re)\n",
+                           mpc_realref (x), mpc_imagref (x),
+                           mpc_realref (y), mpc_imagref (y));
+              mpfr_printf ("   mpc_%s gives (%Re,%Re)\n", BAR,
+                           mpc_realref (z), mpc_imagref (z));
+              mpfr_printf ("      c%s gives (%Re,%Re)\n", BAR,
+                           mpc_realref (t), mpc_imagref (t));
+            }
+          errors ++;
+          if (err_re > max_err_re)
+            max_err_re = err_re;
+          if (err_im > max_err_im)
+            max_err_im = err_im;
+        }
+    }
+  mpc_clear (x);
+  mpc_clear (y);
+  mpc_clear (z);
+  mpc_clear (t);
+  printf ("Number of errors for %s: %lu/%lu (max %lu,%lu)\n", BAR, errors, n,
+          max_err_re, max_err_im);
+}
+
+#undef FOO
+#undef CFOO