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