libgfortran/intrinsics/trigd.inc

   1 /* Implementation of the degree trignometric functions COSD, SIND, TAND.
   2    Copyright (C) 2020-2023 Free Software Foundation, Inc.
   3    Contributed by Steven G. Kargl <kargl@gcc.gnu.org>
   4    and Fritz Reese <foreese@gcc.gnu.org>
   5
   6 This file is part of the GNU Fortran runtime library (libgfortran).
   7
   8 Libgfortran is free software; you can redistribute it and/or
   9 modify it under the terms of the GNU General Public
  10 License as published by the Free Software Foundation; either
  11 version 3 of the License, or (at your option) any later version.
  12
  13 Libgfortran is distributed in the hope that it will be useful,
  14 but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16 GNU General Public License for more details.
  17
  18 Under Section 7 of GPL version 3, you are granted additional
  19 permissions described in the GCC Runtime Library Exception, version
  20 3.1, as published by the Free Software Foundation.
  21
  22 You should have received a copy of the GNU General Public License and
  23 a copy of the GCC Runtime Library Exception along with this program;
  24 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  25 <http://www.gnu.org/licenses/>.  */
  26
  27
  28 /*
  29
  30 This file is included from both the FE and the runtime library code.
  31 Operations are generalized using GMP/MPFR functions. When included from
  32 libgfortran, these should be overridden using macros which will use native
  33 operations conforming to the same API. From the FE, the GMP/MPFR functions can
  34 be used as-is.
  35
  36 The following macros are used and must be defined, unless listed as [optional]:
  37
  38 FTYPE
  39     Type name for the real-valued parameter.
  40     Variables of this type are constructed/destroyed using mpfr_init()
  41     and mpfr_clear.
  42
  43 RETTYPE
  44     Return type of the functions.
  45
  46 RETURN(x)
  47     Insert code to return a value.
  48     The parameter x is the result variable, which was also the input parameter.
  49
  50 ITYPE
  51     Type name for integer types.
  52
  53 SIND, COSD, TRIGD
  54     Names for the degree-valued trig functions defined by this module.
  55
  56 ENABLE_SIND, ENABLE_COSD, ENABLE_TAND
  57     Whether the degree-valued trig functions can be enabled.
  58
  59 ERROR_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
  64 ISFINITE(x)
  65     Whether x is a regular number or zero (not inf or NaN).
  66
  67 D2R(x)
  68     Convert x from radians to degrees.
  69
  70 SET_COSD30(x)
  71     Set x to COSD(30), or equivalently, SIND(60).
  72
  73 TINY_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
  77 COSD_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
  81 SIND_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
  91 RETTYPE
  92 SIND (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
 237 RETTYPE
 238 COSD (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
 385 RETTYPE
 386 TAND (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: */