1/* Implementation of the degree trignometric functions COSD, SIND, TAND. 2 Copyright (C) 2020-2022 Free Software Foundation, Inc. 3 Contributed by Steven G. Kargl <kargl@gcc.gnu.org> 4 and Fritz Reese <foreese@gcc.gnu.org> 5 6This file is part of the GNU Fortran runtime library (libgfortran). 7 8Libgfortran is free software; you can redistribute it and/or 9modify it under the terms of the GNU General Public 10License as published by the Free Software Foundation; either 11version 3 of the License, or (at your option) any later version. 12 13Libgfortran is distributed in the hope that it will be useful, 14but WITHOUT ANY WARRANTY; without even the implied warranty of 15MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16GNU General Public License for more details. 17 18Under Section 7 of GPL version 3, you are granted additional 19permissions described in the GCC Runtime Library Exception, version 203.1, as published by the Free Software Foundation. 21 22You should have received a copy of the GNU General Public License and 23a copy of the GCC Runtime Library Exception along with this program; 24see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 25<http://www.gnu.org/licenses/>. */ 26 27 28/* 29 30This file is included from both the FE and the runtime library code. 31Operations are generalized using GMP/MPFR functions. When included from 32libgfortran, these should be overridden using macros which will use native 33operations conforming to the same API. From the FE, the GMP/MPFR functions can 34be used as-is. 35 36The following macros are used and must be defined, unless listed as [optional]: 37 38FTYPE 39 Type name for the real-valued parameter. 40 Variables of this type are constructed/destroyed using mpfr_init() 41 and mpfr_clear. 42 43RETTYPE 44 Return type of the functions. 45 46RETURN(x) 47 Insert code to return a value. 48 The parameter x is the result variable, which was also the input parameter. 49 50ITYPE 51 Type name for integer types. 52 53SIND, COSD, TRIGD 54 Names for the degree-valued trig functions defined by this module. 55 56ENABLE_SIND, ENABLE_COSD, ENABLE_TAND 57 Whether the degree-valued trig functions can be enabled. 58 59ERROR_RETURN(f, k, x) 60 If ENABLE_<xxx>D is not defined, this is substituted to assert an 61 error condition for function f, kind k, and parameter x. 62 The function argument is one of {sind, cosd, tand}. 63 64ISFINITE(x) 65 Whether x is a regular number or zero (not inf or NaN). 66 67D2R(x) 68 Convert x from radians to degrees. 69 70SET_COSD30(x) 71 Set x to COSD(30), or equivalently, SIND(60). 72 73TINY_LITERAL [optional] 74 Value subtracted from 1 to cause raise INEXACT for COSD(x) for x << 1. 75 If not set, COSD(x) for x <= COSD_SMALL_LITERAL simply returns 1. 76 77COSD_SMALL_LITERAL [optional] 78 Value such that x <= COSD_SMALL_LITERAL implies COSD(x) = 1 to within the 79 precision of FTYPE. If not set, this condition is not checked. 80 81SIND_SMALL_LITERAL [optional] 82 Value such that x <= SIND_SMALL_LITERAL implies SIND(x) = D2R(x) to within 83 the precision of FTYPE. If not set, this condition is not checked. 84 85*/ 86 87 88#ifdef SIND 89/* Compute sind(x) = sin(x * pi / 180). */ 90 91RETTYPE 92SIND (FTYPE x) 93{ 94#ifdef ENABLE_SIND 95 if (ISFINITE (x)) 96 { 97 FTYPE s, one; 98 99 /* sin(-x) = - sin(x). */ 100 mpfr_init (s); 101 mpfr_init_set_ui (one, 1, GFC_RND_MODE); 102 mpfr_copysign (s, one, x, GFC_RND_MODE); 103 mpfr_clear (one); 104 105#ifdef SIND_SMALL_LITERAL 106 /* sin(x) = x as x -> 0; but only for some precisions. */ 107 FTYPE ax; 108 mpfr_init (ax); 109 mpfr_abs (ax, x, GFC_RND_MODE); 110 if (mpfr_cmp_ld (ax, SIND_SMALL_LITERAL) < 0) 111 { 112 D2R (x); 113 mpfr_clear (ax); 114 return x; 115 } 116 117 mpfr_swap (x, ax); 118 mpfr_clear (ax); 119 120#else 121 mpfr_abs (x, x, GFC_RND_MODE); 122#endif /* SIND_SMALL_LITERAL */ 123 124 /* Reduce angle to x in [0,360]. */ 125 FTYPE period; 126 mpfr_init_set_ui (period, 360, GFC_RND_MODE); 127 mpfr_fmod (x, x, period, GFC_RND_MODE); 128 mpfr_clear (period); 129 130 /* Special cases with exact results. */ 131 ITYPE n; 132 mpz_init (n); 133 if (mpfr_get_z (n, x, GFC_RND_MODE) == 0 && mpz_divisible_ui_p (n, 30)) 134 { 135 /* Flip sign for odd n*pi (x is % 360 so this is only for 180). 136 This respects sgn(sin(x)) = sgn(d/dx sin(x)) = sgn(cos(x)). */ 137 if (mpz_divisible_ui_p (n, 180)) 138 { 139 mpfr_set_ui (x, 0, GFC_RND_MODE); 140 if (mpz_cmp_ui (n, 180) == 0) 141 mpfr_neg (s, s, GFC_RND_MODE); 142 } 143 else if (mpz_divisible_ui_p (n, 90)) 144 mpfr_set_si (x, (mpz_cmp_ui (n, 90) == 0 ? 1 : -1), GFC_RND_MODE); 145 else if (mpz_divisible_ui_p (n, 60)) 146 { 147 SET_COSD30 (x); 148 if (mpz_cmp_ui (n, 180) >= 0) 149 mpfr_neg (x, x, GFC_RND_MODE); 150 } 151 else 152 mpfr_set_ld (x, (mpz_cmp_ui (n, 180) < 0 ? 0.5L : -0.5L), 153 GFC_RND_MODE); 154 } 155 156 /* Fold [0,360] into the range [0,45], and compute either SIN() or 157 COS() depending on symmetry of shifting into the [0,45] range. */ 158 else 159 { 160 bool fold_cos = false; 161 if (mpfr_cmp_ui (x, 180) <= 0) 162 { 163 if (mpfr_cmp_ui (x, 90) <= 0) 164 { 165 if (mpfr_cmp_ui (x, 45) > 0) 166 { 167 /* x = COS(D2R(90 - x)) */ 168 mpfr_ui_sub (x, 90, x, GFC_RND_MODE); 169 fold_cos = true; 170 } 171 } 172 else 173 { 174 if (mpfr_cmp_ui (x, 135) <= 0) 175 { 176 mpfr_sub_ui (x, x, 90, GFC_RND_MODE); 177 fold_cos = true; 178 } 179 else 180 mpfr_ui_sub (x, 180, x, GFC_RND_MODE); 181 } 182 } 183 184 else if (mpfr_cmp_ui (x, 270) <= 0) 185 { 186 if (mpfr_cmp_ui (x, 225) <= 0) 187 mpfr_sub_ui (x, x, 180, GFC_RND_MODE); 188 else 189 { 190 mpfr_ui_sub (x, 270, x, GFC_RND_MODE); 191 fold_cos = true; 192 } 193 mpfr_neg (s, s, GFC_RND_MODE); 194 } 195 196 else 197 { 198 if (mpfr_cmp_ui (x, 315) <= 0) 199 { 200 mpfr_sub_ui (x, x, 270, GFC_RND_MODE); 201 fold_cos = true; 202 } 203 else 204 mpfr_ui_sub (x, 360, x, GFC_RND_MODE); 205 mpfr_neg (s, s, GFC_RND_MODE); 206 } 207 208 D2R (x); 209 210 if (fold_cos) 211 mpfr_cos (x, x, GFC_RND_MODE); 212 else 213 mpfr_sin (x, x, GFC_RND_MODE); 214 } 215 216 mpfr_mul (x, x, s, GFC_RND_MODE); 217 mpz_clear (n); 218 mpfr_clear (s); 219 } 220 221 /* Return NaN for +-Inf and NaN and raise exception. */ 222 else 223 mpfr_sub (x, x, x, GFC_RND_MODE); 224 225 RETURN (x); 226 227#else 228 ERROR_RETURN(sind, KIND, x); 229#endif // ENABLE_SIND 230} 231#endif // SIND 232 233 234#ifdef COSD 235/* Compute cosd(x) = cos(x * pi / 180). */ 236 237RETTYPE 238COSD (FTYPE x) 239{ 240#ifdef ENABLE_COSD 241#if defined(TINY_LITERAL) && defined(COSD_SMALL_LITERAL) 242 static const volatile FTYPE tiny = TINY_LITERAL; 243#endif 244 245 if (ISFINITE (x)) 246 { 247#ifdef COSD_SMALL_LITERAL 248 FTYPE ax; 249 mpfr_init (ax); 250 251 mpfr_abs (ax, x, GFC_RND_MODE); 252 /* No spurious underflows!. In radians, cos(x) = 1-x*x/2 as x -> 0. */ 253 if (mpfr_cmp_ld (ax, COSD_SMALL_LITERAL) <= 0) 254 { 255 mpfr_set_ui (x, 1, GFC_RND_MODE); 256#ifdef TINY_LITERAL 257 /* Cause INEXACT. */ 258 if (!mpfr_zero_p (ax)) 259 mpfr_sub_d (x, x, tiny, GFC_RND_MODE); 260#endif 261 262 mpfr_clear (ax); 263 return x; 264 } 265 266 mpfr_swap (x, ax); 267 mpfr_clear (ax); 268#else 269 mpfr_abs (x, x, GFC_RND_MODE); 270#endif /* COSD_SMALL_LITERAL */ 271 272 /* Reduce angle to ax in [0,360]. */ 273 FTYPE period; 274 mpfr_init_set_ui (period, 360, GFC_RND_MODE); 275 mpfr_fmod (x, x, period, GFC_RND_MODE); 276 mpfr_clear (period); 277 278 /* Special cases with exact results. 279 Return negative zero for cosd(270) for consistency with libm cos(). */ 280 ITYPE n; 281 mpz_init (n); 282 if (mpfr_get_z (n, x, GFC_RND_MODE) == 0 && mpz_divisible_ui_p (n, 30)) 283 { 284 if (mpz_divisible_ui_p (n, 180)) 285 mpfr_set_si (x, (mpz_cmp_ui (n, 180) == 0 ? -1 : 1), 286 GFC_RND_MODE); 287 else if (mpz_divisible_ui_p (n, 90)) 288 mpfr_set_zero (x, 0); 289 else if (mpz_divisible_ui_p (n, 60)) 290 { 291 mpfr_set_ld (x, 0.5, GFC_RND_MODE); 292 if (mpz_cmp_ui (n, 60) != 0 && mpz_cmp_ui (n, 300) != 0) 293 mpfr_neg (x, x, GFC_RND_MODE); 294 } 295 else 296 { 297 SET_COSD30 (x); 298 if (mpz_cmp_ui (n, 30) != 0 && mpz_cmp_ui (n, 330) != 0) 299 mpfr_neg (x, x, GFC_RND_MODE); 300 } 301 } 302 303 /* Fold [0,360] into the range [0,45], and compute either SIN() or 304 COS() depending on symmetry of shifting into the [0,45] range. */ 305 else 306 { 307 bool neg = false; 308 bool fold_sin = false; 309 if (mpfr_cmp_ui (x, 180) <= 0) 310 { 311 if (mpfr_cmp_ui (x, 90) <= 0) 312 { 313 if (mpfr_cmp_ui (x, 45) > 0) 314 { 315 mpfr_ui_sub (x, 90, x, GFC_RND_MODE); 316 fold_sin = true; 317 } 318 } 319 else 320 { 321 if (mpfr_cmp_ui (x, 135) <= 0) 322 { 323 mpfr_sub_ui (x, x, 90, GFC_RND_MODE); 324 fold_sin = true; 325 } 326 else 327 mpfr_ui_sub (x, 180, x, GFC_RND_MODE); 328 neg = true; 329 } 330 } 331 332 else if (mpfr_cmp_ui (x, 270) <= 0) 333 { 334 if (mpfr_cmp_ui (x, 225) <= 0) 335 mpfr_sub_ui (x, x, 180, GFC_RND_MODE); 336 else 337 { 338 mpfr_ui_sub (x, 270, x, GFC_RND_MODE); 339 fold_sin = true; 340 } 341 neg = true; 342 } 343 344 else 345 { 346 if (mpfr_cmp_ui (x, 315) <= 0) 347 { 348 mpfr_sub_ui (x, x, 270, GFC_RND_MODE); 349 fold_sin = true; 350 } 351 else 352 mpfr_ui_sub (x, 360, x, GFC_RND_MODE); 353 } 354 355 D2R (x); 356 357 if (fold_sin) 358 mpfr_sin (x, x, GFC_RND_MODE); 359 else 360 mpfr_cos (x, x, GFC_RND_MODE); 361 362 if (neg) 363 mpfr_neg (x, x, GFC_RND_MODE); 364 } 365 366 mpz_clear (n); 367 } 368 369 /* Return NaN for +-Inf and NaN and raise exception. */ 370 else 371 mpfr_sub (x, x, x, GFC_RND_MODE); 372 373 RETURN (x); 374 375#else 376 ERROR_RETURN(cosd, KIND, x); 377#endif // ENABLE_COSD 378} 379#endif // COSD 380 381 382#ifdef TAND 383/* Compute tand(x) = tan(x * pi / 180). */ 384 385RETTYPE 386TAND (FTYPE x) 387{ 388#ifdef ENABLE_TAND 389 if (ISFINITE (x)) 390 { 391 FTYPE s, one; 392 393 /* tan(-x) = - tan(x). */ 394 mpfr_init (s); 395 mpfr_init_set_ui (one, 1, GFC_RND_MODE); 396 mpfr_copysign (s, one, x, GFC_RND_MODE); 397 mpfr_clear (one); 398 399#ifdef SIND_SMALL_LITERAL 400 /* tan(x) = x as x -> 0; but only for some precisions. */ 401 FTYPE ax; 402 mpfr_init (ax); 403 mpfr_abs (ax, x, GFC_RND_MODE); 404 if (mpfr_cmp_ld (ax, SIND_SMALL_LITERAL) < 0) 405 { 406 D2R (x); 407 mpfr_clear (ax); 408 return x; 409 } 410 411 mpfr_swap (x, ax); 412 mpfr_clear (ax); 413 414#else 415 mpfr_abs (x, x, GFC_RND_MODE); 416#endif /* SIND_SMALL_LITERAL */ 417 418 /* Reduce angle to x in [0,360]. */ 419 FTYPE period; 420 mpfr_init_set_ui (period, 360, GFC_RND_MODE); 421 mpfr_fmod (x, x, period, GFC_RND_MODE); 422 mpfr_clear (period); 423 424 /* Special cases with exact results. */ 425 ITYPE n; 426 mpz_init (n); 427 if (mpfr_get_z (n, x, GFC_RND_MODE) == 0 && mpz_divisible_ui_p (n, 45)) 428 { 429 if (mpz_divisible_ui_p (n, 180)) 430 mpfr_set_zero (x, 0); 431 432 /* Though mathematically NaN is more appropriate for tan(n*90), 433 returning +/-Inf offers the advantage that 1/tan(n*90) returns 0, 434 which is mathematically sound. In fact we rely on this behavior 435 to implement COTAND(x) = 1 / TAND(x). 436 */ 437 else if (mpz_divisible_ui_p (n, 90)) 438 mpfr_set_inf (x, mpz_cmp_ui (n, 90) == 0 ? 0 : 1); 439 440 else 441 { 442 mpfr_set_ui (x, 1, GFC_RND_MODE); 443 if (mpz_cmp_ui (n, 45) != 0 && mpz_cmp_ui (n, 225) != 0) 444 mpfr_neg (x, x, GFC_RND_MODE); 445 } 446 } 447 448 else 449 { 450 /* Fold [0,360] into the range [0,90], and compute TAN(). */ 451 if (mpfr_cmp_ui (x, 180) <= 0) 452 { 453 if (mpfr_cmp_ui (x, 90) > 0) 454 { 455 mpfr_ui_sub (x, 180, x, GFC_RND_MODE); 456 mpfr_neg (s, s, GFC_RND_MODE); 457 } 458 } 459 else 460 { 461 if (mpfr_cmp_ui (x, 270) <= 0) 462 { 463 mpfr_sub_ui (x, x, 180, GFC_RND_MODE); 464 } 465 else 466 { 467 mpfr_ui_sub (x, 360, x, GFC_RND_MODE); 468 mpfr_neg (s, s, GFC_RND_MODE); 469 } 470 } 471 472 D2R (x); 473 mpfr_tan (x, x, GFC_RND_MODE); 474 } 475 476 mpfr_mul (x, x, s, GFC_RND_MODE); 477 mpz_clear (n); 478 mpfr_clear (s); 479 } 480 481 /* Return NaN for +-Inf and NaN and raise exception. */ 482 else 483 mpfr_sub (x, x, x, GFC_RND_MODE); 484 485 RETURN (x); 486 487#else 488 ERROR_RETURN(tand, KIND, x); 489#endif // ENABLE_TAND 490} 491#endif // TAND 492 493/* vim: set ft=c: */ 494