[tests/] Add generic tests using templates.

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

2 files changed, 560 insertions, 7 deletions

diff --git a/tests/Makefile.am b/tests/Makefile.am
index bee7180..73d71a1 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -50,12 +50,12 @@ libmpc_tests_la_SOURCES=mpc-tests.h random.c tgeneric.c read_data.c	\
   print_parameter.c read_description.c read_line.c tpl_gmp.c		\
   tpl_mpc.c tpl_mpfr.c tpl_native.c
 
-TEMPLATES = data_check.tpl abs.dsc acos.dsc acosh.dsc asin.dsc		\
-  asinh.dsc atan.dsc atanh.dsc add.dsc add_fr.dsc add_si.dsc add_ui.dsc	\
-  arg.dsc conj.dsc cos.dsc cosh.dsc div.dsc div_fr.dsc exp.dsc fma.dsc	\
-  fr_div.dsc fr_sub.dsc log.dsc log10.dsc mul.dsc mul_fr.dsc neg.dsc	\
-  norm.dsc pow.dsc pow_fr.dsc pow_si.dsc pow_ui.dsc pow_z.dsc proj.dsc	\
-  sin.dsc sinh.dsc sub.dsc sub_fr.dsc tan.dsc tanh.dsc
+DESCRIPTIONS = abs.dsc acos.dsc acosh.dsc asin.dsc asinh.dsc atan.dsc	\
+  atanh.dsc add.dsc add_fr.dsc add_si.dsc add_ui.dsc arg.dsc conj.dsc	\
+  cos.dsc cosh.dsc div.dsc div_fr.dsc exp.dsc fma.dsc fr_div.dsc	\
+  fr_sub.dsc log.dsc log10.dsc mul.dsc mul_fr.dsc neg.dsc norm.dsc	\
+  pow.dsc pow_fr.dsc pow_si.dsc pow_ui.dsc pow_z.dsc proj.dsc sin.dsc	\
+  sinh.dsc sub.dsc sub_fr.dsc tan.dsc tanh.dsc
 DATA_SETS = abs.dat acos.dat acosh.dat asin.dat asinh.dat atan.dat	\
   atanh.dat add.dat add_fr.dat arg.dat conj.dat cos.dat cosh.dat	\
   div.dat div_fr.dat exp.dat fma.dat fr_div.dat fr_sub.dat		\
@@ -63,7 +63,7 @@ DATA_SETS = abs.dat acos.dat acosh.dat asin.dat asinh.dat atan.dat	\
   pow.dat pow_fr.dat pow_si.dat pow_ui.dat pow_z.dat proj.dat sin.dat	\
   sinh.dat sqr.dat sqrt.dat strtoc.dat sub.dat sub_fr.dat tan.dat	\
   tanh.dat
-EXTRA_DIST = $(DATA_SETS) $(TEMPLATES)
+EXTRA_DIST = data_check.tpl tgeneric.tpl $(DATA_SETS) $(DESCRIPTIONS)
 
 TESTS_ENVIRONMENT = $(VALGRIND)
 TESTS = $(check_PROGRAMS)
