1 /* Test file for mpfr_sqr.
2
3 Copyright 2004-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 #define TEST_FUNCTION mpfr_sqr
26 #include "tgeneric.c"
27
28 static int
29 inexact_sign (int x)
30 {
31 return (x < 0) ? -1 : (x > 0);
32 }
33
34 static void
35 error1 (mpfr_rnd_t rnd, mpfr_prec_t prec,
36 mpfr_t in, mpfr_ptr outmul, mpfr_ptr outsqr)
37 {
38 printf("ERROR: for %s and prec=%lu\nINPUT=", mpfr_print_rnd_mode(rnd),
39 (unsigned long) prec);
40 mpfr_dump(in);
41 printf("OutputMul="); mpfr_dump(outmul);
42 printf("OutputSqr="); mpfr_dump(outsqr);
43 exit(1);
44 }
45
46 static void
47 error2 (mpfr_rnd_t rnd, mpfr_prec_t prec, mpfr_ptr in, mpfr_ptr out,
48 int inexactmul, int inexactsqr)
49 {
50 printf("ERROR: for %s and prec=%lu\nINPUT=", mpfr_print_rnd_mode(rnd),
51 (unsigned long) prec);
52 mpfr_dump(in);
53 printf("Output="); mpfr_dump(out);
54 printf("InexactMul= %d InexactSqr= %d\n", inexactmul, inexactsqr);
55 exit(1);
56 }
57
58 static void
59 check_random (mpfr_prec_t p)
60 {
61 mpfr_t x,y,z;
62 int r;
63 int i, inexact1, inexact2;
64
65 mpfr_inits2 (p, x, y, z, (mpfr_ptr) 0);
66 for(i = 0 ; i < 500 ; i++)
67 {
68 mpfr_urandomb (x, RANDS);
69 if (MPFR_IS_PURE_FP(x))
70 RND_LOOP_NO_RNDF (r)
71 {
72 /* the following call to mpfr_mul with identical arguments is
73 intentional (to compare with mpfr_sqr) */
74 inexact1 = mpfr_mul (y, x, x, (mpfr_rnd_t) r);
75 inexact2 = mpfr_sqr (z, x, (mpfr_rnd_t) r);
76 if (mpfr_cmp (y, z))
77 error1 ((mpfr_rnd_t) r,p,x,y,z);
78 if (inexact_sign (inexact1) != inexact_sign (inexact2))
79 error2 ((mpfr_rnd_t) r,p,x,y,inexact1,inexact2);
80 }
81 }
82 mpfr_clears (x, y, z, (mpfr_ptr) 0);
83 }
84
85 static void
86 check_special (void)
87 {
88 mpfr_t x, y;
89 mpfr_exp_t emin;
90
91 mpfr_init (x);
92 mpfr_init (y);
93
94 mpfr_set_nan (x);
95 mpfr_sqr (y, x, MPFR_RNDN);
96 MPFR_ASSERTN (mpfr_nan_p (y));
97
98 mpfr_set_inf (x, 1);
99 mpfr_sqr (y, x, MPFR_RNDN);
100 MPFR_ASSERTN (mpfr_inf_p (y) && mpfr_sgn (y) > 0);
101
102 mpfr_set_inf (x, -1);
103 mpfr_sqr (y, x, MPFR_RNDN);
104 MPFR_ASSERTN (mpfr_inf_p (y) && mpfr_sgn (y) > 0);
105
106 mpfr_set_ui (x, 0, MPFR_RNDN);
107 mpfr_sqr (y, x, MPFR_RNDN);
108 MPFR_ASSERTN (mpfr_zero_p (y));
109
110 emin = mpfr_get_emin ();
111 set_emin (0);
112 mpfr_set_ui (x, 1, MPFR_RNDN);
113 mpfr_div_2ui (x, x, 1, MPFR_RNDN);
114 MPFR_ASSERTN (!mpfr_zero_p (x));
115 mpfr_sqr (y, x, MPFR_RNDN);
116 MPFR_ASSERTN (mpfr_zero_p (y));
117 set_emin (emin);
118
119 mpfr_clear (y);
120 mpfr_clear (x);
121 }
122
123 /* check ax < __gmpfr_emin with rnd_mode == MPFR_RNDN, rb = 0 and sb <> 0 */
124 static void
125 test_underflow (void)
126 {
127 mpfr_t x, y;
128 mpfr_exp_t emin;
129
130 emin = mpfr_get_emin ();
131 set_emin (0);
132
133 mpfr_init2 (x, 24);
134 mpfr_init2 (y, 24);
135
136 mpfr_set_ui_2exp (x, 11863283, -24, MPFR_RNDN);
137 /* x^2 = 0.011111111111111111111111101101100111111010101001*2^(-48)
138 thus we have an underflow */
139 mpfr_clear_underflow ();
140 mpfr_sqr (y, x, MPFR_RNDN);
141 MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 1, -1) == 0);
142 MPFR_ASSERTN(mpfr_underflow_p ());
143
144 mpfr_set_prec (x, 69);
145 mpfr_set_prec (y, 69);
146 mpfr_set_str_binary (x, "0.101101010000010011110011001100111111100111011110011001001000010001011");
147 mpfr_clear_underflow ();
148 mpfr_sqr (y, x, MPFR_RNDN);
149 MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 1, -1) == 0);
150 MPFR_ASSERTN(mpfr_underflow_p ());
151
152 mpfr_clear (y);
153 mpfr_clear (x);
154
155 set_emin (emin);
156 }
157
158 /* Test of a bug seen with GCC 4.5.2 and GMP 5.0.1 on m68k (m68000 target).
159 https://sympa.inria.fr/sympa/arc/mpfr/2011-05/msg00003.html
160 https://sympa.inria.fr/sympa/arc/mpfr/2011-05/msg00041.html
161 */
162 static void
163 check_mpn_sqr (void)
164 {
165 #if GMP_NUMB_BITS == 32 && __GNU_MP_VERSION >= 5
166 /* Note: since we test a low-level bug, src is initialized
167 without using a GMP function, just in case. */
168 mp_limb_t src[5] =
169 { 0x90000000, 0xbaa55f4f, 0x2cbec4d9, 0xfcef3242, 0xda827999 };
170 mp_limb_t exd[10] =
171 { 0x00000000, 0x31000000, 0xbd4bc59a, 0x41fbe2b5, 0x33471e7e,
172 0x90e826a7, 0xbaa55f4f, 0x2cbec4d9, 0xfcef3242, 0xba827999 };
173 mp_limb_t dst[10];
174 int i;
175
176 mpn_sqr (dst, src, 5); /* new in GMP 5 */
177 for (i = 0; i < 10; i++)
178 {
179 if (dst[i] != exd[i])
180 {
181 printf ("Error in check_mpn_sqr\n");
182 printf ("exd[%d] = 0x%08lx\n", i, (unsigned long) exd[i]);
183 printf ("dst[%d] = 0x%08lx\n", i, (unsigned long) dst[i]);
184 printf ("Note: This is not a bug in MPFR, but a bug in"
185 " either GMP or, more\nprobably, in the compiler."
186 " It may cause other tests to fail.\n");
187 exit (1);
188 }
189 }
190 #endif
191 }
192
193 static void
194 coverage (mpfr_prec_t pmax)
195 {
196 mpfr_prec_t p;
197
198 for (p = MPFR_PREC_MIN; p <= pmax; p++)
199 {
200 mpfr_t a, b, c;
201 int inex;
202 mpfr_exp_t emin;
203
204 mpfr_init2 (a, p);
205 mpfr_init2 (b, p);
206 mpfr_init2 (c, p);
207
208 /* exercise carry in most significant bits of a, with overflow */
209 mpfr_set_ui_2exp (b, 1, mpfr_get_emax (), MPFR_RNDZ);
210 mpfr_sqrt (b, b, MPFR_RNDU);
211 mpfr_div_2ui (c, b, 1, MPFR_RNDN);
212 mpfr_sqr (c, c, MPFR_RNDN);
213 mpfr_clear_flags ();
214 inex = mpfr_sqr (a, b, MPFR_RNDN);
215 /* if EXP(c) > emax-2, there is overflow */
216 if (mpfr_get_exp (c) > mpfr_get_emax () - 2)
217 {
218 MPFR_ASSERTN(inex > 0);
219 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0);
220 MPFR_ASSERTN(mpfr_overflow_p ());
221 }
222 else /* no overflow */
223 {
224 /* 2^p-1 is a square only for p=1 */
225 MPFR_ASSERTN((p == 1 && inex == 0) || (p > 1 && inex < 0));
226 MPFR_ASSERTN(!mpfr_overflow_p ());
227 mpfr_set_ui_2exp (c, 1, mpfr_get_emax (), MPFR_RNDZ);
228 MPFR_ASSERTN(mpfr_equal_p (a, c));
229 }
230
231 /* same as above, with RNDU */
232 mpfr_set_ui_2exp (b, 1, mpfr_get_emax (), MPFR_RNDZ);
233 mpfr_sqrt (b, b, MPFR_RNDU);
234 mpfr_div_2ui (c, b, 1, MPFR_RNDN);
235 mpfr_sqr (c, c, MPFR_RNDU);
236 mpfr_clear_flags ();
237 inex = mpfr_sqr (a, b, MPFR_RNDU);
238 /* if EXP(c) > emax-2, there is overflow */
239 if (mpfr_get_exp (c) > mpfr_get_emax () - 2)
240 {
241 MPFR_ASSERTN(inex > 0);
242 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0);
243 MPFR_ASSERTN(mpfr_overflow_p ());
244 }
245 else /* no overflow */
246 {
247 /* 2^p-1 is a square only for p=1 */
248 MPFR_ASSERTN((p == 1 && inex == 0) || (p > 1 && inex < 0));
249 MPFR_ASSERTN(!mpfr_overflow_p ());
250 mpfr_set_ui_2exp (c, 1, mpfr_get_emax (), MPFR_RNDZ);
251 MPFR_ASSERTN(mpfr_equal_p (a, c));
252 }
253
254 /* exercise trivial overflow */
255 mpfr_set_ui_2exp (b, 1, mpfr_get_emax (), MPFR_RNDZ);
256 mpfr_sqrt (b, b, MPFR_RNDU);
257 mpfr_mul_2ui (b, b, 1, MPFR_RNDN);
258 mpfr_clear_flags ();
259 inex = mpfr_sqr (a, b, MPFR_RNDN);
260 MPFR_ASSERTN(inex > 0);
261 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0);
262 MPFR_ASSERTN(mpfr_overflow_p ());
263
264 /* exercise trivial underflow */
265 mpfr_set_ui_2exp (b, 1, mpfr_get_emin () - 1, MPFR_RNDZ);
266 mpfr_sqrt (b, b, MPFR_RNDU);
267 mpfr_div_2ui (b, b, 1, MPFR_RNDN);
268 mpfr_clear_flags ();
269 inex = mpfr_sqr (a, b, MPFR_RNDN);
270 MPFR_ASSERTN(inex < 0);
271 MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0);
272 MPFR_ASSERTN(mpfr_underflow_p ());
273
274 /* exercise square between 0.5*2^emin and its predecessor (emin even) */
275 emin = mpfr_get_emin ();
276 set_emin (emin + (emin & 1)); /* now emin is even */
277 mpfr_set_ui_2exp (b, 1, mpfr_get_emin () - 1, MPFR_RNDN);
278 inex = mpfr_sqrt (b, b, MPFR_RNDZ);
279 MPFR_ASSERTN(inex != 0); /* sqrt(2) is not exact */
280 mpfr_mul_2ui (c, b, 1, MPFR_RNDN);
281 mpfr_sqr (c, c, MPFR_RNDN);
282 mpfr_clear_flags ();
283 inex = mpfr_sqr (a, b, MPFR_RNDN);
284 if (mpfr_get_exp (c) < mpfr_get_emin () + 2) /* underflow */
285 {
286 /* if c > 0.5*2^(emin+1), we should round to 0.5*2^emin */
287 if (mpfr_cmp_ui_2exp (c, 1, mpfr_get_emin ()) > 0)
288 {
289 MPFR_ASSERTN(inex > 0);
290 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0);
291 MPFR_ASSERTN(mpfr_underflow_p ());
292 }
293 else /* we should round to 0 */
294 {
295 MPFR_ASSERTN(inex < 0);
296 MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0);
297 MPFR_ASSERTN(mpfr_underflow_p ());
298 }
299 }
300 else
301 {
302 MPFR_ASSERTN(inex > 0);
303 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0);
304 MPFR_ASSERTN(!mpfr_underflow_p ());
305 }
306 set_emin (emin);
307
308 /* exercise exact square root 2^(emin-2) for emin even */
309 emin = mpfr_get_emin ();
310 set_emin (emin + (emin & 1)); /* now emin is even */
311 mpfr_set_ui_2exp (b, 1, (mpfr_get_emin () - 2) / 2, MPFR_RNDN);
312 inex = mpfr_sqr (a, b, MPFR_RNDN);
313 MPFR_ASSERTN(inex < 0);
314 MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0);
315 MPFR_ASSERTN(mpfr_underflow_p ());
316 set_emin (emin);
317
318 /* same as above, for RNDU */
319 emin = mpfr_get_emin ();
320 set_emin (emin + (emin & 1)); /* now emin is even */
321 mpfr_set_ui_2exp (b, 1, mpfr_get_emin () - 1, MPFR_RNDN);
322 inex = mpfr_sqrt (b, b, MPFR_RNDZ);
323 MPFR_ASSERTN(inex != 0); /* sqrt(2) is not exact */
324 mpfr_mul_2ui (c, b, 1, MPFR_RNDN);
325 mpfr_sqr (c, c, MPFR_RNDU);
326 mpfr_clear_flags ();
327 inex = mpfr_sqr (a, b, MPFR_RNDU);
328 MPFR_ASSERTN(inex > 0);
329 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0);
330 /* we have underflow if c < 2^(emin+1) */
331 if (mpfr_cmp_ui_2exp (c, 1, mpfr_get_emin () + 1) < 0)
332 MPFR_ASSERTN(mpfr_underflow_p ());
333 else
334 MPFR_ASSERTN(!mpfr_underflow_p ());
335 set_emin (emin);
336
337 mpfr_clear (a);
338 mpfr_clear (b);
339 mpfr_clear (c);
340 }
341 }
342
343 int
344 main (void)
345 {
346 mpfr_prec_t p;
347
348 tests_start_mpfr ();
349
350 coverage (1024);
351 check_mpn_sqr ();
352 check_special ();
353 test_underflow ();
354
355 for (p = MPFR_PREC_MIN; p < 200; p++)
356 check_random (p);
357
358 test_generic (MPFR_PREC_MIN, 200, 15);
359 data_check ("data/sqr", mpfr_sqr, "mpfr_sqr");
360 bad_cases (mpfr_sqr, mpfr_sqrt, "mpfr_sqr", 8, -256, 255, 4, 128, 800, 50);
361
362 tests_end_mpfr ();
363 return 0;
364 }