1/* Test file for mpfr_sin. 2 3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao 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 20http://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 <stdlib.h> 24 25#include "mpfr-test.h" 26 27#ifdef CHECK_EXTERNAL 28static int 29test_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 50static void 51check53 (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 72static void 73check53b (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 94static void 95test_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 162static void 163check_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 202const char xs[] = "0.111011111110110000111000001100000111110E-1"; 203 204static void 205check_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. */ 228static void 229check_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 273int 274main (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_clear (s2); 352 mpfr_clear (c2); 353 mpfr_clear (s); 354 mpfr_clear (c); 355 mpfr_clear (x); 356 357 test_generic (2, 100, 15); 358 test_sign (); 359 check_tiny (); 360 361 data_check ("data/sin", mpfr_sin, "mpfr_sin"); 362 bad_cases (mpfr_sin, mpfr_asin, "mpfr_sin", 256, -40, 0, 4, 128, 800, 50); 363 364 tests_end_mpfr (); 365 return 0; 366} 367