path: root/tests/tgeneric.tpl

/* 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

/* helper functions, defined after tgeneric */
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);
    }
}


/* special cases */

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;
}