1 /* Test file for mpfr_sub1sp.
2
3 Copyright 2003-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #include "mpfr-test.h"
24
25 static void check_special (void);
26 static void check_random (mpfr_prec_t p);
27 static void check_underflow (mpfr_prec_t p);
28 static void check_corner (mpfr_prec_t p);
29
30 static void
31 bug20170109 (void)
32 {
33 mpfr_t a, b, c;
34
35 mpfr_init2 (a, 111);
36 mpfr_init2 (b, 111);
37 mpfr_init2 (c, 111);
38 mpfr_set_str_binary (b, "0.110010010000111111011010101000100010000101101000110000100011010011000100110001100110001010001011100000001101110E1");
39 mpfr_set_str_binary (c, "0.111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111E-63");
40 mpfr_sub (a, b, c, MPFR_RNDN);
41 mpfr_set_str_binary (b, "0.110010010000111111011010101000100010000101101000110000100011001111000100110001100110001010001011100000001101110E1");
42 MPFR_ASSERTN(mpfr_equal_p (a, b));
43 mpfr_clear (a);
44 mpfr_clear (b);
45 mpfr_clear (c);
46 }
47
48 /* check mpfr_sub1sp1 when:
49 (1) p = GMP_NUMB_BITS-1, d = GMP_NUMB_BITS and bp[0] = MPFR_LIMB_HIGHBIT
50 (2) p = 2*GMP_NUMB_BITS-1, d = 2*GMP_NUMB_BITS and b = 1000...000
51 (3) p = 3*GMP_NUMB_BITS-1, d = 3*GMP_NUMB_BITS and b = 1000...000
52 */
53 static void
54 test20170208 (void)
55 {
56 mpfr_t a, b, c;
57 int inex;
58
59 mpfr_inits2 (GMP_NUMB_BITS - 1, a, b, c, (mpfr_ptr) 0);
60
61 /* test (1) */
62 mpfr_set_ui_2exp (b, 1, GMP_NUMB_BITS, MPFR_RNDN);
63 mpfr_set_ui_2exp (c, 1, 0, MPFR_RNDN);
64 inex = mpfr_sub (a, b, c, MPFR_RNDN);
65 /* b-c = 2^GMP_NUMB_BITS-1 which has GMP_NUMB_BITS bits, thus we should
66 round to 2^GMP_NUMB_BITS (even rule) */
67 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, GMP_NUMB_BITS) == 0);
68 MPFR_ASSERTN(inex > 0);
69 inex = mpfr_sub1sp (a, b, c, MPFR_RNDN);
70 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, GMP_NUMB_BITS) == 0);
71 MPFR_ASSERTN(inex > 0);
72
73 mpfr_set_ui_2exp (c, 2, 0, MPFR_RNDN);
74 mpfr_nextbelow (c);
75 /* now c = 2 - 2^(1-GMP_NUMB_BITS) */
76 inex = mpfr_sub (a, b, c, MPFR_RNDN);
77 /* b-c = 2^GMP_NUMB_BITS-2+2^(1-GMP_NUMB_BITS), which should
78 round to 2^GMP_NUMB_BITS-2. We check by directly inspecting the bit
79 field of a, since mpfr_cmp_ui might not work if unsigned long is shorter
80 than mp_limb_t, and we don't want to use mpfr_add_ui or mpfr_sub_ui
81 to construct the expected result. */
82 MPFR_ASSERTN(MPFR_MANT(a)[0] == (mp_limb_t) -2);
83 MPFR_ASSERTN(MPFR_EXP(a) == GMP_NUMB_BITS);
84 MPFR_ASSERTN(inex < 0);
85 inex = mpfr_sub1sp (a, b, c, MPFR_RNDN);
86 MPFR_ASSERTN(MPFR_MANT(a)[0] == (mp_limb_t) -2);
87 MPFR_ASSERTN(MPFR_EXP(a) == GMP_NUMB_BITS);
88 MPFR_ASSERTN(inex < 0);
89
90 /* test (2) */
91 mpfr_set_prec (a, 2 * GMP_NUMB_BITS - 1);
92 mpfr_set_prec (b, 2 * GMP_NUMB_BITS - 1);
93 mpfr_set_prec (c, 2 * GMP_NUMB_BITS - 1);
94 mpfr_set_ui_2exp (b, 1, 2 * GMP_NUMB_BITS, MPFR_RNDN);
95 mpfr_set_ui_2exp (c, 1, 0, MPFR_RNDN);
96 inex = mpfr_sub (a, b, c, MPFR_RNDN);
97 /* b-c = 2^(2*GMP_NUMB_BITS)-1 which has 2*GMP_NUMB_BITS bits, thus we should
98 round to 2^(2*GMP_NUMB_BITS) (even rule) */
99 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, 2 * GMP_NUMB_BITS) == 0);
100 MPFR_ASSERTN(inex > 0);
101 inex = mpfr_sub1sp (a, b, c, MPFR_RNDN);
102 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, 2 * GMP_NUMB_BITS) == 0);
103 MPFR_ASSERTN(inex > 0);
104
105 mpfr_set_ui_2exp (c, 2, 0, MPFR_RNDN);
106 mpfr_nextbelow (c);
107 /* now c = 2 - 2^(1-2*GMP_NUMB_BITS) */
108 inex = mpfr_sub (a, b, c, MPFR_RNDN);
109 /* b-c = 2^(2*GMP_NUMB_BITS)-2+2^(1-2*GMP_NUMB_BITS), which should
110 round to 2^(2*GMP_NUMB_BITS)-2. We check by directly inspecting the bit
111 field of a, since mpfr_cmp_ui might not work if unsigned long is shorter
112 than mp_limb_t, and we don't want to use mpfr_add_ui or mpfr_sub_ui
113 to construct the expected result. */
114 MPFR_ASSERTN(MPFR_MANT(a)[1] == (mp_limb_t) -1);
115 MPFR_ASSERTN(MPFR_MANT(a)[0] == (mp_limb_t) -2);
116 MPFR_ASSERTN(MPFR_EXP(a) == 2 * GMP_NUMB_BITS);
117 MPFR_ASSERTN(inex < 0);
118 inex = mpfr_sub1sp (a, b, c, MPFR_RNDN);
119 MPFR_ASSERTN(MPFR_MANT(a)[1] == (mp_limb_t) -1);
120 MPFR_ASSERTN(MPFR_MANT(a)[0] == (mp_limb_t) -2);
121 MPFR_ASSERTN(MPFR_EXP(a) == 2 * GMP_NUMB_BITS);
122 MPFR_ASSERTN(inex < 0);
123
124 /* test (3) */
125 mpfr_set_prec (a, 3 * GMP_NUMB_BITS - 1);
126 mpfr_set_prec (b, 3 * GMP_NUMB_BITS - 1);
127 mpfr_set_prec (c, 3 * GMP_NUMB_BITS - 1);
128 mpfr_set_ui_2exp (b, 1, 3 * GMP_NUMB_BITS, MPFR_RNDN);
129 mpfr_set_ui_2exp (c, 1, 0, MPFR_RNDN);
130 inex = mpfr_sub (a, b, c, MPFR_RNDN);
131 /* b-c = 2^(3*GMP_NUMB_BITS)-1 which has 3*GMP_NUMB_BITS bits, thus we should
132 round to 3^(2*GMP_NUMB_BITS) (even rule) */
133 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, 3 * GMP_NUMB_BITS) == 0);
134 MPFR_ASSERTN(inex > 0);
135 inex = mpfr_sub1sp (a, b, c, MPFR_RNDN);
136 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, 3 * GMP_NUMB_BITS) == 0);
137 MPFR_ASSERTN(inex > 0);
138
139 mpfr_set_ui_2exp (c, 2, 0, MPFR_RNDN);
140 mpfr_nextbelow (c);
141 /* now c = 2 - 2^(1-3*GMP_NUMB_BITS) */
142 inex = mpfr_sub (a, b, c, MPFR_RNDN);
143 /* b-c = 2^(3*GMP_NUMB_BITS)-2+2^(1-3*GMP_NUMB_BITS), which should
144 round to 2^(3*GMP_NUMB_BITS)-2. We check by directly inspecting the bit
145 field of a, since mpfr_cmp_ui might not work if unsigned long is shorter
146 than mp_limb_t, and we don't want to use mpfr_add_ui or mpfr_sub_ui
147 to construct the expected result. */
148 MPFR_ASSERTN(MPFR_MANT(a)[2] == (mp_limb_t) -1);
149 MPFR_ASSERTN(MPFR_MANT(a)[1] == (mp_limb_t) -1);
150 MPFR_ASSERTN(MPFR_MANT(a)[0] == (mp_limb_t) -2);
151 MPFR_ASSERTN(MPFR_EXP(a) == 3 * GMP_NUMB_BITS);
152 MPFR_ASSERTN(inex < 0);
153 inex = mpfr_sub1sp (a, b, c, MPFR_RNDN);
154 MPFR_ASSERTN(MPFR_MANT(a)[2] == (mp_limb_t) -1);
155 MPFR_ASSERTN(MPFR_MANT(a)[1] == (mp_limb_t) -1);
156 MPFR_ASSERTN(MPFR_MANT(a)[0] == (mp_limb_t) -2);
157 MPFR_ASSERTN(MPFR_EXP(a) == 3 * GMP_NUMB_BITS);
158 MPFR_ASSERTN(inex < 0);
159
160 mpfr_clears (a, b, c, (mpfr_ptr) 0);
161 }
162
163 static void
164 compare_sub_sub1sp (void)
165 {
166 mpfr_t a, b, c, a_ref;
167 mpfr_prec_t p;
168 unsigned long d;
169 int i, inex_ref, inex;
170 int r;
171
172 for (p = 1; p <= 3*GMP_NUMB_BITS; p++)
173 {
174 mpfr_inits2 (p, a, b, c, a_ref, (mpfr_ptr) 0);
175 for (d = 0; d <= p + 2; d++)
176 {
177 /* EXP(b) - EXP(c) = d */
178 for (i = 0; i < 4; i++)
179 {
180 /* for i even, b is the smallest number, for b odd the largest */
181 mpfr_set_ui_2exp (b, 1, d, MPFR_RNDN);
182 if (i & 1)
183 {
184 mpfr_mul_2ui (b, b, 1, MPFR_RNDN);
185 mpfr_nextbelow (b);
186 }
187 mpfr_set_ui_2exp (c, 1, 0, MPFR_RNDN);
188 if (i & 2)
189 {
190 mpfr_mul_2ui (c, c, 1, MPFR_RNDN);
191 mpfr_nextbelow (c);
192 }
193 RND_LOOP_NO_RNDF (r)
194 {
195 /* increase the precision of b to ensure sub1sp is not used */
196 mpfr_prec_round (b, p + 1, MPFR_RNDN);
197 inex_ref = mpfr_sub (a_ref, b, c, (mpfr_rnd_t) r);
198 inex = mpfr_prec_round (b, p, MPFR_RNDN);
199 MPFR_ASSERTN(inex == 0);
200 inex = mpfr_sub1sp (a, b, c, (mpfr_rnd_t) r);
201 if (inex != inex_ref)
202 {
203 printf ("mpfr_sub and mpfr_sub1sp differ for r=%s\n",
204 mpfr_print_rnd_mode ((mpfr_rnd_t) r));
205 printf ("b="); mpfr_dump (b);
206 printf ("c="); mpfr_dump (c);
207 printf ("expected inex=%d and ", inex_ref);
208 mpfr_dump (a_ref);
209 printf ("got inex=%d and ", inex);
210 mpfr_dump (a);
211 exit (1);
212 }
213 MPFR_ASSERTN(mpfr_equal_p (a, a_ref));
214 }
215 }
216 }
217 mpfr_clears (a, b, c, a_ref, (mpfr_ptr) 0);
218 }
219 }
220
221 static void
222 bug20171213 (void)
223 {
224 mpfr_t a, b, c;
225
226 mpfr_init2 (a, 127);
227 mpfr_init2 (b, 127);
228 mpfr_init2 (c, 127);
229 mpfr_set_str_binary (b, "0.1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E1");
230 mpfr_set_str_binary (c, "0.1000011010111101100101100110101111111001011001010000110000000000000000000000000000000000000000000000000000000000000000000000000E-74");
231 mpfr_sub (a, b, c, MPFR_RNDN);
232 mpfr_set_str_binary (b, "0.1111111111111111111111111111111111111111111111111111111111111111111111111101111001010000100110100110010100000001101001101011110E0");
233 MPFR_ASSERTN(mpfr_equal_p (a, b));
234 mpfr_clear (a);
235 mpfr_clear (b);
236 mpfr_clear (c);
237 }
238
239 /* generic test for bug20171213:
240 b = 1.0 with precision p
241 c = 0.1xxx110...0E-e with precision p, with e >= 1, such that the part 1xxx1 has
242 exactly p+1-e bits, thus b-c = 0.111..01... is exact on p+1 bits.
243 Due to the round-to-even rule, b-c should be rounded to 0.111..0.
244 */
245 static void
246 bug20171213_gen (mpfr_prec_t pmax)
247 {
248 mpfr_prec_t p;
249 mpfr_exp_t e;
250 mpfr_t a, b, c, d;
251
252 for (p = MPFR_PREC_MIN; p <= pmax; p++)
253 {
254 for (e = 1; e < p; e++)
255 {
256 mpfr_init2 (a, p);
257 mpfr_init2 (b, p);
258 mpfr_init2 (c, p);
259 mpfr_init2 (d, p);
260 mpfr_set_ui (b, 1, MPFR_RNDN);
261 mpfr_set_ui_2exp (c, 1, p + 1 - e, MPFR_RNDN); /* c = 2^(p + 1 - e) */
262 mpfr_sub_ui (c, c, 1, MPFR_RNDN); /* c = 2^(p + 1 - e) - 1 */
263 mpfr_div_2ui (c, c, p + 1, MPFR_RNDN); /* c = 2^(-e) - 2^(-p-1) */
264 /* the exact difference is 1 - 2^(-e) + 2^(-p-1) */
265 mpfr_sub (a, b, c, MPFR_RNDN);
266 /* check that a = 1 - 2^(-e) */
267 mpfr_set_ui_2exp (d, 1, e, MPFR_RNDN); /* b = 2^e */
268 mpfr_sub_ui (d, d, 1, MPFR_RNDN); /* b = 2^e - 1 */
269 mpfr_div_2ui (d, d, e, MPFR_RNDN); /* b = 1 - 2^(-e) */
270 if (! mpfr_equal_p (a, d))
271 {
272 printf ("bug20171213_gen failed for p=%ld, e=%ld\n",
273 (long) p, (long) e);
274 printf ("b="); mpfr_dump (b);
275 printf ("c="); mpfr_dump (c);
276 printf ("got a="); mpfr_dump (a);
277 printf ("expected d="); mpfr_dump (d);
278 exit (1);
279 }
280 mpfr_clear (a);
281 mpfr_clear (b);
282 mpfr_clear (c);
283 mpfr_clear (d);
284 }
285 }
286 }
287
288 static void
289 coverage (void)
290 {
291 mpfr_t a, b, c, d, u;
292 int inex;
293
294 /* coverage test in mpfr_sub1sp: case d=1, limb > MPFR_LIMB_HIGHBIT, RNDF
295 and also RNDZ */
296 mpfr_init2 (a, 3 * GMP_NUMB_BITS);
297 mpfr_init2 (b, 3 * GMP_NUMB_BITS);
298 mpfr_init2 (c, 3 * GMP_NUMB_BITS);
299 mpfr_init2 (d, 3 * GMP_NUMB_BITS);
300 mpfr_init2 (u, 3 * GMP_NUMB_BITS);
301 mpfr_set_ui (b, 1, MPFR_RNDN);
302 mpfr_nextbelow (b); /* b = 1 - 2^(-p) */
303 mpfr_set_prec (c, GMP_NUMB_BITS);
304 mpfr_set_ui_2exp (c, 1, -1, MPFR_RNDN);
305 mpfr_nextbelow (c);
306 mpfr_nextbelow (c); /* c = 1/2 - 2*2^(-GMP_NUMB_BITS-1) */
307 mpfr_prec_round (c, 3 * GMP_NUMB_BITS, MPFR_RNDN);
308 mpfr_nextbelow (c); /* c = 1/2 - 2*2^(-GMP_NUMB_BITS-1) - 2^(-p-1) */
309 /* b-c = c */
310 mpfr_sub (a, b, c, MPFR_RNDF);
311 mpfr_sub (d, b, c, MPFR_RNDD);
312 mpfr_sub (u, b, c, MPFR_RNDU);
313 /* check a = d or u */
314 MPFR_ASSERTN(mpfr_equal_p (a, d) || mpfr_equal_p (a, u));
315
316 /* coverage test in mpfr_sub1sp: case d=p, RNDN, sb = 0, significand of b
317 is even but b<>2^e, (case 1e) */
318 mpfr_set_prec (a, 3 * GMP_NUMB_BITS);
319 mpfr_set_prec (b, 3 * GMP_NUMB_BITS);
320 mpfr_set_prec (c, 3 * GMP_NUMB_BITS);
321 mpfr_set_ui (b, 1, MPFR_RNDN);
322 mpfr_nextabove (b);
323 mpfr_nextabove (b);
324 mpfr_set_ui_2exp (c, 1, -3 * GMP_NUMB_BITS, MPFR_RNDN);
325 inex = mpfr_sub (a, b, c, MPFR_RNDN);
326 MPFR_ASSERTN(inex > 0);
327 MPFR_ASSERTN(mpfr_equal_p (a, b));
328
329 mpfr_clear (a);
330 mpfr_clear (b);
331 mpfr_clear (c);
332 mpfr_clear (d);
333 mpfr_clear (u);
334 }
335
336 /* bug in mpfr_sub1sp1n, made generic */
337 static void
338 bug20180217 (mpfr_prec_t pmax)
339 {
340 mpfr_t a, b, c;
341 int inex;
342 mpfr_prec_t p;
343
344 for (p = MPFR_PREC_MIN; p <= pmax; p++)
345 {
346 mpfr_init2 (a, p);
347 mpfr_init2 (b, p);
348 mpfr_init2 (c, p);
349 mpfr_set_ui (b, 1, MPFR_RNDN); /* b = 1 */
350 mpfr_set_ui_2exp (c, 1, -p-1, MPFR_RNDN); /* c = 2^(-p-1) */
351 /* a - b = 1 - 2^(-p-1) and should be rounded to 1 (case 2f of
352 mpfr_sub1sp) */
353 inex = mpfr_sub (a, b, c, MPFR_RNDN);
354 MPFR_ASSERTN(inex > 0);
355 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
356 /* check also when a=b */
357 mpfr_set_ui (a, 1, MPFR_RNDN);
358 inex = mpfr_sub (a, a, c, MPFR_RNDN);
359 MPFR_ASSERTN(inex > 0);
360 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
361 /* and when a=c */
362 mpfr_set_ui (b, 1, MPFR_RNDN); /* b = 1 */
363 mpfr_set_ui_2exp (a, 1, -p-1, MPFR_RNDN);
364 inex = mpfr_sub (a, b, a, MPFR_RNDN);
365 MPFR_ASSERTN(inex > 0);
366 MPFR_ASSERTN(mpfr_cmp_ui (a, 1) == 0);
367 mpfr_clear (a);
368 mpfr_clear (b);
369 mpfr_clear (c);
370 }
371 }
372
373 /* bug in revision 12985 with tlog and GMP_CHECK_RANDOMIZE=1534111552615050
374 (introduced in revision 12242, does not affect the 4.0 branch) */
375 static void
376 bug20180813 (void)
377 {
378 mpfr_t a, b, c;
379
380 mpfr_init2 (a, 194);
381 mpfr_init2 (b, 194);
382 mpfr_init2 (c, 194);
383 mpfr_set_str_binary (b, "0.10000111101000100000010000100010110111011100110100000101100111000010101000110110010101011101101011110110001000111001000010110010111010010100011011010100001010001110000101000010101110100110001000E7");
384 mpfr_set_str_binary (c, "0.10000000000000000100001111010001000000100001000101101110111001101000001011001110000101010001101100101010111011010111101100010001110010000101100101110100101000110110101000010100011100001010000101E24");
385 mpfr_sub (a, b, c, MPFR_RNDN);
386 mpfr_set_str_binary (b, "-0.11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111E23");
387 MPFR_ASSERTN(mpfr_equal_p (a, b));
388 mpfr_clear (a);
389 mpfr_clear (b);
390 mpfr_clear (c);
391 }
392
393 /* bug in revision 13599 with tatan and GMP_CHECK_RANDOMIZE=1567609230659336:
394 the values are equal, but the ternary value differs between sub1 and sub1sp
395 (bug introduced with mpfr_sub1sp2n, does not affect the 4.0 branch) */
396 static void
397 bug20190904 (void)
398 {
399 mpfr_t a, b, c;
400 int ret;
401
402 mpfr_init2 (a, 128);
403 mpfr_init2 (b, 128);
404 mpfr_init2 (c, 128);
405 mpfr_set_str_binary (b, "0.11001001000011111101101010100010001000010110100011000010001101001100010011000110011000101000101110000000110111000001110011010001E1");
406 mpfr_set_str_binary (c, "0.10010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010010000000000000000000000E-102");
407 ret = mpfr_sub (a, b, c, MPFR_RNDN);
408 mpfr_set_str_binary (b, "0.11001001000011111101101010100010001000010110100011000010001101001100010011000110011000101000101101111111101111000001110011010001E1");
409 MPFR_ASSERTN(mpfr_equal_p (a, b));
410 MPFR_ASSERTN(ret > 0);
411 mpfr_clear (a);
412 mpfr_clear (b);
413 mpfr_clear (c);
414 }
415
416 int
417 main (void)
418 {
419 mpfr_prec_t p;
420
421 tests_start_mpfr ();
422
423 bug20190904 ();
424 bug20180813 ();
425 bug20180217 (1024);
426 coverage ();
427 compare_sub_sub1sp ();
428 test20170208 ();
429 bug20170109 ();
430 bug20171213 ();
431 bug20171213_gen (256);
432 check_special ();
433 for (p = MPFR_PREC_MIN ; p < 200 ; p++)
434 {
435 check_underflow (p);
436 check_random (p);
437 check_corner (p);
438 }
439
440 tests_end_mpfr ();
441 return 0;
442 }
443
444 #define STD_ERROR \
445 do \
446 { \
447 printf("ERROR: for %s and p=%lu and i=%d:\nY=", \
448 mpfr_print_rnd_mode ((mpfr_rnd_t) r), (unsigned long) p, i); \
449 mpfr_dump (y); \
450 printf ("Z="); mpfr_dump (z); \
451 printf ("Expected: "); mpfr_dump (x2); \
452 printf ("Got : "); mpfr_dump (x); \
453 abort(); \
454 } \
455 while (0)
456
457 #define STD_ERROR2 \
458 do \
459 { \
460 printf("ERROR: for %s and p=%lu and i=%d:\nY=", \
461 mpfr_print_rnd_mode ((mpfr_rnd_t) r), (unsigned long) p, i); \
462 mpfr_dump (y); \
463 printf ("Z="); mpfr_dump (z); \
464 printf ("Expected: "); mpfr_dump (x2); \
465 printf ("Got : "); mpfr_dump (x); \
466 printf ("Wrong inexact flag. Expected %d. Got %d\n", \
467 inexact1, inexact2); \
468 exit(1); \
469 } \
470 while (0)
471
472 static void
473 check_random (mpfr_prec_t p)
474 {
475 mpfr_t x,y,z,x2;
476 int r;
477 int i, inexact1, inexact2;
478
479 mpfr_inits2 (p, x, y, z, x2, (mpfr_ptr) 0);
480
481 for (i = 0 ; i < 500 ; i++)
482 {
483 mpfr_urandomb (y, RANDS);
484 mpfr_urandomb (z, RANDS);
485 if (MPFR_IS_PURE_FP(y) && MPFR_IS_PURE_FP(z))
486 RND_LOOP (r)
487 {
488 if (r == MPFR_RNDF)
489 continue; /* mpfr_sub1 and mpfr_sub1sp could differ,
490 and inexact makes no sense */
491 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
492 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
493 if (mpfr_cmp(x, x2))
494 STD_ERROR;
495 if (inexact1 != inexact2)
496 STD_ERROR2;
497 }
498 }
499
500 mpfr_clears (x, y, z, x2, (mpfr_ptr) 0);
501 }
502
503 static void
504 check_special (void)
505 {
506 mpfr_t x,y,z,x2;
507 int r;
508 mpfr_prec_t p;
509 int i = -1, inexact1, inexact2;
510 mpfr_exp_t es;
511
512 mpfr_inits (x, y, z, x2, (mpfr_ptr) 0);
513
514 RND_LOOP (r)
515 {
516 if (r == MPFR_RNDF)
517 continue; /* comparison between sub1 and sub1sp makes no sense here */
518
519 p = 53;
520 mpfr_set_prec(x, 53);
521 mpfr_set_prec(x2, 53);
522 mpfr_set_prec(y, 53);
523 mpfr_set_prec(z, 53);
524
525 mpfr_set_str_binary (y,
526 "0.10110111101101110010010010011011000001101101011011001E31");
527
528 mpfr_sub1sp (x, y, y, (mpfr_rnd_t) r);
529 if (mpfr_cmp_ui(x, 0))
530 {
531 printf("Error for x-x with p=%lu. Expected 0. Got:",
532 (unsigned long) p);
533 mpfr_dump (x);
534 exit(1);
535 }
536
537 mpfr_set(z, y, (mpfr_rnd_t) r);
538 mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
539 if (mpfr_cmp_ui(x, 0))
540 {
541 printf("Error for x-y with y=x and p=%lu. Expected 0. Got:",
542 (unsigned long) p);
543 mpfr_dump (x);
544 exit(1);
545 }
546 /* diff = 0 */
547 mpfr_set_str_binary (y,
548 "0.10110111101101110010010010011011001001101101011011001E31");
549 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
550 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
551 if (mpfr_cmp(x, x2))
552 STD_ERROR;
553 if (inexact1 != inexact2)
554 STD_ERROR2;
555
556 /* Diff = 1 */
557 mpfr_set_str_binary (y,
558 "0.10110111101101110010010010011011000001101101011011001E30");
559 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
560 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
561 if (mpfr_cmp(x, x2))
562 STD_ERROR;
563 if (inexact1 != inexact2)
564 STD_ERROR2;
565
566 /* Diff = 2 */
567 mpfr_set_str_binary (y,
568 "0.10110111101101110010010010011011000101101101011011001E32");
569 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
570 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
571 if (mpfr_cmp(x, x2))
572 STD_ERROR;
573 if (inexact1 != inexact2)
574 STD_ERROR2;
575
576 /* Diff = 32 */
577 mpfr_set_str_binary (y,
578 "0.10110111101101110010010010011011000001101101011011001E63");
579 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
580 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
581 if (mpfr_cmp(x, x2))
582 STD_ERROR;
583 if (inexact1 != inexact2)
584 STD_ERROR2;
585
586 /* Diff = 52 */
587 mpfr_set_str_binary (y,
588 "0.10110111101101110010010010011011010001101101011011001E83");
589 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
590 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
591 if (mpfr_cmp(x, x2))
592 STD_ERROR;
593 if (inexact1 != inexact2)
594 STD_ERROR2;
595
596 /* Diff = 53 */
597 mpfr_set_str_binary (y,
598 "0.10110111101101110010010010011111000001101101011011001E31");
599 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
600 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
601 if (mpfr_cmp(x, x2))
602 STD_ERROR;
603 if (inexact1 != inexact2)
604 STD_ERROR2;
605
606 /* Diff > 200 */
607 mpfr_set_str_binary (y,
608 "0.10110111101101110010010010011011000001101101011011001E331");
609 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
610 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
611 if (mpfr_cmp(x, x2))
612 STD_ERROR;
613 if (inexact1 != inexact2)
614 STD_ERROR2;
615
616 mpfr_set_str_binary (y,
617 "0.10000000000000000000000000000000000000000000000000000E31");
618 mpfr_set_str_binary (z,
619 "0.11111111111111111111111111111111111111111111111111111E30");
620 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
621 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
622 if (mpfr_cmp(x, x2))
623 STD_ERROR;
624 if (inexact1 != inexact2)
625 STD_ERROR2;
626
627 mpfr_set_str_binary (y,
628 "0.10000000000000000000000000000000000000000000000000000E31");
629 mpfr_set_str_binary (z,
630 "0.11111111111111111111111111111111111111111111111111111E29");
631 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
632 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
633 if (mpfr_cmp(x, x2))
634 STD_ERROR;
635 if (inexact1 != inexact2)
636 STD_ERROR2;
637
638 mpfr_set_str_binary (y,
639 "0.10000000000000000000000000000000000000000000000000000E52");
640 mpfr_set_str_binary (z,
641 "0.10000000000010000000000000000000000000000000000000000E00");
642 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
643 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
644 if (mpfr_cmp(x, x2))
645 STD_ERROR;
646 if (inexact1 != inexact2)
647 STD_ERROR2;
648
649 mpfr_set_str_binary (y,
650 "0.11100000000000000000000000000000000000000000000000000E53");
651 mpfr_set_str_binary (z,
652 "0.10000000000000000000000000000000000000000000000000000E00");
653 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
654 inexact2 = mpfr_sub1sp(z, y, z, (mpfr_rnd_t) r);
655 mpfr_set(x, z, (mpfr_rnd_t) r);
656 if (mpfr_cmp(x, x2))
657 STD_ERROR;
658 if (inexact1 != inexact2)
659 STD_ERROR2;
660
661 mpfr_set_str_binary (y,
662 "0.10000000000000000000000000000000000000000000000000000E53");
663 mpfr_set_str_binary (z,
664 "0.10100000000000000000000000000000000000000000000000000E00");
665 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
666 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
667 if (mpfr_cmp(x, x2))
668 STD_ERROR;
669 if (inexact1 != inexact2)
670 STD_ERROR2;
671
672 mpfr_set_str_binary (y,
673 "0.10000000000000000000000000000000000000000000000000000E54");
674 mpfr_set_str_binary (z,
675 "0.10100000000000000000000000000000000000000000000000000E00");
676 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
677 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
678 if (mpfr_cmp(x, x2))
679 STD_ERROR;
680 if (inexact1 != inexact2)
681 STD_ERROR2;
682
683 p = 63;
684 mpfr_set_prec(x, p);
685 mpfr_set_prec(x2, p);
686 mpfr_set_prec(y, p);
687 mpfr_set_prec(z, p);
688 mpfr_set_str_binary (y,
689 "0.100000000000000000000000000000000000000000000000000000000000000E62");
690 mpfr_set_str_binary (z,
691 "0.110000000000000000000000000000000000000000000000000000000000000E00");
692 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
693 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
694 if (mpfr_cmp(x, x2))
695 STD_ERROR;
696 if (inexact1 != inexact2)
697 STD_ERROR2;
698
699 p = 64;
700 mpfr_set_prec(x, 64);
701 mpfr_set_prec(x2, 64);
702 mpfr_set_prec(y, 64);
703 mpfr_set_prec(z, 64);
704
705 mpfr_set_str_binary (y,
706 "0.1100000000000000000000000000000000000000000000000000000000000000E31");
707 mpfr_set_str_binary (z,
708 "0.1111111111111111111111111110000000000000000000000000011111111111E29");
709 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
710 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
711 if (mpfr_cmp(x, x2))
712 STD_ERROR;
713 if (inexact1 != inexact2)
714 STD_ERROR2;
715
716 mpfr_set_str_binary (y,
717 "0.1000000000000000000000000000000000000000000000000000000000000000E63");
718 mpfr_set_str_binary (z,
719 "0.1011000000000000000000000000000000000000000000000000000000000000E00");
720 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
721 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
722 if (mpfr_cmp(x, x2))
723 STD_ERROR;
724 if (inexact1 != inexact2)
725 STD_ERROR2;
726
727 mpfr_set_str_binary (y,
728 "0.1000000000000000000000000000000000000000000000000000000000000000E63");
729 mpfr_set_str_binary (z,
730 "0.1110000000000000000000000000000000000000000000000000000000000000E00");
731 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
732 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
733 if (mpfr_cmp(x, x2))
734 STD_ERROR;
735 if (inexact1 != inexact2)
736 STD_ERROR2;
737
738 mpfr_set_str_binary (y,
739 "0.10000000000000000000000000000000000000000000000000000000000000E63");
740 mpfr_set_str_binary (z,
741 "0.10000000000000000000000000000000000000000000000000000000000000E00");
742 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
743 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
744 if (mpfr_cmp(x, x2))
745 STD_ERROR;
746 if (inexact1 != inexact2)
747 STD_ERROR2;
748
749 mpfr_set_str_binary (y,
750 "0.1000000000000000000000000000000000000000000000000000000000000000E64");
751 mpfr_set_str_binary (z,
752 "0.1010000000000000000000000000000000000000000000000000000000000000E00");
753 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
754 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
755 if (mpfr_cmp(x, x2))
756 STD_ERROR;
757 if (inexact1 != inexact2)
758 STD_ERROR2;
759
760 MPFR_SET_NAN(x);
761 MPFR_SET_NAN(x2);
762 mpfr_set_str_binary (y,
763 "0.1000000000000000000000000000000000000000000000000000000000000000"
764 "E-1073741823");
765 mpfr_set_str_binary (z,
766 "0.1100000000000000000000000000000000000000000000000000000000000000"
767 "E-1073741823");
768 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
769 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
770 if (mpfr_cmp(x, x2))
771 STD_ERROR;
772 if (inexact1 != inexact2)
773 STD_ERROR2;
774
775 p = 9;
776 mpfr_set_prec(x, p);
777 mpfr_set_prec(x2, p);
778 mpfr_set_prec(y, p);
779 mpfr_set_prec(z, p);
780
781 mpfr_set_str_binary (y, "0.100000000E1");
782 mpfr_set_str_binary (z, "0.100000000E-8");
783 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
784 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
785 if (mpfr_cmp(x, x2))
786 STD_ERROR;
787 if (inexact1 != inexact2)
788 STD_ERROR2;
789
790 p = 34;
791 mpfr_set_prec(x, p);
792 mpfr_set_prec(x2, p);
793 mpfr_set_prec(y, p);
794 mpfr_set_prec(z, p);
795
796 mpfr_set_str_binary (y, "-0.1011110000111100010111011100110100E-18");
797 mpfr_set_str_binary (z, "0.1000101010110011010101011110000000E-14");
798 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
799 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
800 if (mpfr_cmp(x, x2))
801 STD_ERROR;
802 if (inexact1 != inexact2)
803 STD_ERROR2;
804
805 p = 124;
806 mpfr_set_prec(x, p);
807 mpfr_set_prec(x2, p);
808 mpfr_set_prec(y, p);
809 mpfr_set_prec(z, p);
810
811 mpfr_set_str_binary (y,
812 "0.1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E1");
813 mpfr_set_str_binary (z,
814 "0.1011111000100111000011001000011101010101101100101010101001000001110100001101110110001110111010000011101001100010111110001100E-31");
815 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
816 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
817 if (mpfr_cmp(x, x2))
818 STD_ERROR;
819 if (inexact1 != inexact2)
820 STD_ERROR2;
821
822 p = 288;
823 mpfr_set_prec(x, p);
824 mpfr_set_prec(x2, p);
825 mpfr_set_prec(y, p);
826 mpfr_set_prec(z, p);
827
828 mpfr_set_str_binary (y,
829 "0.111000110011000001000111101010111011110011101001101111111110000011100101000001001010110010101010011001010100000001110011110001010101101010001011101110100100001011110100110000101101100011010001001011011010101010000010001101001000110010010111111011110001111101001000101101001100101100101000E80");
830 mpfr_set_str_binary (z,
831 "-0.100001111111101001011010001100110010100111001110000110011101001011010100001000000100111011010110110010000000000010101101011000010000110001110010100001100101011100100100001011000100011110000001010101000100011101001000010111100000111000111011001000100100011000100000010010111000000100100111E-258");
832 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
833 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
834 if (mpfr_cmp(x, x2))
835 STD_ERROR;
836 if (inexact1 != inexact2)
837 STD_ERROR2;
838
839 p = 85;
840 mpfr_set_prec(x, p);
841 mpfr_set_prec(x2, p);
842 mpfr_set_prec(y, p);
843 mpfr_set_prec(z, p);
844
845 mpfr_set_str_binary (y,
846 "0.1111101110100110110110100010101011101001100010100011110110110010010011101100101111100E-4");
847 mpfr_set_str_binary (z,
848 "0.1111101110100110110110100010101001001000011000111000011101100101110100001110101010110E-4");
849 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
850 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
851 if (mpfr_cmp(x, x2))
852 STD_ERROR;
853 if (inexact1 != inexact2)
854 STD_ERROR2;
855
856 p = 64;
857 mpfr_set_prec(x, p); mpfr_set_prec(x2, p);
858 mpfr_set_prec(y, p); mpfr_set_prec(z, p);
859
860 mpfr_set_str_binary (y,
861 "0.11000000000000000000000000000000"
862 "00000000000000000000000000000000E1");
863 mpfr_set_str_binary (z,
864 "0.10000000000000000000000000000000"
865 "00000000000000000000000000000001E0");
866 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
867 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
868 if (mpfr_cmp(x, x2))
869 STD_ERROR;
870 if (inexact1 != inexact2)
871 STD_ERROR2;
872
873 mpfr_set_str_binary (y,
874 "0.11000000000000000000000000000000"
875 "000000000000000000000000000001E1");
876 mpfr_set_str_binary (z,
877 "0.10000000000000000000000000000000"
878 "00000000000000000000000000000001E0");
879 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
880 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
881 if (mpfr_cmp(x, x2))
882 STD_ERROR;
883 if (inexact1 != inexact2)
884 STD_ERROR2;
885
886 es = mpfr_get_emin ();
887 set_emin (-1024);
888
889 mpfr_set_str_binary (y,
890 "0.10000000000000000000000000000000"
891 "000000000000000000000000000000E-1023");
892 mpfr_set_str_binary (z,
893 "0.10000000000000000000000000000000"
894 "00000000000000000000000000000001E-1023");
895 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
896 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
897 if (mpfr_cmp(x, x2))
898 STD_ERROR;
899 if (inexact1 != inexact2)
900 STD_ERROR2;
901
902 mpfr_set_str_binary (y,
903 "0.10000000000000000000000000000000"
904 "000000000000000000000000000000E-1023");
905 mpfr_set_str_binary (z,
906 "0.1000000000000000000000000000000"
907 "000000000000000000000000000000E-1023");
908 inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
909 inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
910 if (mpfr_cmp(x, x2))
911 STD_ERROR;
912 if (inexact1 != inexact2)
913 STD_ERROR2;
914
915 set_emin (es);
916 }
917
918 mpfr_clears (x, y, z, x2, (mpfr_ptr) 0);
919 }
920
921 static void
922 check_underflow (mpfr_prec_t p)
923 {
924 mpfr_t x, y, z;
925 int inexact;
926
927 mpfr_inits2 (p, x, y, z, (mpfr_ptr) 0);
928
929 if (p >= 2) /* we need p >= 2 so that 3 is exact */
930 {
931 mpfr_set_ui_2exp (y, 4, mpfr_get_emin () - 2, MPFR_RNDN);
932 mpfr_set_ui_2exp (z, 3, mpfr_get_emin () - 2, MPFR_RNDN);
933 inexact = mpfr_sub (x, y, z, MPFR_RNDN);
934 if (inexact >= 0 || (mpfr_cmp_ui (x, 0) != 0))
935 {
936 printf ("4*2^(emin-2) - 3*2^(emin-2) with RNDN failed for p=%ld\n",
937 (long) p);
938 printf ("Expected inexact < 0 with x=0\n");
939 printf ("Got inexact=%d with x=", inexact);
940 mpfr_dump (x);
941 exit (1);
942 }
943 }
944
945 if (p >= 3) /* we need p >= 3 so that 5 is exact */
946 {
947 mpfr_set_ui_2exp (y, 5, mpfr_get_emin () - 2, MPFR_RNDN);
948 mpfr_set_ui_2exp (z, 4, mpfr_get_emin () - 2, MPFR_RNDN);
949 inexact = mpfr_sub (x, y, z, MPFR_RNDN);
950 if (inexact >= 0 || (mpfr_cmp_ui (x, 0) != 0))
951 {
952 printf ("5*2^(emin-2) - 4*2^(emin-2) with RNDN failed for p=%ld\n",
953 (long) p);
954 printf ("Expected inexact < 0 with x=0\n");
955 printf ("Got inexact=%d with x=", inexact);
956 mpfr_dump (x);
957 exit (1);
958 }
959 }
960
961 mpfr_clears (x, y, z, (mpfr_ptr) 0);
962 }
963
964 /* check corner cases of mpfr_sub1sp in case d = 1 and limb = MPFR_LIMB_HIGHBIT */
965 static void
966 check_corner (mpfr_prec_t p)
967 {
968 mpfr_t x, y, z;
969 mpfr_exp_t e;
970 int inex, odd;
971
972 if (p < 4) /* ensures that the initial value of z is > 1 below */
973 return;
974
975 mpfr_inits2 (p, x, y, z, (mpfr_ptr) 0);
976 mpfr_const_pi (y, MPFR_RNDN);
977 mpfr_set_ui (z, 2, MPFR_RNDN);
978 inex = mpfr_sub (z, y, z, MPFR_RNDN); /* z is near pi-2, thus y-z is near 2 */
979 MPFR_ASSERTN(inex == 0);
980 for (e = 0; e < p; e++)
981 {
982 /* add 2^(-e) to z */
983 mpfr_mul_2ui (z, z, e, MPFR_RNDN);
984 inex = mpfr_add_ui (z, z, 1, MPFR_RNDN);
985 MPFR_ASSERTN(inex == 0);
986 mpfr_div_2ui (z, z, e, MPFR_RNDN);
987
988 /* compute x = y - z which should be exact, near 2-2^(-e) */
989 inex = mpfr_sub (x, y, z, MPFR_RNDN);
990 MPFR_ASSERTN(inex == 0);
991 MPFR_ASSERTN(mpfr_get_exp (x) == 1);
992
993 /* restore initial z */
994 mpfr_mul_2ui (z, z, e, MPFR_RNDN);
995 inex = mpfr_sub_ui (z, z, 1, MPFR_RNDN);
996 MPFR_ASSERTN(inex == 0);
997 mpfr_div_2ui (z, z, e, MPFR_RNDN);
998
999 /* subtract 2^(-e) to z */
1000 mpfr_mul_2ui (z, z, e, MPFR_RNDN);
1001 inex = mpfr_sub_ui (z, z, 1, MPFR_RNDN);
1002 MPFR_ASSERTN(inex == 0);
1003 mpfr_div_2ui (z, z, e, MPFR_RNDN);
1004
1005 /* ensure last significant bit of z is 0 so that y-z is exact */
1006 odd = mpfr_min_prec (z) == p;
1007 if (odd) /* add one ulp to z */
1008 mpfr_nextabove (z);
1009
1010 /* compute x = y - z which should be exact, near 2+2^(-e) */
1011 inex = mpfr_sub (x, y, z, MPFR_RNDN);
1012 MPFR_ASSERTN(inex == 0);
1013 MPFR_ASSERTN(mpfr_get_exp (x) == 2);
1014
1015 /* restore initial z */
1016 if (odd)
1017 mpfr_nextbelow (z);
1018 mpfr_mul_2ui (z, z, e, MPFR_RNDN);
1019 inex = mpfr_add_ui (z, z, 1, MPFR_RNDN);
1020 MPFR_ASSERTN(inex == 0);
1021 mpfr_div_2ui (z, z, e, MPFR_RNDN);
1022 }
1023 mpfr_clears (x, y, z, (mpfr_ptr) 0);
1024 }