1 /* Test file for mpfr_rec_sqrt.
2
3 Copyright 2008-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 <time.h>
24 #include "mpfr-test.h"
25
26 #define TEST_FUNCTION mpfr_rec_sqrt
27 #define TEST_RANDOM_POS 8 /* 8/512 = 1/64 of the tested numbers are negative */
28 #include "tgeneric.c"
29
30 static void
31 special (void)
32 {
33 mpfr_t x, y;
34 int inex;
35
36 mpfr_init (x);
37 mpfr_init (y);
38
39 /* rec_sqrt(NaN) = NaN */
40 mpfr_set_nan (x);
41 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN);
42 MPFR_ASSERTN(mpfr_nan_p (x) && inex == 0);
43
44 /* rec_sqrt(+Inf) = +0 */
45 mpfr_set_inf (x, 1);
46 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN);
47 MPFR_ASSERTN(mpfr_zero_p (x) && MPFR_IS_POS(x) && inex == 0);
48
49 /* rec_sqrt(-Inf) = NaN */
50 mpfr_set_inf (x, -1);
51 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN);
52 MPFR_ASSERTN(mpfr_nan_p (x) && inex == 0);
53
54 /* rec_sqrt(+0) = +Inf */
55 mpfr_set_ui (x, 0, MPFR_RNDN);
56 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN);
57 MPFR_ASSERTN(mpfr_inf_p (x) && MPFR_IS_POS(x) && inex == 0);
58
59 /* rec_sqrt(-0) = +Inf */
60 mpfr_set_ui (x, 0, MPFR_RNDN);
61 mpfr_neg (x, x, MPFR_RNDN);
62 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN);
63 MPFR_ASSERTN(mpfr_inf_p (x) && MPFR_IS_POS(x) && inex == 0);
64
65 /* rec_sqrt(-1) = NaN */
66 mpfr_set_si (x, -1, MPFR_RNDN);
67 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN);
68 MPFR_ASSERTN(mpfr_nan_p (x) && inex == 0);
69
70 /* rec_sqrt(1) = 1 */
71 mpfr_set_ui (x, 1, MPFR_RNDN);
72 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN);
73 MPFR_ASSERTN((mpfr_cmp_ui (x, 1) == 0) && (inex == 0));
74
75 mpfr_set_prec (x, 23);
76 mpfr_set_prec (y, 33);
77 mpfr_set_str_binary (x, "1.0001110110101001010100e-1");
78 inex = mpfr_rec_sqrt (y, x, MPFR_RNDU);
79 mpfr_set_prec (x, 33);
80 mpfr_set_str_binary (x, "1.01010110101110100100100101011");
81 MPFR_ASSERTN (inex > 0 && mpfr_cmp (x, y) == 0);
82
83 mpfr_clear (x);
84 mpfr_clear (y);
85 }
86
87 /* Worst case incorrectly rounded in r5573, found with the bad_cases test */
88 static void
89 bad_case1 (void)
90 {
91 mpfr_t x, y, z;
92
93 mpfr_init2 (x, 72);
94 mpfr_inits2 (6, y, z, (mpfr_ptr) 0);
95 mpfr_set_str (x, "1.08310518720928b30e@-120", 16, MPFR_RNDN);
96 mpfr_set_str (z, "f.8@59", 16, MPFR_RNDN);
97 /* z = rec_sqrt(x) rounded on 6 bits toward 0, the exact value
98 being ~= f.bffffffffffffffffa11@59. */
99 mpfr_rec_sqrt (y, x, MPFR_RNDZ);
100 if (mpfr_cmp0 (y, z) != 0)
101 {
102 printf ("Error in bad_case1\nexpected ");
103 mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN);
104 printf ("\ngot ");
105 mpfr_out_str (stdout, 16, 0, y, MPFR_RNDN);
106 printf ("\n");
107 exit (1);
108 }
109 mpfr_clears (x, y, z, (mpfr_ptr) 0);
110 }
111
112 static int
113 pm2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
114 {
115 return mpfr_pow_si (y, x, -2, rnd_mode);
116 }
117
118 /* exercises corner cases with inputs around 1 or 2 */
119 static void
120 bad_case2 (void)
121 {
122 mpfr_t r, u;
123 mpfr_prec_t pr, pu;
124 int rnd;
125
126 for (pr = MPFR_PREC_MIN; pr <= 192; pr++)
127 for (pu = MPFR_PREC_MIN; pu <= 192; pu++)
128 {
129 mpfr_init2 (r, pr);
130 mpfr_init2 (u, pu);
131
132 mpfr_set_ui (u, 1, MPFR_RNDN);
133 RND_LOOP (rnd)
134 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd);
135
136 mpfr_nextbelow (u);
137 RND_LOOP (rnd)
138 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd);
139
140 mpfr_nextbelow (u);
141 RND_LOOP (rnd)
142 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd);
143
144 mpfr_set_ui (u, 1, MPFR_RNDN);
145 mpfr_nextabove (u);
146 RND_LOOP (rnd)
147 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd);
148
149 mpfr_nextabove (u);
150 RND_LOOP (rnd)
151 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd);
152
153 mpfr_set_ui (u, 2, MPFR_RNDN);
154 RND_LOOP (rnd)
155 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd);
156
157 mpfr_nextbelow (u);
158 RND_LOOP (rnd)
159 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd);
160
161 mpfr_nextbelow (u);
162 RND_LOOP (rnd)
163 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd);
164
165 mpfr_set_ui (u, 2, MPFR_RNDN);
166 mpfr_nextabove (u);
167 RND_LOOP (rnd)
168 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd);
169
170 mpfr_nextabove (u);
171 RND_LOOP (rnd)
172 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd);
173
174 mpfr_clear (r);
175 mpfr_clear (u);
176 }
177 }
178
179 /* Before commits 270f4df6b3a49caae1cf564dcdc1c55b1c5989eb (master) and
180 934dd8842b4bdeb919a73123203bc8ce56db38d1 (4.2 branch) on 2023-04-17,
181 this was giving a stack overflow in mpfr_rec_sqrt due to a Ziv loop
182 where the working precision was increased additively instead of the
183 standard Ziv loop using the MPFR_ZIV_* macros.
184 */
185 static void
186 bad_case3 (void)
187 {
188 mpfr_t x, y;
189 int inex;
190
191 mpfr_init2 (x, 123456);
192 mpfr_init2 (y, 4);
193 mpfr_set_ui (y, 9, MPFR_RNDN);
194 mpfr_ui_div (x, 1, y, MPFR_RNDN);
195 inex = mpfr_rec_sqrt (y, x, MPFR_RNDN);
196 /* Let's also check the result, though this is not the real purpose
197 of this test (a stack overflow just makes the program crash).
198 1/9 = 0.111000111000111000111000111000111000...E-3 and since the
199 precision 123456 is divisible by 6, x > 1/9. Thus 1/sqrt(x) < 3. */
200 if (mpfr_cmp_ui0 (y, 3) != 0 || inex <= 0)
201 {
202 printf ("Error in bad_case3: expected 3 with inex > 0, got ");
203 mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN);
204 printf (" with inex=%d\n", inex);
205 exit (1);
206 }
207 mpfr_clear (x);
208 mpfr_clear (y);
209 }
210
211 /* timing test for n limbs (so that we can compare with GMP speed -s n) */
212 static void
213 test (unsigned long n)
214 {
215 mpfr_prec_t p = n * GMP_NUMB_BITS;
216 mpfr_t x, y, z;
217 gmp_randstate_t state;
218 double t;
219
220 gmp_randinit_default (state);
221 mpfr_init2 (x, p);
222 mpfr_init2 (y, p);
223 mpfr_init2 (z, p);
224 mpfr_urandom (x, state, MPFR_RNDN);
225 mpfr_urandom (y, state, MPFR_RNDN);
226
227 /* multiplication */
228 t = clock ();
229 mpfr_mul (z, x, y, MPFR_RNDN);
230 t = clock () - t;
231 printf ("mpfr_mul: %.3gs\n", t / (double) CLOCKS_PER_SEC);
232
233 /* squaring */
234 t = clock ();
235 mpfr_sqr (z, x, MPFR_RNDN);
236 t = clock () - t;
237 printf ("mpfr_sqr: %.3gs\n", t / (double) CLOCKS_PER_SEC);
238
239 /* square root */
240 t = clock ();
241 mpfr_sqrt (z, x, MPFR_RNDN);
242 t = clock () - t;
243 printf ("mpfr_sqrt: %.3gs\n", t / (double) CLOCKS_PER_SEC);
244
245 /* inverse square root */
246 t = clock ();
247 mpfr_rec_sqrt (z, x, MPFR_RNDN);
248 t = clock () - t;
249 printf ("mpfr_rec_sqrt: %.3gs\n", t / (double) CLOCKS_PER_SEC);
250
251 mpfr_clear (x);
252 mpfr_clear (y);
253 mpfr_clear (z);
254 gmp_randclear (state);
255 }
256
257 int
258 main (int argc, char *argv[])
259 {
260 tests_start_mpfr ();
261
262 if (argc == 2) /* trec_sqrt n */
263 {
264 unsigned long n = strtoul (argv[1], NULL, 10);
265 test (n);
266 goto end;
267 }
268
269 special ();
270 bad_case1 ();
271 bad_case2 ();
272 bad_case3 ();
273 test_generic (MPFR_PREC_MIN, 300, 15);
274
275 data_check ("data/rec_sqrt", mpfr_rec_sqrt, "mpfr_rec_sqrt");
276 bad_cases (mpfr_rec_sqrt, pm2, "mpfr_rec_sqrt", 0, -256, 255, 4, 128,
277 800, 50);
278 bad_cases (mpfr_rec_sqrt, pm2, "mpfr_rec_sqrt", 0, -256, 255, 9999, 9999,
279 120000, 1);
280
281 end:
282 tests_end_mpfr ();
283 return 0;
284 }