On Fri, Mar 6, 2020 at 6:18 PM Steve Kargl
<s...@troutmask.apl.washington.edu> wrote:
[...]
> TL;DR version.
>
>   Fix the simplification and handling of the degree trigonometric functions.
>   This includes fixing a number of ICEs.  See PR 93871.

An updated version of the patch is attached. Regression tests (and new
test cases) are pending.

Changes since Steve's patch of note:

* libgfortran/intrinsics/trigd.c now indirectly includes trigd.inc
(formerly trigd_inc.c) through trigd_lib.inc. trigd.inc is now written
using GMP/MPFR functions, so that the same file can be included from
the front-end for simplification. The GMP/MPFR functions are
pre-processed into native code by trigd_lib.inc. This ensures that
both the FE and the library are using the same code for resolving
these functions, preventing future maintenance woes.

* TAND(90, 270) returns +Infinity and -Infinity, instead of NaN. This
is compatible with behavior of (at least some) DEC compilers.

* COTAN[D] are implemented as -TAN[D](x + 90 degrees) rather than 1 /
COTAN[D] for speed and to avoid singularities (though the new
implementation of TAND allows 1 / TAND to return zero when TAND
returns infinity).

* SMALL thresholds and signs of some specific values are corrected.
For REAL(10) and REAL(16), the SMALL threshold for SIND(x) = D2R(x) is
eliminated, since there would be loss of precision. The other
thresholds achieve COSD(x) = 1 and SIND(x) = D2R(x) without loss of
precision.

ChangeLogs are below. I will post an update soon after I perform
regression tests, which include some new tests which I will add.

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.

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.
diff --git a/gcc/fortran/gfortran.h b/gcc/fortran/gfortran.h
index 96037629f5f..371f77be64d 100644
--- a/gcc/fortran/gfortran.h
+++ b/gcc/fortran/gfortran.h
@@ -357,6 +357,7 @@ enum gfc_isym_id
   GFC_ISYM_ACCESS,
   GFC_ISYM_ACHAR,
   GFC_ISYM_ACOS,
+  GFC_ISYM_ACOSD,
   GFC_ISYM_ACOSH,
   GFC_ISYM_ADJUSTL,
   GFC_ISYM_ADJUSTR,
@@ -369,10 +370,13 @@ enum gfc_isym_id
   GFC_ISYM_ANINT,
   GFC_ISYM_ANY,
   GFC_ISYM_ASIN,
+  GFC_ISYM_ASIND,
   GFC_ISYM_ASINH,
   GFC_ISYM_ASSOCIATED,
   GFC_ISYM_ATAN,
   GFC_ISYM_ATAN2,
+  GFC_ISYM_ATAN2D,
+  GFC_ISYM_ATAND,
   GFC_ISYM_ATANH,
   GFC_ISYM_ATOMIC_ADD,
   GFC_ISYM_ATOMIC_AND,
@@ -410,8 +414,10 @@ enum gfc_isym_id
   GFC_ISYM_CONJG,
   GFC_ISYM_CONVERSION,
   GFC_ISYM_COS,
+  GFC_ISYM_COSD,
   GFC_ISYM_COSH,
   GFC_ISYM_COTAN,
+  GFC_ISYM_COTAND,
   GFC_ISYM_COUNT,
   GFC_ISYM_CPU_TIME,
   GFC_ISYM_CSHIFT,
@@ -598,6 +604,7 @@ enum gfc_isym_id
   GFC_ISYM_SIGNAL,
   GFC_ISYM_SI_KIND,
   GFC_ISYM_SIN,
+  GFC_ISYM_SIND,
   GFC_ISYM_SINH,
   GFC_ISYM_SIZE,
   GFC_ISYM_SLEEP,
@@ -618,6 +625,7 @@ enum gfc_isym_id
   GFC_ISYM_SYSTEM,
   GFC_ISYM_SYSTEM_CLOCK,
   GFC_ISYM_TAN,
+  GFC_ISYM_TAND,
   GFC_ISYM_TANH,
   GFC_ISYM_TEAM_NUMBER,
   GFC_ISYM_THIS_IMAGE,
diff --git a/gcc/fortran/intrinsic.c b/gcc/fortran/intrinsic.c
index 3012187ddae..17f5efc6566 100644
--- a/gcc/fortran/intrinsic.c
+++ b/gcc/fortran/intrinsic.c
@@ -3281,116 +3281,130 @@ add_functions (void)
 
   make_generic ("loc", GFC_ISYM_LOC, GFC_STD_GNU);
 
