tgmpop.c revision 1.1 1 1.1 mrg /* Test file for mpfr_add_[q,z], mpfr_sub_[q,z], mpfr_div_[q,z], mpfr_mul_[q,z]
2 1.1 mrg and mpfr_cmp_[q,z]
3 1.1 mrg
4 1.1 mrg Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
5 1.1 mrg Contributed by the Arenaire and Cacao projects, INRIA.
6 1.1 mrg
7 1.1 mrg This file is part of the GNU MPFR Library.
8 1.1 mrg
9 1.1 mrg The GNU MPFR Library is free software; you can redistribute it and/or modify
10 1.1 mrg it under the terms of the GNU Lesser General Public License as published by
11 1.1 mrg the Free Software Foundation; either version 3 of the License, or (at your
12 1.1 mrg option) any later version.
13 1.1 mrg
14 1.1 mrg The GNU MPFR Library is distributed in the hope that it will be useful, but
15 1.1 mrg WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 1.1 mrg or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 1.1 mrg License for more details.
18 1.1 mrg
19 1.1 mrg You should have received a copy of the GNU Lesser General Public License
20 1.1 mrg along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
21 1.1 mrg http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 1.1 mrg 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 1.1 mrg
24 1.1 mrg #include <stdio.h>
25 1.1 mrg #include <stdlib.h>
26 1.1 mrg #include "mpfr-test.h"
27 1.1 mrg
28 1.1 mrg #define CHECK_FOR(str, cond) \
29 1.1 mrg if ((cond) == 0) { \
30 1.1 mrg printf ("Special case error " str ". Inexact value=%d\n", res); \
31 1.1 mrg printf ("Get "); mpfr_dump (y); \
32 1.1 mrg printf ("X= "); mpfr_dump (x); \
33 1.1 mrg printf ("Z= "); mpz_dump (mpq_numref(z)); \
34 1.1 mrg printf (" /"); mpz_dump (mpq_denref(z)); \
35 1.1 mrg exit (1); \
36 1.1 mrg }
37 1.1 mrg
38 1.1 mrg static void
39 1.1 mrg special (void)
40 1.1 mrg {
41 1.1 mrg mpfr_t x, y;
42 1.1 mrg mpq_t z;
43 1.1 mrg int res = 0;
44 1.1 mrg
45 1.1 mrg mpfr_init (x);
46 1.1 mrg mpfr_init (y);
47 1.1 mrg mpq_init (z);
48 1.1 mrg
49 1.1 mrg /* cancellation in mpfr_add_q */
50 1.1 mrg mpfr_set_prec (x, 60);
51 1.1 mrg mpfr_set_prec (y, 20);
52 1.1 mrg mpz_set_str (mpq_numref (z), "-187207494", 10);
53 1.1 mrg mpz_set_str (mpq_denref (z), "5721", 10);
54 1.1 mrg mpfr_set_str_binary (x, "11111111101001011011100101100011011110010011100010000100001E-44");
55 1.1 mrg mpfr_add_q (y, x, z, MPFR_RNDN);
56 1.1 mrg CHECK_FOR ("cancelation in add_q", mpfr_cmp_ui_2exp (y, 256783, -64) == 0);
57 1.1 mrg
58 1.1 mrg mpfr_set_prec (x, 19);
59 1.1 mrg mpfr_set_str_binary (x, "0.1011110101110011100E0");
60 1.1 mrg mpz_set_str (mpq_numref (z), "187207494", 10);
61 1.1 mrg mpz_set_str (mpq_denref (z), "5721", 10);
62 1.1 mrg mpfr_set_prec (y, 29);
63 1.1 mrg mpfr_add_q (y, x, z, MPFR_RNDD);
64 1.1 mrg mpfr_set_prec (x, 29);
65 1.1 mrg mpfr_set_str_binary (x, "11111111101001110011010001001E-14");
66 1.1 mrg CHECK_FOR ("cancelation in add_q", mpfr_cmp (x,y) == 0);
67 1.1 mrg
68 1.1 mrg /* Inf */
69 1.1 mrg mpfr_set_inf (x, 1);
70 1.1 mrg mpz_set_str (mpq_numref (z), "395877315", 10);
71 1.1 mrg mpz_set_str (mpq_denref (z), "3508975966", 10);
72 1.1 mrg mpfr_set_prec (y, 118);
73 1.1 mrg mpfr_add_q (y, x, z, MPFR_RNDU);
74 1.1 mrg CHECK_FOR ("inf", mpfr_inf_p (y) && mpfr_sgn (y) > 0);
75 1.1 mrg mpfr_sub_q (y, x, z, MPFR_RNDU);
76 1.1 mrg CHECK_FOR ("inf", mpfr_inf_p (y) && mpfr_sgn (y) > 0);
77 1.1 mrg
78 1.1 mrg /* Nan */
79 1.1 mrg MPFR_SET_NAN (x);
80 1.1 mrg mpfr_add_q (y, x, z, MPFR_RNDU);
81 1.1 mrg CHECK_FOR ("nan", mpfr_nan_p (y));
82 1.1 mrg mpfr_sub_q (y, x, z, MPFR_RNDU);
83 1.1 mrg CHECK_FOR ("nan", mpfr_nan_p (y));
84 1.1 mrg
85 1.1 mrg /* Exact value */
86 1.1 mrg mpfr_set_prec (x, 60);
87 1.1 mrg mpfr_set_prec (y, 60);
88 1.1 mrg mpfr_set_str1 (x, "0.5");
89 1.1 mrg mpz_set_str (mpq_numref (z), "3", 10);
90 1.1 mrg mpz_set_str (mpq_denref (z), "2", 10);
91 1.1 mrg res = mpfr_add_q (y, x, z, MPFR_RNDU);
92 1.1 mrg CHECK_FOR ("0.5+3/2", mpfr_cmp_ui(y, 2)==0 && res==0);
93 1.1 mrg res = mpfr_sub_q (y, x, z, MPFR_RNDU);
94 1.1 mrg CHECK_FOR ("0.5-3/2", mpfr_cmp_si(y, -1)==0 && res==0);
95 1.1 mrg
96 1.1 mrg /* Inf Rationnal */
97 1.1 mrg mpq_set_ui (z, 1, 0);
98 1.1 mrg mpfr_set_str1 (x, "0.5");
99 1.1 mrg res = mpfr_add_q (y, x, z, MPFR_RNDN);
100 1.1 mrg CHECK_FOR ("0.5+1/0", mpfr_inf_p (y) && MPFR_SIGN (y) > 0 && res == 0);
101 1.1 mrg res = mpfr_sub_q (y, x, z, MPFR_RNDN);
102 1.1 mrg CHECK_FOR ("0.5-1/0", mpfr_inf_p (y) && MPFR_SIGN (y) < 0 && res == 0);
103 1.1 mrg mpq_set_si (z, -1, 0);
104 1.1 mrg res = mpfr_add_q (y, x, z, MPFR_RNDN);
105 1.1 mrg CHECK_FOR ("0.5+ -1/0", mpfr_inf_p (y) && MPFR_SIGN (y) < 0 && res == 0);
106 1.1 mrg res = mpfr_sub_q (y, x, z, MPFR_RNDN);
107 1.1 mrg CHECK_FOR ("0.5- -1/0", mpfr_inf_p (y) && MPFR_SIGN (y) > 0 && res == 0);
108 1.1 mrg res = mpfr_div_q (y, x, z, MPFR_RNDN);
109 1.1 mrg CHECK_FOR ("0.5 / (-1/0)", mpfr_zero_p (y) && MPFR_SIGN (y) < 0 && res == 0);
110 1.1 mrg
111 1.1 mrg /* 0 */
112 1.1 mrg mpq_set_ui (z, 0, 1);
113 1.1 mrg mpfr_set_ui (x, 42, MPFR_RNDN);
114 1.1 mrg res = mpfr_add_q (y, x, z, MPFR_RNDN);
115 1.1 mrg CHECK_FOR ("42+0/1", mpfr_cmp_ui (y, 42) == 0 && res == 0);
116 1.1 mrg res = mpfr_sub_q (y, x, z, MPFR_RNDN);
117 1.1 mrg CHECK_FOR ("42-0/1", mpfr_cmp_ui (y, 42) == 0 && res == 0);
118 1.1 mrg res = mpfr_mul_q (y, x, z, MPFR_RNDN);
119 1.1 mrg CHECK_FOR ("42*0/1", mpfr_zero_p (y) && MPFR_SIGN (y) > 0 && res == 0);
120 1.1 mrg res = mpfr_div_q (y, x, z, MPFR_RNDN);
121 1.1 mrg CHECK_FOR ("42/(0/1)", mpfr_inf_p (y) && MPFR_SIGN (y) > 0 && res == 0);
122 1.1 mrg
123 1.1 mrg
124 1.1 mrg mpq_clear (z);
125 1.1 mrg mpfr_clear (x);
126 1.1 mrg mpfr_clear (y);
127 1.1 mrg }
128 1.1 mrg
129 1.1 mrg static void
130 1.1 mrg check_for_zero (void)
131 1.1 mrg {
132 1.1 mrg /* Check that 0 is unsigned! */
133 1.1 mrg mpq_t q;
134 1.1 mrg mpz_t z;
135 1.1 mrg mpfr_t x;
136 1.1 mrg int r;
137 1.1 mrg mpfr_sign_t i;
138 1.1 mrg
139 1.1 mrg mpfr_init (x);
140 1.1 mrg mpz_init (z);
141 1.1 mrg mpq_init (q);
142 1.1 mrg
143 1.1 mrg mpz_set_ui (z, 0);
144 1.1 mrg mpq_set_ui (q, 0, 1);
145 1.1 mrg
146 1.1 mrg MPFR_SET_ZERO (x);
147 1.1 mrg RND_LOOP (r)
148 1.1 mrg {
149 1.1 mrg for (i = MPFR_SIGN_NEG ; i <= MPFR_SIGN_POS ;
150 1.1 mrg i+=MPFR_SIGN_POS-MPFR_SIGN_NEG)
151 1.1 mrg {
152 1.1 mrg MPFR_SET_SIGN(x, i);
153 1.1 mrg mpfr_add_z (x, x, z, (mpfr_rnd_t) r);
154 1.1 mrg if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
155 1.1 mrg {
156 1.1 mrg printf("GMP Zero errors for add_z & rnd=%s & s=%d\n",
157 1.1 mrg mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
158 1.1 mrg mpfr_dump (x);
159 1.1 mrg exit (1);
160 1.1 mrg }
161 1.1 mrg mpfr_sub_z (x, x, z, (mpfr_rnd_t) r);
162 1.1 mrg if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
163 1.1 mrg {
164 1.1 mrg printf("GMP Zero errors for sub_z & rnd=%s & s=%d\n",
165 1.1 mrg mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
166 1.1 mrg mpfr_dump (x);
167 1.1 mrg exit (1);
168 1.1 mrg }
169 1.1 mrg mpfr_mul_z (x, x, z, (mpfr_rnd_t) r);
170 1.1 mrg if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
171 1.1 mrg {
172 1.1 mrg printf("GMP Zero errors for mul_z & rnd=%s & s=%d\n",
173 1.1 mrg mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
174 1.1 mrg mpfr_dump (x);
175 1.1 mrg exit (1);
176 1.1 mrg }
177 1.1 mrg mpfr_add_q (x, x, q, (mpfr_rnd_t) r);
178 1.1 mrg if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
179 1.1 mrg {
180 1.1 mrg printf("GMP Zero errors for add_q & rnd=%s & s=%d\n",
181 1.1 mrg mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
182 1.1 mrg mpfr_dump (x);
183 1.1 mrg exit (1);
184 1.1 mrg }
185 1.1 mrg mpfr_sub_q (x, x, q, (mpfr_rnd_t) r);
186 1.1 mrg if (!MPFR_IS_ZERO(x) || MPFR_SIGN(x)!=i)
187 1.1 mrg {
188 1.1 mrg printf("GMP Zero errors for sub_q & rnd=%s & s=%d\n",
189 1.1 mrg mpfr_print_rnd_mode ((mpfr_rnd_t) r), i);
190 1.1 mrg mpfr_dump (x);
191 1.1 mrg exit (1);
192 1.1 mrg }
193 1.1 mrg }
194 1.1 mrg }
195 1.1 mrg
196 1.1 mrg mpq_clear (q);
197 1.1 mrg mpz_clear (z);
198 1.1 mrg mpfr_clear (x);
199 1.1 mrg }
200 1.1 mrg
201 1.1 mrg static void
202 1.1 mrg test_cmp_z (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
203 1.1 mrg {
204 1.1 mrg mpfr_t x, z;
205 1.1 mrg mpz_t y;
206 1.1 mrg mpfr_prec_t p;
207 1.1 mrg int res1, res2;
208 1.1 mrg int n;
209 1.1 mrg
210 1.1 mrg mpfr_init (x);
211 1.1 mrg mpfr_init2 (z, MPFR_PREC_MIN);
212 1.1 mrg mpz_init (y);
213 1.1 mrg for(p=pmin ; p < pmax ; p++)
214 1.1 mrg {
215 1.1 mrg mpfr_set_prec (x, p);
216 1.1 mrg for ( n = 0; n < nmax ; n++)
217 1.1 mrg {
218 1.1 mrg mpfr_urandomb (x, RANDS);
219 1.1 mrg mpz_urandomb (y, RANDS, 1024);
220 1.1 mrg if (!MPFR_IS_SINGULAR (x))
221 1.1 mrg {
222 1.1 mrg mpfr_sub_z (z, x, y, MPFR_RNDN);
223 1.1 mrg res1 = mpfr_sgn (z);
224 1.1 mrg res2 = mpfr_cmp_z (x, y);
225 1.1 mrg if (res1 != res2)
226 1.1 mrg {
227 1.1 mrg printf("Error for mpfr_cmp_z: res=%d sub_z gives %d\n",
228 1.1 mrg res2, res1);
229 1.1 mrg exit (1);
230 1.1 mrg }
231 1.1 mrg }
232 1.1 mrg }
233 1.1 mrg }
234 1.1 mrg mpz_clear (y);
235 1.1 mrg mpfr_clear (x);
236 1.1 mrg mpfr_clear (z);
237 1.1 mrg }
238 1.1 mrg
239 1.1 mrg static void
240 1.1 mrg test_cmp_q (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
241 1.1 mrg {
242 1.1 mrg mpfr_t x, z;
243 1.1 mrg mpq_t y;
244 1.1 mrg mpfr_prec_t p;
245 1.1 mrg int res1, res2;
246 1.1 mrg int n;
247 1.1 mrg
248 1.1 mrg mpfr_init (x);
249 1.1 mrg mpfr_init2 (z, MPFR_PREC_MIN);
250 1.1 mrg mpq_init (y);
251 1.1 mrg for(p=pmin ; p < pmax ; p++)
252 1.1 mrg {
253 1.1 mrg mpfr_set_prec (x, p);
254 1.1 mrg for (n = 0 ; n < nmax ; n++)
255 1.1 mrg {
256 1.1 mrg mpfr_urandomb (x, RANDS);
257 1.1 mrg mpq_set_ui (y, randlimb (), randlimb() );
258 1.1 mrg if (!MPFR_IS_SINGULAR (x))
259 1.1 mrg {
260 1.1 mrg mpfr_sub_q (z, x, y, MPFR_RNDN);
261 1.1 mrg res1 = mpfr_sgn (z);
262 1.1 mrg res2 = mpfr_cmp_q (x, y);
263 1.1 mrg if (res1 != res2)
264 1.1 mrg {
265 1.1 mrg printf("Error for mpfr_cmp_q: res=%d sub_z gives %d\n",
266 1.1 mrg res2, res1);
267 1.1 mrg exit (1);
268 1.1 mrg }
269 1.1 mrg }
270 1.1 mrg }
271 1.1 mrg }
272 1.1 mrg mpq_clear (y);
273 1.1 mrg mpfr_clear (x);
274 1.1 mrg mpfr_clear (z);
275 1.1 mrg }
276 1.1 mrg
277 1.1 mrg static void
278 1.1 mrg test_cmp_f (mpfr_prec_t pmin, mpfr_prec_t pmax, int nmax)
279 1.1 mrg {
280 1.1 mrg mpfr_t x, z;
281 1.1 mrg mpf_t y;
282 1.1 mrg mpfr_prec_t p;
283 1.1 mrg int res1, res2;
284 1.1 mrg int n;
285 1.1 mrg
286 1.1 mrg mpfr_init (x);
287 1.1 mrg mpfr_init2 (z, pmax+GMP_NUMB_BITS);
288 1.1 mrg mpf_init2 (y, MPFR_PREC_MIN);
289 1.1 mrg for(p=pmin ; p < pmax ; p+=3)
290 1.1 mrg {
291 1.1 mrg mpfr_set_prec (x, p);
292 1.1 mrg mpf_set_prec (y, p);
293 1.1 mrg for ( n = 0; n < nmax ; n++)
294 1.1 mrg {
295 1.1 mrg mpfr_urandomb (x, RANDS);
296 1.1 mrg mpf_urandomb (y, RANDS, p);
297 1.1 mrg if (!MPFR_IS_SINGULAR (x))
298 1.1 mrg {
299 1.1 mrg mpfr_set_f (z, y, MPFR_RNDN);
300 1.1 mrg mpfr_sub (z, x, z, MPFR_RNDN);
301 1.1 mrg res1 = mpfr_sgn (z);
302 1.1 mrg res2 = mpfr_cmp_f (x, y);
303 1.1 mrg if (res1 != res2)
304 1.1 mrg {
305 1.1 mrg printf("Error for mpfr_cmp_f: res=%d sub gives %d\n",
306 1.1 mrg res2, res1);
307 1.1 mrg exit (1);
308 1.1 mrg }
309 1.1 mrg }
310 1.1 mrg }
311 1.1 mrg }
312 1.1 mrg mpf_clear (y);
313 1.1 mrg mpfr_clear (x);
314 1.1 mrg mpfr_clear (z);
315 1.1 mrg }
316 1.1 mrg
317 1.1 mrg static void
318 1.1 mrg test_specialz (int (*mpfr_func)(mpfr_ptr, mpfr_srcptr, mpz_srcptr, mpfr_rnd_t),
319 1.1 mrg void (*mpz_func)(mpz_ptr, mpz_srcptr, mpz_srcptr),
320 1.1 mrg const char *op)
321 1.1 mrg {
322 1.1 mrg mpfr_t x1, x2;
323 1.1 mrg mpz_t z1, z2;
324 1.1 mrg int res;
325 1.1 mrg
326 1.1 mrg mpfr_inits2 (128, x1, x2, (mpfr_ptr) 0);
327 1.1 mrg mpz_init (z1); mpz_init(z2);
328 1.1 mrg mpz_fac_ui (z1, 19); /* 19!+1 fits perfectly in a 128 bits mantissa */
329 1.1 mrg mpz_add_ui (z1, z1, 1);
330 1.1 mrg mpz_fac_ui (z2, 20); /* 20!+1 fits perfectly in a 128 bits mantissa */
331 1.1 mrg mpz_add_ui (z2, z2, 1);
332 1.1 mrg
333 1.1 mrg res = mpfr_set_z(x1, z1, MPFR_RNDN);
334 1.1 mrg if (res)
335 1.1 mrg {
336 1.1 mrg printf("Specialz %s: set_z1 error\n", op);
337 1.1 mrg exit(1);
338 1.1 mrg }
339 1.1 mrg mpfr_set_z (x2, z2, MPFR_RNDN);
340 1.1 mrg if (res)
341 1.1 mrg {
342 1.1 mrg printf("Specialz %s: set_z2 error\n", op);
343 1.1 mrg exit(1);
344 1.1 mrg }
345 1.1 mrg
346 1.1 mrg /* (19!+1) * (20!+1) fits in a 128 bits number */
347 1.1 mrg res = mpfr_func(x1, x1, z2, MPFR_RNDN);
348 1.1 mrg if (res)
349 1.1 mrg {
350 1.1 mrg printf("Specialz %s: wrong inexact flag.\n", op);
351 1.1 mrg exit(1);
352 1.1 mrg }
353 1.1 mrg mpz_func(z1, z1, z2);
354 1.1 mrg res = mpfr_set_z (x2, z1, MPFR_RNDN);
355 1.1 mrg if (res)
356 1.1 mrg {
357 1.1 mrg printf("Specialz %s: set_z2 error\n", op);
358 1.1 mrg exit(1);
359 1.1 mrg }
360 1.1 mrg if (mpfr_cmp(x1, x2))
361 1.1 mrg {
362 1.1 mrg printf("Specialz %s: results differ.\nx1=", op);
363 1.1 mrg mpfr_print_binary(x1);
364 1.1 mrg printf("\nx2=");
365 1.1 mrg mpfr_print_binary(x2);
366 1.1 mrg printf ("\nZ2=");
367 1.1 mrg mpz_out_str (stdout, 2, z1);
368 1.1 mrg putchar('\n');
369 1.1 mrg exit(1);
370 1.1 mrg }
371 1.1 mrg
372 1.1 mrg mpz_set_ui (z1, 1);
373 1.1 mrg mpz_set_ui (z2, 0);
374 1.1 mrg mpfr_set_ui (x1, 1, MPFR_RNDN);
375 1.1 mrg mpz_func (z1, z1, z2);
376 1.1 mrg res = mpfr_func(x1, x1, z2, MPFR_RNDN);
377 1.1 mrg mpfr_set_z (x2, z1, MPFR_RNDN);
378 1.1 mrg if (mpfr_cmp(x1, x2))
379 1.1 mrg {
380 1.1 mrg printf("Specialz %s: results differ(2).\nx1=", op);
381 1.1 mrg mpfr_print_binary(x1);
382 1.1 mrg printf("\nx2=");
383 1.1 mrg mpfr_print_binary(x2);
384 1.1 mrg putchar('\n');
385 1.1 mrg exit(1);
386 1.1 mrg }
387 1.1 mrg
388 1.1 mrg mpz_clear (z1); mpz_clear(z2);
389 1.1 mrg mpfr_clears (x1, x2, (mpfr_ptr) 0);
390 1.1 mrg }
391 1.1 mrg
392 1.1 mrg static void
393 1.1 mrg test_genericz (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N,
394 1.1 mrg int (*func)(mpfr_ptr, mpfr_srcptr, mpz_srcptr, mpfr_rnd_t),
395 1.1 mrg const char *op)
396 1.1 mrg {
397 1.1 mrg mpfr_prec_t prec;
398 1.1 mrg mpfr_t arg1, dst_big, dst_small, tmp;
399 1.1 mrg mpz_t arg2;
400 1.1 mrg mpfr_rnd_t rnd;
401 1.1 mrg int inexact, compare, compare2;
402 1.1 mrg unsigned int n;
403 1.1 mrg
404 1.1 mrg mpfr_inits (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
405 1.1 mrg mpz_init (arg2);
406 1.1 mrg
407 1.1 mrg for (prec = p0; prec <= p1; prec++)
408 1.1 mrg {
409 1.1 mrg mpfr_set_prec (arg1, prec);
410 1.1 mrg mpfr_set_prec (tmp, prec);
411 1.1 mrg mpfr_set_prec (dst_small, prec);
412 1.1 mrg
413 1.1 mrg for (n=0; n<N; n++)
414 1.1 mrg {
415 1.1 mrg mpfr_urandomb (arg1, RANDS);
416 1.1 mrg mpz_urandomb (arg2, RANDS, 1024);
417 1.1 mrg rnd = RND_RAND ();
418 1.1 mrg mpfr_set_prec (dst_big, 2*prec);
419 1.1 mrg compare = func(dst_big, arg1, arg2, rnd);
420 1.1 mrg if (mpfr_can_round (dst_big, 2*prec, rnd, rnd, prec))
421 1.1 mrg {
422 1.1 mrg mpfr_set (tmp, dst_big, rnd);
423 1.1 mrg inexact = func(dst_small, arg1, arg2, rnd);
424 1.1 mrg if (mpfr_cmp (tmp, dst_small))
425 1.1 mrg {
426 1.1 mrg printf ("Results differ for prec=%u rnd_mode=%s and %s_z:\n"
427 1.1 mrg "arg1=",
428 1.1 mrg (unsigned) prec, mpfr_print_rnd_mode (rnd), op);
429 1.1 mrg mpfr_print_binary (arg1);
430 1.1 mrg printf("\narg2=");
431 1.1 mrg mpz_out_str (stdout, 2, arg2);
432 1.1 mrg printf ("\ngot ");
433 1.1 mrg mpfr_dump (dst_small);
434 1.1 mrg printf ("expected ");
435 1.1 mrg mpfr_dump (tmp);
436 1.1 mrg printf ("approx ");
437 1.1 mrg mpfr_dump (dst_big);
438 1.1 mrg exit (1);
439 1.1 mrg }
440 1.1 mrg compare2 = mpfr_cmp (tmp, dst_big);
441 1.1 mrg /* if rounding to nearest, cannot know the sign of t - f(x)
442 1.1 mrg because of composed rounding: y = o(f(x)) and t = o(y) */
443 1.1 mrg if (compare * compare2 >= 0)
444 1.1 mrg compare = compare + compare2;
445 1.1 mrg else
446 1.1 mrg compare = inexact; /* cannot determine sign(t-f(x)) */
447 1.1 mrg if (((inexact == 0) && (compare != 0)) ||
448 1.1 mrg ((inexact > 0) && (compare <= 0)) ||
449 1.1 mrg ((inexact < 0) && (compare >= 0)))
450 1.1 mrg {
451 1.1 mrg printf ("Wrong inexact flag for rnd=%s and %s_z:\n"
452 1.1 mrg "expected %d, got %d\n",
453 1.1 mrg mpfr_print_rnd_mode (rnd), op, compare, inexact);
454 1.1 mrg printf ("\narg1="); mpfr_print_binary (arg1);
455 1.1 mrg printf ("\narg2="); mpz_out_str(stdout, 2, arg2);
456 1.1 mrg printf ("\ndstl="); mpfr_print_binary (dst_big);
457 1.1 mrg printf ("\ndsts="); mpfr_print_binary (dst_small);
458 1.1 mrg printf ("\ntmp ="); mpfr_dump (tmp);
459 1.1 mrg exit (1);
460 1.1 mrg }
461 1.1 mrg }
462 1.1 mrg }
463 1.1 mrg }
464 1.1 mrg
465 1.1 mrg mpz_clear (arg2);
466 1.1 mrg mpfr_clears (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
467 1.1 mrg }
468 1.1 mrg
469 1.1 mrg static void
470 1.1 mrg test_genericq (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N,
471 1.1 mrg int (*func)(mpfr_ptr, mpfr_srcptr, mpq_srcptr, mpfr_rnd_t),
472 1.1 mrg const char *op)
473 1.1 mrg {
474 1.1 mrg mpfr_prec_t prec;
475 1.1 mrg mpfr_t arg1, dst_big, dst_small, tmp;
476 1.1 mrg mpq_t arg2;
477 1.1 mrg mpfr_rnd_t rnd;
478 1.1 mrg int inexact, compare, compare2;
479 1.1 mrg unsigned int n;
480 1.1 mrg
481 1.1 mrg mpfr_inits (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
482 1.1 mrg mpq_init (arg2);
483 1.1 mrg
484 1.1 mrg for (prec = p0; prec <= p1; prec++)
485 1.1 mrg {
486 1.1 mrg mpfr_set_prec (arg1, prec);
487 1.1 mrg mpfr_set_prec (tmp, prec);
488 1.1 mrg mpfr_set_prec (dst_small, prec);
489 1.1 mrg
490 1.1 mrg for (n=0; n<N; n++)
491 1.1 mrg {
492 1.1 mrg mpfr_urandomb (arg1, RANDS);
493 1.1 mrg mpq_set_ui (arg2, randlimb (), randlimb() );
494 1.1 mrg mpq_canonicalize (arg2);
495 1.1 mrg rnd = RND_RAND ();
496 1.1 mrg mpfr_set_prec (dst_big, prec+10);
497 1.1 mrg compare = func(dst_big, arg1, arg2, rnd);
498 1.1 mrg if (mpfr_can_round (dst_big, prec+10, rnd, rnd, prec))
499 1.1 mrg {
500 1.1 mrg mpfr_set (tmp, dst_big, rnd);
501 1.1 mrg inexact = func(dst_small, arg1, arg2, rnd);
502 1.1 mrg if (mpfr_cmp (tmp, dst_small))
503 1.1 mrg {
504 1.1 mrg printf ("Results differ for prec=%u rnd_mode=%s and %s_q:\n"
505 1.1 mrg "arg1=",
506 1.1 mrg (unsigned) prec, mpfr_print_rnd_mode (rnd), op);
507 1.1 mrg mpfr_print_binary (arg1);
508 1.1 mrg printf("\narg2=");
509 1.1 mrg mpq_out_str(stdout, 2, arg2);
510 1.1 mrg printf ("\ngot ");
511 1.1 mrg mpfr_print_binary (dst_small);
512 1.1 mrg printf ("\nexpected ");
513 1.1 mrg mpfr_print_binary (tmp);
514 1.1 mrg printf ("\napprox ");
515 1.1 mrg mpfr_print_binary (dst_big);
516 1.1 mrg putchar('\n');
517 1.1 mrg exit (1);
518 1.1 mrg }
519 1.1 mrg compare2 = mpfr_cmp (tmp, dst_big);
520 1.1 mrg /* if rounding to nearest, cannot know the sign of t - f(x)
521 1.1 mrg because of composed rounding: y = o(f(x)) and t = o(y) */
522 1.1 mrg if (compare * compare2 >= 0)
523 1.1 mrg compare = compare + compare2;
524 1.1 mrg else
525 1.1 mrg compare = inexact; /* cannot determine sign(t-f(x)) */
526 1.1 mrg if (((inexact == 0) && (compare != 0)) ||
527 1.1 mrg ((inexact > 0) && (compare <= 0)) ||
528 1.1 mrg ((inexact < 0) && (compare >= 0)))
529 1.1 mrg {
530 1.1 mrg printf ("Wrong inexact flag for rnd=%s and %s_q:\n"
531 1.1 mrg "expected %d, got %d",
532 1.1 mrg mpfr_print_rnd_mode (rnd), op, compare, inexact);
533 1.1 mrg printf ("\narg1="); mpfr_print_binary (arg1);
534 1.1 mrg printf ("\narg2="); mpq_out_str(stdout, 2, arg2);
535 1.1 mrg printf ("\ndstl="); mpfr_print_binary (dst_big);
536 1.1 mrg printf ("\ndsts="); mpfr_print_binary (dst_small);
537 1.1 mrg printf ("\ntmp ="); mpfr_print_binary (tmp);
538 1.1 mrg putchar('\n');
539 1.1 mrg exit (1);
540 1.1 mrg }
541 1.1 mrg }
542 1.1 mrg }
543 1.1 mrg }
544 1.1 mrg
545 1.1 mrg mpq_clear (arg2);
546 1.1 mrg mpfr_clears (arg1, dst_big, dst_small, tmp, (mpfr_ptr) 0);
547 1.1 mrg }
548 1.1 mrg
549 1.1 mrg static void
550 1.1 mrg test_specialq (mpfr_prec_t p0, mpfr_prec_t p1, unsigned int N,
551 1.1 mrg int (*mpfr_func)(mpfr_ptr, mpfr_srcptr, mpq_srcptr, mpfr_rnd_t),
552 1.1 mrg void (*mpq_func)(mpq_ptr, mpq_srcptr, mpq_srcptr),
553 1.1 mrg const char *op)
554 1.1 mrg {
555 1.1 mrg mpfr_t fra, frb, frq;
556 1.1 mrg mpq_t q1, q2, qr;
557 1.1 mrg unsigned int n;
558 1.1 mrg mpfr_prec_t prec;
559 1.1 mrg
560 1.1 mrg for (prec = p0 ; prec < p1 ; prec++)
561 1.1 mrg {
562 1.1 mrg mpfr_inits2 (prec, fra, frb, frq, (mpfr_ptr) 0);
563 1.1 mrg mpq_init (q1); mpq_init(q2); mpq_init (qr);
564 1.1 mrg
565 1.1 mrg for( n = 0 ; n < N ; n++)
566 1.1 mrg {
567 1.1 mrg mpq_set_ui(q1, randlimb(), randlimb() );
568 1.1 mrg mpq_set_ui(q2, randlimb(), randlimb() );
569 1.1 mrg mpq_canonicalize (q1);
570 1.1 mrg mpq_canonicalize (q2);
571 1.1 mrg mpq_func (qr, q1, q2);
572 1.1 mrg mpfr_set_q (fra, q1, MPFR_RNDD);
573 1.1 mrg mpfr_func (fra, fra, q2, MPFR_RNDD);
574 1.1 mrg mpfr_set_q (frb, q1, MPFR_RNDU);
575 1.1 mrg mpfr_func (frb, frb, q2, MPFR_RNDU);
576 1.1 mrg mpfr_set_q (frq, qr, MPFR_RNDN);
577 1.1 mrg /* We should have fra <= qr <= frb */
578 1.1 mrg if ( (mpfr_cmp(fra, frq) > 0) || (mpfr_cmp (frq, frb) > 0))
579 1.1 mrg {
580 1.1 mrg printf("Range error for prec=%lu and %s", prec, op);
581 1.1 mrg printf ("\nq1="); mpq_out_str(stdout, 2, q1);
582 1.1 mrg printf ("\nq2="); mpq_out_str(stdout, 2, q2);
583 1.1 mrg printf ("\nfr_dn="); mpfr_print_binary (fra);
584 1.1 mrg printf ("\nfr_q ="); mpfr_print_binary (frq);
585 1.1 mrg printf ("\nfr_up="); mpfr_print_binary (frb);
586 1.1 mrg putchar('\n');
587 1.1 mrg exit (1);
588 1.1 mrg }
589 1.1 mrg }
590 1.1 mrg
591 1.1 mrg mpq_clear (q1); mpq_clear (q2); mpq_clear (qr);
592 1.1 mrg mpfr_clears (fra, frb, frq, (mpfr_ptr) 0);
593 1.1 mrg }
594 1.1 mrg }
595 1.1 mrg
596 1.1 mrg int
597 1.1 mrg main (int argc, char *argv[])
598 1.1 mrg {
599 1.1 mrg tests_start_mpfr ();
600 1.1 mrg
601 1.1 mrg special ();
602 1.1 mrg
603 1.1 mrg test_specialz (mpfr_add_z, mpz_add, "add");
604 1.1 mrg test_specialz (mpfr_sub_z, mpz_sub, "sub");
605 1.1 mrg test_specialz (mpfr_mul_z, mpz_mul, "mul");
606 1.1 mrg test_genericz (2, 100, 100, mpfr_add_z, "add");
607 1.1 mrg test_genericz (2, 100, 100, mpfr_sub_z, "sub");
608 1.1 mrg test_genericz (2, 100, 100, mpfr_mul_z, "mul");
609 1.1 mrg test_genericz (2, 100, 100, mpfr_div_z, "div");
610 1.1 mrg
611 1.1 mrg test_genericq (2, 100, 100, mpfr_add_q, "add");
612 1.1 mrg test_genericq (2, 100, 100, mpfr_sub_q, "sub");
613 1.1 mrg test_genericq (2, 100, 100, mpfr_mul_q, "mul");
614 1.1 mrg test_genericq (2, 100, 100, mpfr_div_q, "div");
615 1.1 mrg test_specialq (2, 100, 100, mpfr_mul_q, mpq_mul, "mul");
616 1.1 mrg test_specialq (2, 100, 100, mpfr_div_q, mpq_div, "div");
617 1.1 mrg test_specialq (2, 100, 100, mpfr_add_q, mpq_add, "add");
618 1.1 mrg test_specialq (2, 100, 100, mpfr_sub_q, mpq_sub, "sub");
619 1.1 mrg
620 1.1 mrg test_cmp_z (2, 100, 100);
621 1.1 mrg test_cmp_q (2, 100, 100);
622 1.1 mrg test_cmp_f (2, 100, 100);
623 1.1 mrg
624 1.1 mrg check_for_zero ();
625 1.1 mrg
626 1.1 mrg tests_end_mpfr ();
627 1.1 mrg return 0;
628 1.1 mrg }
629 1.1 mrg
630