Fix PR fortran/93871 and re-implement degree-valued trigonometric intrinsics.

2020-04-01  Fritz Reese  <foreese@gcc.gnu.org>
	    Steven G. Kargl  <kargl@gcc.gnu.org>

gcc/fortran/ChangeLog

	PR fortran/93871
	* gfortran.h (GFC_ISYM_ACOSD, GFC_ISYM_ASIND, GFC_ISYM_ATAN2D,
	GFC_ISYM_ATAND, GFC_ISYM_COSD, GFC_ISYM_COTAND, GFC_ISYM_SIND,
	GFC_ISYM_TAND): New.
	* intrinsic.c (add_functions): Remove check for flag_dec_math.
	Give degree trig functions simplification and name resolution
	functions (e.g, gfc_simplify_atrigd () and gfc_resolve_atrigd ()).
	(do_simplify): Remove special casing of degree trig functions.
	* intrinsic.h (gfc_simplify_acosd, gfc_simplify_asind,
	gfc_simplify_atand, gfc_simplify_cosd, gfc_simplify_cotand,
	gfc_simplify_sind, gfc_simplify_tand, gfc_resolve_trigd2): Add new
	prototypes.
	(gfc_simplify_atrigd, gfc_simplify_trigd, gfc_resolve_cotan,
	resolve_atrigd): Remove prototypes of deleted functions.
	* iresolve.c (is_trig_resolved, copy_replace_function_shallow,
	gfc_resolve_cotan, get_radians, get_degrees, resolve_trig_call,
	gfc_resolve_atrigd, gfc_resolve_atan2d): Delete functions.
	(gfc_resolve_trigd, gfc_resolve_trigd2): Resolve to library functions.
	* simplify.c (rad2deg, deg2rad, gfc_simplify_acosd, gfc_simplify_asind,
	gfc_simplify_atand, gfc_simplify_atan2d, gfc_simplify_cosd,
	gfc_simplify_sind, gfc_simplify_tand, gfc_simplify_cotand): New
	functions.
	(gfc_simplify_atan2): Fix error message.
	(simplify_trig_call, gfc_simplify_trigd, gfc_simplify_atrigd,
	radians_f): Delete functions.
	* trans-intrinsic.c: Add LIB_FUNCTION decls for sind, cosd, tand.
	(rad2deg, gfc_conv_intrinsic_atrigd, gfc_conv_intrinsic_cotan,
	gfc_conv_intrinsic_cotand, gfc_conv_intrinsic_atan2d): New functions.
	(gfc_conv_intrinsic_function): Handle ACOSD, ASIND, ATAND, COTAN,
	COTAND, ATAN2D.
	* trigd_fe.inc: New file. Included by simplify.c to implement
	simplify_sind, simplify_cosd, simplify_tand with code common to the
	libgfortran implementation.

gcc/testsuite/ChangeLog

	PR fortran/93871
	* gfortran.dg/dec_math.f90: Extend coverage to real(10) and real(16).
	* gfortran.dg/dec_math_2.f90: New test.
	* gfortran.dg/dec_math_3.f90: Likewise.
	* gfortran.dg/dec_math_4.f90: Likewise.
	* gfortran.dg/dec_math_5.f90: Likewise.

libgfortran/ChangeLog

	PR fortran/93871
	* Makefile.am, Makefile.in: New make rule for intrinsics/trigd.c.
	* gfortran.map: New routines for {sind, cosd, tand}X{r4, r8, r10, r16}.
	* intrinsics/trigd.c, intrinsics/trigd_lib.inc, intrinsics/trigd.inc:
	New files. Defines native degree-valued trig functions.

7 files changed, 852 insertions, 5 deletions

diff --git a/libgfortran/ChangeLog b/libgfortran/ChangeLog
index bef6306..e33d349 100644
--- a/libgfortran/ChangeLog
+++ b/libgfortran/ChangeLog
@@ -1,3 +1,12 @@
+2020-04-01  Fritz Reese  <foreese@gcc.gnu.org>
+	    Steven G. Kargl  <kargl@gcc.gnu.org>
+
+	PR fortran/93871
+	* Makefile.am, Makefile.in: New make rule for intrinsics/trigd.c.
+	* gfortran.map: New routines for {sind, cosd, tand}X{r4, r8, r10, r16}.
+	* intrinsics/trigd.c, intrinsics/trigd_lib.inc, intrinsics/trigd.inc:
+	New files. Defines native degree-valued trig functions.
+
 2020-02-18  Thomas Koenig  <tkoenig@gcc.gnu.org>
 
 	PR fortran/93599
