1/* Test file for mpfr_rec_sqrt.
2
3Copyright 2008-2023 Free Software Foundation, Inc.
4Contributed by the AriC and Caramba projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 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
30static void
31special (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 */
88static void
89bad_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
112static int
113pm2 (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 */
119static void
120bad_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*/
185static void
186bad_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) */
212static void
213test (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
257int
258main (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}
285