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
double_rounding_mpfr(mpfr_ptr lowprec,mpfr_ptr hiprec,int hiprec_inex,mpfr_rnd_t hiprec_rnd)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
double_rounding_mpc(mpc_ptr lowprec,mpc_ptr hiprec,int hiprec_inex,mpc_rnd_t hiprec_rnd)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
double_rounding(mpc_fun_param_t * params)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