1 /* double_rounding.c -- Functions for checking double rounding. 2 3 Copyright (C) 2013, 2014 INRIA 4 5 This file is part of GNU MPC. 6 7 GNU MPC is free software; you can redistribute it and/or modify it under 8 the terms of the GNU Lesser General Public License as published by the 9 Free Software Foundation; either version 3 of the License, or (at your 10 option) any later version. 11 12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY 13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for 15 more details. 16 17 You should have received a copy of the GNU Lesser General Public License 18 along with this program. If not, see http://www.gnu.org/licenses/ . 19 */ 20 21 #include "mpc-tests.h" 22 23 /* return 1 if double rounding occurs; 24 return 0 otherwise */ 25 static int 26 double_rounding_mpfr (mpfr_ptr lowprec, 27 mpfr_ptr hiprec, int hiprec_inex, mpfr_rnd_t hiprec_rnd) 28 { 29 mpfr_exp_t hiprec_err; 30 mpfr_rnd_t lowprec_rnd = hiprec_rnd; 31 mpfr_prec_t lowprec_prec = mpfr_get_prec (lowprec); 32 33 /* hiprec error is bounded by one ulp */ 34 hiprec_err = mpfr_get_prec (hiprec) - 1; 35 36 if (hiprec_rnd == MPFR_RNDN) 37 /* when rounding to nearest, use the trick for determining the 38 correct ternary value which is described in the MPFR 39 documentation */ 40 { 41 hiprec_err++; /* error is bounded by one half-ulp */ 42 lowprec_rnd = MPFR_RNDZ; 43 lowprec_prec++; 44 } 45 46 return (hiprec_inex == 0 47 || mpfr_can_round (hiprec, hiprec_err, hiprec_rnd, 48 lowprec_rnd, lowprec_prec)); 49 } 50 51 /* return 1 if double rounding occurs; 52 return 0 otherwise */ 53 static int 54 double_rounding_mpc (mpc_ptr lowprec, 55 mpc_ptr hiprec, int hiprec_inex, mpc_rnd_t hiprec_rnd) 56 { 57 mpfr_ptr lowprec_re = mpc_realref (lowprec); 58 mpfr_ptr lowprec_im = mpc_imagref (lowprec); 59 mpfr_ptr hiprec_re = mpc_realref (hiprec); 60 mpfr_ptr hiprec_im = mpc_imagref (hiprec); 61 int inex_re = MPC_INEX_RE (hiprec_inex); 62 int inex_im = MPC_INEX_IM (hiprec_inex); 63 mpfr_rnd_t rnd_re = MPC_RND_RE (hiprec_rnd); 64 mpfr_rnd_t rnd_im = MPC_RND_IM (hiprec_rnd); 65 66 return (double_rounding_mpfr (lowprec_re, hiprec_re, inex_re, rnd_re) 67 && double_rounding_mpfr (lowprec_im, hiprec_im, inex_im, rnd_im)); 68 } 69 70 /* check whether double rounding occurs; if not, round extra precise output 71 value and set reference parameter */ 72 int 73 double_rounding (mpc_fun_param_t *params) 74 { 75 int out; 76 const int offset = params->nbout + params->nbin; 77 int rnd_index = offset - params->nbrnd; 78 79 for (out = 0; out < params->nbout; out++) { 80 if (params->T[out] == MPC) 81 { 82 int inex; 83 84 MPC_ASSERT ((params->T[0] == MPC_INEX) 85 || (params->T[0] == MPCC_INEX)); 86 MPC_ASSERT ((params->T[offset] == MPC_INEX) 87 || (params->T[offset] == MPCC_INEX)); 88 MPC_ASSERT (params->T[out + offset] == MPC); 89 MPC_ASSERT (params->T[rnd_index] == MPC_RND); 90 91 /* 92 For the time being, there may be either one or two rounding modes; 93 in the latter case, we assume that there are three outputs: 94 the inexact value and two complex numbers. 95 */ 96 inex = (params->nbrnd == 1 ? params->P[0].mpc_inex 97 : (out == 1 ? MPC_INEX1 (params->P[0].mpcc_inex) 98 : MPC_INEX2 (params->P[0].mpcc_inex))); 99 100 if (double_rounding_mpc (params->P[out + offset].mpc_data.mpc, 101 params->P[out].mpc, 102 inex, 103 params->P[rnd_index].mpc_rnd)) 104 /* the high-precision value and the exact value round to the same 105 low-precision value */ 106 { 107 int inex_hp, inex_re, inex_im; 108 inex = mpc_set (params->P[out + offset].mpc_data.mpc, 109 params->P[out].mpc, 110 params->P[rnd_index].mpc_rnd); 111 params->P[out + offset].mpc_data.known_sign_real = -1; 112 params->P[out + offset].mpc_data.known_sign_imag = -1; 113 114 /* no double rounding means that the ternary value may come from 115 the high-precision calculation or from the rounding */ 116 if (params->nbrnd == 1) 117 inex_hp = params->P[0].mpc_inex; 118 else /* nbrnd == 2 */ 119 if (out == 1) 120 inex_hp = MPC_INEX1 (params->P[0].mpcc_inex); 121 else /* out == 2 */ 122 inex_hp = MPC_INEX2 (params->P[0].mpcc_inex); 123 124 if (MPC_INEX_RE (inex) == 0) 125 inex_re = MPC_INEX_RE (inex_hp); 126 else 127 inex_re = MPC_INEX_RE (inex); 128 if (MPC_INEX_IM (inex) == 0) 129 inex_im = MPC_INEX_IM (inex_hp); 130 else 131 inex_im = MPC_INEX_IM (inex); 132 133 if (params->nbrnd == 1) { 134 params->P[offset].mpc_inex_data.real = inex_re; 135 params->P[offset].mpc_inex_data.imag = inex_im; 136 } 137 else /* nbrnd == 2 */ 138 if (out == 1) 139 params->P[offset].mpcc_inex = MPC_INEX (inex_re, inex_im); 140 else /* out == 2 */ 141 params->P[offset].mpcc_inex 142 = MPC_INEX12 (params->P[offset].mpcc_inex, 143 MPC_INEX (inex_re, inex_im)); 144 145 rnd_index++; 146 } 147 else 148 /* double rounding occurs */ 149 return 1; 150 } 151 else if (params->T[out] == MPFR) 152 { 153 MPC_ASSERT (params->T[0] == MPFR_INEX); 154 MPC_ASSERT (params->T[offset] == MPFR_INEX); 155 MPC_ASSERT (params->T[out + offset] == MPFR); 156 MPC_ASSERT (params->T[rnd_index] == MPFR_RND); 157 158 if (double_rounding_mpfr (params->P[out + offset].mpfr_data.mpfr, 159 params->P[out].mpfr, 160 params->P[0].mpfr_inex, 161 params->P[rnd_index].mpfr_rnd)) 162 /* the hight-precision value and the exact value round to the same 163 low-precision value */ 164 { 165 int inex; 166 inex = mpfr_set (params->P[out + offset].mpfr_data.mpfr, 167 params->P[out].mpfr, 168 params->P[rnd_index].mpfr_rnd); 169 params->P[out + offset].mpfr_data.known_sign = -1; 170 171 /* no double rounding means that the ternary value may comes from 172 the high-precision calculation or from the rounding */ 173 if (inex == 0) 174 params->P[offset].mpfr_inex = params->P[0].mpfr_inex; 175 else 176 params->P[offset].mpfr_inex = inex; 177 178 rnd_index++; 179 } 180 else 181 /* double rounding occurs */ 182 return 1; 183 } 184 } 185 return 0; 186 } 187