123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493 |
- /* Implementation of the degree trignometric functions COSD, SIND, TAND.
- Copyright (C) 2020-2022 Free Software Foundation, Inc.
- Contributed by Steven G. Kargl <kargl@gcc.gnu.org>
- and Fritz Reese <foreese@gcc.gnu.org>
- This file is part of the GNU Fortran runtime library (libgfortran).
- Libgfortran is free software; you can redistribute it and/or
- modify it under the terms of the GNU General Public
- License as published by the Free Software Foundation; either
- version 3 of the License, or (at your option) any later version.
- Libgfortran is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
- Under Section 7 of GPL version 3, you are granted additional
- permissions described in the GCC Runtime Library Exception, version
- 3.1, as published by the Free Software Foundation.
- You should have received a copy of the GNU General Public License and
- a copy of the GCC Runtime Library Exception along with this program;
- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
- <http://www.gnu.org/licenses/>. */
- /*
- This file is included from both the FE and the runtime library code.
- Operations are generalized using GMP/MPFR functions. When included from
- libgfortran, these should be overridden using macros which will use native
- operations conforming to the same API. From the FE, the GMP/MPFR functions can
- be used as-is.
- The following macros are used and must be defined, unless listed as [optional]:
- FTYPE
- Type name for the real-valued parameter.
- Variables of this type are constructed/destroyed using mpfr_init()
- and mpfr_clear.
- RETTYPE
- Return type of the functions.
- RETURN(x)
- Insert code to return a value.
- The parameter x is the result variable, which was also the input parameter.
- ITYPE
- Type name for integer types.
- SIND, COSD, TRIGD
- Names for the degree-valued trig functions defined by this module.
- ENABLE_SIND, ENABLE_COSD, ENABLE_TAND
- Whether the degree-valued trig functions can be enabled.
- ERROR_RETURN(f, k, x)
- If ENABLE_<xxx>D is not defined, this is substituted to assert an
- error condition for function f, kind k, and parameter x.
- The function argument is one of {sind, cosd, tand}.
- ISFINITE(x)
- Whether x is a regular number or zero (not inf or NaN).
- D2R(x)
- Convert x from radians to degrees.
- SET_COSD30(x)
- Set x to COSD(30), or equivalently, SIND(60).
- TINY_LITERAL [optional]
- Value subtracted from 1 to cause raise INEXACT for COSD(x) for x << 1.
- If not set, COSD(x) for x <= COSD_SMALL_LITERAL simply returns 1.
- COSD_SMALL_LITERAL [optional]
- Value such that x <= COSD_SMALL_LITERAL implies COSD(x) = 1 to within the
- precision of FTYPE. If not set, this condition is not checked.
- SIND_SMALL_LITERAL [optional]
- Value such that x <= SIND_SMALL_LITERAL implies SIND(x) = D2R(x) to within
- the precision of FTYPE. If not set, this condition is not checked.
- */
- #ifdef SIND
- /* Compute sind(x) = sin(x * pi / 180). */
- RETTYPE
- SIND (FTYPE x)
- {
- #ifdef ENABLE_SIND
- if (ISFINITE (x))
- {
- FTYPE s, one;
- /* sin(-x) = - sin(x). */
- mpfr_init (s);
- mpfr_init_set_ui (one, 1, GFC_RND_MODE);
- mpfr_copysign (s, one, x, GFC_RND_MODE);
- mpfr_clear (one);
- #ifdef SIND_SMALL_LITERAL
- /* sin(x) = x as x -> 0; but only for some precisions. */
- FTYPE ax;
- mpfr_init (ax);
- mpfr_abs (ax, x, GFC_RND_MODE);
- if (mpfr_cmp_ld (ax, SIND_SMALL_LITERAL) < 0)
- {
- D2R (x);
- mpfr_clear (ax);
- return x;
- }
- mpfr_swap (x, ax);
- mpfr_clear (ax);
- #else
- mpfr_abs (x, x, GFC_RND_MODE);
- #endif /* SIND_SMALL_LITERAL */
- /* Reduce angle to x in [0,360]. */
- FTYPE period;
- mpfr_init_set_ui (period, 360, GFC_RND_MODE);
- mpfr_fmod (x, x, period, GFC_RND_MODE);
- mpfr_clear (period);
- /* Special cases with exact results. */
- ITYPE n;
- mpz_init (n);
- if (mpfr_get_z (n, x, GFC_RND_MODE) == 0 && mpz_divisible_ui_p (n, 30))
- {
- /* Flip sign for odd n*pi (x is % 360 so this is only for 180).
- This respects sgn(sin(x)) = sgn(d/dx sin(x)) = sgn(cos(x)). */
- if (mpz_divisible_ui_p (n, 180))
- {
- mpfr_set_ui (x, 0, GFC_RND_MODE);
- if (mpz_cmp_ui (n, 180) == 0)
- mpfr_neg (s, s, GFC_RND_MODE);
- }
- else if (mpz_divisible_ui_p (n, 90))
- mpfr_set_si (x, (mpz_cmp_ui (n, 90) == 0 ? 1 : -1), GFC_RND_MODE);
- else if (mpz_divisible_ui_p (n, 60))
- {
- SET_COSD30 (x);
- if (mpz_cmp_ui (n, 180) >= 0)
- mpfr_neg (x, x, GFC_RND_MODE);
- }
- else
- mpfr_set_ld (x, (mpz_cmp_ui (n, 180) < 0 ? 0.5L : -0.5L),
- GFC_RND_MODE);
- }
- /* Fold [0,360] into the range [0,45], and compute either SIN() or
- COS() depending on symmetry of shifting into the [0,45] range. */
- else
- {
- bool fold_cos = false;
- if (mpfr_cmp_ui (x, 180) <= 0)
- {
- if (mpfr_cmp_ui (x, 90) <= 0)
- {
- if (mpfr_cmp_ui (x, 45) > 0)
- {
- /* x = COS(D2R(90 - x)) */
- mpfr_ui_sub (x, 90, x, GFC_RND_MODE);
- fold_cos = true;
- }
- }
- else
- {
- if (mpfr_cmp_ui (x, 135) <= 0)
- {
- mpfr_sub_ui (x, x, 90, GFC_RND_MODE);
- fold_cos = true;
- }
- else
- mpfr_ui_sub (x, 180, x, GFC_RND_MODE);
- }
- }
- else if (mpfr_cmp_ui (x, 270) <= 0)
- {
- if (mpfr_cmp_ui (x, 225) <= 0)
- mpfr_sub_ui (x, x, 180, GFC_RND_MODE);
- else
- {
- mpfr_ui_sub (x, 270, x, GFC_RND_MODE);
- fold_cos = true;
- }
- mpfr_neg (s, s, GFC_RND_MODE);
- }
- else
- {
- if (mpfr_cmp_ui (x, 315) <= 0)
- {
- mpfr_sub_ui (x, x, 270, GFC_RND_MODE);
- fold_cos = true;
- }
- else
- mpfr_ui_sub (x, 360, x, GFC_RND_MODE);
- mpfr_neg (s, s, GFC_RND_MODE);
- }
- D2R (x);
- if (fold_cos)
- mpfr_cos (x, x, GFC_RND_MODE);
- else
- mpfr_sin (x, x, GFC_RND_MODE);
- }
- mpfr_mul (x, x, s, GFC_RND_MODE);
- mpz_clear (n);
- mpfr_clear (s);
- }
- /* Return NaN for +-Inf and NaN and raise exception. */
- else
- mpfr_sub (x, x, x, GFC_RND_MODE);
- RETURN (x);
- #else
- ERROR_RETURN(sind, KIND, x);
- #endif // ENABLE_SIND
- }
- #endif // SIND
- #ifdef COSD
- /* Compute cosd(x) = cos(x * pi / 180). */
- RETTYPE
- COSD (FTYPE x)
- {
- #ifdef ENABLE_COSD
- #if defined(TINY_LITERAL) && defined(COSD_SMALL_LITERAL)
- static const volatile FTYPE tiny = TINY_LITERAL;
- #endif
- if (ISFINITE (x))
- {
- #ifdef COSD_SMALL_LITERAL
- FTYPE ax;
- mpfr_init (ax);
- mpfr_abs (ax, x, GFC_RND_MODE);
- /* No spurious underflows!. In radians, cos(x) = 1-x*x/2 as x -> 0. */
- if (mpfr_cmp_ld (ax, COSD_SMALL_LITERAL) <= 0)
- {
- mpfr_set_ui (x, 1, GFC_RND_MODE);
- #ifdef TINY_LITERAL
- /* Cause INEXACT. */
- if (!mpfr_zero_p (ax))
- mpfr_sub_d (x, x, tiny, GFC_RND_MODE);
- #endif
- mpfr_clear (ax);
- return x;
- }
- mpfr_swap (x, ax);
- mpfr_clear (ax);
- #else
- mpfr_abs (x, x, GFC_RND_MODE);
- #endif /* COSD_SMALL_LITERAL */
- /* Reduce angle to ax in [0,360]. */
- FTYPE period;
- mpfr_init_set_ui (period, 360, GFC_RND_MODE);
- mpfr_fmod (x, x, period, GFC_RND_MODE);
- mpfr_clear (period);
- /* Special cases with exact results.
- Return negative zero for cosd(270) for consistency with libm cos(). */
- ITYPE n;
- mpz_init (n);
- if (mpfr_get_z (n, x, GFC_RND_MODE) == 0 && mpz_divisible_ui_p (n, 30))
- {
- if (mpz_divisible_ui_p (n, 180))
- mpfr_set_si (x, (mpz_cmp_ui (n, 180) == 0 ? -1 : 1),
- GFC_RND_MODE);
- else if (mpz_divisible_ui_p (n, 90))
- mpfr_set_zero (x, 0);
- else if (mpz_divisible_ui_p (n, 60))
- {
- mpfr_set_ld (x, 0.5, GFC_RND_MODE);
- if (mpz_cmp_ui (n, 60) != 0 && mpz_cmp_ui (n, 300) != 0)
- mpfr_neg (x, x, GFC_RND_MODE);
- }
- else
- {
- SET_COSD30 (x);
- if (mpz_cmp_ui (n, 30) != 0 && mpz_cmp_ui (n, 330) != 0)
- mpfr_neg (x, x, GFC_RND_MODE);
- }
- }
- /* Fold [0,360] into the range [0,45], and compute either SIN() or
- COS() depending on symmetry of shifting into the [0,45] range. */
- else
- {
- bool neg = false;
- bool fold_sin = false;
- if (mpfr_cmp_ui (x, 180) <= 0)
- {
- if (mpfr_cmp_ui (x, 90) <= 0)
- {
- if (mpfr_cmp_ui (x, 45) > 0)
- {
- mpfr_ui_sub (x, 90, x, GFC_RND_MODE);
- fold_sin = true;
- }
- }
- else
- {
- if (mpfr_cmp_ui (x, 135) <= 0)
- {
- mpfr_sub_ui (x, x, 90, GFC_RND_MODE);
- fold_sin = true;
- }
- else
- mpfr_ui_sub (x, 180, x, GFC_RND_MODE);
- neg = true;
- }
- }
- else if (mpfr_cmp_ui (x, 270) <= 0)
- {
- if (mpfr_cmp_ui (x, 225) <= 0)
- mpfr_sub_ui (x, x, 180, GFC_RND_MODE);
- else
- {
- mpfr_ui_sub (x, 270, x, GFC_RND_MODE);
- fold_sin = true;
- }
- neg = true;
- }
- else
- {
- if (mpfr_cmp_ui (x, 315) <= 0)
- {
- mpfr_sub_ui (x, x, 270, GFC_RND_MODE);
- fold_sin = true;
- }
- else
- mpfr_ui_sub (x, 360, x, GFC_RND_MODE);
- }
- D2R (x);
- if (fold_sin)
- mpfr_sin (x, x, GFC_RND_MODE);
- else
- mpfr_cos (x, x, GFC_RND_MODE);
- if (neg)
- mpfr_neg (x, x, GFC_RND_MODE);
- }
- mpz_clear (n);
- }
- /* Return NaN for +-Inf and NaN and raise exception. */
- else
- mpfr_sub (x, x, x, GFC_RND_MODE);
- RETURN (x);
- #else
- ERROR_RETURN(cosd, KIND, x);
- #endif // ENABLE_COSD
- }
- #endif // COSD
- #ifdef TAND
- /* Compute tand(x) = tan(x * pi / 180). */
- RETTYPE
- TAND (FTYPE x)
- {
- #ifdef ENABLE_TAND
- if (ISFINITE (x))
- {
- FTYPE s, one;
- /* tan(-x) = - tan(x). */
- mpfr_init (s);
- mpfr_init_set_ui (one, 1, GFC_RND_MODE);
- mpfr_copysign (s, one, x, GFC_RND_MODE);
- mpfr_clear (one);
- #ifdef SIND_SMALL_LITERAL
- /* tan(x) = x as x -> 0; but only for some precisions. */
- FTYPE ax;
- mpfr_init (ax);
- mpfr_abs (ax, x, GFC_RND_MODE);
- if (mpfr_cmp_ld (ax, SIND_SMALL_LITERAL) < 0)
- {
- D2R (x);
- mpfr_clear (ax);
- return x;
- }
- mpfr_swap (x, ax);
- mpfr_clear (ax);
- #else
- mpfr_abs (x, x, GFC_RND_MODE);
- #endif /* SIND_SMALL_LITERAL */
- /* Reduce angle to x in [0,360]. */
- FTYPE period;
- mpfr_init_set_ui (period, 360, GFC_RND_MODE);
- mpfr_fmod (x, x, period, GFC_RND_MODE);
- mpfr_clear (period);
- /* Special cases with exact results. */
- ITYPE n;
- mpz_init (n);
- if (mpfr_get_z (n, x, GFC_RND_MODE) == 0 && mpz_divisible_ui_p (n, 45))
- {
- if (mpz_divisible_ui_p (n, 180))
- mpfr_set_zero (x, 0);
- /* Though mathematically NaN is more appropriate for tan(n*90),
- returning +/-Inf offers the advantage that 1/tan(n*90) returns 0,
- which is mathematically sound. In fact we rely on this behavior
- to implement COTAND(x) = 1 / TAND(x).
- */
- else if (mpz_divisible_ui_p (n, 90))
- mpfr_set_inf (x, mpz_cmp_ui (n, 90) == 0 ? 0 : 1);
- else
- {
- mpfr_set_ui (x, 1, GFC_RND_MODE);
- if (mpz_cmp_ui (n, 45) != 0 && mpz_cmp_ui (n, 225) != 0)
- mpfr_neg (x, x, GFC_RND_MODE);
- }
- }
- else
- {
- /* Fold [0,360] into the range [0,90], and compute TAN(). */
- if (mpfr_cmp_ui (x, 180) <= 0)
- {
- if (mpfr_cmp_ui (x, 90) > 0)
- {
- mpfr_ui_sub (x, 180, x, GFC_RND_MODE);
- mpfr_neg (s, s, GFC_RND_MODE);
- }
- }
- else
- {
- if (mpfr_cmp_ui (x, 270) <= 0)
- {
- mpfr_sub_ui (x, x, 180, GFC_RND_MODE);
- }
- else
- {
- mpfr_ui_sub (x, 360, x, GFC_RND_MODE);
- mpfr_neg (s, s, GFC_RND_MODE);
- }
- }
- D2R (x);
- mpfr_tan (x, x, GFC_RND_MODE);
- }
- mpfr_mul (x, x, s, GFC_RND_MODE);
- mpz_clear (n);
- mpfr_clear (s);
- }
- /* Return NaN for +-Inf and NaN and raise exception. */
- else
- mpfr_sub (x, x, x, GFC_RND_MODE);
- RETURN (x);
- #else
- ERROR_RETURN(tand, KIND, x);
- #endif // ENABLE_TAND
- }
- #endif // TAND
- /* vim: set ft=c: */
|