Merge from trunk at:

	commit 56638b9b1853666f575928f8baf17f70e4ed3517
	Author: GCC Administrator <gccadmin@gcc.gnu.org>
	Date:   Wed Jun 17 00:16:36 2020 +0000

	    Daily bump.

4 files changed, 1086 insertions, 0 deletions

diff --git a/libgfortran/intrinsics/c99_functions.c b/libgfortran/intrinsics/c99_functions.c
index 3ec5aeb..b75d177 100644
--- a/libgfortran/intrinsics/c99_functions.c
+++ b/libgfortran/intrinsics/c99_functions.c
@@ -229,6 +229,17 @@ ceilf (float x)
 }
 #endif
 
+#if !defined(HAVE_COPYSIGN) && defined(HAVE_INLINE_BUILTIN_COPYSIGN)
+#define HAVE_COPYSIGN 1
+double copysign (double x, double y);
+
+double
+copysign (double x, double y)
+{
+  return __builtin_copysign (x, y);
+}
+#endif
+
 #ifndef HAVE_COPYSIGNF
 #define HAVE_COPYSIGNF 1
 float copysignf (float x, float y);
@@ -240,6 +251,17 @@ copysignf (float x, float y)
 }
 #endif
 
+#if !defined(HAVE_COPYSIGNL) && defined(HAVE_INLINE_BUILTIN_COPYSIGNL)
+#define HAVE_COPYSIGNL 1
+long double copysignl (long double x, long double y);
+
+long double
+copysignl (long double x, long double y)
+{
+  return __builtin_copysignl (x, y);
+}
+#endif
+
 #ifndef HAVE_COSF
 #define HAVE_COSF 1
 float cosf (float x);
@@ -273,6 +295,17 @@ expf (float x)
 }
 #endif
 
+#if !defined(HAVE_FABS) && defined(HAVE_INLINE_BUILTIN_FABS)
+#define HAVE_FABS 1
+double fabs (double x);
+
+double
+fabs (double x)
+{
+  return __builtin_fabs (x);
+}
+#endif
+
 #ifndef HAVE_FABSF
 #define HAVE_FABSF 1
 float fabsf (float x);
@@ -284,6 +317,17 @@ fabsf (float x)
 }
 #endif
 
+#if !defined(HAVE_FABSL) && defined(HAVE_INLINE_BUILTIN_FABSL)
+#define HAVE_FABSL 1
+long double fabsl (long double x);
+
+long double
+fabsl (long double x)
+{
+  return __builtin_fabsl (x);
+}
+#endif
+
 #ifndef HAVE_FLOORF
 #define HAVE_FLOORF 1
 float floorf (float x);
@@ -2112,3 +2156,36 @@ lgammaf (float x)
   return (float) lgamma ((double) x);
 }
 #endif
