1 /* Test file for mpfr_sin. 2 3 Copyright 2001-2018 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 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 "mpfr-test.h" 24 25 #ifdef CHECK_EXTERNAL 26 static int 27 test_sin (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode) 28 { 29 int res; 30 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_get_prec (a)>=53; 31 if (ok) 32 { 33 mpfr_print_raw (b); 34 } 35 res = mpfr_sin (a, b, rnd_mode); 36 if (ok) 37 { 38 printf (" "); 39 mpfr_print_raw (a); 40 printf ("\n"); 41 } 42 return res; 43 } 44 #else 45 #define test_sin mpfr_sin 46 #endif 47 48 static void 49 check53 (const char *xs, const char *sin_xs, mpfr_rnd_t rnd_mode) 50 { 51 mpfr_t xx, s; 52 53 mpfr_init2 (xx, 53); 54 mpfr_init2 (s, 53); 55 mpfr_set_str1 (xx, xs); /* should be exact */ 56 test_sin (s, xx, rnd_mode); 57 if (mpfr_cmp_str1 (s, sin_xs)) 58 { 59 printf ("mpfr_sin failed for x=%s, rnd=%s\n", 60 xs, mpfr_print_rnd_mode (rnd_mode)); 61 printf ("mpfr_sin gives sin(x)="); 62 mpfr_out_str (stdout, 10, 0, s, MPFR_RNDN); 63 printf (", expected %s\n", sin_xs); 64 exit (1); 65 } 66 mpfr_clear (xx); 67 mpfr_clear (s); 68 } 69 70 static void 71 check53b (const char *xs, const char *sin_xs, mpfr_rnd_t rnd_mode) 72 { 73 mpfr_t xx, s; 74 75 mpfr_init2 (xx, 53); 76 mpfr_init2 (s, 53); 77 mpfr_set_str (xx, xs, 2, MPFR_RNDN); /* should be exact */ 78 test_sin (s, xx, rnd_mode); 79 if (mpfr_cmp_str (s, sin_xs, 2, MPFR_RNDN)) 80 { 81 printf ("mpfr_sin failed in rounding mode %s for\n x = %s\n", 82 mpfr_print_rnd_mode (rnd_mode), xs); 83 printf (" got "); 84 mpfr_out_str (stdout, 2, 0, s, MPFR_RNDN); 85 printf ("\nexpected %s\n", sin_xs); 86 exit (1); 87 } 88 mpfr_clear (xx); 89 mpfr_clear (s); 90 } 91 92 static void 93 test_sign (void) 94 { 95 mpfr_t pid, piu, x, y; 96 int p, k; 97 98 mpfr_init2 (pid, 4096); 99 mpfr_const_pi (pid, MPFR_RNDD); 100 mpfr_init2 (piu, 4096); 101 mpfr_const_pi (piu, MPFR_RNDU); 102 mpfr_init (x); 103 mpfr_init2 (y, 2); 104 for (p = 8; p <= 128; p++) 105 for (k = 2; k <= 6; k += 2) 106 { 107 mpfr_set_prec (x, p); 108 mpfr_mul_ui (x, pid, k, MPFR_RNDD); 109 test_sin (y, x, MPFR_RNDN); 110 if (MPFR_IS_POS (y)) 111 { 112 printf ("Error in test_sign for sin(%dpi-epsilon), prec = %d" 113 " for argument.\nResult should have been negative.\n", 114 k, p); 115 exit (1); 116 } 117 mpfr_mul_ui (x, piu, k, MPFR_RNDU); 118 test_sin (y, x, MPFR_RNDN); 119 if (MPFR_IS_NEG (y)) 120 { 121 printf ("Error in test_sign for sin(%dpi+epsilon), prec = %d" 122 " for argument.\nResult should have been positive.\n", 123 k, p); 124 exit (1); 125 } 126 } 127 128 /* worst case on 53 bits */ 129 mpfr_set_prec (x, 53); 130 mpfr_set_prec (y, 53); 131 mpfr_set_str (x, "6134899525417045", 10, MPFR_RNDN); 132 test_sin (y, x, MPFR_RNDN); 133 mpfr_set_str_binary (x, "11011010111101011110111100010101010101110000000001011E-106"); 134 MPFR_ASSERTN(mpfr_cmp (x, y) == 0); 135 136 /* Bug on Special cases */ 137 mpfr_set_str_binary (x, "0.100011011010111101E-32"); 138 test_sin (y, x, MPFR_RNDN); 139 if (mpfr_cmp_str (y, "0.10001101101011110100000000000000000000000000000000000E-32", 2, MPFR_RNDN)) 140 { 141 printf("sin special 97 error:\nx="); 142 mpfr_dump (x); printf("y="); 143 mpfr_dump (y); 144 exit (1); 145 } 146 147 mpfr_set_prec (x, 53); 148 mpfr_set_prec (y, 53); 149 mpfr_set_str_binary (x, "1.1001001000011111101101010100010001000010110100010011"); 150 mpfr_set_str_binary (y, "1.1111111111111111111111111111111111111111111111111111e-1"); 151 test_sin (x, x, MPFR_RNDZ); 152 MPFR_ASSERTN(mpfr_cmp (x, y) == 0); 153 154 mpfr_clear (pid); 155 mpfr_clear (piu); 156 mpfr_clear (x); 157 mpfr_clear (y); 158 } 159 160 static void 161 check_nans (void) 162 { 163 mpfr_t x, y; 164 165 mpfr_init2 (x, 123L); 166 mpfr_init2 (y, 123L); 167 168 mpfr_set_nan (x); 169 test_sin (y, x, MPFR_RNDN); 170 if (! mpfr_nan_p (y)) 171 { 172 printf ("Error: sin(NaN) != NaN\n"); 173 exit (1); 174 } 175 176 mpfr_set_inf (x, 1); 177 test_sin (y, x, MPFR_RNDN); 178 if (! mpfr_nan_p (y)) 179 { 180 printf ("Error: sin(Inf) != NaN\n"); 181 exit (1); 182 } 183 184 mpfr_set_inf (x, -1); 185 test_sin (y, x, MPFR_RNDN); 186 if (! mpfr_nan_p (y)) 187 { 188 printf ("Error: sin(-Inf) != NaN\n"); 189 exit (1); 190 } 191 192 mpfr_clear (x); 193 mpfr_clear (y); 194 } 195 196 #define TEST_FUNCTION test_sin 197 #define REDUCE_EMAX 262143 /* otherwise arg. reduction is too expensive */ 198 #include "tgeneric.c" 199 200 const char xs[] = "0.111011111110110000111000001100000111110E-1"; 201 202 static void 203 check_regression (void) 204 { 205 mpfr_t x, y; 206 mpfr_prec_t p; 207 int i; 208 209 p = strlen (xs) - 2 - 3; 210 mpfr_inits2 (p, x, y, (mpfr_ptr) 0); 211 212 mpfr_set_str (x, xs, 2, MPFR_RNDN); 213 i = mpfr_sin (y, x, MPFR_RNDN); 214 if (i >= 0 215 || mpfr_cmp_str (y, "0.111001110011110011110001010110011101110E-1", 216 2, MPFR_RNDN)) 217 { 218 printf ("Regression test failed (1) i=%d\ny=", i); 219 mpfr_dump (y); 220 exit (1); 221 } 222 mpfr_clears (x, y, (mpfr_ptr) 0); 223 } 224 225 /* Test provided by Christopher Creutzig, 2007-05-21. */ 226 static void 227 check_tiny (void) 228 { 229 mpfr_t x, y; 230 231 mpfr_init2 (x, 53); 232 mpfr_init2 (y, 53); 233 234 mpfr_set_ui (x, 1, MPFR_RNDN); 235 mpfr_set_exp (x, mpfr_get_emin ()); 236 mpfr_sin (y, x, MPFR_RNDD); 237 if (mpfr_cmp (x, y) < 0) 238 { 239 printf ("Error in check_tiny: got sin(x) > x for x = 2^(emin-1)\n"); 240 exit (1); 241 } 242 243 mpfr_sin (y, x, MPFR_RNDU); 244 mpfr_mul_2ui (y, y, 1, MPFR_RNDU); 245 if (mpfr_cmp (x, y) > 0) 246 { 247 printf ("Error in check_tiny: got sin(x) < x/2 for x = 2^(emin-1)\n"); 248 exit (1); 249 } 250 251 mpfr_neg (x, x, MPFR_RNDN); 252 mpfr_sin (y, x, MPFR_RNDU); 253 if (mpfr_cmp (x, y) > 0) 254 { 255 printf ("Error in check_tiny: got sin(x) < x for x = -2^(emin-1)\n"); 256 exit (1); 257 } 258 259 mpfr_sin (y, x, MPFR_RNDD); 260 mpfr_mul_2ui (y, y, 1, MPFR_RNDD); 261 if (mpfr_cmp (x, y) < 0) 262 { 263 printf ("Error in check_tiny: got sin(x) > x/2 for x = -2^(emin-1)\n"); 264 exit (1); 265 } 266 267 mpfr_clear (y); 268 mpfr_clear (x); 269 } 270 271 int 272 main (int argc, char *argv[]) 273 { 274 mpfr_t x, c, s, c2, s2; 275 276 tests_start_mpfr (); 277 278 check_regression (); 279 check_nans (); 280 281 /* worst case from PhD thesis of Vincent Lefe`vre: x=8980155785351021/2^54 */ 282 check53 ("4.984987858808754279e-1", "4.781075595393330379e-1", MPFR_RNDN); 283 check53 ("4.984987858808754279e-1", "4.781075595393329824e-1", MPFR_RNDD); 284 check53 ("4.984987858808754279e-1", "4.781075595393329824e-1", MPFR_RNDZ); 285 check53 ("4.984987858808754279e-1", "4.781075595393330379e-1", MPFR_RNDU); 286 check53 ("1.00031274099908640274", "8.416399183372403892e-1", MPFR_RNDN); 287 check53 ("1.00229256850978698523", "8.427074524447979442e-1", MPFR_RNDZ); 288 check53 ("1.00288304857059840103", "8.430252033025980029e-1", MPFR_RNDZ); 289 check53 ("1.00591265847407274059", "8.446508805292128885e-1", MPFR_RNDN); 290 291 /* Other worst cases showing a bug introduced on 2005-01-29 in rev 3248 */ 292 check53b ("1.0111001111010111010111111000010011010001110001111011e-21", 293 "1.0111001111010111010111111000010011010001101001110001e-21", 294 MPFR_RNDU); 295 check53b ("1.1011101111111010000001010111000010000111100100101101", 296 "1.1111100100101100001111100000110011110011010001010101e-1", 297 MPFR_RNDU); 298 299 mpfr_init2 (x, 2); 300 301 mpfr_set_str (x, "0.5", 10, MPFR_RNDN); 302 test_sin (x, x, MPFR_RNDD); 303 if (mpfr_cmp_ui_2exp (x, 3, -3)) /* x != 0.375 = 3/8 */ 304 { 305 printf ("mpfr_sin(0.5, MPFR_RNDD) failed with precision=2\n"); 306 exit (1); 307 } 308 309 /* bug found by Kevin Ryde */ 310 mpfr_const_pi (x, MPFR_RNDN); 311 mpfr_mul_ui (x, x, 3L, MPFR_RNDN); 312 mpfr_div_ui (x, x, 2L, MPFR_RNDN); 313 test_sin (x, x, MPFR_RNDN); 314 if (mpfr_cmp_ui (x, 0) >= 0) 315 { 316 printf ("Error: wrong sign for sin(3*Pi/2)\n"); 317 exit (1); 318 } 319 320 /* Can fail on an assert */ 321 mpfr_set_prec (x, 53); 322 mpfr_set_str (x, "77291789194529019661184401408", 10, MPFR_RNDN); 323 mpfr_init2 (c, 4); mpfr_init2 (s, 42); 324 mpfr_init2 (c2, 4); mpfr_init2 (s2, 42); 325 326 test_sin (s, x, MPFR_RNDN); 327 mpfr_cos (c, x, MPFR_RNDN); 328 mpfr_sin_cos (s2, c2, x, MPFR_RNDN); 329 if (mpfr_cmp (c2, c)) 330 { 331 printf("cos differs for x=77291789194529019661184401408"); 332 exit (1); 333 } 334 if (mpfr_cmp (s2, s)) 335 { 336 printf("sin differs for x=77291789194529019661184401408"); 337 exit (1); 338 } 339 340 mpfr_set_str_binary (x, "1.1001001000011111101101010100010001000010110100010011"); 341 test_sin (x, x, MPFR_RNDZ); 342 if (mpfr_cmp_str (x, "1.1111111111111111111111111111111111111111111111111111e-1", 2, MPFR_RNDN)) 343 { 344 printf ("Error for x= 1.1001001000011111101101010100010001000010110100010011\nGot "); 345 mpfr_dump (x); 346 exit (1); 347 } 348 349 mpfr_set_prec (s, 9); 350 mpfr_set_prec (x, 190); 351 mpfr_const_pi (x, MPFR_RNDN); 352 mpfr_sin (s, x, MPFR_RNDZ); 353 if (mpfr_cmp_str (s, "0.100000101e-196", 2, MPFR_RNDN)) 354 { 355 printf ("Error for x ~= pi\n"); 356 mpfr_dump (s); 357 exit (1); 358 } 359 360 mpfr_clear (s2); 361 mpfr_clear (c2); 362 mpfr_clear (s); 363 mpfr_clear (c); 364 mpfr_clear (x); 365 366 test_generic (MPFR_PREC_MIN, 100, 15); 367 test_generic (MPFR_SINCOS_THRESHOLD-1, MPFR_SINCOS_THRESHOLD+1, 2); 368 test_sign (); 369 check_tiny (); 370 371 data_check ("data/sin", mpfr_sin, "mpfr_sin"); 372 bad_cases (mpfr_sin, mpfr_asin, "mpfr_sin", 256, -40, 0, 4, 128, 800, 50); 373 374 tests_end_mpfr (); 375 return 0; 376 } 377