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