-  if (flag_dec_math)
-    {
-      add_sym_1 ("acosd", GFC_ISYM_ACOS, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dr, GFC_STD_GNU,
-		 gfc_check_fn_r, gfc_simplify_atrigd, gfc_resolve_atrigd,
-		 x, BT_REAL, dr, REQUIRED);
-
-      add_sym_1 ("dacosd", GFC_ISYM_ACOS, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dd, GFC_STD_GNU,
-		 gfc_check_fn_d, gfc_simplify_atrigd, gfc_resolve_atrigd,
-		 x, BT_REAL, dd, REQUIRED);
 
-      make_generic ("acosd", GFC_ISYM_ACOS, GFC_STD_GNU);
+  /* The next of intrinsic subprogram are the degree trignometric functions.
+     These were hidden behind the -fdec-math option, but are now simply
+     included as extensions to the set of intrinsic subprograms.  */
 
-      add_sym_1 ("asind", GFC_ISYM_ASIN, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dr, GFC_STD_GNU,
-		 gfc_check_fn_r, gfc_simplify_atrigd, gfc_resolve_atrigd,
-		 x, BT_REAL, dr, REQUIRED);
+  add_sym_1 ("acosd", GFC_ISYM_ACOSD, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dr, GFC_STD_GNU,
+	     gfc_check_fn_r, gfc_simplify_acosd, gfc_resolve_trigd,
+	     x, BT_REAL, dr, REQUIRED);
 
-      add_sym_1 ("dasind", GFC_ISYM_ASIN, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dd, GFC_STD_GNU,
-		 gfc_check_fn_d, gfc_simplify_atrigd, gfc_resolve_atrigd,
-		 x, BT_REAL, dd, REQUIRED);
+  add_sym_1 ("dacosd", GFC_ISYM_ACOSD, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dd, GFC_STD_GNU,
+	     gfc_check_fn_d, gfc_simplify_acosd, gfc_resolve_trigd,
+	     x, BT_REAL, dd, REQUIRED);
 
-      make_generic ("asind", GFC_ISYM_ASIN, GFC_STD_GNU);
+  make_generic ("acosd", GFC_ISYM_ACOSD, GFC_STD_GNU);
 
-      add_sym_1 ("atand", GFC_ISYM_ATAN, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dr, GFC_STD_GNU,
-		 gfc_check_fn_r, gfc_simplify_atrigd, gfc_resolve_atrigd,
-		 x, BT_REAL, dr, REQUIRED);
+  add_sym_1 ("asind", GFC_ISYM_ASIND, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dr, GFC_STD_GNU,
+	     gfc_check_fn_r, gfc_simplify_asind, gfc_resolve_trigd,
+	     x, BT_REAL, dr, REQUIRED);
 
-      add_sym_1 ("datand", GFC_ISYM_ATAN, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dd, GFC_STD_GNU,
-		 gfc_check_fn_d, gfc_simplify_atrigd, gfc_resolve_atrigd,
-		 x, BT_REAL, dd, REQUIRED);
+  add_sym_1 ("dasind", GFC_ISYM_ASIND, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dd, GFC_STD_GNU,
+	     gfc_check_fn_d, gfc_simplify_asind, gfc_resolve_trigd,
+	     x, BT_REAL, dd, REQUIRED);
 
-      make_generic ("atand", GFC_ISYM_ATAN, GFC_STD_GNU);
+  make_generic ("asind", GFC_ISYM_ASIND, GFC_STD_GNU);
 
-      add_sym_2 ("atan2d",GFC_ISYM_ATAN2,CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dr, GFC_STD_GNU,
-		 gfc_check_atan2, gfc_simplify_atan2d, gfc_resolve_atan2d,
-		 y, BT_REAL, dr, REQUIRED, x, BT_REAL, dr, REQUIRED);
+  add_sym_1 ("atand", GFC_ISYM_ATAND, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dr, GFC_STD_GNU,
+	     gfc_check_fn_r, gfc_simplify_atand, gfc_resolve_trigd,
+	     x, BT_REAL, dr, REQUIRED);
 
-      add_sym_2 ("datan2d",GFC_ISYM_ATAN2,CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dd, GFC_STD_GNU,
-		 gfc_check_datan2, gfc_simplify_atan2d, gfc_resolve_atan2d,
-		 y, BT_REAL, dd, REQUIRED, x, BT_REAL, dd, REQUIRED);
+  add_sym_1 ("datand", GFC_ISYM_ATAND, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dd, GFC_STD_GNU,
+	     gfc_check_fn_d, gfc_simplify_atand, gfc_resolve_trigd,
+	     x, BT_REAL, dd, REQUIRED);
 
-      make_generic ("atan2d", GFC_ISYM_ATAN2, GFC_STD_GNU);
+  make_generic ("atand", GFC_ISYM_ATAND, GFC_STD_GNU);
 
-      add_sym_1 ("cosd", GFC_ISYM_COS, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dr, GFC_STD_GNU,
-		 gfc_check_fn_r, gfc_simplify_trigd, gfc_resolve_trigd,
-		 x, BT_REAL, dr, REQUIRED);
+  add_sym_2 ("atan2d", GFC_ISYM_ATAN2D, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dr, GFC_STD_GNU,
+	     gfc_check_atan2, gfc_simplify_atan2d, gfc_resolve_trigd2,
+	     y, BT_REAL, dr, REQUIRED,
+	     x, BT_REAL, dr, REQUIRED);
 
-      add_sym_1 ("dcosd", GFC_ISYM_COS, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dd, GFC_STD_GNU,
-		 gfc_check_fn_d, gfc_simplify_trigd, gfc_resolve_trigd,
-		 x, BT_REAL, dd, REQUIRED);
+  add_sym_2 ("datan2d", GFC_ISYM_ATAN2D, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dd, GFC_STD_GNU,
+	     gfc_check_datan2, gfc_simplify_atan2d, gfc_resolve_trigd2,
+	     y, BT_REAL, dd, REQUIRED,
+	     x, BT_REAL, dd, REQUIRED);
 
-      make_generic ("cosd", GFC_ISYM_COS, GFC_STD_GNU);
+  make_generic ("atan2d", GFC_ISYM_ATAN2D, GFC_STD_GNU);
 
-      add_sym_1 ("cotan", GFC_ISYM_COTAN, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dr, GFC_STD_GNU,
-		 gfc_check_fn_rc2008, gfc_simplify_cotan, gfc_resolve_cotan,
-		 x, BT_REAL, dr, REQUIRED);
+  add_sym_1 ("cosd", GFC_ISYM_COSD, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dr, GFC_STD_GNU,
+	     gfc_check_fn_r, gfc_simplify_cosd, gfc_resolve_trigd,
+	     x, BT_REAL, dr, REQUIRED);
 
-      add_sym_1 ("dcotan", GFC_ISYM_COTAN, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dd, GFC_STD_GNU,
-		 gfc_check_fn_d, gfc_simplify_cotan, gfc_resolve_cotan,
-		 x, BT_REAL, dd, REQUIRED);
+  add_sym_1 ("dcosd", GFC_ISYM_COSD, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dd, GFC_STD_GNU,
+	     gfc_check_fn_d, gfc_simplify_cosd, gfc_resolve_trigd,
+	     x, BT_REAL, dd, REQUIRED);
 
-      make_generic ("cotan", GFC_ISYM_COTAN, GFC_STD_GNU);
+  make_generic ("cosd", GFC_ISYM_COSD, GFC_STD_GNU);
 
-      add_sym_1 ("cotand", GFC_ISYM_COTAN, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dr, GFC_STD_GNU,
-		 gfc_check_fn_r, gfc_simplify_trigd, gfc_resolve_trigd,
-		 x, BT_REAL, dr, REQUIRED);
+  add_sym_1 ("cotan", GFC_ISYM_COTAN, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dr, GFC_STD_GNU,
+	     gfc_check_fn_rc2008, gfc_simplify_cotan, gfc_resolve_trigd,
+	     x, BT_REAL, dr, REQUIRED);
 
-      add_sym_1 ("dcotand",GFC_ISYM_COTAN, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dd, GFC_STD_GNU,
-		 gfc_check_fn_d, gfc_simplify_trigd, gfc_resolve_trigd,
-		 x, BT_REAL, dd, REQUIRED);
+  add_sym_1 ("dcotan", GFC_ISYM_COTAN, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dd, GFC_STD_GNU,
+	     gfc_check_fn_d, gfc_simplify_cotan, gfc_resolve_trigd,
+	     x, BT_REAL, dd, REQUIRED);
 
-      make_generic ("cotand", GFC_ISYM_COTAN, GFC_STD_GNU);
+  add_sym_1 ("ccotan", GFC_ISYM_COTAN, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_COMPLEX, dz, GFC_STD_GNU,
+	     NULL, gfc_simplify_cotan, gfc_resolve_trigd,
+	     x, BT_COMPLEX, dz, REQUIRED);
 
-      add_sym_1 ("sind", GFC_ISYM_SIN, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dr, GFC_STD_GNU,
-		 gfc_check_fn_r, gfc_simplify_trigd, gfc_resolve_trigd,
-		 x, BT_REAL, dr, REQUIRED);
+  add_sym_1 ("zcotan", GFC_ISYM_COTAN, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_COMPLEX, dd, GFC_STD_GNU,
+	     NULL, gfc_simplify_cotan, gfc_resolve_trigd,
+	     x, BT_COMPLEX, dd, REQUIRED);
 
-      add_sym_1 ("dsind", GFC_ISYM_SIN, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dd, GFC_STD_GNU,
-		 gfc_check_fn_d, gfc_simplify_trigd, gfc_resolve_trigd,
-		 x, BT_REAL, dd, REQUIRED);
+  make_generic ("cotan", GFC_ISYM_COTAN, GFC_STD_GNU);
 
-      make_generic ("sind", GFC_ISYM_SIN, GFC_STD_GNU);
+  add_sym_1 ("cotand", GFC_ISYM_COTAND, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dr, GFC_STD_GNU,
+	     gfc_check_fn_r, gfc_simplify_cotand, gfc_resolve_trigd,
+	     x, BT_REAL, dr, REQUIRED);
 
-      add_sym_1 ("tand", GFC_ISYM_TAN, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dr, GFC_STD_GNU,
-		 gfc_check_fn_r, gfc_simplify_trigd, gfc_resolve_trigd,
-		 x, BT_REAL, dr, REQUIRED);
+  add_sym_1 ("dcotand", GFC_ISYM_COTAND, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dd, GFC_STD_GNU,
+	     gfc_check_fn_d, gfc_simplify_cotand, gfc_resolve_trigd,
+	     x, BT_REAL, dd, REQUIRED);
 
-      add_sym_1 ("dtand", GFC_ISYM_TAN, CLASS_ELEMENTAL, ACTUAL_YES, BT_REAL,
-		 dd, GFC_STD_GNU,
-		 gfc_check_fn_d, gfc_simplify_trigd, gfc_resolve_trigd,
-		 x, BT_REAL, dd, REQUIRED);
+  make_generic ("cotand", GFC_ISYM_COTAND, GFC_STD_GNU);
 
-      make_generic ("tand", GFC_ISYM_TAN, GFC_STD_GNU);
-    }
+  add_sym_1 ("sind", GFC_ISYM_SIND, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dr, GFC_STD_GNU,
+	     gfc_check_fn_r, gfc_simplify_sind, gfc_resolve_trigd,
+	     x, BT_REAL, dr, REQUIRED);
+
+  add_sym_1 ("dsind", GFC_ISYM_SIND, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dd, GFC_STD_GNU,
+	     gfc_check_fn_d, gfc_simplify_sind, gfc_resolve_trigd,
+	     x, BT_REAL, dd, REQUIRED);
+
+  make_generic ("sind", GFC_ISYM_SIND, GFC_STD_GNU);
+
+  add_sym_1 ("tand", GFC_ISYM_TAND, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dr, GFC_STD_GNU,
+	     gfc_check_fn_r, gfc_simplify_tand, gfc_resolve_trigd,
+	     x, BT_REAL, dr, REQUIRED);
+
+  add_sym_1 ("dtand", GFC_ISYM_TAND, CLASS_ELEMENTAL, ACTUAL_YES,
+	     BT_REAL, dd, GFC_STD_GNU,
+	     gfc_check_fn_d, gfc_simplify_tand, gfc_resolve_trigd,
+	     x, BT_REAL, dd, REQUIRED);
+
+  make_generic ("tand", GFC_ISYM_TAND, GFC_STD_GNU);
 
   /* The following function is internally used for coarray libray functions.
      "make_from_module" makes it inaccessible for external users.  */
