/* double_rounding.c -- Functions for checking double rounding.
Copyright (C) 2013, 2014 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/ .
*/
#include "mpc-tests.h"
/* return 1 if double rounding occurs;
return 0 otherwise */
static int
double_rounding_mpfr (mpfr_ptr lowprec,
mpfr_ptr hiprec, int hiprec_inex, mpfr_rnd_t hiprec_rnd)
{
mpfr_exp_t hiprec_err;
mpfr_rnd_t lowprec_rnd = hiprec_rnd;
mpfr_prec_t lowprec_prec = mpfr_get_prec (lowprec);
/* hiprec error is bounded by one ulp */
hiprec_err = mpfr_get_prec (hiprec) - 1;
if (hiprec_rnd == MPFR_RNDN)
/* when rounding to nearest, use the trick for determining the
correct ternary value which is described in the MPFR
documentation */
{
hiprec_err++; /* error is bounded by one half-ulp */
lowprec_rnd = MPFR_RNDZ;
lowprec_prec++;
}
return (hiprec_inex == 0
|| mpfr_can_round (hiprec, hiprec_err, hiprec_rnd,
lowprec_rnd, lowprec_prec));
}
/* return 1 if double rounding occurs;
return 0 otherwise */
static int
double_rounding_mpc (mpc_ptr lowprec,
mpc_ptr hiprec, int hiprec_inex, mpc_rnd_t hiprec_rnd)
{
mpfr_ptr lowprec_re = mpc_realref (lowprec);
mpfr_ptr lowprec_im = mpc_imagref (lowprec);
mpfr_ptr hiprec_re = mpc_realref (hiprec);
mpfr_ptr hiprec_im = mpc_imagref (hiprec);
int inex_re = MPC_INEX_RE (hiprec_inex);
int inex_im = MPC_INEX_IM (hiprec_inex);
mpfr_rnd_t rnd_re = MPC_RND_RE (hiprec_rnd);
mpfr_rnd_t rnd_im = MPC_RND_IM (hiprec_rnd);
return (double_rounding_mpfr (lowprec_re, hiprec_re, inex_re, rnd_re)
&& double_rounding_mpfr (lowprec_im, hiprec_im, inex_im, rnd_im));
}
/* check whether double rounding occurs; if not, round extra precise output
value and set reference parameter */
int
double_rounding (mpc_fun_param_t *params)
{
int out;
const int offset = params->nbout + params->nbin;
int rnd_index = offset - params->nbrnd;
for (out = 0; out < params->nbout; out++) {
if (params->T[out] == MPC)
{
int inex;
MPC_ASSERT ((params->T[0] == MPC_INEX)
|| (params->T[0] == MPCC_INEX));
MPC_ASSERT ((params->T[offset] == MPC_INEX)
|| (params->T[offset] == MPCC_INEX));
MPC_ASSERT (params->T[out + offset] == MPC);
MPC_ASSERT (params->T[rnd_index] == MPC_RND);
/*
For the time being, there may be either one or two rounding modes;
in the latter case, we assume that there are three outputs:
the inexact value and two complex numbers.
*/
inex = (params->nbrnd == 1 ? params->P[0].mpc_inex
: (out == 1 ? MPC_INEX1 (params->P[0].mpcc_inex)
: MPC_INEX2 (params->P[0].mpcc_inex)));
if (double_rounding_mpc (params->P[out + offset].mpc_data.mpc,
params->P[out].mpc,
inex,
params->P[rnd_index].mpc_rnd))
/* the high-precision value and the exact value round to the same
low-precision value */
{
int inex_hp, inex_re, inex_im;
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;
/* no double rounding means that the ternary value may come from
the high-precision calculation or from the rounding */
if (params->nbrnd == 1)
inex_hp = params->P[0].mpc_inex;
else /* nbrnd == 2 */
if (out == 1)
inex_hp = MPC_INEX1 (params->P[0].mpcc_inex);
else /* out == 2 */
inex_hp = MPC_INEX2 (params->P[0].mpcc_inex);
if (MPC_INEX_RE (inex) == 0)
inex_re = MPC_INEX_RE (inex_hp);
else
inex_re = MPC_INEX_RE (inex);
if (MPC_INEX_IM (inex) == 0)
inex_im = MPC_INEX_IM (inex_hp);
else
inex_im = MPC_INEX_IM (inex);
if (params->nbrnd == 1) {
params->P[offset].mpc_inex_data.real = inex_re;
params->P[offset].mpc_inex_data.imag = inex_im;
}
else /* nbrnd == 2 */
if (out == 1)
params->P[offset].mpcc_inex = MPC_INEX (inex_re, inex_im);
else /* out == 2 */
params->P[offset].mpcc_inex
= MPC_INEX12 (params->P[offset].mpcc_inex,
MPC_INEX (inex_re, inex_im));
rnd_index++;
}
else
/* double rounding occurs */
return 1;
}
else if (params->T[out] == MPFR)
{
MPC_ASSERT (params->T[0] == MPFR_INEX);
MPC_ASSERT (params->T[offset] == MPFR_INEX);
MPC_ASSERT (params->T[out + offset] == MPFR);
MPC_ASSERT (params->T[rnd_index] == MPFR_RND);
if (double_rounding_mpfr (params->P[out + offset].mpfr_data.mpfr,
params->P[out].mpfr,
params->P[0].mpfr_inex,
params->P[rnd_index].mpfr_rnd))
/* the hight-precision value and the exact value round to the same
low-precision value */
{
int inex;
inex = mpfr_set (params->P[out + offset].mpfr_data.mpfr,
params->P[out].mpfr,
params->P[rnd_index].mpfr_rnd);
params->P[out + offset].mpfr_data.known_sign = -1;
/* no double rounding means that the ternary value may comes from
the high-precision calculation or from the rounding */
if (inex == 0)
params->P[offset].mpfr_inex = params->P[0].mpfr_inex;
else
params->P[offset].mpfr_inex = inex;
rnd_index++;
}
else
/* double rounding occurs */
return 1;
}
}
return 0;
}