diff --git a/libgfortran/Makefile.am b/libgfortran/Makefile.am
index 295a2d4..8ca0f6c 100644
--- a/libgfortran/Makefile.am
+++ b/libgfortran/Makefile.am
@@ -141,6 +141,7 @@ intrinsics/reshape_generic.c \
 intrinsics/reshape_packed.c \
 intrinsics/selected_int_kind.f90 \
 intrinsics/selected_real_kind.f90 \
+intrinsics/trigd.c \
 intrinsics/unpack_generic.c \
 runtime/in_pack_generic.c \
 runtime/in_unpack_generic.c
diff --git a/libgfortran/Makefile.in b/libgfortran/Makefile.in
index a6804c1..97a978a 100644
--- a/libgfortran/Makefile.in
+++ b/libgfortran/Makefile.in
@@ -422,8 +422,9 @@ am__objects_58 = associated.lo abort.lo args.lo cshift0.lo eoshift0.lo \
 	pack_generic.lo selected_char_kind.lo size.lo \
 	spread_generic.lo string_intrinsics.lo rand.lo random.lo \
 	reshape_generic.lo reshape_packed.lo selected_int_kind.lo \
-	selected_real_kind.lo unpack_generic.lo in_pack_generic.lo \
-	in_unpack_generic.lo $(am__objects_56) $(am__objects_57)
+	selected_real_kind.lo trigd.lo unpack_generic.lo \
+	in_pack_generic.lo in_unpack_generic.lo $(am__objects_56) \
+	$(am__objects_57)
 @IEEE_SUPPORT_TRUE@am__objects_59 = ieee_arithmetic.lo \
 @IEEE_SUPPORT_TRUE@	ieee_exceptions.lo ieee_features.lo
 am__objects_60 =
@@ -771,9 +772,9 @@ gfor_helper_src = intrinsics/associated.c intrinsics/abort.c \
 	intrinsics/rand.c intrinsics/random.c \
 	intrinsics/reshape_generic.c intrinsics/reshape_packed.c \
 	intrinsics/selected_int_kind.f90 \
-	intrinsics/selected_real_kind.f90 intrinsics/unpack_generic.c \
-	runtime/in_pack_generic.c runtime/in_unpack_generic.c \
-	$(am__append_3) $(am__append_4)
+	intrinsics/selected_real_kind.f90 intrinsics/trigd.c \
+	intrinsics/unpack_generic.c runtime/in_pack_generic.c \
+	runtime/in_unpack_generic.c $(am__append_3) $(am__append_4)
 @IEEE_SUPPORT_FALSE@gfor_ieee_src = 
 @IEEE_SUPPORT_TRUE@gfor_ieee_src = \
 @IEEE_SUPPORT_TRUE@ieee/ieee_arithmetic.F90 \
