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 /* timing test for n limbs (so that we can compare with GMP speed -s n) */ 180 static void 181 test (unsigned long n) 182 { 183 mpfr_prec_t p = n * GMP_NUMB_BITS; 184 mpfr_t x, y, z; 185 gmp_randstate_t state; 186 double t; 187 188 gmp_randinit_default (state); 189 mpfr_init2 (x, p); 190 mpfr_init2 (y, p); 191 mpfr_init2 (z, p); 192 mpfr_urandom (x, state, MPFR_RNDN); 193 mpfr_urandom (y, state, MPFR_RNDN); 194 195 /* multiplication */ 196 t = clock (); 197 mpfr_mul (z, x, y, MPFR_RNDN); 198 t = clock () - t; 199 printf ("mpfr_mul: %.3gs\n", t / (double) CLOCKS_PER_SEC); 200 201 /* squaring */ 202 t = clock (); 203 mpfr_sqr (z, x, MPFR_RNDN); 204 t = clock () - t; 205 printf ("mpfr_sqr: %.3gs\n", t / (double) CLOCKS_PER_SEC); 206 207 /* square root */ 208 t = clock (); 209 mpfr_sqrt (z, x, MPFR_RNDN); 210 t = clock () - t; 211 printf ("mpfr_sqrt: %.3gs\n", t / (double) CLOCKS_PER_SEC); 212 213 /* inverse square root */ 214 t = clock (); 215 mpfr_rec_sqrt (z, x, MPFR_RNDN); 216 t = clock () - t; 217 printf ("mpfr_rec_sqrt: %.3gs\n", t / (double) CLOCKS_PER_SEC); 218 219 mpfr_clear (x); 220 mpfr_clear (y); 221 mpfr_clear (z); 222 gmp_randclear (state); 223 } 224 225 int 226 main (int argc, char *argv[]) 227 { 228 tests_start_mpfr (); 229 230 if (argc == 2) /* trec_sqrt n */ 231 { 232 unsigned long n = strtoul (argv[1], NULL, 10); 233 test (n); 234 goto end; 235 } 236 237 special (); 238 bad_case1 (); 239 bad_case2 (); 240 test_generic (MPFR_PREC_MIN, 300, 15); 241 242 data_check ("data/rec_sqrt", mpfr_rec_sqrt, "mpfr_rec_sqrt"); 243 bad_cases (mpfr_rec_sqrt, pm2, "mpfr_rec_sqrt", 0, -256, 255, 4, 128, 244 800, 50); 245 246 end: 247 tests_end_mpfr (); 248 return 0; 249 } 250