1 /* Test file for mpfr_add and mpfr_sub.
2
3 Copyright 1999-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 #define N 30000
24
25 #include <float.h>
26
27 #include "mpfr-test.h"
28
29 /* If the precisions are the same, we want to test both mpfr_add1sp
30 and mpfr_add1. */
31
32 /* FIXME: modify check() to test the ternary value and the flags. */
33
34 static int usesp;
35
36 static int
37 test_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
38 {
39 int res;
40 #ifdef CHECK_EXTERNAL
41 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c);
42 if (ok)
43 {
44 mpfr_print_raw (b);
45 printf (" ");
46 mpfr_print_raw (c);
47 }
48 #endif
49 if (usesp || MPFR_ARE_SINGULAR(b,c) || MPFR_SIGN(b) != MPFR_SIGN(c))
50 res = mpfr_add (a, b, c, rnd_mode);
51 else
52 {
53 if (MPFR_GET_EXP(b) < MPFR_GET_EXP(c))
54 res = mpfr_add1(a, c, b, rnd_mode);
55 else
56 res = mpfr_add1(a, b, c, rnd_mode);
57 }
58 #ifdef CHECK_EXTERNAL
59 if (ok)
60 {
61 printf (" ");
62 mpfr_print_raw (a);
63 printf ("\n");
64 }
65 #endif
66 return res;
67 }
68
69 /* checks that xs+ys gives the expected result zs */
70 static void
71 check (const char *xs, const char *ys, mpfr_rnd_t rnd_mode,
72 unsigned int px, unsigned int py, unsigned int pz, const char *zs)
73 {
74 mpfr_t xx,yy,zz;
75
76 mpfr_init2 (xx, px);
77 mpfr_init2 (yy, py);
78 mpfr_init2 (zz, pz);
79
80 mpfr_set_str1 (xx, xs);
81 mpfr_set_str1 (yy, ys);
82 test_add (zz, xx, yy, rnd_mode);
83 if (mpfr_cmp_str1 (zz, zs) )
84 {
85 printf ("expected sum is %s, got ", zs);
86 mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
87 printf ("\nmpfr_add failed for x=%s y=%s with rnd_mode=%s\n",
88 xs, ys, mpfr_print_rnd_mode (rnd_mode));
89 exit (1);
90 }
91 mpfr_clears (xx, yy, zz, (mpfr_ptr) 0);
92 }
93
94 static void
95 check2b (const char *xs, int px,
96 const char *ys, int py,
97 const char *rs, int pz,
98 mpfr_rnd_t rnd_mode)
99 {
100 mpfr_t xx, yy, zz;
101
102 mpfr_init2 (xx,px);
103 mpfr_init2 (yy,py);
104 mpfr_init2 (zz,pz);
105 mpfr_set_str_binary (xx, xs);
106 mpfr_set_str_binary (yy, ys);
107 test_add (zz, xx, yy, rnd_mode);
108 if (mpfr_cmp_str (zz, rs, 2, MPFR_RNDN))
109 {
110 printf ("(2) x=%s,%d y=%s,%d pz=%d,rnd=%s\n",
111 xs, px, ys, py, pz, mpfr_print_rnd_mode (rnd_mode));
112 printf ("got "); mpfr_dump (zz);
113 mpfr_set_str(zz, rs, 2, MPFR_RNDN);
114 printf ("instead of "); mpfr_dump (zz);
115 exit (1);
116 }
117 mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
118 }
119
120 static void
121 check64 (void)
122 {
123 mpfr_t x, t, u;
124
125 mpfr_init (x);
126 mpfr_init (t);
127 mpfr_init (u);
128
129 mpfr_set_prec (x, 29);
130 mpfr_set_str_binary (x, "1.1101001000101111011010010110e-3");
131 mpfr_set_prec (t, 58);
132 mpfr_set_str_binary (t, "0.11100010011111001001100110010111110110011000000100101E-1");
133 mpfr_set_prec (u, 29);
134 test_add (u, x, t, MPFR_RNDD);
135 mpfr_set_str_binary (t, "1.0101011100001000011100111110e-1");
136 if (mpfr_cmp (u, t))
137 {
138 printf ("mpfr_add(u, x, t) failed for prec(x)=29, prec(t)=58\n");
139 printf ("expected "); mpfr_out_str (stdout, 2, 29, t, MPFR_RNDN);
140 puts ("");
141 printf ("got "); mpfr_out_str (stdout, 2, 29, u, MPFR_RNDN);
142 puts ("");
143 exit(1);
144 }
145
146 mpfr_set_prec (x, 4);
147 mpfr_set_str_binary (x, "-1.0E-2");
148 mpfr_set_prec (t, 2);
149 mpfr_set_str_binary (t, "-1.1e-2");
150 mpfr_set_prec (u, 2);
151 test_add (u, x, t, MPFR_RNDN);
152 if ((mp_limb_t) (MPFR_MANT(u)[0] << 2))
153 {
154 printf ("result not normalized for prec=2\n");
155 mpfr_dump (u);
156 exit (1);
157 }
158 mpfr_set_str_binary (t, "-1.0e-1");
159 if (mpfr_cmp (u, t))
160 {
161 printf ("mpfr_add(u, x, t) failed for prec(x)=4, prec(t)=2\n");
162 printf ("expected -1.0e-1\n");
163 printf ("got "); mpfr_out_str (stdout, 2, 4, u, MPFR_RNDN);
164 puts ("");
165 exit (1);
166 }
167
168 mpfr_set_prec (x, 8);
169 mpfr_set_str_binary (x, "-0.10011010"); /* -77/128 */
170 mpfr_set_prec (t, 4);
171 mpfr_set_str_binary (t, "-1.110e-5"); /* -7/128 */
172 mpfr_set_prec (u, 4);
173 test_add (u, x, t, MPFR_RNDN); /* should give -5/8 */
174 mpfr_set_str_binary (t, "-1.010e-1");
175 if (mpfr_cmp (u, t)) {
176 printf ("mpfr_add(u, x, t) failed for prec(x)=8, prec(t)=4\n");
177 printf ("expected -1.010e-1\n");
178 printf ("got "); mpfr_out_str (stdout, 2, 4, u, MPFR_RNDN);
179 puts ("");
180 exit (1);
181 }
182
183 mpfr_set_prec (x, 112); mpfr_set_prec (t, 98); mpfr_set_prec (u, 54);
184 mpfr_set_str_binary (x, "-0.11111100100000000011000011100000101101010001000111E-401");
185 mpfr_set_str_binary (t, "0.10110000100100000101101100011111111011101000111000101E-464");
186 test_add (u, x, t, MPFR_RNDN);
187 if (mpfr_cmp (u, x))
188 {
189 printf ("mpfr_add(u, x, t) failed for prec(x)=112, prec(t)=98\n");
190 exit (1);
191 }
192
193 mpfr_set_prec (x, 92); mpfr_set_prec (t, 86); mpfr_set_prec (u, 53);
194 mpfr_set_str (x, "-5.03525136761487735093e-74", 10, MPFR_RNDN);
195 mpfr_set_str (t, "8.51539046314262304109e-91", 10, MPFR_RNDN);
196 test_add (u, x, t, MPFR_RNDN);
197 if (mpfr_cmp_str1 (u, "-5.0352513676148773509283672e-74") )
198 {
199 printf ("mpfr_add(u, x, t) failed for prec(x)=92, prec(t)=86\n");
200 exit (1);
201 }
202
203 mpfr_set_prec(x, 53); mpfr_set_prec(t, 76); mpfr_set_prec(u, 76);
204 mpfr_set_str_binary(x, "-0.10010010001001011011110000000000001010011011011110001E-32");
205 mpfr_set_str_binary(t, "-0.1011000101110010000101111111011111010001110011110111100110101011110010011111");
206 mpfr_sub(u, x, t, MPFR_RNDU);
207 mpfr_set_str_binary(t, "0.1011000101110010000101111111011100111111101010011011110110101011101000000100");
208 if (mpfr_cmp(u,t))
209 {
210 printf ("expect "); mpfr_dump (t);
211 printf ("mpfr_add failed for precisions 53-76\n");
212 exit (1);
213 }
214 mpfr_set_prec(x, 53); mpfr_set_prec(t, 108); mpfr_set_prec(u, 108);
215 mpfr_set_str_binary(x, "-0.10010010001001011011110000000000001010011011011110001E-32");
216 mpfr_set_str_binary(t, "-0.101100010111001000010111111101111101000111001111011110011010101111001001111000111011001110011000000000111111");
217 mpfr_sub(u, x, t, MPFR_RNDU);
218 mpfr_set_str_binary(t, "0.101100010111001000010111111101110011111110101001101111011010101110100000001011000010101110011000000000111111");
219 if (mpfr_cmp(u,t))
220 {
221 printf ("expect "); mpfr_dump (t);
222 printf ("mpfr_add failed for precisions 53-108\n");
223 exit (1);
224 }
225 mpfr_set_prec(x, 97); mpfr_set_prec(t, 97); mpfr_set_prec(u, 97);
226 mpfr_set_str_binary(x, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010000000000000000E-39");
227 mpfr_set_ui(t, 1, MPFR_RNDN);
228 test_add (u, x, t, MPFR_RNDN);
229 mpfr_set_str_binary(x, "0.1000000000000000000000000000000000000000111110110000100000000101100011011110100000101111100010001E1");
230 if (mpfr_cmp(u,x))
231 {
232 printf ("mpfr_add failed for precision 97\n");
233 exit (1);
234 }
235 mpfr_set_prec(x, 128); mpfr_set_prec(t, 128); mpfr_set_prec(u, 128);
236 mpfr_set_str_binary(x, "0.10101011111001001010111011001000101100111101000000111111111011010100001100011101010001010111111101111010100110111111100101100010E-4");
237 mpfr_set(t, x, MPFR_RNDN);
238 mpfr_sub(u, x, t, MPFR_RNDN);
239 mpfr_set_prec(x, 96); mpfr_set_prec(t, 96); mpfr_set_prec(u, 96);
240 mpfr_set_str_binary(x, "0.111000000001110100111100110101101001001010010011010011100111100011010100011001010011011011000010E-4");
241 mpfr_set(t, x, MPFR_RNDN);
242 mpfr_sub(u, x, t, MPFR_RNDN);
243 mpfr_set_prec(x, 85); mpfr_set_prec(t, 85); mpfr_set_prec(u, 85);
244 mpfr_set_str_binary(x, "0.1111101110100110110110100010101011101001100010100011110110110010010011101100101111100E-4");
245 mpfr_set_str_binary(t, "0.1111101110100110110110100010101001001000011000111000011101100101110100001110101010110E-4");
246 mpfr_sub(u, x, t, MPFR_RNDU);
247 mpfr_sub(x, x, t, MPFR_RNDU);
248 if (mpfr_cmp(x, u) != 0)
249 {
250 printf ("Error in mpfr_sub: u=x-t and x=x-t give different results\n");
251 exit (1);
252 }
253 if (! MPFR_IS_NORMALIZED (u))
254 {
255 printf ("Error in mpfr_sub: result is not msb-normalized (1)\n");
256 exit (1);
257 }
258 mpfr_set_prec(x, 65); mpfr_set_prec(t, 65); mpfr_set_prec(u, 65);
259 mpfr_set_str_binary(x, "0.10011010101000110101010000000011001001001110001011101011111011101E623");
260 mpfr_set_str_binary(t, "0.10011010101000110101010000000011001001001110001011101011111011100E623");
261 mpfr_sub(u, x, t, MPFR_RNDU);
262 if (mpfr_cmp_ui_2exp(u, 1, 558))
263 { /* 2^558 */
264 printf ("Error (1) in mpfr_sub\n");
265 exit (1);
266 }
267
268 mpfr_set_prec(x, 64); mpfr_set_prec(t, 64); mpfr_set_prec(u, 64);
269 mpfr_set_str_binary(x, "0.1000011110101111011110111111000011101011101111101101101100000100E-220");
270 mpfr_set_str_binary(t, "0.1000011110101111011110111111000011101011101111101101010011111101E-220");
271 test_add (u, x, t, MPFR_RNDU);
272 if ((MPFR_MANT(u)[0] & 1) != 1)
273 {
274 printf ("error in mpfr_add with rnd_mode=MPFR_RNDU\n");
275 printf ("b= "); mpfr_dump (x);
276 printf ("c= "); mpfr_dump (t);
277 printf ("b+c="); mpfr_dump (u);
278 exit (1);
279 }
280
281 /* bug found by Norbert Mueller, 14 Sep 2000 */
282 mpfr_set_prec(x, 56); mpfr_set_prec(t, 83); mpfr_set_prec(u, 10);
283 mpfr_set_str_binary(x, "0.10001001011011001111101100110100000101111010010111010111E-7");
284 mpfr_set_str_binary(t, "0.10001001011011001111101100110100000101111010010111010111000000000111110110110000100E-7");
285 mpfr_sub(u, x, t, MPFR_RNDU);
286
287 /* array bound write found by Norbert Mueller, 26 Sep 2000 */
288 mpfr_set_prec(x, 109); mpfr_set_prec(t, 153); mpfr_set_prec(u, 95);
289 mpfr_set_str_binary(x,"0.1001010000101011101100111000110001111111111111111111111111111111111111111111111111111111111111100000000000000E33");
290 mpfr_set_str_binary(t,"-0.100101000010101110110011100011000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011100101101000000100100001100110111E33");
291 test_add (u, x, t, MPFR_RNDN);
292
293 /* array bound writes found by Norbert Mueller, 27 Sep 2000 */
294 mpfr_set_prec(x, 106); mpfr_set_prec(t, 53); mpfr_set_prec(u, 23);
295 mpfr_set_str_binary(x, "-0.1000011110101111111001010001000100001011000000000000000000000000000000000000000000000000000000000000000000E-59");
296 mpfr_set_str_binary(t, "-0.10000111101011111110010100010001101100011100110100000E-59");
297 mpfr_sub(u, x, t, MPFR_RNDN);
298 mpfr_set_prec(x, 177); mpfr_set_prec(t, 217); mpfr_set_prec(u, 160);
299 mpfr_set_str_binary(x, "-0.111010001011010000111001001010010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E35");
300 mpfr_set_str_binary(t, "0.1110100010110100001110010010100100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111011010011100001111001E35");
301 test_add (u, x, t, MPFR_RNDN);
302 mpfr_set_prec(x, 214); mpfr_set_prec(t, 278); mpfr_set_prec(u, 207);
303 mpfr_set_str_binary(x, "0.1000100110100110101101101101000000010000100111000001001110001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E66");
304 mpfr_set_str_binary(t, "-0.10001001101001101011011011010000000100001001110000010011100010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001111011111001001100011E66");
305 test_add (u, x, t, MPFR_RNDN);
306 mpfr_set_prec(x, 32); mpfr_set_prec(t, 247); mpfr_set_prec(u, 223);
307 mpfr_set_str_binary(x, "0.10000000000000000000000000000000E1");
308 mpfr_set_str_binary(t, "0.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100000110001110100000100011110000101110110011101110100110110111111011010111100100000000000000000000000000E0");
309 mpfr_sub(u, x, t, MPFR_RNDN);
310 if (! MPFR_IS_NORMALIZED (u))
311 {
312 printf ("Error in mpfr_sub: result is not msb-normalized (2)\n");
313 exit (1);
314 }
315
316 /* bug found by Nathalie Revol, 21 March 2001 */
317 mpfr_set_prec (x, 65);
318 mpfr_set_prec (t, 65);
319 mpfr_set_prec (u, 65);
320 mpfr_set_str_binary (x, "0.11100100101101001100111011111111110001101001000011101001001010010E-35");
321 mpfr_set_str_binary (t, "0.10000000000000000000000000000000000001110010010110100110011110000E1");
322 mpfr_sub (u, t, x, MPFR_RNDU);
323 if (! MPFR_IS_NORMALIZED (u))
324 {
325 printf ("Error in mpfr_sub: result is not msb-normalized (3)\n");
326 exit (1);
327 }
328
329 /* bug found by Fabrice Rouillier, 27 Mar 2001 */
330 mpfr_set_prec (x, 107);
331 mpfr_set_prec (t, 107);
332 mpfr_set_prec (u, 107);
333 mpfr_set_str_binary (x, "0.10111001001111010010001000000010111111011011011101000001001000101000000000000000000000000000000000000000000E315");
334 mpfr_set_str_binary (t, "0.10000000000000000000000000000000000101110100100101110110000001100101011111001000011101111100100100111011000E350");
335 mpfr_sub (u, x, t, MPFR_RNDU);
336 if (! MPFR_IS_NORMALIZED (u))
337 {
338 printf ("Error in mpfr_sub: result is not msb-normalized (4)\n");
339 exit (1);
340 }
341
342 /* checks that NaN flag is correctly reset */
343 mpfr_set_ui (t, 1, MPFR_RNDN);
344 mpfr_set_ui (u, 1, MPFR_RNDN);
345 mpfr_set_nan (x);
346 test_add (x, t, u, MPFR_RNDN);
347 if (mpfr_cmp_ui (x, 2))
348 {
349 printf ("Error in mpfr_add: 1+1 gives ");
350 mpfr_out_str(stdout, 10, 0, x, MPFR_RNDN);
351 exit (1);
352 }
353
354 mpfr_clear(x); mpfr_clear(t); mpfr_clear(u);
355 }
356
357 /* check case when c does not overlap with a, but both b and c count
358 for rounding */
359 static void
360 check_case_1b (void)
361 {
362 mpfr_t a, b, c;
363 unsigned int prec_a, prec_b, prec_c, dif;
364
365 mpfr_init (a);
366 mpfr_init (b);
367 mpfr_init (c);
368
369 {
370 prec_a = MPFR_PREC_MIN + (randlimb () % 63);
371 mpfr_set_prec (a, prec_a);
372 for (prec_b = prec_a + 2; prec_b <= 64; prec_b++)
373 {
374 dif = prec_b - prec_a;
375 mpfr_set_prec (b, prec_b);
376 /* b = 1 - 2^(-prec_a) + 2^(-prec_b) */
377 mpfr_set_ui (b, 1, MPFR_RNDN);
378 mpfr_div_2ui (b, b, dif, MPFR_RNDN);
379 mpfr_sub_ui (b, b, 1, MPFR_RNDN);
380 mpfr_div_2ui (b, b, prec_a, MPFR_RNDN);
381 mpfr_add_ui (b, b, 1, MPFR_RNDN);
382 for (prec_c = dif; prec_c <= 64; prec_c++)
383 {
384 /* c = 2^(-prec_a) - 2^(-prec_b) */
385 mpfr_set_prec (c, prec_c);
386 mpfr_set_si (c, -1, MPFR_RNDN);
387 mpfr_div_2ui (c, c, dif, MPFR_RNDN);
388 mpfr_add_ui (c, c, 1, MPFR_RNDN);
389 mpfr_div_2ui (c, c, prec_a, MPFR_RNDN);
390 test_add (a, b, c, MPFR_RNDN);
391 if (mpfr_cmp_ui (a, 1) != 0)
392 {
393 printf ("case (1b) failed for prec_a=%u, prec_b=%u,"
394 " prec_c=%u\n", prec_a, prec_b, prec_c);
395 printf ("b="); mpfr_dump (b);
396 printf ("c="); mpfr_dump (c);
397 printf ("a="); mpfr_dump (a);
398 exit (1);
399 }
400 }
401 }
402 }
403
404 mpfr_clear (a);
405 mpfr_clear (b);
406 mpfr_clear (c);
407 }
408
409 /* check case when c overlaps with a */
410 static void
411 check_case_2 (void)
412 {
413 mpfr_t a, b, c, d;
414
415 mpfr_init2 (a, 300);
416 mpfr_init2 (b, 800);
417 mpfr_init2 (c, 500);
418 mpfr_init2 (d, 800);
419
420 mpfr_set_str_binary(a, "1E110"); /* a = 2^110 */
421 mpfr_set_str_binary(b, "1E900"); /* b = 2^900 */
422 mpfr_set_str_binary(c, "1E500"); /* c = 2^500 */
423 test_add (c, c, a, MPFR_RNDZ); /* c = 2^500 + 2^110 */
424 mpfr_sub (d, b, c, MPFR_RNDZ); /* d = 2^900 - 2^500 - 2^110 */
425 test_add (b, b, c, MPFR_RNDZ); /* b = 2^900 + 2^500 + 2^110 */
426 test_add (a, b, d, MPFR_RNDZ); /* a = 2^901 */
427 if (mpfr_cmp_ui_2exp (a, 1, 901))
428 {
429 printf ("b + d fails for b=2^900+2^500+2^110, d=2^900-2^500-2^110\n");
430 printf ("expected 1.0e901, got ");
431 mpfr_out_str (stdout, 2, 0, a, MPFR_RNDN);
432 printf ("\n");
433 exit (1);
434 }
435
436 mpfr_clear (a);
437 mpfr_clear (b);
438 mpfr_clear (c);
439 mpfr_clear (d);
440 }
441
442 /* checks when source and destination are equal */
443 static void
444 check_same (void)
445 {
446 mpfr_t x;
447
448 mpfr_init(x); mpfr_set_ui(x, 1, MPFR_RNDZ);
449 test_add (x, x, x, MPFR_RNDZ);
450 if (mpfr_cmp_ui (x, 2))
451 {
452 printf ("Error when all 3 operands are equal\n");
453 exit (1);
454 }
455 mpfr_clear(x);
456 }
457
458 #define check53(x, y, r, z) check(x, y, r, 53, 53, 53, z)
459
460 #define MAX_PREC 256
461
462 static void
463 check_inexact (void)
464 {
465 mpfr_t x, y, z, u;
466 mpfr_prec_t px, py, pu, pz;
467 int inexact, cmp;
468 mpfr_rnd_t rnd;
469
470 mpfr_init (x);
471 mpfr_init (y);
472 mpfr_init (z);
473 mpfr_init (u);
474
475 mpfr_set_prec (x, 2);
476 mpfr_set_str_binary (x, "0.1E-4");
477 mpfr_set_prec (u, 33);
478 mpfr_set_str_binary (u, "0.101110100101101100000000111100000E-1");
479 mpfr_set_prec (y, 31);
480 inexact = test_add (y, x, u, MPFR_RNDN);
481
482 if (inexact != 0)
483 {
484 printf ("Wrong ternary value (2): expected 0, got %d\n", inexact);
485 exit (1);
486 }
487
488 mpfr_set_prec (x, 2);
489 mpfr_set_str_binary (x, "0.1E-4");
490 mpfr_set_prec (u, 33);
491 mpfr_set_str_binary (u, "0.101110100101101100000000111100000E-1");
492 mpfr_set_prec (y, 28);
493 inexact = test_add (y, x, u, MPFR_RNDN);
494
495 if (inexact != 0)
496 {
497 printf ("Wrong ternary value (1): expected 0, got %d\n", inexact);
498 exit (1);
499 }
500
501 for (px = 2; px < MAX_PREC; px++)
502 {
503 mpfr_set_prec (x, px);
504
505 do
506 {
507 mpfr_urandomb (x, RANDS);
508 }
509 while (mpfr_cmp_ui (x, 0) == 0);
510
511 for (pu = 2; pu < MAX_PREC; pu++)
512 {
513 mpfr_set_prec (u, pu);
514
515 do
516 {
517 mpfr_urandomb (u, RANDS);
518 }
519 while (mpfr_cmp_ui (u, 0) == 0);
520
521 py = MPFR_PREC_MIN + (randlimb () % (MAX_PREC - 1));
522 mpfr_set_prec (y, py);
523 pz = mpfr_cmpabs (x, u) >= 0 ?
524 MPFR_EXP(x) - MPFR_EXP(u) :
525 MPFR_EXP(u) - MPFR_EXP(x);
526 /* x + u is exactly representable with precision
527 abs(EXP(x)-EXP(u)) + max(prec(x), prec(u)) + 1 */
528 pz = pz + MAX(MPFR_PREC(x), MPFR_PREC(u)) + 1;
529 mpfr_set_prec (z, pz);
530
531 rnd = RND_RAND_NO_RNDF ();
532 inexact = test_add (z, x, u, rnd);
533 if (inexact != 0)
534 {
535 printf ("z <- x + u should be exact\n");
536 printf ("x="); mpfr_dump (x);
537 printf ("u="); mpfr_dump (u);
538 printf ("z="); mpfr_dump (z);
539 exit (1);
540 }
541
542 rnd = RND_RAND_NO_RNDF ();
543 inexact = test_add (y, x, u, rnd);
544 cmp = mpfr_cmp (y, z);
545 if ((inexact == 0 && cmp != 0) ||
546 (inexact > 0 && cmp <= 0) ||
547 (inexact < 0 && cmp >= 0))
548 {
549 printf ("Wrong ternary value for rnd=%s\n",
550 mpfr_print_rnd_mode (rnd));
551 printf ("expected %d, got %d\n", cmp, inexact);
552 printf ("x="); mpfr_dump (x);
553 printf ("u="); mpfr_dump (u);
554 printf ("y= "); mpfr_dump (y);
555 printf ("x+u="); mpfr_dump (z);
556 exit (1);
557 }
558 }
559 }
560
561 mpfr_clear (x);
562 mpfr_clear (y);
563 mpfr_clear (z);
564 mpfr_clear (u);
565 }
566
567 static void
568 check_nans (void)
569 {
570 mpfr_t s, x, y;
571
572 mpfr_init2 (x, 8L);
573 mpfr_init2 (y, 8L);
574 mpfr_init2 (s, 8L);
575
576 /* +inf + -inf == nan */
577 mpfr_set_inf (x, 1);
578 mpfr_set_inf (y, -1);
579 test_add (s, x, y, MPFR_RNDN);
580 MPFR_ASSERTN (mpfr_nan_p (s));
581
582 /* +inf + 1 == +inf */
583 mpfr_set_inf (x, 1);
584 mpfr_set_ui (y, 1L, MPFR_RNDN);
585 test_add (s, x, y, MPFR_RNDN);
586 MPFR_ASSERTN (mpfr_inf_p (s));
587 MPFR_ASSERTN (mpfr_sgn (s) > 0);
588
589 /* -inf + 1 == -inf */
590 mpfr_set_inf (x, -1);
591 mpfr_set_ui (y, 1L, MPFR_RNDN);
592 test_add (s, x, y, MPFR_RNDN);
593 MPFR_ASSERTN (mpfr_inf_p (s));
594 MPFR_ASSERTN (mpfr_sgn (s) < 0);
595
596 /* 1 + +inf == +inf */
597 mpfr_set_ui (x, 1L, MPFR_RNDN);
598 mpfr_set_inf (y, 1);
599 test_add (s, x, y, MPFR_RNDN);
600 MPFR_ASSERTN (mpfr_inf_p (s));
601 MPFR_ASSERTN (mpfr_sgn (s) > 0);
602
603 /* 1 + -inf == -inf */
604 mpfr_set_ui (x, 1L, MPFR_RNDN);
605 mpfr_set_inf (y, -1);
606 test_add (s, x, y, MPFR_RNDN);
607 MPFR_ASSERTN (mpfr_inf_p (s));
608 MPFR_ASSERTN (mpfr_sgn (s) < 0);
609
610 mpfr_clear (x);
611 mpfr_clear (y);
612 mpfr_clear (s);
613 }
614
615 static void
616 check_alloc (void)
617 {
618 mpfr_t a;
619
620 mpfr_init2 (a, 10000);
621 mpfr_set_prec (a, 53);
622 mpfr_set_ui (a, 15236, MPFR_RNDN);
623 test_add (a, a, a, MPFR_RNDN);
624 mpfr_mul (a, a, a, MPFR_RNDN);
625 mpfr_div (a, a, a, MPFR_RNDN);
626 mpfr_sub (a, a, a, MPFR_RNDN);
627 mpfr_clear (a);
628 }
629
630 static void
631 check_overflow (void)
632 {
633 mpfr_t a, b, c;
634 mpfr_prec_t prec_a, prec_b, prec_c;
635 int r, up;
636
637 mpfr_init (a);
638 mpfr_init (b);
639 mpfr_init (c);
640
641 RND_LOOP_NO_RNDF (r)
642 for (prec_a = 2; prec_a <= 128; prec_a += 2)
643 for (prec_b = 2; prec_b <= 128; prec_b += 2)
644 for (prec_c = 2; prec_c <= 128; prec_c += 2)
645 {
646 mpfr_set_prec (a, prec_a);
647 mpfr_set_prec (b, prec_b);
648 mpfr_set_prec (c, prec_c);
649
650 mpfr_setmax (b, mpfr_get_emax ());
651
652 up = r == MPFR_RNDA || r == MPFR_RNDU || r == MPFR_RNDN;
653
654 /* set c with overlap with bits of b: will always overflow */
655 mpfr_set_ui_2exp (c, 1, mpfr_get_emax () - prec_b / 2, MPFR_RNDN);
656 mpfr_nextbelow (c);
657 mpfr_clear_overflow ();
658 test_add (a, b, c, (mpfr_rnd_t) r);
659 if (!mpfr_overflow_p () || (up && !mpfr_inf_p (a)))
660 {
661 printf ("No overflow (1) in check_overflow for rnd=%s\n",
662 mpfr_print_rnd_mode ((mpfr_rnd_t) r));
663 printf ("b="); mpfr_dump (b);
664 printf ("c="); mpfr_dump (c);
665 printf ("a="); mpfr_dump (a);
666 exit (1);
667 }
668
669 if (r == MPFR_RNDZ || r == MPFR_RNDD || prec_a >= prec_b + prec_c)
670 continue;
671
672 /* set c to 111...111 so that ufp(c) = 1/2 ulp(b): will only overflow
673 when prec_a < prec_b + prec_c, and rounding up or to nearest */
674 mpfr_set_ui_2exp (c, 1, mpfr_get_emax () - prec_b, MPFR_RNDN);
675 mpfr_nextbelow (c);
676 mpfr_clear_overflow ();
677 test_add (a, b, c, (mpfr_rnd_t) r);
678 /* b + c is exactly representable iff prec_a >= prec_b + prec_c */
679 if (!mpfr_overflow_p () || !mpfr_inf_p (a))
680 {
681 printf ("No overflow (2) in check_overflow for rnd=%s\n",
682 mpfr_print_rnd_mode ((mpfr_rnd_t) r));
683 printf ("b="); mpfr_dump (b);
684 printf ("c="); mpfr_dump (c);
685 printf ("a="); mpfr_dump (a);
686 exit (1);
687 }
688 }
689
690 mpfr_set_prec (b, 256);
691 mpfr_setmax (b, mpfr_get_emax ());
692 mpfr_set_prec (c, 256);
693 mpfr_set_ui (c, 1, MPFR_RNDN);
694 mpfr_set_exp (c, mpfr_get_emax () - 512);
695 mpfr_set_prec (a, 256);
696 mpfr_clear_overflow ();
697 mpfr_add (a, b, c, MPFR_RNDU);
698 if (!mpfr_overflow_p ())
699 {
700 printf ("No overflow (3) in check_overflow\n");
701 printf ("b="); mpfr_dump (b);
702 printf ("c="); mpfr_dump (c);
703 printf ("a="); mpfr_dump (a);
704 exit (1);
705 }
706
707 mpfr_clear (a);
708 mpfr_clear (b);
709 mpfr_clear (c);
710 }
711
712 static void
713 check_1111 (void)
714 {
715 mpfr_t one;
716 long n;
717
718 mpfr_init2 (one, MPFR_PREC_MIN);
719 mpfr_set_ui (one, 1, MPFR_RNDN);
720 for (n = 0; n < N; n++)
721 {
722 mpfr_prec_t prec_a, prec_b, prec_c;
723 mpfr_exp_t tb=0, tc, diff;
724 mpfr_t a, b, c, s;
725 int m = 512;
726 int sb, sc;
727 int inex_a, inex_s;
728 mpfr_rnd_t rnd_mode;
729
730 prec_a = MPFR_PREC_MIN + (randlimb () % m);
731 prec_b = MPFR_PREC_MIN + (randlimb () % m);
732 /* we need prec_c > 1 so that % (prec_c - 1) is well defined below */
733 do prec_c = MPFR_PREC_MIN + (randlimb () % m); while (prec_c == 1);
734 mpfr_init2 (a, prec_a);
735 mpfr_init2 (b, prec_b);
736 mpfr_init2 (c, prec_c);
737 /* we need prec_b - (sb != 2) > 0 below */
738 do sb = randlimb () % 3; while (prec_b - (sb != 2) == 0);
739 if (sb != 0)
740 {
741 tb = 1 + (randlimb () % (prec_b - (sb != 2)));
742 mpfr_div_2ui (b, one, tb, MPFR_RNDN);
743 if (sb == 2)
744 mpfr_neg (b, b, MPFR_RNDN);
745 test_add (b, b, one, MPFR_RNDN);
746 }
747 else
748 mpfr_set (b, one, MPFR_RNDN);
749 tc = 1 + (randlimb () % (prec_c - 1));
750 mpfr_div_2ui (c, one, tc, MPFR_RNDN);
751 sc = RAND_BOOL ();
752 if (sc)
753 mpfr_neg (c, c, MPFR_RNDN);
754 test_add (c, c, one, MPFR_RNDN);
755 diff = (randlimb () % (2*m)) - m;
756 mpfr_mul_2si (c, c, diff, MPFR_RNDN);
757 rnd_mode = RND_RAND_NO_RNDF ();
758 inex_a = test_add (a, b, c, rnd_mode);
759 mpfr_init2 (s, MPFR_PREC_MIN + 2*m);
760 inex_s = test_add (s, b, c, MPFR_RNDN); /* exact */
761 if (inex_s)
762 {
763 printf ("check_1111: result should have been exact.\n");
764 exit (1);
765 }
766 inex_s = mpfr_prec_round (s, prec_a, rnd_mode);
767 if ((inex_a < 0 && inex_s >= 0) ||
768 (inex_a == 0 && inex_s != 0) ||
769 (inex_a > 0 && inex_s <= 0) ||
770 !mpfr_equal_p (a, s))
771 {
772 printf ("check_1111: results are different.\n");
773 printf ("prec_a = %d, prec_b = %d, prec_c = %d\n",
774 (int) prec_a, (int) prec_b, (int) prec_c);
775 printf ("tb = %d, tc = %d, diff = %d, rnd = %s\n",
776 (int) tb, (int) tc, (int) diff,
777 mpfr_print_rnd_mode (rnd_mode));
778 printf ("sb = %d, sc = %d\n", sb, sc);
779 printf ("b = "); mpfr_dump (b);
780 printf ("c = "); mpfr_dump (c);
781 printf ("a = "); mpfr_dump (a);
782 printf ("s = "); mpfr_dump (s);
783 printf ("inex_a = %d, inex_s = %d\n", inex_a, inex_s);
784 exit (1);
785 }
786 mpfr_clear (a);
787 mpfr_clear (b);
788 mpfr_clear (c);
789 mpfr_clear (s);
790 }
791 mpfr_clear (one);
792 }
793
794 static void
795 check_1minuseps (void)
796 {
797 static mpfr_prec_t prec_a[] = {
798 MPFR_PREC_MIN, 30, 31, 32, 33, 62, 63, 64, 65, 126, 127, 128, 129
799 };
800 static int supp_b[] = {
801 0, 1, 2, 3, 4, 29, 30, 31, 32, 33, 34, 35, 61, 62, 63, 64, 65, 66, 67
802 };
803 mpfr_t a, b, c;
804 unsigned int ia, ib, ic;
805
806 mpfr_init2 (c, MPFR_PREC_MIN);
807
808 for (ia = 0; ia < numberof (prec_a); ia++)
809 for (ib = 0; ib < numberof (supp_b); ib++)
810 {
811 mpfr_prec_t prec_b;
812 int rnd_mode;
813
814 prec_b = prec_a[ia] + supp_b[ib];
815
816 mpfr_init2 (a, prec_a[ia]);
817 mpfr_init2 (b, prec_b);
818
819 mpfr_set_ui (c, 1, MPFR_RNDN);
820 mpfr_div_ui (b, c, prec_a[ia], MPFR_RNDN);
821 mpfr_sub (b, c, b, MPFR_RNDN); /* b = 1 - 2^(-prec_a) */
822
823 for (ic = 0; ic < numberof (supp_b); ic++)
824 RND_LOOP (rnd_mode)
825 {
826 mpfr_t s;
827 int inex_a, inex_s;
828
829 if (rnd_mode == MPFR_RNDF)
830 continue; /* inex_a, inex_s make no sense */
831
832 mpfr_set_ui (c, 1, MPFR_RNDN);
833 mpfr_div_ui (c, c, prec_a[ia] + supp_b[ic], MPFR_RNDN);
834 inex_a = test_add (a, b, c, (mpfr_rnd_t) rnd_mode);
835 mpfr_init2 (s, 256);
836 inex_s = test_add (s, b, c, MPFR_RNDN); /* exact */
837 if (inex_s)
838 {
839 printf ("check_1minuseps: result should have been exact "
840 "(ia = %u, ib = %u, ic = %u)\n", ia, ib, ic);
841 exit (1);
842 }
843 inex_s = mpfr_prec_round (s, prec_a[ia], (mpfr_rnd_t) rnd_mode);
844 if ((inex_a < 0 && inex_s >= 0) ||
845 (inex_a == 0 && inex_s != 0) ||
846 (inex_a > 0 && inex_s <= 0) ||
847 !mpfr_equal_p (a, s))
848 {
849 printf ("check_1minuseps: results are different.\n");
850 printf ("ia = %u, ib = %u, ic = %u\n", ia, ib, ic);
851 exit (1);
852 }
853 mpfr_clear (s);
854 }
855
856 mpfr_clear (a);
857 mpfr_clear (b);
858 }
859
860 mpfr_clear (c);
861 }
862
863 /* Test case bk == 0 in add1.c (b has entirely been read and
864 c hasn't been taken into account). */
865 static void
866 coverage_bk_eq_0 (void)
867 {
868 mpfr_t a, b, c;
869 int inex;
870
871 mpfr_init2 (a, GMP_NUMB_BITS);
872 mpfr_init2 (b, 2 * GMP_NUMB_BITS);
873 mpfr_init2 (c, GMP_NUMB_BITS);
874
875 mpfr_set_ui_2exp (b, 1, 2 * GMP_NUMB_BITS, MPFR_RNDN);
876 mpfr_sub_ui (b, b, 1, MPFR_RNDN);
877 /* b = 111...111 (in base 2) where the 1's fit 2 whole limbs */
878
879 mpfr_set_ui_2exp (c, 1, -1, MPFR_RNDN); /* c = 1/2 */
880
881 inex = mpfr_add (a, b, c, MPFR_RNDU);
882 mpfr_set_ui_2exp (c, 1, 2 * GMP_NUMB_BITS, MPFR_RNDN);
883 if (! mpfr_equal_p (a, c))
884 {
885 printf ("Error in coverage_bk_eq_0\n");
886 printf ("Expected ");
887 mpfr_dump (c);
888 printf ("Got ");
889 mpfr_dump (a);
890 exit (1);
891 }
892 MPFR_ASSERTN (inex > 0);
893
894 mpfr_clear (a);
895 mpfr_clear (b);
896 mpfr_clear (c);
897 }
898
899 static void
900 tests (void)
901 {
902 check_alloc ();
903 check_nans ();
904 check_inexact ();
905 check_case_1b ();
906 check_case_2 ();
907 check64();
908 coverage_bk_eq_0 ();
909
910 check("293607738.0", "1.9967571564050541e-5", MPFR_RNDU, 64, 53, 53,
911 "2.9360773800002003e8");
912 check("880524.0", "-2.0769715792901673e-5", MPFR_RNDN, 64, 53, 53,
913 "8.8052399997923023e5");
914 check("1196426492.0", "-1.4218093058435347e-3", MPFR_RNDN, 64, 53, 53,
915 "1.1964264919985781e9");
916 check("982013018.0", "-8.941829477291838e-7", MPFR_RNDN, 64, 53, 53,
917 "9.8201301799999905e8");
918 check("1092583421.0", "1.0880649218158844e9", MPFR_RNDN, 64, 53, 53,
919 "2.1806483428158846e9");
920 check("1.8476886419022969e-6", "961494401.0", MPFR_RNDN, 53, 64, 53,
921 "9.6149440100000179e8");
922 check("-2.3222118418069868e5", "1229318102.0", MPFR_RNDN, 53, 64, 53,
923 "1.2290858808158193e9");
924 check("-3.0399171300395734e-6", "874924868.0", MPFR_RNDN, 53, 64, 53,
925 "8.749248679999969e8");
926 check("9.064246624706179e1", "663787413.0", MPFR_RNDN, 53, 64, 53,
927 "6.6378750364246619e8");
928 check("-1.0954322421551264e2", "281806592.0", MPFR_RNDD, 53, 64, 53,
929 "2.8180648245677572e8");
930 check("5.9836930386056659e-8", "1016217213.0", MPFR_RNDN, 53, 64, 53,
931 "1.0162172130000001e9");
932 check("-1.2772161928500301e-7", "1237734238.0", MPFR_RNDN, 53, 64, 53,
933 "1.2377342379999998e9");
934 check("-4.567291988483277e8", "1262857194.0", MPFR_RNDN, 53, 64, 53,
935 "8.0612799515167236e8");
936 check("4.7719471752925262e7", "196089880.0", MPFR_RNDN, 53, 53, 53,
937 "2.4380935175292528e8");
938 check("4.7719471752925262e7", "196089880.0", MPFR_RNDN, 53, 64, 53,
939 "2.4380935175292528e8");
940 check("-1.716113812768534e-140", "1271212614.0", MPFR_RNDZ, 53, 64, 53,
941 "1.2712126139999998e9");
942 check("-1.2927455200185474e-50", "1675676122.0", MPFR_RNDD, 53, 64, 53,
943 "1.6756761219999998e9");
944
945 check53("1.22191250737771397120e+20", "948002822.0", MPFR_RNDN,
946 "122191250738719408128.0");
947 check53("9966027674114492.0", "1780341389094537.0", MPFR_RNDN,
948 "11746369063209028.0");
949 check53("2.99280481918991653800e+272", "5.34637717585790933424e+271",
950 MPFR_RNDN, "3.5274425367757071711e272");
951 check_same();
952 check53("6.14384195492641560499e-02", "-6.14384195401037683237e-02",
953 MPFR_RNDU, "9.1603877261370314499e-12");
954 check53("1.16809465359248765399e+196", "7.92883212101990665259e+196",
955 MPFR_RNDU, "9.0969267746123943065e196");
956 check53("3.14553393112021279444e-67", "3.14553401015952024126e-67", MPFR_RNDU,
957 "6.2910679412797336946e-67");
958
959 check53("5.43885304644369509058e+185","-1.87427265794105342763e-57",MPFR_RNDN,
960 "5.4388530464436950905e185");
961 check53("5.43885304644369509058e+185","-1.87427265794105342763e-57",MPFR_RNDZ,
962 "5.4388530464436944867e185");
963 check53("5.43885304644369509058e+185","-1.87427265794105342763e-57",MPFR_RNDU,
964 "5.4388530464436950905e185");
965 check53("5.43885304644369509058e+185","-1.87427265794105342763e-57",MPFR_RNDD,
966 "5.4388530464436944867e185");
967
968 check2b("1.001010101110011000000010100101110010111001010000001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e358",187,
969 "-1.11100111001101100010001111111110101101110001000000000000000000000000000000000000000000e160",87,
970 "1.001010101110011000000010100101110010111001010000000111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111e358",178,
971 MPFR_RNDD);
972 check2b("-1.111100100011100111010101010101001010100100000111001000000000000000000e481",70,
973 "1.1111000110100011110101111110110010010000000110101000000000000000e481",65,
974 "-1.001010111111101011010000001100011101100101000000000000000000e472",61,
975 MPFR_RNDD);
976 check2b("1.0100010111010000100101000000111110011100011001011010000000000000000000000000000000e516",83,
977 "-1.1001111000100001011100000001001100110011110010111111000000e541",59,
978 "-1.1001111000100001011011110111000001001011100000011110100000110001110011010011000000000000000000000000000000000000000000000000e541",125,
979 MPFR_RNDZ);
980 check2b("-1.0010111100000100110001011011010000000011000111101000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e261",155,
981 "-1.00111110100011e239",15,
982 "-1.00101111000001001100101010101110001100110001111010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e261",159,
983 MPFR_RNDD);
984 check2b("-1.110111000011111011000000001001111101101001010100111000000000000000000000000e880",76,
985 "-1.1010010e-634",8,
986 "-1.11011100001111101100000000100111110110100101010011100000000000000000000000e880",75,
987 MPFR_RNDZ);
988 check2b("1.00100100110110101001010010101111000001011100100101010000000000000000000000000000e-530",81,
989 "-1.101101111100000111000011001010110011001011101001110100000e-908",58,
990 "1.00100100110110101001010010101111000001011100100101010e-530",54,
991 MPFR_RNDN);
992 check2b("1.0101100010010111101000000001000010010010011000111011000000000000000000000000000000000000000000000000000000000000000000e374",119,
993 "1.11100101100101e358",15,
994 "1.01011000100110011000010110100100100100100110001110110000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e374",150,
995 MPFR_RNDZ);
996 check2b("-1.10011001000010100000010100100110110010011111101111000000000000000000000000000000000000000000000000000000000000000000e-172",117,
997 "1.111011100000101010110000100100110100100001001000011100000000e-173",61,
998 "-1.0100010000001001010110011011101001001011101011110001000000000000000e-173",68,
999 MPFR_RNDZ);
1000 check2b("-1.011110000111101011100001100110100011100101000011011000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-189",175,
1001 "1.1e631",2,
1002 "1.011111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111e631",115,
1003 MPFR_RNDZ);
1004 check2b("-1.101011101001011101100011001000001100010100001101011000000000000000000000000000000000000000000e-449",94,
1005 "-1.01111101010111000011000110011101000111001100110111100000000000000e-429",66,
1006 "-1.01111101010111000100110010000110100100101111111111101100010100001101011000000000000000000000000000000000000000e-429",111,
1007 MPFR_RNDU);
1008 check2b("-1.101011101001011101100011001000001100010100001101011000000000000000000000000000000000000000000e-449",94,
1009 "-1.01111101010111000011000110011101000111001100110111100000000000000e-429",66,
1010 "-1.01111101010111000100110010000110100100101111111111101100010100001101011000000000000000000000000000000000000000e-429",111,
1011 MPFR_RNDD);
1012 check2b("-1.1001000011101000110000111110010100100101110101111100000000000000000000000000000000000000000000000000000000e-72",107,
1013 "-1.001100011101100100010101101010101011010010111111010000000000000000000000000000e521",79,
1014 "-1.00110001110110010001010110101010101101001011111101000000000000000000000000000000000000000000000001e521",99,
1015 MPFR_RNDD);
1016 check2b("-1.01010001111000000101010100100100110101011011100001110000000000e498",63,
1017 "1.010000011010101111000100111100011100010101011110010100000000000e243",64,
1018 "-1.010100011110000001010101001001001101010110111000011100000000000e498",64,
1019 MPFR_RNDN);
1020 check2b("1.00101100010101000011010000011000111111011110010111000000000000000000000000000000000000000000000000000000000e178",108,
1021 "-1.10101101010101000110011011111001001101111111110000100000000e160",60,
1022 "1.00101100010100111100100011000011111001000010011101110010000000001111100000000000000000000000000000000000e178",105,
1023 MPFR_RNDN);
1024 check2b("1.00110011010100111110011010110100111101110101100100110000000000000000000000000000000000000000000000e559",99,
1025 "-1.011010110100111011100110100110011100000000111010011000000000000000e559",67,
1026 "-1.101111111101011111111111001001100100011100001001100000000000000000000000000000000000000000000e556",94,
1027 MPFR_RNDU);
1028 check2b("-1.100000111100101001100111011100011011000001101001111100000000000000000000000000e843",79,
1029 "-1.1101101010110000001001000100001100110011000110110111000000000000000000000000000000000000000000e414",95,
1030 "-1.1000001111001010011001110111000110110000011010100000e843",53,
1031 MPFR_RNDD);
1032 check2b("-1.110110010110100010100011000110111001010000010111110000000000e-415",61,
1033 "-1.0000100101100001111100110011111111110100011101101011000000000000000000e751",71,
1034 "-1.00001001011000011111001100111111111101000111011010110e751",54,
1035 MPFR_RNDN);
1036 check2b("-1.1011011011110001001101010101001000010100010110111101000000000000000000000e258",74,
1037 "-1.00011100010110110101001011000100100000100010101000010000000000000000000000000000000000000000000000e268",99,
1038 "-1.0001110011001001000011110001000111010110101011110010011011110100000000000000000000000000000000000000e268",101,
1039 MPFR_RNDD);
1040 check2b("-1.1011101010011101011000000100100110101101101110000001000000000e629",62,
1041 "1.111111100000011100100011100000011101100110111110111000000000000000000000000000000000000000000e525",94,
1042 "-1.101110101001110101100000010010011010110110111000000011111111111111111111111111111111111111111111111111101e629",106,
1043 MPFR_RNDD);
1044 check2b("1.111001000010001100010000001100000110001011110111011000000000000000000000000000000000000e152",88,
1045 "1.111110111001100100000100111111010111000100111111001000000000000000e152",67,
1046 "1.1110111111011110000010101001011011101010000110110100e153",53,
1047 MPFR_RNDN);
1048 check2b("1.000001100011110010110000110100001010101101111011110100e696",55,
1049 "-1.1011001111011100100001011110100101010101110111010101000000000000000000000000000000000000000000000000000000000000e730",113,
1050 "-1.1011001111011100100001011110100100010100010011100010e730",53,
1051 MPFR_RNDN);
1052 check2b("-1.11010111100001001111000001110101010010001111111001100000000000000000000000000000000000000000000000000000000000e530",111,
1053 "1.01110100010010000000010110111101011101000001111101100000000000000000000000000000000000000000000000e530",99,
1054 "-1.1000110011110011101010101101111101010011011111000000000000000e528",62,
1055 MPFR_RNDD);
1056 check2b("-1.0001100010010100111101101011101000100100010011100011000000000000000000000000000000000000000000000000000000000e733",110,
1057 "-1.001000000111110010100101010100110111001111011011001000000000000000000000000000000000000000000000000000000000e710",109,
1058 "-1.000110001001010011111000111110110001110110011000110110e733",55,
1059 MPFR_RNDN);
1060 check2b("-1.1101011110000100111100000111010101001000111111100110000000000000000000000e530",74,
1061 "1.01110100010010000000010110111101011101000001111101100000000000000000000000000000000000000000000000000000000000e530",111,
1062 "-1.10001100111100111010101011011111010100110111110000000000000000000000000000e528",75,
1063 MPFR_RNDU);
1064 check2b("1.00110011010100111110011010110100111101110101100100110000000000000000000000000000000000000000000000e559",99,
1065 "-1.011010110100111011100110100110011100000000111010011000000000000000e559",67,
1066 "-1.101111111101011111111111001001100100011100001001100000000000000000000000000000000000000000000e556",94,
1067 MPFR_RNDU);
1068 check2b("-1.100101111110110000000110111111011010011101101111100100000000000000e-624",67,
1069 "1.10111010101110100000010110101000000000010011100000100000000e-587",60,
1070 "1.1011101010111010000001011010011111110100011110001011111111001000000100101100010010000011100000000000000000000e-587",110,
1071 MPFR_RNDU);
1072 check2b("-1.10011001000010100000010100100110110010011111101111000000000000000000000000000000000000000000000000000000000000000000e-172",117,
1073 "1.111011100000101010110000100100110100100001001000011100000000e-173",61,
1074 "-1.0100010000001001010110011011101001001011101011110001000000000000000e-173",68,
1075 MPFR_RNDZ);
1076 check2b("1.1000111000110010101001010011010011101100010110001001000000000000000000000000000000000000000000000000e167",101,
1077 "1.0011110010000110000000101100100111000001110110110000000000000000000000000e167",74,
1078 "1.01100101010111000101001111111111010101110001100111001000000000000000000000000000000000000000000000000000e168",105,
1079 MPFR_RNDZ);
1080 check2b("1.100101111111110010100101110111100001110000100001010000000000000000000000000000000000000000000000e808",97,
1081 "-1.1110011001100000100000111111110000110010100111001011000000000000000000000000000000e807",83,
1082 "1.01001001100110001100011111000000000001011010010111010000000000000000000000000000000000000000000e807",96,
1083 MPFR_RNDN);
1084 check2b("1e128",128,
1085 "1e0",128,
1086 "100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001e0",256,
1087 MPFR_RNDN);
1088 check2b("0.10000000000000000011101111111110110110001100010110100011010000110101101101010010000110110011100110001101110010111011011110010101e1",128,
1089 "0.11110011001100100001111011100010010100011011011111111110101001101001101011100101011110101111111100010010100010100111110110011001e-63",128,
1090 "0.10000000000000000011101111111110110110001100010110100011010001000100111010000100001110100001101111011111100000111011011000111100e1",128,
1091 MPFR_RNDN);
1092 check2b("0.101000111001110100111100000010110010011101001011101001001101101000010100010000011100110110110011011000001110110101100111010110011101000111010011110001010001110110110011111011011100101011111100e-2",192,
1093 "0.101101100000111110011110000110001000110110010010110000111101110010111010100110000000000110010100010001000001011011101010111010101011010110011011011110110111001000111001101000010101111010010010e-130",192,
1094 "0.101000111001110100111100000010110010011101001011101001001101101000010100010000011100110110110011011000001110110101100111010110101000011111100011011000110011011001000001100000001000111011011001e-2",192,
1095 MPFR_RNDN);
1096
1097 /* Checking double precision (53 bits) */
1098 check53("-8.22183238641455905806e-19", "7.42227178769761587878e-19",MPFR_RNDD,
1099 "-7.9956059871694317927e-20");
1100 check53("5.82106394662028628236e+234","-5.21514064202368477230e+89",MPFR_RNDD,
1101 "5.8210639466202855763e234");
1102 check53("5.72931679569871602371e+122","-5.72886070363264321230e+122",
1103 MPFR_RNDN, "4.5609206607281141508e118");
1104 check53("-5.09937369394650450820e+238", "2.70203299854862982387e+250",
1105 MPFR_RNDD, "2.7020329985435301323e250");
1106 check53("-2.96695924472363684394e+27", "1.22842938251111500000e+16",MPFR_RNDD,
1107 "-2.96695924471135255027e27");
1108 check53("1.74693641655743793422e-227", "-7.71776956366861843469e-229",
1109 MPFR_RNDN, "1.669758720920751867e-227");
1110 /* x = -7883040437021647.0; for (i=0; i<468; i++) x = x / 2.0;*/
1111 check53("-1.03432206392780011159e-125", "1.30127034799251347548e-133",
1112 MPFR_RNDN,
1113 "-1.0343220509150965661100887242027378881805094180354e-125");
1114 check53("1.05824655795525779205e+71", "-1.06022698059744327881e+71",MPFR_RNDZ,
1115 "-1.9804226421854867632e68");
1116 check53("-5.84204911040921732219e+240", "7.26658169050749590763e+240",
1117 MPFR_RNDD, "1.4245325800982785854e240");
1118 check53("1.00944884131046636376e+221","2.33809162651471520268e+215",MPFR_RNDN,
1119 "1.0094511794020929787e221");
1120 /*x = 7045852550057985.0; for (i=0; i<986; i++) x = x / 2.0;*/
1121 check53("4.29232078932667367325e-278",
1122 "1.0773525047389793833221116707010783793203080117586e-281"
1123 , MPFR_RNDU, "4.2933981418314132787e-278");
1124 check53("5.27584773801377058681e-80", "8.91207657803547196421e-91", MPFR_RNDN,
1125 "5.2758477381028917269e-80");
1126 check53("2.99280481918991653800e+272", "5.34637717585790933424e+271",
1127 MPFR_RNDN, "3.5274425367757071711e272");
1128 check53("4.67302514390488041733e-184", "2.18321376145645689945e-190",
1129 MPFR_RNDN, "4.6730273271186420541e-184");
1130 check53("5.57294120336300389254e+71", "2.60596167942024924040e+65", MPFR_RNDZ,
1131 "5.5729438093246831053e71");
1132 check53("6.6052588496951015469e24", "4938448004894539.0", MPFR_RNDU,
1133 "6.6052588546335505068e24");
1134 check53("1.23056185051606761523e-190", "1.64589756643433857138e-181",
1135 MPFR_RNDU, "1.6458975676649006598e-181");
1136 check53("2.93231171510175981584e-280", "3.26266919161341483877e-273",
1137 MPFR_RNDU, "3.2626694848445867288e-273");
1138 check53("5.76707395945001907217e-58", "4.74752971449827687074e-51", MPFR_RNDD,
1139 "4.747530291205672325e-51");
1140 check53("277363943109.0", "11.0", MPFR_RNDN, "277363943120.0");
1141 check53("1.44791789689198883921e-140", "-1.90982880222349071284e-121",
1142 MPFR_RNDN, "-1.90982880222349071e-121");
1143
1144
1145 /* tests for particular cases (Vincent Lefevre, 22 Aug 2001) */
1146 check53("9007199254740992.0", "1.0", MPFR_RNDN, "9007199254740992.0");
1147 check53("9007199254740994.0", "1.0", MPFR_RNDN, "9007199254740996.0");
1148 check53("9007199254740992.0", "-1.0", MPFR_RNDN, "9007199254740991.0");
1149 check53("9007199254740994.0", "-1.0", MPFR_RNDN, "9007199254740992.0");
1150 check53("9007199254740996.0", "-1.0", MPFR_RNDN, "9007199254740996.0");
1151
1152 check_overflow ();
1153 check_1111 ();
1154 check_1minuseps ();
1155 }
1156
1157 static void
1158 check_extreme (void)
1159 {
1160 mpfr_t u, v, w, x, y;
1161 int i, inex, r;
1162
1163 mpfr_inits2 (32, u, v, w, x, y, (mpfr_ptr) 0);
1164 mpfr_setmin (u, mpfr_get_emax ());
1165 mpfr_setmax (v, mpfr_get_emin ());
1166 mpfr_setmin (w, mpfr_get_emax () - 40);
1167 RND_LOOP (r)
1168 for (i = 0; i < 2; i++)
1169 {
1170 mpfr_add (x, u, v, (mpfr_rnd_t) r);
1171 mpfr_set_prec (y, 64);
1172 inex = mpfr_add (y, u, w, MPFR_RNDN);
1173 MPFR_ASSERTN (inex == 0);
1174 mpfr_prec_round (y, 32, (mpfr_rnd_t) r);
1175 /* for RNDF, the rounding directions are not necessarily consistent */
1176 if (! mpfr_equal_p (x, y) && r != MPFR_RNDF)
1177 {
1178 printf ("Error in u + v (%s)\n",
1179 mpfr_print_rnd_mode ((mpfr_rnd_t) r));
1180 printf ("u = ");
1181 mpfr_dump (u);
1182 printf ("v = ");
1183 mpfr_dump (v);
1184 printf ("Expected ");
1185 mpfr_dump (y);
1186 printf ("Got ");
1187 mpfr_dump (x);
1188 exit (1);
1189 }
1190 mpfr_neg (v, v, MPFR_RNDN);
1191 mpfr_neg (w, w, MPFR_RNDN);
1192 }
1193 mpfr_clears (u, v, w, x, y, (mpfr_ptr) 0);
1194 }
1195
1196 static void
1197 testall_rndf (mpfr_prec_t pmax)
1198 {
1199 mpfr_t a, b, c, d;
1200 mpfr_prec_t pa, pb, pc;
1201 mpfr_exp_t eb;
1202
1203 for (pa = MPFR_PREC_MIN; pa <= pmax; pa++)
1204 {
1205 mpfr_init2 (a, pa);
1206 mpfr_init2 (d, pa);
1207 for (pb = MPFR_PREC_MIN; pb <= pmax; pb++)
1208 {
1209 mpfr_init2 (b, pb);
1210 for (eb = 0; eb <= pmax + 3; eb ++)
1211 {
1212 mpfr_set_ui_2exp (b, 1, eb, MPFR_RNDN);
1213 while (mpfr_cmp_ui_2exp (b, 1, eb + 1) < 0)
1214 {
1215 for (pc = MPFR_PREC_MIN; pc <= pmax; pc++)
1216 {
1217 mpfr_init2 (c, pc);
1218 mpfr_set_ui (c, 1, MPFR_RNDN);
1219 while (mpfr_cmp_ui (c, 2) < 0)
1220 {
1221 mpfr_add (a, b, c, MPFR_RNDF);
1222 mpfr_add (d, b, c, MPFR_RNDD);
1223 if (!mpfr_equal_p (a, d))
1224 {
1225 mpfr_add (d, b, c, MPFR_RNDU);
1226 if (!mpfr_equal_p (a, d))
1227 {
1228 printf ("Error: mpfr_add(a,b,c,RNDF) does not "
1229 "match RNDD/RNDU\n");
1230 printf ("b="); mpfr_dump (b);
1231 printf ("c="); mpfr_dump (c);
1232 printf ("a="); mpfr_dump (a);
1233 exit (1);
1234 }
1235 }
1236 mpfr_nextabove (c);
1237 }
1238 mpfr_clear (c);
1239 }
1240 mpfr_nextabove (b);
1241 }
1242 }
1243 mpfr_clear (b);
1244 }
1245 mpfr_clear (a);
1246 mpfr_clear (d);
1247 }
1248 }
1249
1250 static void
1251 test_rndf_exact (mp_size_t pmax)
1252 {
1253 mpfr_t a, b, c, d;
1254 mpfr_prec_t pa, pb, pc;
1255 mpfr_exp_t eb;
1256
1257 for (pa = MPFR_PREC_MIN; pa <= pmax; pa++)
1258 {
1259 /* only check pa mod GMP_NUMB_BITS = -2, -1, 0, 1, 2 */
1260 if ((pa + 2) % GMP_NUMB_BITS > 4)
1261 continue;
1262 mpfr_init2 (a, pa);
1263 mpfr_init2 (d, pa);
1264 for (pb = MPFR_PREC_MIN; pb <= pmax; pb++)
1265 {
1266 if ((pb + 2) % GMP_NUMB_BITS > 4)
1267 continue;
1268 mpfr_init2 (b, pb);
1269 for (eb = 0; eb <= pmax + 3; eb ++)
1270 {
1271 mpfr_urandomb (b, RANDS);
1272 mpfr_mul_2ui (b, b, eb, MPFR_RNDN);
1273 for (pc = MPFR_PREC_MIN; pc <= pmax; pc++)
1274 {
1275 if ((pc + 2) % GMP_NUMB_BITS > 4)
1276 continue;
1277 mpfr_init2 (c, pc);
1278 mpfr_urandomb (c, RANDS);
1279 mpfr_add (a, b, c, MPFR_RNDF);
1280 mpfr_add (d, b, c, MPFR_RNDD);
1281 if (!mpfr_equal_p (a, d))
1282 {
1283 mpfr_add (d, b, c, MPFR_RNDU);
1284 if (!mpfr_equal_p (a, d))
1285 {
1286 printf ("Error: mpfr_add(a,b,c,RNDF) does not "
1287 "match RNDD/RNDU\n");
1288 printf ("b="); mpfr_dump (b);
1289 printf ("c="); mpfr_dump (c);
1290 printf ("a="); mpfr_dump (a);
1291 exit (1);
1292 }
1293 }
1294
1295 /* now make the low bits from c match those from b */
1296 mpfr_sub (c, b, d, MPFR_RNDN);
1297 mpfr_add (a, b, c, MPFR_RNDF);
1298 mpfr_add (d, b, c, MPFR_RNDD);
1299 if (!mpfr_equal_p (a, d))
1300 {
1301 mpfr_add (d, b, c, MPFR_RNDU);
1302 if (!mpfr_equal_p (a, d))
1303 {
1304 printf ("Error: mpfr_add(a,b,c,RNDF) does not "
1305 "match RNDD/RNDU\n");
1306 printf ("b="); mpfr_dump (b);
1307 printf ("c="); mpfr_dump (c);
1308 printf ("a="); mpfr_dump (a);
1309 exit (1);
1310 }
1311 }
1312
1313 mpfr_clear (c);
1314 }
1315 }
1316 mpfr_clear (b);
1317 }
1318 mpfr_clear (a);
1319 mpfr_clear (d);
1320 }
1321 }
1322
1323 #define TEST_FUNCTION test_add
1324 #define TWO_ARGS
1325 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
1326 #include "tgeneric.c"
1327
1328 int
1329 main (int argc, char *argv[])
1330 {
1331 tests_start_mpfr ();
1332
1333 usesp = 0;
1334 tests ();
1335
1336 #ifndef CHECK_EXTERNAL /* no need to check twice */
1337 usesp = 1;
1338 tests ();
1339 #endif
1340
1341 test_rndf_exact (200);
1342 testall_rndf (7);
1343 check_extreme ();
1344
1345 test_generic (MPFR_PREC_MIN, 1000, 100);
1346
1347 tests_end_mpfr ();
1348 return 0;
1349 }