@@ -2252,6 +2253,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/time.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/transfer.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/transfer128.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/trigd.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/umask.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/unit.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/unix.Plo@am__quote@
@@ -6404,6 +6406,13 @@ reshape_packed.lo: intrinsics/reshape_packed.c
 @AMDEP_TRUE@@am__fastdepCC_FALSE@	DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCC_FALSE@	$(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o reshape_packed.lo `test -f 'intrinsics/reshape_packed.c' || echo '$(srcdir)/'`intrinsics/reshape_packed.c
 
+trigd.lo: intrinsics/trigd.c
+@am__fastdepCC_TRUE@	$(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT trigd.lo -MD -MP -MF $(DEPDIR)/trigd.Tpo -c -o trigd.lo `test -f 'intrinsics/trigd.c' || echo '$(srcdir)/'`intrinsics/trigd.c
+@am__fastdepCC_TRUE@	$(AM_V_at)$(am__mv) $(DEPDIR)/trigd.Tpo $(DEPDIR)/trigd.Plo
+@AMDEP_TRUE@@am__fastdepCC_FALSE@	$(AM_V_CC)source='intrinsics/trigd.c' object='trigd.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@	DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@	$(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o trigd.lo `test -f 'intrinsics/trigd.c' || echo '$(srcdir)/'`intrinsics/trigd.c
+
 unpack_generic.lo: intrinsics/unpack_generic.c
 @am__fastdepCC_TRUE@	$(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT unpack_generic.lo -MD -MP -MF $(DEPDIR)/unpack_generic.Tpo -c -o unpack_generic.lo `test -f 'intrinsics/unpack_generic.c' || echo '$(srcdir)/'`intrinsics/unpack_generic.c
 @am__fastdepCC_TRUE@	$(AM_V_at)$(am__mv) $(DEPDIR)/unpack_generic.Tpo $(DEPDIR)/unpack_generic.Plo
diff --git a/libgfortran/gfortran.map b/libgfortran/gfortran.map
index 3601bc2..ebf1a6f 100644
--- a/libgfortran/gfortran.map
+++ b/libgfortran/gfortran.map
@@ -1606,4 +1606,16 @@ GFORTRAN_9.2 {
 GFORTRAN_10 {
   global:
   _gfortran_os_error_at;
+  _gfortran_sind_r4;
+  _gfortran_sind_r8;
+  _gfortran_sind_r10;
+  _gfortran_sind_r16;
+  _gfortran_cosd_r4;
+  _gfortran_cosd_r8;
+  _gfortran_cosd_r10;
+  _gfortran_cosd_r16;
+  _gfortran_tand_r4;
+  _gfortran_tand_r8;
+  _gfortran_tand_r10;
+  _gfortran_tand_r16;
 } GFORTRAN_9.2;
diff --git a/libgfortran/intrinsics/trigd.c b/libgfortran/intrinsics/trigd.c
new file mode 100644
index 0000000..8169906
--- /dev/null
+++ b/libgfortran/intrinsics/trigd.c
@@ -0,0 +1,205 @@
+/* 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>
+
+
+/*
+   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)).
+
+ */
+
+/* Build _gfortran_sind_r4, _gfortran_cosd_r4, and _gfortran_tand_r4  */
+
+#define FTYPE       GFC_REAL_4
+#define SIND        sind_r4
+#define COSD        cosd_r4
+#define TAND        tand_r4
+#define SUFFIX(x)   x ## f
+
+#define TINY        0x1.p-100f	/* ~= 7.889e-31 */
+#define COSD_SMALL  0x1.p-7f	/*  = 7.8125e-3 */
+#define SIND_SMALL  0x1.p-5f	/*  = 3.125e-2 */
+#define COSD30      8.66025388e-01f
+
+#define PIO180H     1.74560547e-02f	/* high 12 bits.  */
+#define PIO180L    -2.76216747e-06f	/* Next 24 bits.  */
+
+#include "trigd_lib.inc"
+
+#undef FTYPE
+#undef TINY
+#undef COSD_SMALL
+#undef SIND_SMALL
+#undef COSD30
+#undef PIO180H
+#undef PIO180L
+#undef SIND
+#undef COSD
+#undef TAND
+#undef SUFFIX
+
+
+/* Build _gfortran_sind_r8, _gfortran_cosd_r8, and _gfortran_tand_r8.  */
+
+#define FTYPE       GFC_REAL_8
+#define SIND        sind_r8
+#define COSD        cosd_r8
+#define TAND        tand_r8
+#define SUFFIX(x)   x
+
+#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.  */
+
+#include "trigd_lib.inc"
+
+#undef FTYPE
+#undef TINY
+#undef COSD_SMALL
+#undef SIND_SMALL
+#undef COSD30
+#undef PIO180H
+#undef PIO180L
+#undef SIND
+#undef COSD
+#undef TAND
+#undef SUFFIX
+
+
+/* Build _gfortran_sind_r10, _gfortran_cosd_r10, and _gfortran_tand_r10.  */
+
+#ifdef HAVE_GFC_REAL_10
+
+#define FTYPE       GFC_REAL_10
+#define SIND        sind_r10
+#define COSD        cosd_r10
+#define TAND        tand_r10
+#define SUFFIX(x)   x ## l	/* L */
+
+#define TINY        0x1.p-16400L	/* ~= 1.28e-4937 (min exp -16494) */
+#define COSD_SMALL  0x1.p-26L	/* ~= 1.490e-8 */
+#undef  SIND_SMALL		/* not precise */
+#define COSD30       8.66025403784438646787e-01L
+
+#define PIO180H     1.74532925229868851602e-02L	/* high 32 bits */
+#define PIO180L    -3.04358939097084072823e-12L	/* Next 64 bits */
+
+#include "trigd_lib.inc"
+#undef FTYPE
+#undef TINY
+#undef COSD_SMALL
+#undef SIND_SMALL
+#undef COSD30
+#undef PIO180H
+#undef PIO180L
+#undef SIND
+#undef COSD
+#undef TAND
+#undef SUFFIX
+#endif /* HAVE_GFC_REAL_10 */
+
+
+/* Build _gfortran_sind_r16, _gfortran_cosd_r16, and _gfortran_tand_r16.  */
+
+#ifdef HAVE_GFC_REAL_16
+
+#define FTYPE       GFC_REAL_16
+#define SIND        sind_r16
+#define COSD        cosd_r16
+#define TAND        tand_r16
+
+#ifdef GFC_REAL_16_IS_FLOAT128	/* libquadmath.  */
+#define SUFFIX(x) x ## q
+#else
+#define SUFFIX(x) x ## l
+#endif /* GFC_REAL_16_IS_FLOAT128  */
+
+#define TINY        SUFFIX(0x1.p-16400)	/* ~= 1.28e-4937 */
+#define COSD_SMALL  SUFFIX(0x1.p-51)	/* ~= 4.441e-16 */
+#undef  SIND_SMALL		/* not precise */
+#define COSD30      SUFFIX(8.66025403784438646763723170752936183e-01)
+#define PIO180H     SUFFIX(1.74532925199433197605003442731685936e-02)
+#define PIO180L     SUFFIX(-2.39912634365882824665106671063098954e-17)
+
+#include "trigd_lib.inc"
+
+#undef FTYPE
+#undef COSD_SMALL
+#undef SIND_SMALL
+#undef COSD30
+#undef PIO180H
+#undef PIO180L
+#undef PIO180
+#undef D2R
+#undef CPYSGN
+#undef FABS
+#undef FMOD
+#undef SIN
+#undef COS
+#undef TAN
+#undef SIND
+#undef COSD
+#undef TAND
+#undef SUFFIX
+#endif /* HAVE_GFC_REAL_16 */
diff --git a/libgfortran/intrinsics/trigd.inc b/libgfortran/intrinsics/trigd.inc
new file mode 100644
index 0000000..98bfae7
--- /dev/null
+++ b/libgfortran/intrinsics/trigd.inc
@@ -0,0 +1,464 @@
+/* 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 and GMP/FMPR functions are used and must be defined.
+
+
+Types and names:
+
+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.
+
+
+Literal values:
+
+TINY [optional]
+    Value subtracted from 1 to cause rase INEXACT for COSD(x)
+    for x << 1. If not set, COSD(x) for x <= COSD_SMALL simply returns 1.
+
+COSD_SMALL [optional]
+    Value such that x <= COSD_SMALL implies COSD(x) = 1 to within the
+    precision of FTYPE. If not set, this condition is not checked.
+
+SIND_SMALL [optional]
+    Value such that x <= SIND_SMALL implies SIND(x) = D2R(x) to within
+    the precision of FTYPE. If not set, this condition is not checked.
+
+COSD30
+    Value of SIND(60) and COSD(30).
+
+*/
+
+
+/* Compute sind(x) = sin(x * pi / 180). */
+
+RETTYPE
+SIND (FTYPE x)
+{
+  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
+      /* 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) < 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 */
+
+      /* 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);
+}
+
+
+/* Compute cosd(x) = cos(x * pi / 180).  */
+
+RETTYPE
+COSD (FTYPE x)
+{
+#if defined(TINY) && defined(COSD_SMALL)
+  static const volatile FTYPE tiny = TINY;
+#endif
+
+  if (ISFINITE (x))
+    {
+#ifdef COSD_SMALL
+      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) <= 0)
+	{
+	  mpfr_set_ui (x, 1, GFC_RND_MODE);
+#ifdef TINY
+	  /* 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 */
+
+      /* 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);
+}
+
+
+/* Compute tand(x) = tan(x * pi / 180).  */
+
+RETTYPE
+TAND (FTYPE x)
+{
+  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
+      /* 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) < 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 */
+
+      /* 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);
+}
+
+/* vim: set ft=c: */
diff --git a/libgfortran/intrinsics/trigd_lib.inc b/libgfortran/intrinsics/trigd_lib.inc
new file mode 100644
index 0000000..b6d4145
--- /dev/null
+++ b/libgfortran/intrinsics/trigd_lib.inc
@@ -0,0 +1,147 @@
+/* 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:
+
+FTYPE             -- floating point type
+SIND, COSD, TAND  -- names of the functions to define
+SUFFIX(x)         -- add a literal suffix for floating point constants (f, ...)
+
+COSD_SMALL [optional] -- for x <= COSD_SMALL, COSD(x) = 1 if set
+TINY       [optional] -- subtract from 1 under the above condition 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
+
+ */
+
+#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 = SUFFIX(copysign)((op1), (op2))
+#define mpfr_fmod(rop, x, d, rnd) (rop = SUFFIX(fmod)((x), (d)))
+#define mpfr_abs(rop, op, rnd) (rop = SUFFIX(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 = SUFFIX(copysign)(0, (s)))
+#define mpfr_set_inf(rop, s) (rop = ((s)*-2 + 1) * INFINITY)
+#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 = SUFFIX(sin)(x))
+#define mpfr_cos(rop, x, rnd) (rop = SUFFIX(cos)(x))
+#define mpfr_tan(rop, x, rnd) (rop = SUFFIX(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 FMA(x,y,z)  SUFFIX(fma)((x), (y), (z))
+#define D2R(x) (x = FMA((x), PIO180H, (x) * PIO180L))
+
+#define SET_COSD30(x) (x = COSD30)
+
+
+extern FTYPE SIND (FTYPE);
+export_proto (SIND);
+
+extern FTYPE COSD (FTYPE);
+export_proto (COSD);
+
+extern FTYPE TAND (FTYPE);
+export_proto (TAND);
+
+#include "trigd.inc"
+
+#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: */