1 /* Test file for mpfr_rec_sqrt. 2 3 Copyright 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramel 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 http://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 <stdio.h> 24 #include <stdlib.h> 25 26 #include "mpfr-test.h" 27 28 #if MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0) 29 30 #define TEST_FUNCTION mpfr_rec_sqrt 31 #define TEST_RANDOM_POS 8 /* 8/512 = 1/64 of the tested numbers are negative */ 32 #include "tgeneric.c" 33 34 static void 35 special (void) 36 { 37 mpfr_t x, y; 38 int inex; 39 40 mpfr_init (x); 41 mpfr_init (y); 42 43 /* rec_sqrt(NaN) = NaN */ 44 mpfr_set_nan (x); 45 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN); 46 MPFR_ASSERTN(mpfr_nan_p (x) && inex == 0); 47 48 /* rec_sqrt(+Inf) = +0 */ 49 mpfr_set_inf (x, 1); 50 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN); 51 MPFR_ASSERTN(mpfr_zero_p (x) && MPFR_IS_POS(x) && inex == 0); 52 53 /* rec_sqrt(-Inf) = NaN */ 54 mpfr_set_inf (x, -1); 55 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN); 56 MPFR_ASSERTN(mpfr_nan_p (x) && inex == 0); 57 58 /* rec_sqrt(+0) = +Inf */ 59 mpfr_set_ui (x, 0, MPFR_RNDN); 60 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN); 61 MPFR_ASSERTN(mpfr_inf_p (x) && MPFR_IS_POS(x) && inex == 0); 62 63 /* rec_sqrt(-0) = +Inf */ 64 mpfr_set_ui (x, 0, MPFR_RNDN); 65 mpfr_neg (x, x, MPFR_RNDN); 66 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN); 67 MPFR_ASSERTN(mpfr_inf_p (x) && MPFR_IS_POS(x) && inex == 0); 68 69 /* rec_sqrt(-1) = NaN */ 70 mpfr_set_si (x, -1, MPFR_RNDN); 71 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN); 72 MPFR_ASSERTN(mpfr_nan_p (x) && inex == 0); 73 74 /* rec_sqrt(1) = 1 */ 75 mpfr_set_ui (x, 1, MPFR_RNDN); 76 inex = mpfr_rec_sqrt (x, x, MPFR_RNDN); 77 MPFR_ASSERTN((mpfr_cmp_ui (x, 1) == 0) && (inex == 0)); 78 79 mpfr_set_prec (x, 23); 80 mpfr_set_prec (y, 33); 81 mpfr_set_str_binary (x, "1.0001110110101001010100e-1"); 82 inex = mpfr_rec_sqrt (y, x, MPFR_RNDU); 83 mpfr_set_prec (x, 33); 84 mpfr_set_str_binary (x, "1.01010110101110100100100101011"); 85 MPFR_ASSERTN (inex > 0 && mpfr_cmp (x, y) == 0); 86 87 mpfr_clear (x); 88 mpfr_clear (y); 89 } 90 91 /* Worst case incorrectly rounded in r5573, found with the bad_cases test */ 92 static void 93 bad_case1 (void) 94 { 95 mpfr_t x, y, z; 96 97 mpfr_init2 (x, 72); 98 mpfr_inits2 (6, y, z, (mpfr_ptr) 0); 99 mpfr_set_str (x, "1.08310518720928b30e@-120", 16, MPFR_RNDN); 100 mpfr_set_str (z, "f.8@59", 16, MPFR_RNDN); 101 /* z = rec_sqrt(x) rounded on 6 bits toward 0, the exact value 102 being ~= f.bffffffffffffffffa11@59. */ 103 mpfr_rec_sqrt (y, x, MPFR_RNDZ); 104 if (mpfr_cmp0 (y, z) != 0) 105 { 106 printf ("Error in bad_case1\nexpected "); 107 mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN); 108 printf ("\ngot "); 109 mpfr_out_str (stdout, 16, 0, y, MPFR_RNDN); 110 printf ("\n"); 111 exit (1); 112 } 113 mpfr_clears (x, y, z, (mpfr_ptr) 0); 114 } 115 116 static int 117 pm2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 118 { 119 return mpfr_pow_si (y, x, -2, rnd_mode); 120 } 121 122 /* exercises corner cases with inputs around 1 or 2 */ 123 static void 124 bad_case2 (void) 125 { 126 mpfr_t r, u; 127 mpfr_prec_t pr, pu; 128 int rnd; 129 130 for (pr = MPFR_PREC_MIN; pr <= 192; pr++) 131 for (pu = MPFR_PREC_MIN; pu <= 192; pu++) 132 { 133 mpfr_init2 (r, pr); 134 mpfr_init2 (u, pu); 135 136 mpfr_set_ui (u, 1, MPFR_RNDN); 137 for (rnd = 0; rnd < MPFR_RND_MAX; rnd++) 138 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd); 139 140 mpfr_nextbelow (u); 141 for (rnd = 0; rnd < MPFR_RND_MAX; rnd++) 142 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd); 143 144 mpfr_nextbelow (u); 145 for (rnd = 0; rnd < MPFR_RND_MAX; rnd++) 146 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd); 147 148 mpfr_set_ui (u, 1, MPFR_RNDN); 149 mpfr_nextabove (u); 150 for (rnd = 0; rnd < MPFR_RND_MAX; rnd++) 151 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd); 152 153 mpfr_nextabove (u); 154 for (rnd = 0; rnd < MPFR_RND_MAX; rnd++) 155 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd); 156 157 mpfr_set_ui (u, 2, MPFR_RNDN); 158 for (rnd = 0; rnd < MPFR_RND_MAX; rnd++) 159 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd); 160 161 mpfr_nextbelow (u); 162 for (rnd = 0; rnd < MPFR_RND_MAX; rnd++) 163 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd); 164 165 mpfr_nextbelow (u); 166 for (rnd = 0; rnd < MPFR_RND_MAX; rnd++) 167 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd); 168 169 mpfr_set_ui (u, 2, MPFR_RNDN); 170 mpfr_nextabove (u); 171 for (rnd = 0; rnd < MPFR_RND_MAX; rnd++) 172 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd); 173 174 mpfr_nextabove (u); 175 for (rnd = 0; rnd < MPFR_RND_MAX; rnd++) 176 mpfr_rec_sqrt (r, u, (mpfr_rnd_t) rnd); 177 178 mpfr_clear (r); 179 mpfr_clear (u); 180 } 181 } 182 183 int 184 main (void) 185 { 186 tests_start_mpfr (); 187 188 special (); 189 bad_case1 (); 190 bad_case2 (); 191 test_generic (2, 300, 15); 192 193 data_check ("data/rec_sqrt", mpfr_rec_sqrt, "mpfr_rec_sqrt"); 194 bad_cases (mpfr_rec_sqrt, pm2, "mpfr_rec_sqrt", 8, -256, 255, 4, 128, 195 800, 50); 196 197 tests_end_mpfr (); 198 return 0; 199 } 200 201 #else /* MPFR_VERSION */ 202 203 int 204 main (void) 205 { 206 printf ("Warning! Test disabled for this MPFR version.\n"); 207 return 0; 208 } 209 210 #endif /* MPFR_VERSION */ 211