+
+#ifndef HAVE_FMA
+#define HAVE_FMA 1
+double fma (double, double, double);
+
+double
+fma (double x, double y, double z)
+{
+  return x * y + z;
+}
+#endif
+
+#ifndef HAVE_FMAF
+#define HAVE_FMAF 1
+float fmaf (float, float, float);
+
+float
+fmaf (float x, float y, float z)
+{
+  return fma (x, y, z);
+}
+#endif
+
+#ifndef HAVE_FMAL
+#define HAVE_FMAL 1
+long double fmal (long double, long double, long double);
+
+long double
+fmal (long double x, long double y, long double z)
+{
+  return x * y + z;
+}
+#endif
diff --git a/libgfortran/intrinsics/trigd.c b/libgfortran/intrinsics/trigd.c
new file mode 100644
index 0000000..e1c51c7
--- /dev/null
+++ b/libgfortran/intrinsics/trigd.c
@@ -0,0 +1,291 @@
+/* Implementation of the degree trignometric functions COSD, SIND, TAND.
+   Copyright (C) 2020 Free Software Foundation, Inc.
+   Contributed by Steven G. Kargl <kargl@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/>.  */
+
+#include "libgfortran.h"
+
+#include <math.h>
+
+/* Body of library functions which are cannot be implemented on the current
+ * platform because it lacks a capability, such as an underlying trigonometric
+ * function (sin, cos, tan) or C99 floating-point function (fabs, fmod). */
+#define STRINGIFY_EXPAND(x) #x
+#define ERROR_RETURN(f, k, x) runtime_error (#f " is unavailable for" \
+    " REAL(KIND=" STRINGIFY_EXPAND(k) ") because the system math library" \
+    " lacks support for it"); \
+    RETURN(x)
+
+/*
+   For real x, let {x}_P or x_P be the closest representible number in the
+   floating point representation which uses P binary bits of fractional
+   precision (with IEEE rounding semantics).
+
+   Similarly, let f_P(x) be shorthand for {f(x)}_P.
+
+   Let ulp_P(x) be the unit of least precision for x: in other words the
+   maximal value of |a_P - b_P| where a_P <= x <= b_P and a_P != b_P.
+
+   Let x  ~= y <-> | x - y | <  ulp_P(x - y).
+
+   Let deg(x) be the value of x radians in degrees.
+
+   Values for each precision P were selected as follows.
+
+
+   COSD_SMALL = 2**{-N} such that for all x <= COSD_SMALL:
+
+     * cos(deg(x)) ~= 1, or equivalently:
+
+       |      1 - cos(deg(x))  | < ulp_P(1).
+
+   Unfortunately for SIND (and therefore TAND) a similar relation is only
+   possible for REAL(4) and REAL(8). With REAL(10) and REAL(16), enough
+   precision is available such that sin_P(x) != x_P for some x less than any
+   value. (There are values where this equality holds, but the distance has
+   inflection points.)
+
+   For REAL(4) and REAL(8), we can select SIND_SMALL such that:
+
+     * sin(deg(x)) ~= deg(x), or equivalently:
+
+       | deg(x) - sin(deg(x)) | < ulp_P(deg(x)).
+
+ */
+
+#ifdef HAVE_GFC_REAL_4
+
+/* Build _gfortran_sind_r4, _gfortran_cosd_r4, and _gfortran_tand_r4  */
+
+#define KIND	4
+#define TINY	0x1.p-100	/* ~= 7.889e-31 */
+#define COSD_SMALL  0x1.p-7	/*  = 7.8125e-3 */
+#define SIND_SMALL  0x1.p-5	/*  = 3.125e-2 */
+#define COSD30      8.66025388e-01
+#define PIO180H     1.74560547e-02	/* high 12 bits.  */
+#define PIO180L    -2.76216747e-06	/* Next 24 bits.  */
+
+#if defined(HAVE_FABSF) && defined(HAVE_FMODF) && defined(HAVE_COPYSIGNF)
+
+#ifdef HAVE_SINF
+#define ENABLE_SIND
+#endif
+
+#ifdef HAVE_COSF
+#define ENABLE_COSD
+#endif
+
+#ifdef HAVE_TANF
+#define ENABLE_TAND
+#endif
+
+#endif /* HAVE_FABSF && HAVE_FMODF && HAVE_COPYSIGNF */
+
+#ifdef GFC_REAL_4_INFINITY
+#define HAVE_INFINITY_KIND
+#endif
+
+#include "trigd_lib.inc"
+
+#undef KIND
+#undef TINY
+#undef COSD_SMALL
+#undef SIND_SMALL
+#undef COSD30
+#undef PIO180H
+#undef PIO180L
+#undef ENABLE_SIND
+#undef ENABLE_COSD
+#undef ENABLE_TAND
+#undef HAVE_INFINITY_KIND
+
+#endif /* HAVE_GFC_REAL_4... */
+
+
+#ifdef HAVE_GFC_REAL_8
+
+/* Build _gfortran_sind_r8, _gfortran_cosd_r8, and _gfortran_tand_r8  */
+
+#define KIND	8
+#define TINY	0x1.p-1000	/* ~= 9.33e-302 (min exp -1074) */
+#define COSD_SMALL  0x1.p-21	/* ~= 4.768e-7 */
+#define SIND_SMALL  0x1.p-19	/* ~= 9.537e-7 */
+#define COSD30      8.6602540378443860e-01
+#define PIO180H     1.7453283071517944e-02	/* high 21 bits.  */
+#define PIO180L     9.4484253514332993e-09	/* Next 53 bits.  */
+
+#if defined(HAVE_FABS) && defined(HAVE_FMOD) && defined(HAVE_COPYSIGN)
+
+#ifdef HAVE_SIN
+#define ENABLE_SIND
+#endif
+
+#ifdef HAVE_COS
+#define ENABLE_COSD
+#endif
+
+#ifdef HAVE_TAN
+#define ENABLE_TAND
+#endif
+
+#endif /* HAVE_FABS && HAVE_FMOD && HAVE_COPYSIGN */
+
+#ifdef GFC_REAL_8_INFINITY
+#define HAVE_INFINITY_KIND
+#endif
+
+#include "trigd_lib.inc"
+
+#undef KIND
+#undef TINY
+#undef COSD_SMALL
+#undef SIND_SMALL
+#undef COSD30
+#undef PIO180H
+#undef PIO180L
+#undef ENABLE_SIND
+#undef ENABLE_COSD
+#undef ENABLE_TAND
+#undef HAVE_INFINITY_KIND
+
+#endif /* HAVE_GFC_REAL_8... */
+
+
+#ifdef HAVE_GFC_REAL_10
+
+/* Build _gfortran_sind_r10, _gfortran_cosd_r10, and _gfortran_tand_r10  */
+
+#define KIND	10
+#define TINY	0x1.p-16400	/* ~= 1.28e-4937 (min exp -16494) */
+#define COSD_SMALL  0x1.p-26	/* ~= 1.490e-8 */
+#undef  SIND_SMALL		/* not precise */
+#define COSD30      8.66025403784438646787e-01
+#define PIO180H     1.74532925229868851602e-02	/* high 32 bits */
+#define PIO180L    -3.04358939097084072823e-12	/* Next 64 bits */
+
+#if defined(HAVE_FABSL) && defined(HAVE_FMODL) && defined(HAVE_COPYSIGNL)
+
+#ifdef HAVE_SINL
+#define ENABLE_SIND
+#endif
+
+#ifdef HAVE_COSL
+#define ENABLE_COSD
+#endif
+
+#ifdef HAVE_TANL
+#define ENABLE_TAND
+#endif
+
+#endif /* HAVE_FABSL && HAVE_FMODL && HAVE_COPYSIGNL */
+
+#ifdef GFC_REAL_10_INFINITY
+#define HAVE_INFINITY_KIND
+#endif
+
+#include "trigd_lib.inc"
+
+#undef KIND
+#undef TINY
+#undef COSD_SMALL
+#undef SIND_SMALL
+#undef COSD30
+#undef PIO180H
+#undef PIO180L
+#undef ENABLE_SIND
+#undef ENABLE_COSD
+#undef ENABLE_TAND
+#undef HAVE_INFINITY_KIND
+
+#endif /* HAVE_GFC_REAL_10 */
+
+
+#ifdef HAVE_GFC_REAL_16
+
+/* Build _gfortran_sind_r16, _gfortran_cosd_r16, and _gfortran_tand_r16  */
+
+#define KIND	16
+#define TINY	0x1.p-16400	/* ~= 1.28e-4937 */
+#undef  SIND_SMALL		/* not precise */
+
+#if GFC_REAL_16_DIGITS == 64
+/* 80 bit precision, use constants from REAL(10).  */
+#define COSD_SMALL  0x1.p-26	/* ~= 1.490e-8 */
+#define COSD30      8.66025403784438646787e-01
+#define PIO180H     1.74532925229868851602e-02	/* high 32 bits */
+#define PIO180L    -3.04358939097084072823e-12	/* Next 64 bits */
+
+#else
+/* Proper float128 precision.  */
+#define COSD_SMALL  0x1.p-51	/* ~= 4.441e-16 */
+#define COSD30      8.66025403784438646763723170752936183e-01
+#define PIO180H     1.74532925199433197605003442731685936e-02
+#define PIO180L     -2.39912634365882824665106671063098954e-17
+#endif
+
+#ifdef GFC_REAL_16_IS_LONG_DOUBLE
+
+#if defined(HAVE_FABSL) && defined(HAVE_FMODL) && defined(HAVE_COPYSIGNL)
+
+#ifdef HAVE_SINL
+#define ENABLE_SIND
+#endif
+
+#ifdef HAVE_COSL
+#define ENABLE_COSD
+#endif
+
+#ifdef HAVE_TANL
+#define ENABLE_TAND
+#endif
+
+#endif /* HAVE_FABSL && HAVE_FMODL && HAVE_COPYSIGNL */
+
+#else
+
+/* libquadmath: HAVE_*Q are never defined.  They must be available.  */
+#define ENABLE_SIND
+#define ENABLE_COSD
+#define ENABLE_TAND
+
+#endif /* GFC_REAL_16_IS_LONG_DOUBLE */
+
+#ifdef GFC_REAL_16_INFINITY
+#define HAVE_INFINITY_KIND
+#endif
+
+#include "trigd_lib.inc"
+
+#undef KIND
+#undef TINY
+#undef COSD_SMALL
+#undef SIND_SMALL
+#undef COSD30
+#undef PIO180H
+#undef PIO180L
+#undef ENABLE_SIND
+#undef ENABLE_COSD
+#undef ENABLE_TAND
+#undef HAVE_INFINITY_KIND
+
+#endif /* HAVE_GFC_REAL_16 */
diff --git a/libgfortran/intrinsics/trigd.inc b/libgfortran/intrinsics/trigd.inc
new file mode 100644
index 0000000..ed228e8
--- /dev/null
+++ b/libgfortran/intrinsics/trigd.inc
@@ -0,0 +1,493 @@
+/* Implementation of the degree trignometric functions COSD, SIND, TAND.
+   Copyright (C) 2020 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: */
diff --git a/libgfortran/intrinsics/trigd_lib.inc b/libgfortran/intrinsics/trigd_lib.inc
new file mode 100644
index 0000000..e90f9de
--- /dev/null
+++ b/libgfortran/intrinsics/trigd_lib.inc
@@ -0,0 +1,225 @@
+/* Stub for defining degree-valued trigonometric functions in libgfortran.
+   Copyright (C) 2020 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 replaces all GMP/MPFR functions used by trigd.inc with native versions.
+The precision is defined by FTYPE defined before including this file.
+The module which includes this file must define the following:
+
+KIND               -- floating point kind (4, 8, 10, 16)
+HAVE_INFINITY_KIND -- defined iff the platform has GFC_REAL_<KIND>_INFINITY
+
+TINY       [optional] -- subtract from 1 under the above condition if set
+COSD_SMALL [optional] -- for x <= COSD_SMALL, COSD(x) = 1 if set
+SIND_SMALL [optional] -- for x <= SIND_SMALL, SIND(x) = D2R(x) if set
+COSD30                -- literal value of COSD(30) to the precision of FTYPE
+PIO180H               -- upper bits of pi/180 for FMA
+PIO180L               -- lower bits of pi/180 for FMA
+
+ */
+
+/* FTYPE := GFC_REAL_<K> */
+#define FTYPE CONCAT_EXPAND(GFC_REAL_,KIND)
+
+/* LITERAL_SUFFIX := GFC_REAL_<K>_LITERAL_SUFFIX */
+#define LITERAL_SUFFIX CONCAT_EXPAND(FTYPE,_LITERAL_SUFFIX)
+
+/* LITERAL(X) := GFC_REAL_<K>_LITERAL(X) */
+#define LITERAL(x) CONCAT_EXPAND(x,LITERAL_SUFFIX)
+
+#define SIND CONCAT_EXPAND(sind_r, KIND)
+#define COSD CONCAT_EXPAND(cosd_r, KIND)
+#define TAND CONCAT_EXPAND(tand_r, KIND)
+
+#ifdef HAVE_INFINITY_KIND
+/* GFC_REAL_X_INFINITY */
+#define INFINITY_KIND CONCAT_EXPAND(FTYPE, _INFINITY)
+#else
+/* GFC_REAL_X_HUGE */
+#define INFINITY_KIND CONCAT_EXPAND(FTYPE, _HUGE)
+#endif
+
+#define CONCAT(x,y) x ## y
+#define CONCAT_EXPAND(x,y) CONCAT(x,y)
+
+#define COPYSIGN LITERAL(copysign)
+#define FMOD LITERAL(fmod)
+#define FABS LITERAL(fabs)
+#define FMA LITERAL(fma)
+#define SIN LITERAL(sin)
+#define COS LITERAL(cos)
+#define TAN LITERAL(tan)
+
+#ifdef TINY
+#define TINY_LITERAL LITERAL(TINY)
+#endif
+
+#ifdef COSD_SMALL
+#define COSD_SMALL_LITERAL LITERAL(COSD_SMALL)
+#endif
+
+#ifdef SIND_SMALL
+#define SIND_SMALL_LITERAL LITERAL(SIND_SMALL)
+#endif
+
+#define COSD30_LITERAL LITERAL(COSD30)
+#define PIO180H_LITERAL LITERAL(PIO180H)
+#define PIO180L_LITERAL LITERAL(PIO180L)
+
+#define ITYPE int
+#define GFC_RND_MODE 0
+#define RETTYPE FTYPE
+#define RETURN(x) return (x)
+
+#define ISFINITE(x) isfinite(x)
+#define mpfr_init(x) do { } while (0)
+#define mpfr_init_set_ui(x, v, rnd) (x = (v))
+#define mpfr_clear(x) do { } while (0)
+#define mpfr_swap(x, y) do { FTYPE z = y; y = x; x = z; } while (0)
+#define mpfr_copysign(rop, op1, op2, rnd) rop = COPYSIGN((op1), (op2))
+#define mpfr_fmod(rop, x, d, rnd) (rop = FMOD((x), (d)))
+#define mpfr_abs(rop, op, rnd) (rop = FABS(op))
+#define mpfr_cmp_ld(x, y) ((x) - (y))
+#define mpfr_cmp_ui(x, n) ((x) - (n))
+#define mpfr_zero_p(x) ((x) == 0)
+#define mpfr_set(rop, x, rnd) (rop = (x))
+#define mpfr_set_zero(rop, s) (rop = COPYSIGN(0, (s)))
+#define mpfr_set_inf(rop, s) (rop = ((s)*-2 + 1) * INFINITY_KIND)
+#define mpfr_set_ui(rop, n, rnd) (rop = (n))
+#define mpfr_set_si(rop, n, rnd) (rop = (n))
+#define mpfr_set_ld(rop, x, rnd) (rop = (x))
+#define mpfr_set_si_2exp(rop, op, exp, rnd) (rop = (0x1.p##exp))
+#define mpfr_get_z(rop, x, rnd) ((rop = (int)(x)), (rop - (x)))
+#define mpfr_mul(rop, op1, op2, rnd)    (rop = ((op1) * (op2)))
+#define mpfr_sub_d(rop, op1, op2, rnd)  (rop = ((op1) - (op2)))
+#define mpfr_sub_ui(rop, op1, op2, rnd) (rop = ((op1) - (op2)))
+#define mpfr_sub(rop, op1, op2, rnd)    (rop = ((op1) - (op2)))
+#define mpfr_ui_sub(rop, op1, op2, rnd) (rop = ((op1) - (op2)))
+#define mpfr_neg(rop, op, rnd) (rop = -(op))
+#define mpfr_sin(rop, x, rnd) (rop = SIN(x))
+#define mpfr_cos(rop, x, rnd) (rop = COS(x))
+#define mpfr_tan(rop, x, rnd) (rop = TAN(x))
+
+#define mpz_init(n) do { } while (0)
+#define mpz_clear(x) do { } while (0)
+#define mpz_cmp_ui(x, y) ((x) - (y))
+#define mpz_divisible_ui_p(n, d) ((n) % (d) == 0)
+
+#define D2R(x) (x = FMA((x), PIO180H_LITERAL, (x) * PIO180L_LITERAL))
+
+#define SET_COSD30(x) (x = COSD30_LITERAL)
+
+#ifdef SIND
+extern FTYPE SIND (FTYPE);
+export_proto (SIND);
+#endif
+
+#ifdef COSD
+extern FTYPE COSD (FTYPE);
+export_proto (COSD);
+#endif
+
+#ifdef TAND
+extern FTYPE TAND (FTYPE);
+export_proto (TAND);
+#endif
+
+#include "trigd.inc"
+
+#undef FTYPE
+#undef LITERAL_SUFFIX
+#undef LITERAL
+#undef CONCAT3
+#undef CONCAT3_EXPAND
+#undef CONCAT
+#undef CONCAT_EXPAND
+#undef SIND
+#undef COSD
+#undef TAND
+#undef INFINITY_KIND
+
+#undef COPYSIGN
+#undef FMOD
+#undef FABS
+#undef FMA
+#undef SIN
+#undef COS
+#undef TAN
+
+#undef TINY_LITERAL
+#undef COSD_SMALL_LITERAL
+#undef SIND_SMALL_LITERAL
+#undef COSD30_LITERAL
+#undef PIO180H_LITERAL
+#undef PIO180L_LITERAL
+
+#undef ITYPE
+#undef GFC_RND_MODE
+#undef RETTYPE
+#undef RETURN
+
+#undef ISFINITE
+#undef mpfr_signbit
+
+#undef mpfr_init
+#undef mpfr_init_set_ui
+#undef mpfr_clear
+#undef mpfr_swap
+#undef mpfr_fmod
+#undef mpfr_abs
+#undef mpfr_cmp_ld
+#undef mpfr_cmp_ui
+#undef mpfr_zero_p
+#undef mpfr_set
+#undef mpfr_set_zero
+#undef mpfr_set_inf
+#undef mpfr_set_ui
+#undef mpfr_set_si
+#undef mpfr_set_ld
+#undef mpfr_set_si_2exp
+#undef mpfr_get_z
+#undef mpfr_mul_si
+#undef mpfr_sub_d
+#undef mpfr_sub_ui
+#undef mpfr_sub
+#undef mpfr_ui_sub
+#undef mpfr_neg
+#undef mpfr_sin
+#undef mpfr_cos
+#undef mpfr_tan
+
+#undef mpz_init
+#undef mpz_clear
+#undef mpz_cmp_ui
+#undef mpz_divisible_ui_p
+
+#undef FMA
+#undef D2R
+
+#undef SET_COSD30
+
+
+/* vim: set ft=c: */