diff options
-rw-r--r-- | tests/Makefile.am | 14 | ||||
-rw-r--r-- | tests/tgeneric.tpl | 553 |
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 (¶ms, description_file); + init_parameters (¶ms); + + /* ask for enough memory */ + set_output_precision (¶ms, 4 * prec_max); + set_input_precision (¶ms, prec_max); + set_reference_precision (¶ms, 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 (¶ms, 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 (¶ms); + for (special = 0; special < last_special ; special++) + check_against_quadruple_precision (¶ms, prec_max, exp_min, exp_max, + special); + + clear_parameters (¶ms); +} + +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); + } +} |