diff --git a/tests/tgeneric.tpl b/tests/tgeneric.tpl
new file mode 100644
index 0000000..befcc65
--- /dev/null
+++ b/tests/tgeneric.tpl
@@ -0,0 +1,553 @@
+/* tgeneric.tpl -- template file for generic tests.
+
+Copyright (C) 2008, 2009, 2010, 2011, 2012, 2013 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/ .
+*/
+
+#ifndef MPC_FUNCTION_CALL
+#error Define MPC_FUNCTION_CALL before including 'data_check.tpl'.
+#endif
+
+#include "mpc-tests.h"
+
+/* helper functions, defined after tgeneric */
+static void set_output_precision (mpc_fun_param_t *params, mpfr_prec_t prec);
+static void set_input_precision  (mpc_fun_param_t *params, mpfr_prec_t prec);
+static void set_reference_precision (mpc_fun_param_t *params,
+                                     mpfr_prec_t prec);
+static void first_rnd_mode       (mpc_fun_param_t *params);
+static int  is_valid_rnd_mode    (mpc_fun_param_t *params);
+static void next_rnd_mode        (mpc_fun_param_t *params);
+static int  double_rounding      (mpc_fun_param_t *params);
+static int  count_special_cases  (mpc_fun_param_t *params);
+static void random_params        (mpc_fun_param_t *params,
+                                  mpfr_exp_t exp_min, mpfr_exp_t exp_max,
+                                  int special);
+static void check_against_quadruple_precision (mpc_fun_param_t *params,
+                                               mpfr_prec_t prec,
+                                               mpfr_exp_t exp_min,
+                                               mpfr_exp_t exp_max,
+                                               int special);
+
+
+/* tgeneric(desc, prec_min, prec_max, step, exp_max) checks rounding with
+   random numbers:
+   - with precision ranging from prec_min to prec_max with an increment of
+   step,
+   - with exponent between -exp_max and exp_max.
+   - for pure real, pure imaginary and infinite random parameters.
+
+   It also checks parameter reuse.
+*/
+static void
+tgeneric_template (const char *description_file,
+                   mpfr_prec_t prec_min, mpfr_prec_t prec_max, mpfr_prec_t step,
+                   mpfr_exp_t exp_max)
+{
+  int special = 0;
+  int last_special;
+  mpfr_prec_t prec;
+  mpfr_exp_t exp_min;
+  mpc_fun_param_t params;
+
+  read_description (&params, description_file);
+  init_parameters (&params);
+
+  /* ask for enough memory */
+  set_output_precision (&params, 4 * prec_max);
+  set_input_precision (&params, prec_max);
+  set_reference_precision (&params, prec_max);
+
+  /* sanity checks */
+  exp_min = mpfr_get_emin ();
+  if (exp_max <= 0 || exp_max > mpfr_get_emax ())
+    exp_max = mpfr_get_emax();
+  if (-exp_max > exp_min)
+    exp_min = - exp_max;
+  if (step < 1)
+    step = 1;
+
+  /* check consistency with quadruple precision for random parameters */
+  for (prec = prec_min; prec <= prec_max; prec += step)
+    check_against_quadruple_precision (&params, prec, exp_min, exp_max, -1);
+
+
+  /* check consistency with quadruple precision for special values:
+     pure real, pure imaginary, or infinite arguments */
+  last_special = count_special_cases (&params);
+  for (special = 0; special < last_special ; special++)
+    check_against_quadruple_precision (&params, prec_max, exp_min, exp_max,
+                                       special);
+
+  clear_parameters (&params);
+}
+
+static void
+check_against_quadruple_precision (mpc_fun_param_t *params,
+                                   mpfr_prec_t prec,
+                                   mpfr_exp_t exp_min, mpfr_exp_t exp_max,
+                                   int special)
+{
+  mpc_operand_t *P = params->P; /* developer-friendly alias, used in macros */
+
+  set_input_precision (params, prec);
+  set_reference_precision (params, prec);
+
+  set_output_precision (params, 4 * prec);
+  random_params (params, exp_min, exp_max, special);
+
+  for (first_rnd_mode (params);
+       is_valid_rnd_mode (params);
+       next_rnd_mode (params))
+    {
+      MPC_FUNCTION_CALL;
+      while (double_rounding (params))
+        /* try another input parameters until no double rounding occurs when
+           the extra-precise result is rounded to working precision */
+        {
+          random_params (params, exp_min, exp_max, special);
+          MPC_FUNCTION_CALL;
+        }
+      set_output_precision (params, prec);
+
+      MPC_FUNCTION_CALL;
+      check_data (NULL, params, 0);
+
+#ifdef MPC_FUNCTION_CALL_SYMMETRIC
+      MPC_FUNCTION_CALL_SYMMETRIC;
+      check_data (NULL, params, 0);
+#endif
+
+#ifdef MPC_FUNCTION_CALL_REUSE_OP1
+      if (copy_parameter (params, 1, 2) == 0)
+        {
+          MPC_FUNCTION_CALL_REUSE_OP1;
+          check_data (NULL, params, 2);
+        }
+#endif
+
+#ifdef MPC_FUNCTION_CALL_REUSE_OP2
+      if (copy_parameter (params, 1, 3) == 0)
+        {
+          MPC_FUNCTION_CALL_REUSE_OP2;
+          check_data (NULL, params, 3);
+        }
+#endif
+
+#ifdef MPC_FUNCTION_CALL_REUSE_OP3
+      if (copy_parameter (params, 1, 4) == 0)
+        {
+          MPC_FUNCTION_CALL_REUSE_OP3;
+          check_data (NULL, params, 4);
+        }
+#endif
+
+      set_output_precision (params, 4 * prec);
+    }
+}
+
+static void
+set_output_precision (mpc_fun_param_t *params, mpfr_prec_t prec)
+{
+  int out;
+  const int start = 0;
+  const int end = params->nbout;
+  for (out = start; out < end; out++)
+    {
+      if (params->T[out] == MPFR)
+        mpfr_set_prec (params->P[out].mpfr, prec);
+      else if (params->T[out] == MPC)
+        mpc_set_prec (params->P[out].mpc, prec);
+    }
+}
+
+static void
+set_input_precision  (mpc_fun_param_t *params, mpfr_prec_t prec)
+{
+  int i;
+  const int start = params->nbout;
+  const int end = start + params->nbin;
+  for (i = start; i < end; i++)
+    {
+      if (params->T[i] == MPFR)
+        mpfr_set_prec (params->P[i].mpfr, prec);
+      else if (params->T[i] == MPC)
+        mpc_set_prec (params->P[i].mpc, prec);
+    }
+}
+
+static void
+set_reference_precision (mpc_fun_param_t *params, mpfr_prec_t prec)
+{
+  int i;
+  const int start = params->nbout + params->nbin;
+  const int end = start + params->nbout;
+  for (i = start; i < end; i++)
+    {
+      if (params->T[i] == MPFR)
+        mpfr_set_prec (params->P[i].mpfr_data.mpfr, prec);
+      else if (params->T[i] == MPC)
+        mpc_set_prec (params->P[i].mpc_data.mpc, prec);
+    }
+}
+
+enum {
+  SPECIAL_MINF,
+  SPECIAL_MZERO,
+  SPECIAL_PZERO,
+  SPECIAL_PINF,
+  SPECIAL_COUNT,
+} special_case;
+
+static int
+count_special_cases (mpc_fun_param_t *params)
+{
+  int i;
+  const int start = params->nbout;
+  const int end = start + params->nbin - 1; /* the last input parameter is the
+                                               rounding mode */
+  int count = 0;
+
+  for (i = start; i < end; i++)
+    {
+      if (params->T[i] == MPFR)
+        count += SPECIAL_COUNT;
+      else if (params->T[i] == MPC)
+        /* special + i x random and random + i x special */
+        count += 2 * SPECIAL_COUNT;
+    }
+  return count;
+}
+
+static void
+special_mpfr (mpfr_ptr x, int special)
+{
+  switch (special)
+    {
+    case SPECIAL_MINF:
+      mpfr_set_inf (x, -1);
+      break;
+    case SPECIAL_MZERO:
+      mpfr_set_zero (x, -1);
+      break;
+    case SPECIAL_PZERO:
+      mpfr_set_zero (x, +1);
+      break;
+    case SPECIAL_PINF:
+      mpfr_set_inf (x, +1);
+      break;
+   }
+}
+
+static void
+special_random_mpc (mpc_ptr z, mpfr_exp_t exp_min, mpfr_exp_t exp_max,
+                    int special)
+{
+  mpfr_ptr special_part;
+  mpfr_ptr random_part;
+  int mpfr_special;
+
+  if (special < SPECIAL_COUNT)
+    {
+      mpfr_special = special;
+      special_part = mpc_realref (z);
+      random_part  = mpc_imagref (z);
+    }
+  else
+    {
+      mpfr_special = special - SPECIAL_COUNT;
+      special_part = mpc_imagref (z);
+      random_part  = mpc_realref (z);
+    }
+
+  special_mpfr (special_part, mpfr_special);
+  test_random_mpfr (random_part, exp_min, exp_max, 128);
+}
+
+static void
+random_params (mpc_fun_param_t *params,
+               mpfr_exp_t exp_min, mpfr_exp_t exp_max,
+               int special)
+{
+  int i;
+  int base_index = 0;
+  const int start = params->nbout;
+  const int end = start + params->nbin - 1; /* the last input parameter is the
+                                               rounding mode */
+
+  for (i = start; i < end; i++)
+    {
+      switch (params->T[i])
+        {
+        case NATIVE_INT:
+        case NATIVE_UL:    case NATIVE_L:
+          /* TODO: draw random value */
+          break;
+
+        case NATIVE_D:     case NATIVE_LD:
+        case NATIVE_DC:    case NATIVE_LDC:
+          /* TODO: draw random value */
+          break;
+
+        case NATIVE_IM:    case NATIVE_UIM:
+          /* TODO: draw random value */
+          break;
+
+        case GMP_Z:
+          /* TODO: draw random value */
+          break;
+        case GMP_Q:
+          /* TODO: draw random value */
+          break;
+        case GMP_F:
+          /* TODO: draw random value */
+          break;
+
+        case MPFR:
+          if (base_index <= special
+              && special - base_index < SPECIAL_COUNT)
+            special_mpfr (params->P[i].mpfr, special - base_index);
+          else
+            test_random_mpfr (params->P[i].mpfr, exp_min, exp_max, 128);
+          base_index += SPECIAL_COUNT;
+          break;
+
+        case MPC:
+          if (base_index <= special
+              && special - base_index < 2 * SPECIAL_COUNT)
+            special_random_mpc (params->P[i].mpc, exp_min, exp_max,
+                                special - base_index);
+          else
+            test_random_mpc (params->P[i].mpc, exp_min, exp_max, 128);
+          base_index += 2 * SPECIAL_COUNT;
+          break;
+
+        case NATIVE_STRING:
+        case MPFR_INEX:    case MPFR_RND:
+        case MPC_INEX:     case MPC_RND:
+          /* unsupported types */
+          exit (1);
+          break;
+        }
+    }
+}
+
+/* return 1 if double rounding occurs;
+   return 0 otherwise */
+static int
+double_rounding_mpc (mpc_ptr lowprec, mpc_ptr hiprec, int inex, mpc_rnd_t rnd)
+{
+        /* reference value in high precision */
+        mpfr_ptr hiprec_re = mpc_realref (hiprec);
+        mpfr_ptr hiprec_im = mpc_imagref (hiprec);
+        mpfr_exp_t hiprec_err_re;
+        mpfr_exp_t hiprec_err_im;
+        /* destination precision and rounding mode */
+        mpfr_rnd_t lowprec_rnd_re = MPC_RND_RE (rnd);
+        mpfr_rnd_t lowprec_rnd_im = MPC_RND_IM (rnd);
+        mpfr_prec_t lowprec_prec_re = MPC_PREC_RE (lowprec);
+        mpfr_prec_t lowprec_prec_im = MPC_PREC_IM (lowprec);
+
+        /* real hiprec error is bounded by one ulp */
+        hiprec_err_re = MPC_PREC_RE (hiprec) - 1;
+        hiprec_err_im = MPC_PREC_IM (hiprec) - 1;
+
+        if (MPC_RND_RE (rnd) == MPFR_RNDN)
+          /* when rounding to nearest, use the trick for determining the
+             correct ternary value which is described in MPFR
+             documentation */
+          {
+            hiprec_err_re++; /* error is bounded by one half-ulp */
+            lowprec_rnd_re = MPFR_RNDZ;
+            lowprec_prec_re++;
+          }
+        if (MPC_RND_IM (rnd) == MPFR_RNDN)
+          /* same trick as above */
+          {
+            hiprec_err_im++;
+            lowprec_rnd_im = MPFR_RNDZ;
+            lowprec_prec_im++;
+          }
+
+        return ((mpfr_nan_p (hiprec_re)
+                 || mpfr_zero_p (hiprec_re)
+                 || mpfr_inf_p (hiprec_re)
+                 || MPC_INEX_RE (inex) == 0
+                 || mpfr_can_round (hiprec_re, hiprec_err_re,
+                                    MPC_RND_RE (rnd),
+                                    lowprec_rnd_re, lowprec_prec_re))
+                && (mpfr_nan_p (hiprec_im)
+                    || mpfr_zero_p (hiprec_im)
+                    || mpfr_inf_p (hiprec_im)
+                    || MPC_INEX_IM (inex) == 0
+                    || mpfr_can_round (hiprec_im, hiprec_err_im,
+                                       MPC_RND_IM (rnd),
+                                       lowprec_rnd_im, lowprec_prec_im)));
+}
+
+/* check whether double rounding occurs; if not, round extra precise output
+   value and set reference parameter */
+static int
+double_rounding (mpc_fun_param_t *params)
+{
+  int out;
+  const int offset = params->nbout + params->nbin;
+  const int rnd_index = params->nbout + params->nbin - 1;
+
+  for (out = 0; out < params->nbout; out++) {
+    if (params->T[out] == MPC)
+      {
+        MPC_ASSERT (params->T[0] == MPC_INEX);
+        MPC_ASSERT (params->T[offset] == MPC_INEX);
+        MPC_ASSERT (params->T[out + offset] == MPC);
+        MPC_ASSERT (params->T[rnd_index] == MPC_RND);
+
+        if (double_rounding_mpc (params->P[out + offset].mpc_data.mpc,
+                                 params->P[out].mpc,
+                                 params->P[0].mpc_inex,
+                                 params->P[rnd_index].mpc_rnd))
+          /* the hight-precision value and the exact value round to the same
+             low-precision value */
+          {
+            int inex;
+            inex = mpc_set (params->P[out + offset].mpc_data.mpc,
+                            params->P[out].mpc,
+                            params->P[rnd_index].mpc_rnd);
+            params->P[out + offset].mpc_data.known_sign_real = -1;
+            params->P[out + offset].mpc_data.known_sign_imag = -1;
+            params->P[offset].mpc_inex_data.real = MPC_INEX_RE (inex);
+            params->P[offset].mpc_inex_data.imag = MPC_INEX_IM (inex);
+          }
+        else
+          /* double rounding occurs */
+          return 1;
+      }
+    else if (params->T[out] == MPFR)
+      {
+        /* TODO */
+      }
+  }
+  return 0;
+}
+
+
+/* helper functions for iterating over mpfr rounding modes */
+
+#define FIRST_MPFR_RND_MODE MPFR_RNDN
+
+static mpfr_rnd_t
+next_mpfr_rnd_mode (mpfr_rnd_t curr)
+{
+  switch (curr)
+    {
+    case MPFR_RNDN:
+      return MPFR_RNDZ;
+    case MPFR_RNDZ:
+      return MPFR_RNDU;
+    case MPFR_RNDU:
+      return MPFR_RNDD;
+    default:
+      /* return invalid guard value in mpfr_rnd_t */
+#if MPFR_VERSION_MAJOR < 3
+      return MPFR_RNDNA;
+#else
+      return MPFR_RNDA; /* valid rounding type, but not used in mpc */
+#endif
+    }
+}
+
+static int
+is_valid_mpfr_rnd_mode (mpfr_rnd_t curr)
+/* returns 1 if curr is a valid rounding mode, and 0 otherwise */
+{
+   if (   curr == MPFR_RNDN || curr == MPFR_RNDZ
+       || curr == MPFR_RNDU || curr == MPFR_RNDD)
+      return 1;
+   else
+      return 0;
+}
+
+static void
+first_rnd_mode (mpc_fun_param_t *params)
+{
+  int rnd_mode_index = params->nbout + params->nbin - 1;
+  switch (params->T[rnd_mode_index])
+    {
+    case MPC_RND:
+      params->P[rnd_mode_index].mpc_rnd =
+        MPC_RND(FIRST_MPFR_RND_MODE, FIRST_MPFR_RND_MODE);
+      break;
+    case MPFR_RND:
+      params->P[rnd_mode_index].mpfr_rnd = FIRST_MPFR_RND_MODE;
+      break;
+    default:
+      printf ("The rounding mode is expected to be"
+              " the last input parameter.\n");
+      exit (-1);
+    }
+}
+
+static void
+next_rnd_mode (mpc_fun_param_t *params)
+{
+  mpfr_rnd_t rnd_re, rnd_im;
+  int rnd_mode_index = params->nbout + params->nbin - 1;
+  switch (params->T[rnd_mode_index])
+    {
+    case MPC_RND:
+      rnd_re = MPC_RND_RE (params->P[rnd_mode_index].mpc_rnd);
+      rnd_im = MPC_RND_IM (params->P[rnd_mode_index].mpc_rnd);
+      rnd_im = next_mpfr_rnd_mode (rnd_im);
+      if (!is_valid_mpfr_rnd_mode (rnd_im))
+        {
+          rnd_re = next_mpfr_rnd_mode (rnd_re);
+          rnd_im = FIRST_MPFR_RND_MODE;
+        }
+      params->P[rnd_mode_index].mpc_rnd = MPC_RND(rnd_re, rnd_im);
+      break;
+    case MPFR_RND:
+      params->P[rnd_mode_index].mpfr_rnd =
+        next_mpfr_rnd_mode (params->P[rnd_mode_index].mpfr_rnd);
+      break;
+    default:
+      printf ("The rounding mode is expected to be"
+              " the last input parameter.\n");
+      exit (-1);
+    }
+}
+
+static int
+is_valid_rnd_mode (mpc_fun_param_t *params)
+/* returns 1 if curr is a valid rounding mode, and 0 otherwise */
+{
+  mpfr_rnd_t rnd_re, rnd_im;
+  int rnd_mode_index = params->nbout + params->nbin - 1;
+  switch (params->T[rnd_mode_index])
+    {
+    case MPC_RND:
+      rnd_re = MPC_RND_RE (params->P[rnd_mode_index].mpc_rnd);
+      rnd_im = MPC_RND_IM (params->P[rnd_mode_index].mpc_rnd);
+      return is_valid_mpfr_rnd_mode (rnd_re)
+        && is_valid_mpfr_rnd_mode (rnd_im);
+    case MPFR_RND:
+      return is_valid_mpfr_rnd_mode (params->P[rnd_mode_index].mpfr_rnd);
+    default:
+      printf ("The rounding mode is expected to be"
+              " the last input parameter.\n");
+      exit (-1);
+    }
+}