libgfortran/intrinsics/trigd.inc

/* Implementation of the degree trignometric functions COSD, SIND, TAND.
   Copyright (C) 2020-2021 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: */