@@ -4566,15 +4580,6 @@ do_simplify (gfc_intrinsic_sym *specific, gfc_expr *e)
       goto finish;
     }
 
-  /* Some math intrinsics need to wrap the original expression.  */
-  if (specific->simplify.f1 == gfc_simplify_trigd
-      || specific->simplify.f1 == gfc_simplify_atrigd
-      || specific->simplify.f1 == gfc_simplify_cotan)
-    {
-      result = (*specific->simplify.f1) (e);
-      goto finish;
-    }
-
   if (specific->simplify.f1 == NULL)
     {
       result = NULL;
diff --git a/gcc/fortran/intrinsic.h b/gcc/fortran/intrinsic.h
index d0456742a56..166ae792939 100644
--- a/gcc/fortran/intrinsic.h
+++ b/gcc/fortran/intrinsic.h
@@ -237,13 +237,14 @@ bool gfc_check_unlink_sub (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_abs (gfc_expr *);
 gfc_expr *gfc_simplify_achar (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_acos (gfc_expr *);
+gfc_expr *gfc_simplify_acosd (gfc_expr *);
 gfc_expr *gfc_simplify_acosh (gfc_expr *);
 gfc_expr *gfc_simplify_adjustl (gfc_expr *);
 gfc_expr *gfc_simplify_adjustr (gfc_expr *);
 gfc_expr *gfc_simplify_aimag (gfc_expr *);
 gfc_expr *gfc_simplify_aint (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_all (gfc_expr *, gfc_expr *);
-gfc_expr *gfc_simplify_atrigd (gfc_expr *);
+gfc_expr *gfc_simplify_asind (gfc_expr *);
 gfc_expr *gfc_simplify_dint (gfc_expr *);
 gfc_expr *gfc_simplify_anint (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_dnint (gfc_expr *);
@@ -252,6 +253,7 @@ gfc_expr *gfc_simplify_any (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_asin (gfc_expr *);
 gfc_expr *gfc_simplify_asinh (gfc_expr *);
 gfc_expr *gfc_simplify_atan (gfc_expr *);
+gfc_expr *gfc_simplify_atand (gfc_expr *);
 gfc_expr *gfc_simplify_atanh (gfc_expr *);
 gfc_expr *gfc_simplify_atan2 (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_atan2d (gfc_expr *, gfc_expr *);
@@ -277,8 +279,10 @@ gfc_expr *gfc_simplify_compiler_version (void);
 gfc_expr *gfc_simplify_complex (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_conjg (gfc_expr *);
 gfc_expr *gfc_simplify_cos (gfc_expr *);
+gfc_expr *gfc_simplify_cosd (gfc_expr *);
 gfc_expr *gfc_simplify_cosh (gfc_expr *);
 gfc_expr *gfc_simplify_cotan (gfc_expr *);
+gfc_expr *gfc_simplify_cotand (gfc_expr *);
 gfc_expr *gfc_simplify_count (gfc_expr *, gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_cshift (gfc_expr *, gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_dcmplx (gfc_expr *, gfc_expr *);
@@ -404,6 +408,7 @@ gfc_expr *gfc_simplify_shifta (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_shiftl (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_shiftr (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_sin (gfc_expr *);
+gfc_expr *gfc_simplify_sind (gfc_expr *);
 gfc_expr *gfc_simplify_sinh (gfc_expr *);
 gfc_expr *gfc_simplify_size (gfc_expr *, gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_sizeof (gfc_expr *);
@@ -414,13 +419,13 @@ gfc_expr *gfc_simplify_spread (gfc_expr *, gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_sqrt (gfc_expr *);
 gfc_expr *gfc_simplify_sum (gfc_expr *, gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_tan (gfc_expr *);
+gfc_expr *gfc_simplify_tand (gfc_expr *);
 gfc_expr *gfc_simplify_tanh (gfc_expr *);
 gfc_expr *gfc_simplify_this_image (gfc_expr *, gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_tiny (gfc_expr *);
 gfc_expr *gfc_simplify_trailz (gfc_expr *);
 gfc_expr *gfc_simplify_transfer (gfc_expr *, gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_transpose (gfc_expr *);
-gfc_expr *gfc_simplify_trigd (gfc_expr *);
 gfc_expr *gfc_simplify_trim (gfc_expr *);
 gfc_expr *gfc_simplify_ubound (gfc_expr *, gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_ucobound (gfc_expr *, gfc_expr *, gfc_expr *);
@@ -473,7 +478,6 @@ void gfc_resolve_conjg (gfc_expr *, gfc_expr *);
 void gfc_resolve_cos (gfc_expr *, gfc_expr *);
 void gfc_resolve_cosh (gfc_expr *, gfc_expr *);
 void gfc_resolve_count (gfc_expr *, gfc_expr *, gfc_expr *, gfc_expr *);
-void gfc_resolve_cotan (gfc_expr *, gfc_expr *);
 void gfc_resolve_cshift (gfc_expr *, gfc_expr *, gfc_expr *, gfc_expr *);
 void gfc_resolve_ctime (gfc_expr *, gfc_expr *);
 void gfc_resolve_dble (gfc_expr *, gfc_expr *);
@@ -612,7 +616,7 @@ void gfc_resolve_time8 (gfc_expr *);
 void gfc_resolve_transfer (gfc_expr *, gfc_expr *, gfc_expr *, gfc_expr *);
 void gfc_resolve_transpose (gfc_expr *, gfc_expr *);
 void gfc_resolve_trigd (gfc_expr *, gfc_expr *);
-void gfc_resolve_atrigd (gfc_expr *, gfc_expr *);
+void gfc_resolve_trigd2 (gfc_expr *, gfc_expr *, gfc_expr *);
 void gfc_resolve_trim (gfc_expr *, gfc_expr *);
 void gfc_resolve_ttynam (gfc_expr *, gfc_expr *);
 void gfc_resolve_ubound (gfc_expr *, gfc_expr *, gfc_expr *, gfc_expr *);
diff --git a/gcc/fortran/iresolve.c b/gcc/fortran/iresolve.c
index a991c3a5ab4..7ecb6595f59 100644
--- a/gcc/fortran/iresolve.c
+++ b/gcc/fortran/iresolve.c
@@ -689,86 +689,6 @@ gfc_resolve_cosh (gfc_expr *f, gfc_expr *x)
 }
 
 
-/* Our replacement of elements of a trig call with an EXPR_OP (e.g.
-   multiplying the result or operands by a factor to convert to/from degrees)
-   will cause the resolve_* function to be invoked again when resolving the
-   freshly created EXPR_OP.  See gfc_resolve_trigd, gfc_resolve_atrigd,
-   gfc_resolve_cotan.  We must observe this and avoid recursively creating
-   layers of nested EXPR_OP expressions.  */
-
-static bool
-is_trig_resolved (gfc_expr *f)
-{
-  /* We know we've already resolved the function if we see the lib call
-     starting with '__'.  */
-  return (f->value.function.name != NULL
-	  && gfc_str_startswith (f->value.function.name, "__"));
-}
-
-/* Return a shallow copy of the function expression f.  The original expression
-   has its pointers cleared so that it may be freed without affecting the
-   shallow copy.  This is similar to gfc_copy_expr, but doesn't perform a deep
-   copy of the argument list, allowing it to be reused somewhere else,
-   setting the expression up nicely for gfc_replace_expr.  */
-
-static gfc_expr *
-copy_replace_function_shallow (gfc_expr *f)
-{
-  gfc_expr *fcopy;
-  gfc_actual_arglist *args;
-
-  /* The only thing deep-copied in gfc_copy_expr is args.  */
-  args = f->value.function.actual;
-  f->value.function.actual = NULL;
-  fcopy = gfc_copy_expr (f);
-  fcopy->value.function.actual = args;
-
-  /* Clear the old function so the shallow copy is not affected if the old
-     expression is freed.  */
-  f->value.function.name = NULL;
-  f->value.function.isym = NULL;
-  f->value.function.actual = NULL;
-  f->value.function.esym = NULL;
-  f->shape = NULL;
-  f->ref = NULL;
-
-  return fcopy;
-}
-
-
-/* Resolve cotan = cos / sin.  */
-
-void
-gfc_resolve_cotan (gfc_expr *f, gfc_expr *x)
-{
-  gfc_expr *result, *fcopy, *sin;
-  gfc_actual_arglist *sin_args;
-
-  if (is_trig_resolved (f))
-    return;
-
-  /* Compute cotan (x) = cos (x) / sin (x).  */
-  f->value.function.isym = gfc_intrinsic_function_by_id (GFC_ISYM_COS);
-  gfc_resolve_cos (f, x);
-
-  sin_args = gfc_get_actual_arglist ();
-  sin_args->expr = gfc_copy_expr (x);
-
-  sin = gfc_get_expr ();
-  sin->ts = f->ts;
-  sin->where = f->where;
-  sin->expr_type = EXPR_FUNCTION;
-  sin->value.function.isym = gfc_intrinsic_function_by_id (GFC_ISYM_SIN);
-  sin->value.function.actual = sin_args;
-  gfc_resolve_sin (sin, sin_args->expr);
-
-  /* Replace f with cos/sin - we do this in place in f for the caller.  */
-  fcopy = copy_replace_function_shallow (f);
-  result = gfc_divide (fcopy, sin);
-  gfc_replace_expr (f, result);
-}
-
-
 void
 gfc_resolve_count (gfc_expr *f, gfc_expr *mask, gfc_expr *dim, gfc_expr *kind)
 {
@@ -2912,158 +2832,6 @@ gfc_resolve_tanh (gfc_expr *f, gfc_expr *x)
 }
 
 
-/* Build an expression for converting degrees to radians.  */
-
-static gfc_expr *
-get_radians (gfc_expr *deg)
-{
-  gfc_expr *result, *factor;
-  gfc_actual_arglist *mod_args;
-
-  gcc_assert (deg->ts.type == BT_REAL);
-
-  /* Set deg = deg % 360 to avoid offsets from large angles.  */
-  factor = gfc_get_constant_expr (deg->ts.type, deg->ts.kind, &deg->where);
-  mpfr_set_d (factor->value.real, 360.0, GFC_RND_MODE);
-
-  mod_args = gfc_get_actual_arglist ();
-  mod_args->expr = deg;
-  mod_args->next = gfc_get_actual_arglist ();
-  mod_args->next->expr = factor;
-
-  result = gfc_get_expr ();
-  result->ts = deg->ts;
-  result->where = deg->where;
-  result->expr_type = EXPR_FUNCTION;
-  result->value.function.isym = gfc_intrinsic_function_by_id (GFC_ISYM_MOD);
-  result->value.function.actual = mod_args;
-
-  /* Set factor = pi / 180.  */
-  factor = gfc_get_constant_expr (deg->ts.type, deg->ts.kind, &deg->where);
-  mpfr_const_pi (factor->value.real, GFC_RND_MODE);
-  mpfr_div_ui (factor->value.real, factor->value.real, 180, GFC_RND_MODE);
-
-  /* Result is rad = (deg % 360) * (pi / 180).  */
-  result = gfc_multiply (result, factor);
-  return result;
-}
-
-
-/* Build an expression for converting radians to degrees.  */
-
-static gfc_expr *
-get_degrees (gfc_expr *rad)
-{
-  gfc_expr *result, *factor;
-  gfc_actual_arglist *mod_args;
-  mpfr_t tmp;
-
-  gcc_assert (rad->ts.type == BT_REAL);
-
-  /* Set rad = rad % 2pi to avoid offsets from large angles.  */
-  factor = gfc_get_constant_expr (rad->ts.type, rad->ts.kind, &rad->where);
-  mpfr_const_pi (factor->value.real, GFC_RND_MODE);
-  mpfr_mul_ui (factor->value.real, factor->value.real, 2, GFC_RND_MODE);
-
-  mod_args = gfc_get_actual_arglist ();
-  mod_args->expr = rad;
-  mod_args->next = gfc_get_actual_arglist ();
-  mod_args->next->expr = factor;
-
-  result = gfc_get_expr ();
-  result->ts = rad->ts;
-  result->where = rad->where;
-  result->expr_type = EXPR_FUNCTION;
-  result->value.function.isym = gfc_intrinsic_function_by_id (GFC_ISYM_MOD);
-  result->value.function.actual = mod_args;
-
-  /* Set factor = 180 / pi.  */
-  factor = gfc_get_constant_expr (rad->ts.type, rad->ts.kind, &rad->where);
-  mpfr_set_ui (factor->value.real, 180, GFC_RND_MODE);
-  mpfr_init (tmp);
-  mpfr_const_pi (tmp, GFC_RND_MODE);
-  mpfr_div (factor->value.real, factor->value.real, tmp, GFC_RND_MODE);
-  mpfr_clear (tmp);
-
-  /* Result is deg = (rad % 2pi) * (180 / pi).  */
-  result = gfc_multiply (result, factor);
-  return result;
-}
-
-
-/* Resolve a call to a trig function.  */
-
-static void
-resolve_trig_call (gfc_expr *f, gfc_expr *x)
-{
-  switch (f->value.function.isym->id)
-    {
-    case GFC_ISYM_ACOS:
-      return gfc_resolve_acos (f, x);
-    case GFC_ISYM_ASIN:
-      return gfc_resolve_asin (f, x);
-    case GFC_ISYM_ATAN:
-      return gfc_resolve_atan (f, x);
-    case GFC_ISYM_ATAN2:
-      /* NB. arg3 is unused for atan2 */
-      return gfc_resolve_atan2 (f, x, NULL);
-    case GFC_ISYM_COS:
-      return gfc_resolve_cos (f, x);
-    case GFC_ISYM_COTAN:
-      return gfc_resolve_cotan (f, x);
-    case GFC_ISYM_SIN:
-      return gfc_resolve_sin (f, x);
-    case GFC_ISYM_TAN:
-      return gfc_resolve_tan (f, x);
-    default:
-      gcc_unreachable ();
-    }
-}
-
-/* Resolve degree trig function as trigd (x) = trig (radians (x)).  */
-
-void
-gfc_resolve_trigd (gfc_expr *f, gfc_expr *x)
-{
-  if (is_trig_resolved (f))
-    return;
-
-  x = get_radians (x);
-  f->value.function.actual->expr = x;
-
-  resolve_trig_call (f, x);
-}
-
-
-/* Resolve degree inverse trig function as atrigd (x) = degrees (atrig (x)).  */
-
-void
-gfc_resolve_atrigd (gfc_expr *f, gfc_expr *x)
-{
-  gfc_expr *result, *fcopy;
-
-  if (is_trig_resolved (f))
-    return;
-
-  resolve_trig_call (f, x);
-
-  fcopy = copy_replace_function_shallow (f);
-  result = get_degrees (fcopy);
-  gfc_replace_expr (f, result);
-}
-
-
-/* Resolve atan2d(x) = degrees(atan2(x)).  */
-
-void
-gfc_resolve_atan2d (gfc_expr *f, gfc_expr *x, gfc_expr *y ATTRIBUTE_UNUSED)
-{
-  /* Note that we lose the second arg here - that's okay because it is
-     unused in gfc_resolve_atan2 anyway.  */
-  gfc_resolve_atrigd (f, x);
-}
-
-
 /* Resolve failed_images (team, kind).  */
 
 void
@@ -3298,6 +3066,30 @@ gfc_resolve_trim (gfc_expr *f, gfc_expr *string)
 }
 
 
+/* Resolve the degree trignometric functions.  This amounts to setting
+   the function return type-spec from its argument and building a
+   library function names of the form _gfortran_sind_r4.  */
+
+void
+gfc_resolve_trigd (gfc_expr *f, gfc_expr *x)
+{
+  f->ts = x->ts;
+  f->value.function.name
+    = gfc_get_string (PREFIX ("%s_%c%d"), f->value.function.isym->name,
+		      gfc_type_letter (x->ts.type), x->ts.kind);
+}
+
+
+void
+gfc_resolve_trigd2 (gfc_expr *f, gfc_expr *y, gfc_expr *x)
+{
+  f->ts = y->ts;
+  f->value.function.name
+    = gfc_get_string (PREFIX ("%s_%d"), f->value.function.isym->name,
+		      x->ts.kind);
+}
+
+
 void
 gfc_resolve_ubound (gfc_expr *f, gfc_expr *array, gfc_expr *dim, gfc_expr *kind)
 {
diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c
index 66ed925c10d..f63f63c9ef6 100644
--- a/gcc/fortran/simplify.c
+++ b/gcc/fortran/simplify.c
@@ -1107,6 +1107,91 @@ gfc_simplify_asin (gfc_expr *x)
 }
 
 
+/* Convert radians to degrees, i.e., x * 180 / pi.  */
+
+static void
+rad2deg (mpfr_t x)
+{
+  mpfr_t tmp;
+
+  mpfr_init (tmp);
+  mpfr_const_pi (tmp, GFC_RND_MODE);
+  mpfr_mul_ui (x, x, 180, GFC_RND_MODE);
+  mpfr_div (x, x, tmp, GFC_RND_MODE);
+  mpfr_clear (tmp);
+}
+
+
+/* Simplify ACOSD(X) where the returned value has units of degree.  */
+
+gfc_expr *
+gfc_simplify_acosd (gfc_expr *x)
+{
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  if (mpfr_cmp_si (x->value.real, 1) > 0
+      || mpfr_cmp_si (x->value.real, -1) < 0)
+    {
+      gfc_error ("Argument of ACOSD at %L must be between -1 and 1",
+		 &x->where);
+      return &gfc_bad_expr;
+    }
+
+  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+  mpfr_acos (result->value.real, x->value.real, GFC_RND_MODE);
+  rad2deg (result->value.real);
+
+  return range_check (result, "ACOSD");
+}
+
+
+/* Simplify asind (x) where the returned value has units of degree. */
+
+gfc_expr *
+gfc_simplify_asind (gfc_expr *x)
+{
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  if (mpfr_cmp_si (x->value.real, 1) > 0
+      || mpfr_cmp_si (x->value.real, -1) < 0)
+    {
+      gfc_error ("Argument of ASIND at %L must be between -1 and 1",
+		 &x->where);
+      return &gfc_bad_expr;
+    }
+
+  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+  mpfr_asin (result->value.real, x->value.real, GFC_RND_MODE);
+  rad2deg (result->value.real);
+
+  return range_check (result, "ASIND");
+}
+
+
+/* Simplify atand (x) where the returned value has units of degree. */
+
+gfc_expr *
+gfc_simplify_atand (gfc_expr *x)
+{
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+  mpfr_atan (result->value.real, x->value.real, GFC_RND_MODE);
+  rad2deg (result->value.real);
+
+  return range_check (result, "ATAND");
+}
+
+
 gfc_expr *
 gfc_simplify_asinh (gfc_expr *x)
 {
@@ -1208,8 +1293,8 @@ gfc_simplify_atan2 (gfc_expr *y, gfc_expr *x)
 
   if (mpfr_zero_p (y->value.real) && mpfr_zero_p (x->value.real))
     {
-      gfc_error ("If first argument of ATAN2 %L is zero, then the "
-		 "second argument must not be zero", &x->where);
+      gfc_error ("If first argument of ATAN2 at %L is zero, then the "
+		 "second argument must not be zero", &y->where);
       return &gfc_bad_expr;
     }
 
@@ -1736,172 +1821,153 @@ gfc_simplify_conjg (gfc_expr *e)
   return range_check (result, "CONJG");
 }
 
-/* Return the simplification of the constant expression in icall, or NULL
-   if the expression is not constant.  */
 
-static gfc_expr *
-simplify_trig_call (gfc_expr *icall)
-{
-  gfc_isym_id func = icall->value.function.isym->id;
-  gfc_expr *x = icall->value.function.actual->expr;
-
-  /* The actual simplifiers will return NULL for non-constant x.  */
-  switch (func)
-    {
-    case GFC_ISYM_ACOS:
-      return gfc_simplify_acos (x);
-    case GFC_ISYM_ASIN:
-      return gfc_simplify_asin (x);
-    case GFC_ISYM_ATAN:
-      return gfc_simplify_atan (x);
-    case GFC_ISYM_COS:
-      return gfc_simplify_cos (x);
-    case GFC_ISYM_COTAN:
-      return gfc_simplify_cotan (x);
-    case GFC_ISYM_SIN:
-      return gfc_simplify_sin (x);
-    case GFC_ISYM_TAN:
-      return gfc_simplify_tan (x);
-    default:
-      gfc_internal_error ("in simplify_trig_call(): Bad intrinsic");
+/* Simplify atan2d (x) where the unit is degree.  */
+
+gfc_expr *
+gfc_simplify_atan2d (gfc_expr *y, gfc_expr *x)
+{
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT || y->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  if (mpfr_zero_p (y->value.real) && mpfr_zero_p (x->value.real))
+    {
+      gfc_error ("If first argument of ATAN2D at %L is zero, then the "
+		 "second argument must not be zero", &y->where);
+      return &gfc_bad_expr;
     }
+
+  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+  mpfr_atan2 (result->value.real, y->value.real, x->value.real, GFC_RND_MODE);
+  rad2deg (result->value.real);
+
+  return range_check (result, "ATAN2D");
 }
 
-/* Convert a floating-point number from radians to degrees.  */
 
-static void
-degrees_f (mpfr_t x, mpfr_rnd_t rnd_mode)
+gfc_expr *
+gfc_simplify_cos (gfc_expr *x)
 {
-  mpfr_t tmp;
-  mpfr_init (tmp);
+  gfc_expr *result;
 
-  /* Set x = x * 180.  */
-  mpfr_mul_ui (x, x, 180, rnd_mode);
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
 
-  /* Set x = x / pi.  */
-  mpfr_const_pi (tmp, rnd_mode);
-  mpfr_div (x, x, tmp, rnd_mode);
+  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
 
-  mpfr_clear (tmp);
+  switch (x->ts.type)
+    {
+      case BT_REAL:
+	mpfr_cos (result->value.real, x->value.real, GFC_RND_MODE);
+	break;
+
+      case BT_COMPLEX:
+	gfc_set_model_kind (x->ts.kind);
+	mpc_cos (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
+	break;
+
+      default:
+	gfc_internal_error ("in gfc_simplify_cos(): Bad type");
+    }
+
+  return range_check (result, "COS");
 }
 
-/* Convert a floating-point number from degrees to radians.  */
 
 static void
-radians_f (mpfr_t x, mpfr_rnd_t rnd_mode)
+deg2rad (mpfr_t x)
 {
-  mpfr_t tmp;
-  mpfr_init (tmp);
+  mpfr_t d2r;
 
-  /* Set x = x % 360 to avoid offsets with large angles.  */
-  mpfr_set_ui (tmp, 360, rnd_mode);
-  mpfr_fmod (tmp, x, tmp, rnd_mode);
+  mpfr_init (d2r);
+  mpfr_const_pi (d2r, GFC_RND_MODE);
+  mpfr_div_ui (d2r, d2r, 180, GFC_RND_MODE);
+  mpfr_mul (x, x, d2r, GFC_RND_MODE);
+  mpfr_clear (d2r);
+}
 
-  /* Set x = x * pi.  */
-  mpfr_const_pi (tmp, rnd_mode);
-  mpfr_mul (x, x, tmp, rnd_mode);
 
-  /* Set x = x / 180.  */
-  mpfr_div_ui (x, x, 180, rnd_mode);
-
-  mpfr_clear (tmp);
-}
+/* Simplification routines for SIND, COSD, TAND.  */
+#include "trigd_fe.inc"
 
 
-/* Convert argument to radians before calling a trig function.  */
+/* Simplify COSD(X) where X has the unit of degree.  */
 
 gfc_expr *
-gfc_simplify_trigd (gfc_expr *icall)
+gfc_simplify_cosd (gfc_expr *x)
 {
-  gfc_expr *arg;
-
-  arg = icall->value.function.actual->expr;
+  gfc_expr *result;
 
-  if (arg->ts.type != BT_REAL)
-    gfc_internal_error ("in gfc_simplify_trigd(): Bad type");
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
 
-  if (arg->expr_type == EXPR_CONSTANT)
-    /* Convert constant to radians before passing off to simplifier.  */
-    radians_f (arg->value.real, GFC_RND_MODE);
+  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+  mpfr_set (result->value.real, x->value.real, GFC_RND_MODE);
+  simplify_cosd (result->value.real);
 
-  /* Let the usual simplifier take over - we just simplified the arg.  */
-  return simplify_trig_call (icall);
+  return range_check (result, "COSD");
 }
 
-/* Convert result of an inverse trig function to degrees.  */
+
+/* Simplify SIND(X) where X has the unit of degree.  */
 
 gfc_expr *
-gfc_simplify_atrigd (gfc_expr *icall)
+gfc_simplify_sind (gfc_expr *x)
 {
   gfc_expr *result;
 
-  if (icall->value.function.actual->expr->ts.type != BT_REAL)
-    gfc_internal_error ("in gfc_simplify_atrigd(): Bad type");
-
-  /* See if another simplifier has work to do first.  */
-  result = simplify_trig_call (icall);
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
 
-  if (result && result->expr_type == EXPR_CONSTANT)
-    {
-      /* Convert constant to degrees after passing off to actual simplifier.  */
-      degrees_f (result->value.real, GFC_RND_MODE);
-      return result;
-    }
+  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+  mpfr_set (result->value.real, x->value.real, GFC_RND_MODE);
+  simplify_sind (result->value.real);
 
-  /* Let gfc_resolve_atrigd take care of the non-constant case.  */
-  return NULL;
+  return range_check (result, "SIND");
 }
 
-/* Convert the result of atan2 to degrees.  */
+
+/* Simplify TAND(X) where X has the unit of degree.  */
 
 gfc_expr *
-gfc_simplify_atan2d (gfc_expr *y, gfc_expr *x)
+gfc_simplify_tand (gfc_expr *x)
 {
   gfc_expr *result;
 
-  if (x->ts.type != BT_REAL || y->ts.type != BT_REAL)
-    gfc_internal_error ("in gfc_simplify_atan2d(): Bad type");
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
 
-  if (x->expr_type == EXPR_CONSTANT && y->expr_type == EXPR_CONSTANT)
-    {
-      result = gfc_simplify_atan2 (y, x);
-      if (result != NULL)
-	{
-	  degrees_f (result->value.real, GFC_RND_MODE);
-	  return result;
-	}
-    }
+  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+  mpfr_set (result->value.real, x->value.real, GFC_RND_MODE);
+  simplify_tand (result->value.real);
 
-  /* Let gfc_resolve_atan2d take care of the non-constant case.  */
-  return NULL;
+  return range_check (result, "TAND");
 }
 
+
+/* Simplify COTAND(X) where X has the unit of degree.  */
+
 gfc_expr *
-gfc_simplify_cos (gfc_expr *x)
+gfc_simplify_cotand (gfc_expr *x)
 {
   gfc_expr *result;
 
   if (x->expr_type != EXPR_CONSTANT)
     return NULL;
 
+  /* Implement COTAND = -TAND(x+90).
+     TAND offers correct exact values for multiples of 30 degrees.
+     This implementation is also compatible with the behavior of some legacy
+     compilers.  Keep this consistent with gfc_conv_intrinsic_cotand.  */
   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+  mpfr_set (result->value.real, x->value.real, GFC_RND_MODE);
+  mpfr_add_ui (result->value.real, result->value.real, 90, GFC_RND_MODE);
+  simplify_tand (result->value.real);
+  mpfr_neg (result->value.real, result->value.real, GFC_RND_MODE);
 
-  switch (x->ts.type)
-    {
-      case BT_REAL:
-	mpfr_cos (result->value.real, x->value.real, GFC_RND_MODE);
-	break;
-
-      case BT_COMPLEX:
-	gfc_set_model_kind (x->ts.kind);
-	mpc_cos (result->value.complex, x->value.complex, GFC_MPC_RND_MODE);
-	break;
-
-      default:
-	gfc_internal_error ("in gfc_simplify_cos(): Bad type");
-    }
-
-  return range_check (result, "COS");
+  return range_check (result, "COTAND");
 }
 
 
@@ -7778,6 +7844,8 @@ gfc_simplify_sum (gfc_expr *array, gfc_expr *dim, gfc_expr *mask)
 }
 
 
+/* Simplify COTAN(X) where X has the unit of radian.  */
+
 gfc_expr *
 gfc_simplify_cotan (gfc_expr *x)
 {
@@ -7799,8 +7867,8 @@ gfc_simplify_cotan (gfc_expr *x)
       /* There is no builtin mpc_cot, so compute cot = cos / sin.  */
       val = &result->value.complex;
       mpc_init2 (swp, mpfr_get_default_prec ());
-      mpc_cos (swp, x->value.complex, GFC_MPC_RND_MODE);
-      mpc_sin (*val, x->value.complex, GFC_MPC_RND_MODE);
+      mpc_sin_cos (*val, swp, x->value.complex, GFC_MPC_RND_MODE,
+		   GFC_MPC_RND_MODE);
       mpc_div (*val, swp, *val, GFC_MPC_RND_MODE);
       mpc_clear (swp);
       break;
diff --git a/gcc/fortran/trans-intrinsic.c b/gcc/fortran/trans-intrinsic.c
index 00bec1ec1df..fd8809902b7 100644
--- a/gcc/fortran/trans-intrinsic.c
+++ b/gcc/fortran/trans-intrinsic.c
@@ -120,6 +120,9 @@ static GTY(()) gfc_intrinsic_map_t gfc_intrinsic_map[] =
 
   /* Functions in libgfortran.  */
   LIB_FUNCTION (ERFC_SCALED, "erfc_scaled", false),
+  LIB_FUNCTION (SIND, "sind", false),
+  LIB_FUNCTION (COSD, "cosd", false),
+  LIB_FUNCTION (TAND, "tand", false),
 
   /* End the list.  */
   LIB_FUNCTION (NONE, NULL, false)
@@ -4385,6 +4388,181 @@ gfc_conv_intrinsic_anyall (gfc_se * se, gfc_expr * expr, enum tree_code op)
   se->expr = resvar;
 }
 
+
+/* Generate the constant 180 / pi, which is used in the conversion
+   of acosd(), asind(), atand(), atan2d().  */
+
+static tree
+rad2deg (int kind)
+{
+  tree retval;
+  mpfr_t pi, t0;
+
+  gfc_set_model_kind (kind);
+  mpfr_init (pi);
+  mpfr_init (t0);
+  mpfr_set_si (t0, 180, GFC_RND_MODE);
+  mpfr_const_pi (pi, GFC_RND_MODE);
+  mpfr_div (t0, t0, pi, GFC_RND_MODE);
+  retval = gfc_conv_mpfr_to_tree (t0, kind, 0);
+  mpfr_clear (t0);
+  mpfr_clear (pi);
+  return retval;
+}
+
+
+/* ACOSD(x) is translated into ACOS(x) * 180 / pi.
+   ASIND(x) is translated into ASIN(x) * 180 / pi.
+   ATAND(x) is translated into ATAN(x) * 180 / pi.  */
+
+static void
+gfc_conv_intrinsic_atrigd (gfc_se * se, gfc_expr * expr, gfc_isym_id id)
+{
+  tree arg;
+  tree atrigd;
+  tree type;
+
+  type = gfc_typenode_for_spec (&expr->ts);
+
+  gfc_conv_intrinsic_function_args (se, expr, &arg, 1);
+
+  if (id == GFC_ISYM_ACOSD)
+    atrigd = gfc_builtin_decl_for_float_kind (BUILT_IN_ACOS, expr->ts.kind);
+  else if (id == GFC_ISYM_ASIND)
+    atrigd = gfc_builtin_decl_for_float_kind (BUILT_IN_ASIN, expr->ts.kind);
+  else if (id == GFC_ISYM_ATAND)
+    atrigd = gfc_builtin_decl_for_float_kind (BUILT_IN_ATAN, expr->ts.kind);
+  else
+    gcc_unreachable ();
+
+  atrigd = build_call_expr_loc (input_location, atrigd, 1, arg);
+
+  se->expr = fold_build2_loc (input_location, MULT_EXPR, type, atrigd,
+			      fold_convert (type, rad2deg (expr->ts.kind)));
+}
+
+
+/* COTAN(X) is translated into -TAN(X+PI/2) for REAL argument and
+   COS(X) / SIN(X) for COMPLEX argument.  */
+
+static void
+gfc_conv_intrinsic_cotan (gfc_se *se, gfc_expr *expr)
+{
+  gfc_intrinsic_map_t *m;
+  tree arg;
+  tree type;
+
+  type = gfc_typenode_for_spec (&expr->ts);
+  gfc_conv_intrinsic_function_args (se, expr, &arg, 1);
+
+  if (expr->ts.type == BT_REAL)
+    {
+      tree tan;
+      tree tmp;
+      mpfr_t pio2;
+
+      /* Create pi/2.  */
+      gfc_set_model_kind (expr->ts.kind);
+      mpfr_init (pio2);
+      mpfr_const_pi (pio2, GFC_RND_MODE);
+      mpfr_div_ui (pio2, pio2, 2, GFC_RND_MODE);
+      tmp = gfc_conv_mpfr_to_tree (pio2, expr->ts.kind, 0);
+      mpfr_clear (pio2);
+
+      /* Find tan builtin function.  */
+      m = gfc_intrinsic_map;
+      for (; m->id != GFC_ISYM_NONE || m->double_built_in != END_BUILTINS; m++)
+	if (GFC_ISYM_TAN == m->id)
+	  break;
+
+      tmp = fold_build2_loc (input_location, PLUS_EXPR, type, arg, tmp);
+      tan = gfc_get_intrinsic_lib_fndecl (m, expr);
+      tan = build_call_expr_loc (input_location, tan, 1, tmp);
+      se->expr = fold_build1_loc (input_location, NEGATE_EXPR, type, tan);
+    }
+  else
+    {
+      tree sin;
+      tree cos;
+
+      /* Find cos builtin function.  */
+      m = gfc_intrinsic_map;
+      for (; m->id != GFC_ISYM_NONE || m->double_built_in != END_BUILTINS; m++)
+	if (GFC_ISYM_COS == m->id)
+	  break;
+
+      cos = gfc_get_intrinsic_lib_fndecl (m, expr);
+      cos = build_call_expr_loc (input_location, cos, 1, arg);
+
+      /* Find sin builtin function.  */
+      m = gfc_intrinsic_map;
+      for (; m->id != GFC_ISYM_NONE || m->double_built_in != END_BUILTINS; m++)
+	if (GFC_ISYM_SIN == m->id)
+	  break;
+
+      sin = gfc_get_intrinsic_lib_fndecl (m, expr);
+      sin = build_call_expr_loc (input_location, sin, 1, arg);
+
+      /* Divide cos by sin. */
+      se->expr = fold_build2_loc (input_location, RDIV_EXPR, type, cos, sin);
+   }
+}
+
+
+/* COTAND(X) is translated into -TAND(X+90) for REAL argument.  */
+
+static void
+gfc_conv_intrinsic_cotand (gfc_se *se, gfc_expr *expr)
+{
+  tree arg;
+  tree type;
+  tree ninety_tree;
+  mpfr_t ninety;
+
+  type = gfc_typenode_for_spec (&expr->ts);
+  gfc_conv_intrinsic_function_args (se, expr, &arg, 1);
+
+  gfc_set_model_kind (expr->ts.kind);
+
+  /* Build the tree for x + 90.  */
+  mpfr_init_set_ui (ninety, 90, GFC_RND_MODE);
+  ninety_tree = gfc_conv_mpfr_to_tree (ninety, expr->ts.kind, 0);
+  arg = fold_build2_loc (input_location, PLUS_EXPR, type, arg, ninety_tree);
+  mpfr_clear (ninety);
+
+  /* Find tand.  */
+  gfc_intrinsic_map_t *m = gfc_intrinsic_map;
+  for (; m->id != GFC_ISYM_NONE || m->double_built_in != END_BUILTINS; m++)
+    if (GFC_ISYM_TAND == m->id)
+      break;
+
+  tree tand = gfc_get_intrinsic_lib_fndecl (m, expr);
+  tand = build_call_expr_loc (input_location, tand, 1, arg);
+
+  se->expr = fold_build1_loc (input_location, NEGATE_EXPR, type, tand);
+}
+
+
+/* ATAN2D(Y,X) is translated into ATAN2(Y,X) * 180 / PI. */
+
+static void
+gfc_conv_intrinsic_atan2d (gfc_se *se, gfc_expr *expr)
+{
+  tree args[2];
+  tree atan2d;
+  tree type;
+
+  gfc_conv_intrinsic_function_args (se, expr, args, 2);
+  type = TREE_TYPE (args[0]);
+
+  atan2d = gfc_builtin_decl_for_float_kind (BUILT_IN_ATAN2, expr->ts.kind);
+  atan2d = build_call_expr_loc (input_location, atan2d, 2, args[0], args[1]);
+
+  se->expr = fold_build2_loc (input_location, MULT_EXPR, type, atan2d,
+			      rad2deg (expr->ts.kind));
+}
+
+
 /* COUNT(A) = Number of true elements in A.  */
 static void
 gfc_conv_intrinsic_count (gfc_se * se, gfc_expr * expr)
@@ -9895,6 +10073,24 @@ gfc_conv_intrinsic_function (gfc_se * se, gfc_expr * expr)
       gfc_conv_intrinsic_anyall (se, expr, NE_EXPR);
       break;
 
+    case GFC_ISYM_ACOSD:
+    case GFC_ISYM_ASIND:
+    case GFC_ISYM_ATAND:
+      gfc_conv_intrinsic_atrigd (se, expr, expr->value.function.isym->id);
+      break;
+
+    case GFC_ISYM_COTAN:
+      gfc_conv_intrinsic_cotan (se, expr);
+      break;
+
+    case GFC_ISYM_COTAND:
+      gfc_conv_intrinsic_cotand (se, expr);
+      break;
+
+    case GFC_ISYM_ATAN2D:
+      gfc_conv_intrinsic_atan2d (se, expr);
+      break;
+
     case GFC_ISYM_BTEST:
       gfc_conv_intrinsic_btest (se, expr);
       break;
diff --git a/gcc/fortran/trigd_fe.inc b/gcc/fortran/trigd_fe.inc
new file mode 100644
index 00000000000..78ca4416a21
--- /dev/null
+++ b/gcc/fortran/trigd_fe.inc
@@ -0,0 +1,50 @@
+
+
+/* Stub for defining degree-valued trigonemetric functions using MPFR.
+   Copyright (C) 2000-2020 Free Software Foundation, Inc.
+   Contributed by Fritz Reese <fore...@gcc.gnu.org>
+   and Steven G. Kargl <ka...@gcc.gnu.org>
+
+This file is part of GCC.
+
+GCC 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, or (at your option) any later
+version.
+
+GCC 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.
+
+You should have received a copy of the GNU General Public License
+along with GCC; see the file COPYING3.  If not see
+<http://www.gnu.org/licenses/>.  */
+
+#define FTYPE mpfr_t
+#define RETTYPE void
+#define RETURN(x) do { } while (0)
+#define ITYPE mpz_t
+
+#define ISFINITE(x) mpfr_number_p(x)
+#define D2R(x) deg2rad(x)
+
+#define SIND simplify_sind
+#define COSD simplify_cosd
+#define TAND simplify_tand
+
+#ifdef HAVE_GFC_REAL_16
+#define COSD30 8.66025403784438646763723170752936183e-01Q
+#else
+#define COSD30 8.66025403784438646763723170752936183e-01L
+#endif
+
+#define SET_COSD30(x) mpfr_set_ld((x), COSD30, GFC_RND_MODE)
+
+static RETTYPE SIND (FTYPE);
+static RETTYPE COSD (FTYPE);
+static RETTYPE TAND (FTYPE);
+
+#include "../../libgfortran/intrinsics/trigd.inc"
+
+/* vim: set ft=c: */
diff --git a/libgfortran/Makefile.am b/libgfortran/Makefile.am
index 295a2d43564..8ca0f6c290d 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 a6804c1d4cf..97a978aa80f 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 3601bc24414..ebf1a6ff40b 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 00000000000..81699069545
--- /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 <ka...@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 00000000000..98bfae7e839
--- /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 <ka...@gcc.gnu.org>
+   and Fritz Reese <fore...@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 00000000000..b6d4145b995
--- /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 <ka...@gcc.gnu.org>
+   and Fritz Reese <fore...@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: */

Reply via email to