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.

14 files changed, 1479 insertions, 641 deletions

diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog
index 435f93d..fdbb8da 100644
--- a/gcc/fortran/ChangeLog
+++ b/gcc/fortran/ChangeLog
@@ -1,3 +1,40 @@
+2020-04-07  Fritz Reese  <foreese@gcc.gnu.org>
+	    Steven G. Kargl  <kargl@gcc.gnu.org>
+
+	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.
+
 2020-04-06  Steven G. Kargl  <kargl@gcc.gnu.org>
 
 	PR fortran/93686
diff --git a/gcc/fortran/gfortran.h b/gcc/fortran/gfortran.h
index 88e4d92..70a6405 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 3012187..17f5efc 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 d045674..166ae79 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 a991c3a..7ecb659 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 66ed925..f63f63c 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 00bec1e..fd88099 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 0000000..78ca441
--- /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 <foreese@gcc.gnu.org>
+   and Steven G. Kargl <kargl@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/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog
index ab04cc0..096cfce 100644
--- a/gcc/testsuite/ChangeLog
+++ b/gcc/testsuite/ChangeLog
@@ -1,3 +1,12 @@
+2020-04-01  Fritz Reese  <foreese@gcc.gnu.org>
+
+	 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.
+
 2020-04-07  Andre Vieira  <andre.simoesdiasvieira@arm.com>
 
 	* g++.target/arm/mve.exp: New.
diff --git a/gcc/testsuite/gfortran.dg/dec_math.f90 b/gcc/testsuite/gfortran.dg/dec_math.f90
index 2a50f97..cc141ab 100644
--- a/gcc/testsuite/gfortran.dg/dec_math.f90
+++ b/gcc/testsuite/gfortran.dg/dec_math.f90
@@ -1,289 +1,700 @@
-! { dg-options "-fdec-math" }
+! { dg-options "-cpp -std=gnu" }
 ! { dg-do run }
 !
-! Test extra math intrinsics offered by -fdec-math.
+! Test extra math intrinsics formerly offered by -fdec-math,
+! now included with -std=gnu or -std=legacy.
 !
 
-  subroutine cmpf(f1, f2, tolerance, str)
+module dec_math
+
+  implicit none
+
+  real(4), parameter :: pi_f = 3.14159274_4
+  real(8), parameter :: pi_d = 3.1415926535897931_8
+#ifdef __GFC_REAL_10__
+  real(10), parameter :: pi_l = 3.1415926535897932383_10
+#endif
+#ifdef __GFC_REAL_16__
+  real(16), parameter :: pi_q = 3.1415926535897932384626433832795028_16
+#endif
+
+  real(4), parameter :: r2d_f = 180.0_4 / pi_f
+  real(8), parameter :: r2d_d = 180.0_8 / pi_d
+#ifdef __GFC_REAL_10__
+  real(10), parameter :: r2d_l = 180.0_10 / pi_l
+#endif
+#ifdef __GFC_REAL_16__
+  real(16), parameter :: r2d_q = 180.0_16 / pi_q
+#endif
+
+contains
+
+  function d2rf(x)
+    implicit none
+    real(4), intent(in) :: x
+    real(4) :: d2rf
+    d2rf = (x * pi_f) / 180.0_4
+  endfunction
+
+  subroutine cmpf(x, f1, f2, tolerance, str)
     implicit none
-    real(4), intent(in) :: f1, f2, tolerance
+    real(4), intent(in) :: x, f1, f2, tolerance
     character(len=*), intent(in) :: str
     if ( abs(f2 - f1) .gt. tolerance ) then
-      write (*, '(A,F12.6,F12.6)') str, f1, f2
+      write (*, '(A,A,F12.6,A,F12.6,F12.6)') str, "(", x, ")", f1, f2
       STOP 1
     endif
   endsubroutine
 
-  subroutine cmpd(d1, d2, tolerance, str)
+  function d2rd(x)
     implicit none
-    real(8), intent(in) :: d1, d2, tolerance
+    real(8), intent(in) :: x
+    real(8) :: d2rd
+    d2rd = (x * pi_d) / 180.0_8
+  endfunction
+
+  subroutine cmpd(x, d1, d2, tolerance, str)
+    implicit none
+    real(8), intent(in) :: x, d1, d2, tolerance
     character(len=*), intent(in) :: str
     if ( dabs(d2 - d1) .gt. tolerance ) then
-      write (*, '(A,F12.6,F12.6)') str, d1, d2
+      write (*, '(A,A,F18.14,A,F18.14,F18.14)') str, "(", x, ")", d1, d2
       STOP 2
     endif
   endsubroutine
 
-implicit none
+#ifdef __GFC_REAL_10__
+  function d2rl(x)
+    implicit none
+    real(10), intent(in) :: x
+    real(10) :: d2rl
+    d2rl = (x * pi_l) / 180.0_10
+  endfunction
 
-  real(4), parameter :: pi_f = (4.0_4 *  atan(1.0_4))
-  real(8), parameter :: pi_d = (4.0_8 * datan(1.0_8))
-  real(4), parameter :: r2d_f = 180.0_4 / pi_f
-  real(8), parameter :: r2d_d = 180.0_8 / pi_d
-  real(4), parameter :: d2r_f = pi_f / 180.0_4
-  real(8), parameter :: d2r_d = pi_d / 180.0_8
+  subroutine cmpl(x, f1, f2, tolerance, str)
+    implicit none
+    real(10), intent(in) :: x, f1, f2, tolerance
+    character(len=*), intent(in) :: str
+    if ( abs(f2 - f1) .gt. tolerance ) then
+      write (*, '(A,A,F21.17,A,F21.17,F21.17)') str, "(", x, ")", f1, f2
+      STOP 1
+    endif
+  endsubroutine
+#endif
+
+#ifdef __GFC_REAL_16__
+  function d2rq(x)
+    implicit none
+    real(16), intent(in) :: x
+    real(16) :: d2rq
+    d2rq = (x * pi_q) / 180.0_16
+  endfunction
+
+  subroutine cmpq(x, f1, f2, tolerance, str)
+    implicit none
+    real(16), intent(in) :: x, f1, f2, tolerance
+    character(len=*), intent(in) :: str
+    if ( abs(f2 - f1) .gt. tolerance ) then
+      write (*, '(A,A,F34.30,A,F34.30,F34.30)') str, "(", x, ")", f1, f2
+      STOP 1
+    endif
+  endsubroutine
+#endif
+
+end module
+
+use dec_math
+
+implicit none
 
 ! inputs
 real(4) :: f_i1, f_i2
 real(4), volatile :: xf
 real(8) :: d_i1, d_i2
 real(8), volatile :: xd
+#ifdef __GFC_REAL_10__
+real(10) :: l_i1, l_i2
+real(10), volatile :: xl
+#endif
+#ifdef __GFC_REAL_16__
+real(16) :: q_i1, q_i2
+real(16), volatile :: xq
+#endif
 
 ! expected outputs from (oe) default (oxe) expression
 real(4) :: f_oe, f_oxe
 real(8) :: d_oe, d_oxe
+#ifdef __GFC_REAL_10__
+real(10) :: l_oe, l_oxe
+#endif
+#ifdef __GFC_REAL_16__
+real(16) :: q_oe, q_oxe
+#endif
 
 ! actual outputs from (oa) default (oc) constant (ox) expression
 real(4) :: f_oa, f_oc, f_ox
 real(8) :: d_oa, d_oc, d_ox
+#ifdef __GFC_REAL_10__
+real(10) :: l_oa, l_oc, l_ox
+#endif
+#ifdef __GFC_REAL_16__
+real(16) :: q_oa, q_oc, q_ox
+#endif
 
 ! tolerance of the answer: assert |exp-act| <= tol
-real(4) :: f_tol
-real(8) :: d_tol
-
-! equivalence tolerance
-f_tol = 5e-5_4
-d_tol = 5e-6_8
-
-! multiplication factors to test non-constant expressions
+! accept loss of ~four decimal places
+real(4), parameter :: f_tol  = 5e-3_4
+real(8), parameter :: d_tol  = 5e-10_8
+#ifdef __GFC_REAL_10__
+real(10), parameter :: l_tol = 5e-15_10
+#endif
+#ifdef __GFC_REAL_16__
+real(16), parameter :: q_tol = 5e-20_16
+#endif
+
+! volatile multiplication factors to test non-constant expressions
 xf = 2.0_4
 xd = 2.0_8
-
-! Input
-f_i1 = 0.68032123_4
-d_i1 = 0.68032123_8
-
-! Expected
-f_oe =     r2d_f*acos (f_i1)
-f_oxe = xf*r2d_f*acos (f_i1)
-d_oe =     r2d_d*dacos(d_i1)
-d_oxe = xd*r2d_d*dacos(d_i1)
+#ifdef __GFC_REAL_10__
+xl = 2.0_10
+#endif
+#ifdef __GFC_REAL_16__
+xq = 2.0_16
+#endif
+
+! Input -- cos(pi/4)
+f_i1 = 0.707107_4
+d_i1 = 0.707106781186548_8
+#ifdef __GFC_REAL_10__
+l_i1 = 0.707106781186547573_10
+#endif
+#ifdef __GFC_REAL_16__
+q_i1 = 0.707106781186547572737310929369414_16
+#endif
+
+! Expected -- pi/4
+f_oe =  r2d_f * acos (f_i1)
+f_oxe = r2d_f * acos (xf * f_i1)
+d_oe =  r2d_d * acos (d_i1)
+d_oxe = r2d_d * acos (xd * d_i1)
+#ifdef __GFC_REAL_10__
+l_oe =  r2d_l * acos (l_i1)
+l_oxe = r2d_l * acos (xl * l_i1)
+#endif
+#ifdef __GFC_REAL_16__
+q_oe =  r2d_q * acos (q_i1)
+q_oxe = r2d_q * acos (xq * q_i1)
+#endif
 
 ! Actual
 f_oa =    acosd (f_i1)
-f_oc =    acosd (0.68032123_4)
-f_ox = xf*acosd (f_i1)
-d_oa =    dacosd (d_i1)
-d_oc =    dacosd (0.68032123_8)
-d_ox = xd*dacosd (0.68032123_8)
-
-call cmpf(f_oe,  f_oa, f_tol, "( ) acosd")
-call cmpf(f_oe,  f_oc, f_tol, "(c) acosd")
-call cmpf(f_oxe, f_ox, f_tol, "(x) acosd")
-call cmpd(d_oe,  d_oa, d_tol, "( ) dacosd")
-call cmpd(d_oe,  d_oc, d_tol, "(c) dacosd")
-call cmpd(d_oxe, d_ox, d_tol, "(x) dacosd")
+f_oc =    acosd (0.707107_4)
+f_ox = acosd (xf * f_i1)
+d_oa =     acosd (d_i1)
+d_oc =     acosd (0.707106781186548_8)
+d_ox = acosd (xd * 0.707106781186548_8)
+#ifdef __GFC_REAL_10__
+l_oa =    acosd (l_i1)
+l_oc =    acosd (0.707106781186547573_10)
+l_ox = acosd (xl * l_i1)
+#endif
+#ifdef __GFC_REAL_16__
+q_oa =    acosd (q_i1)
+q_oc =    acosd (0.707106781186547572737310929369414_16)
+q_ox = acosd (xq * 0.707106781186547572737310929369414_16)
+#endif
+
+call cmpf(f_i1, f_oe,  f_oa, f_tol, "( ) facosd")
+call cmpf(f_i1, f_oe,  f_oc, f_tol, "(c) facosd")
+call cmpf(f_i1, f_oxe, f_ox, f_tol, "(x) facosd")
+call cmpd(d_i1, d_oe,  d_oa, d_tol, "( ) dacosd")
+call cmpd(d_i1, d_oe,  d_oc, d_tol, "(c) dacosd")
+call cmpd(d_i1, d_oxe, d_ox, d_tol, "(x) dacosd")
+#ifdef __GFC_REAL_10__
+call cmpl(l_i1, l_oe,  l_oa, l_tol, "( ) lacosd")
+call cmpl(l_i1, l_oe,  l_oc, l_tol, "(c) lacosd")
+call cmpl(l_i1, l_oxe, l_ox, l_tol, "(x) lacosd")
+#endif
+#ifdef __GFC_REAL_16__
+call cmpq(q_i1, q_oe,  q_oa, q_tol, "( ) qacosd")
+call cmpq(q_i1, q_oe,  q_oc, q_tol, "(c) qacosd")
+call cmpq(q_i1, q_oxe, q_ox, q_tol, "(x) qacosd")
+#endif
 
 ! Input
 f_i1 = 60.0_4
 d_i1 = 60.0_8
+#ifdef __GFC_REAL_10__
+l_i1 = 60.0_10
+#endif
+#ifdef __GFC_REAL_16__
+q_i1 = 60.0_16
+#endif
 
 ! Expected
-f_oe  =    cos (d2r_f*f_i1)
-f_oxe = xf*cos (d2r_f*f_i1)
-d_oe  =    cos (d2r_d*d_i1)
-d_oxe = xd*cos (d2r_d*d_i1)
+f_oe  = cos (d2rf(f_i1))
+f_oxe = cos (d2rf(xf * f_i1))
+d_oe  = cos (d2rd(d_i1))
+d_oxe = cos (d2rd(xd * d_i1))
+#ifdef __GFC_REAL_10__
+l_oe  = cos (d2rl(l_i1))
+l_oxe = cos (d2rl(xl * l_i1))
+#endif
+#ifdef __GFC_REAL_16__
+q_oe  = cos (d2rq(q_i1))
+q_oxe = cos (d2rq(xq * q_i1))
+#endif
 
 ! Actual
-f_oa =     cosd (f_i1)
-f_oc =     cosd (60.0_4)
-f_ox = xf* cosd (f_i1)
-d_oa =    dcosd (d_i1)
-d_oc =    dcosd (60.0_8)
-d_ox = xd* cosd (d_i1)
-
-call cmpf(f_oe,  f_oa, f_tol, "( ) cosd")
-call cmpf(f_oe,  f_oc, f_tol, "(c) cosd")
-call cmpf(f_oxe, f_ox, f_tol, "(x) cosd")
-call cmpd(d_oe,  d_oa, d_tol, "( ) dcosd")
-call cmpd(d_oe,  d_oc, d_tol, "(c) dcosd")
-call cmpd(d_oxe, d_ox, d_tol, "(x) cosd")
-
-! Input
-f_i1 = 0.79345021_4
-d_i1 = 0.79345021_8
-
-! Expected
-f_oe  =    r2d_f*asin (f_i1)
-f_oxe = xf*r2d_f*asin (f_i1)
-d_oe  =    r2d_d*asin (d_i1)
-d_oxe = xd*r2d_d*asin (d_i1)
+f_oa = cosd (f_i1)
+f_oc = cosd (60.0_4)
+f_ox = cosd (xf * f_i1)
+d_oa = cosd (d_i1)
+d_oc = cosd (60.0_8)
+d_ox = cosd (xd * d_i1)
+#ifdef __GFC_REAL_10__
+l_oa = cosd (l_i1)
+l_oc = cosd (60.0_10)
+l_ox = cosd (xl * l_i1)
+#endif
+#ifdef __GFC_REAL_16__
+q_oa = cosd (q_i1)
+q_oc = cosd (60.0_16)
+q_ox = cosd (xq * q_i1)
+#endif
+
+call cmpf(f_i1, f_oe,  f_oa, f_tol, "( ) fcosd")
+call cmpf(f_i1, f_oe,  f_oc, f_tol, "(c) fcosd")
+call cmpf(f_i1, f_oxe, f_ox, f_tol, "(x) fcosd")
+call cmpd(d_i1, d_oe,  d_oa, d_tol, "( ) dcosd")
+call cmpd(d_i1, d_oe,  d_oc, d_tol, "(c) dcosd")
+call cmpd(d_i1, d_oxe, d_ox, d_tol, "(x) cosd")
+#ifdef __GFC_REAL_10__
+call cmpl(l_i1, l_oe,  l_oa, l_tol, "( ) lcosd")
+call cmpl(l_i1, l_oe,  l_oc, l_tol, "(c) lcosd")
+call cmpl(l_i1, l_oxe, l_ox, l_tol, "(x) lcosd")
+#endif
+#ifdef __GFC_REAL_16__
+call cmpq(q_i1, q_oe,  q_oa, q_tol, "( ) qcosd")
+call cmpq(q_i1, q_oe,  q_oc, q_tol, "(c) qcosd")
+call cmpq(q_i1, q_oxe, q_ox, q_tol, "(x) qcosd")
+#endif
+
+! Input -- sin(pi/4)
+f_i1 = 0.707107_4
+d_i1 = 0.707106781186548_8
+#ifdef __GFC_REAL_10__
+l_i1 = 0.707106781186547573_10
+#endif
+#ifdef __GFC_REAL_16__
+q_i1 = 0.707106781186547572737310929369414_16
+#endif
+
+! Expected -- pi/4
+f_oe  = r2d_f * asin (f_i1)
+f_oxe = r2d_f * asin (xf * f_i1)
+d_oe  = r2d_d * asin (d_i1)
+d_oxe = r2d_d * asin (xd * d_i1)
+#ifdef __GFC_REAL_10__
+l_oe  = r2d_l * asin (l_i1)
+l_oxe = r2d_l * asin (xl * l_i1)
+#endif
+#ifdef __GFC_REAL_16__
+q_oe  = r2d_q * asin (q_i1)
+q_oxe = r2d_q * asin (xq * q_i1)
+#endif
 
 ! Actual
-f_oa =     asind (f_i1)
-f_oc =     asind (0.79345021_4)
-f_ox = xf* asind (f_i1)
-d_oa =    dasind (d_i1)
-d_oc =    dasind (0.79345021_8)
-d_ox = xd* asind (d_i1)
-
-call cmpf(f_oe,  f_oa, f_tol, "( ) asind")
-call cmpf(f_oe,  f_oc, f_tol, "(c) asind")
-call cmpf(f_oxe, f_ox, f_tol, "(x) asind")
-call cmpd(d_oe,  d_oa, d_tol, "( ) dasind")
-call cmpd(d_oe,  d_oc, d_tol, "(c) dasind")
-call cmpd(d_oxe, d_ox, d_tol, "(x) asind")
+f_oa = asind (f_i1)
+f_oc = asind (0.707107_4)
+f_ox = asind (xf * f_i1)
+d_oa = asind (d_i1)
+d_oc = asind (0.707106781186548_8)
+d_ox = asind (xd * d_i1)
+#ifdef __GFC_REAL_10__
+l_oa = asind (l_i1)
+l_oc = asind (0.707106781186547573_10)
+l_ox = asind (xl * l_i1)
+#endif
+#ifdef __GFC_REAL_16__
+q_oa = asind (q_i1)
+q_oc = asind (0.707106781186547572737310929369414_16)
+q_ox = asind (xq * q_i1)
+#endif
+
+call cmpf(f_i1, f_oe,  f_oa, f_tol, "( ) fasind")
+call cmpf(f_i1, f_oe,  f_oc, f_tol, "(c) fasind")
+call cmpf(f_i1, f_oxe, f_ox, f_tol, "(x) fasind")
+call cmpd(d_i1, d_oe,  d_oa, d_tol, "( ) dasind")
+call cmpd(d_i1, d_oe,  d_oc, d_tol, "(c) dasind")
+call cmpd(d_i1, d_oxe, d_ox, d_tol, "(x) asind")
+#ifdef __GFC_REAL_10__
+call cmpl(l_i1, l_oe,  l_oa, l_tol, "( ) lasind")
+call cmpl(l_i1, l_oe,  l_oc, l_tol, "(c) lasind")
+call cmpl(l_i1, l_oxe, l_ox, l_tol, "(x) lasind")
+#endif
+#ifdef __GFC_REAL_16__
+call cmpq(q_i1, q_oe,  q_oa, q_tol, "( ) qasind")
+call cmpq(q_i1, q_oe,  q_oc, q_tol, "(c) qasind")
+call cmpq(q_i1, q_oxe, q_ox, q_tol, "(x) qasind")
+#endif
 
 ! Input
 f_i1 = 60.0_4
 d_i1 = 60.0_8
+#ifdef __GFC_REAL_10__
+l_i1 = 60.0_10
+#endif
+#ifdef __GFC_REAL_16__
+q_i1 = 60.0_16
+#endif
 
 ! Expected
-f_oe  =    sin (d2r_f*f_i1)
-f_oxe = xf*sin (d2r_f*f_i1)
-d_oe  =    sin (d2r_d*d_i1)
-d_oxe = xd*sin (d2r_d*d_i1)
+f_oe  = sin (d2rf(f_i1))
+f_oxe = sin (d2rf(xf * f_i1))
+d_oe  = sin (d2rd(d_i1))
+d_oxe = sin (d2rd(xd * d_i1))
+#ifdef __GFC_REAL_10__
+l_oe  = sin (d2rl(l_i1))
+l_oxe = sin (d2rl(xl * l_i1))
+#endif
+#ifdef __GFC_REAL_16__
+q_oe  = sin (d2rq(q_i1))
+q_oxe = sin (d2rq(xq * q_i1))
+#endif
 
 ! Actual
-f_oa =     sind (f_i1)
-f_oc =     sind (60.0_4)
-f_ox = xf* sind (f_i1)
-d_oa =    dsind (d_i1)
-d_oc =    dsind (60.0_8)
-d_ox = xd* sind (d_i1)
-
-call cmpf(f_oe,  f_oa, f_tol, "( ) sind")
-call cmpf(f_oe,  f_oc, f_tol, "(c) sind")
-call cmpf(f_oxe, f_ox, f_tol, "(x) sind")
-call cmpd(d_oe,  d_oa, d_tol, "( ) dsind")
-call cmpd(d_oe,  d_oc, d_tol, "(c) dsind")
-call cmpd(d_oxe, d_ox, d_tol, "(x) sind")
+f_oa = sind (f_i1)
+f_oc = sind (60.0_4)
+f_ox = sind (xf * f_i1)
+d_oa = sind (d_i1)
+d_oc = sind (60.0_8)
+d_ox = sind (xd * d_i1)
+#ifdef __GFC_REAL_10__
+l_oa = sind (l_i1)
+l_oc = sind (60.0_10)
+l_ox = sind (xl * l_i1)
+#endif
+#ifdef __GFC_REAL_16__
+q_oa = sind (q_i1)
+q_oc = sind (60.0_16)
+q_ox = sind (xq * q_i1)
+#endif
+
+call cmpf(f_i1, f_oe,  f_oa, f_tol, "( ) fsind")
+call cmpf(f_i1, f_oe,  f_oc, f_tol, "(c) fsind")
+call cmpf(f_i1, f_oxe, f_ox, f_tol, "(x) fsind")
+call cmpd(d_i1, d_oe,  d_oa, d_tol, "( ) dsind")
+call cmpd(d_i1, d_oe,  d_oc, d_tol, "(c) dsind")
+call cmpd(d_i1, d_oxe, d_ox, d_tol, "(x) sind")
+#ifdef __GFC_REAL_10__
+call cmpl(l_i1, l_oe,  l_oa, l_tol, "( ) lsind")
+call cmpl(l_i1, l_oe,  l_oc, l_tol, "(c) lsind")
+call cmpl(l_i1, l_oxe, l_ox, l_tol, "(x) lsind")
+#endif
+#ifdef __GFC_REAL_16__
+call cmpq(q_i1, q_oe,  q_oa, q_tol, "( ) qsind")
+call cmpq(q_i1, q_oe,  q_oc, q_tol, "(c) qsind")
+call cmpq(q_i1, q_oxe, q_ox, q_tol, "(x) qsind")
+#endif
 
 ! Input
-f_i1 = 2.679676_4
-f_i2 = 1.0_4
-d_i1 = 2.679676_8
-d_i2 = 1.0_8
+f_i1 = 1.0_4
+f_i2 = 2.0_4
+d_i1 = 1.0_8
+d_i2 = 2.0_8
+#ifdef __GFC_REAL_10__
+l_i1 = 1.0_10
+l_i2 = 2.0_10
+#endif
+#ifdef __GFC_REAL_16__
+q_i1 = 1.0_16
+q_i2 = 2.0_16
+#endif
 
 ! Expected
-f_oe  =    r2d_f*atan2 (f_i1, f_i2)
-f_oxe = xf*r2d_f*atan2 (f_i1, f_i2)
-d_oe  =    r2d_d*atan2 (d_i1, d_i2)
-d_oxe = xd*r2d_d*atan2 (d_i1, d_i2)
+f_oe  = r2d_f * atan2 (f_i1, f_i2)
+f_oxe = r2d_f * atan2 (xf * f_i1, f_i2)
+d_oe  = r2d_d * atan2 (d_i1, d_i2)
+d_oxe = r2d_d * atan2 (xd * d_i1, d_i2)
+#ifdef __GFC_REAL_10__
+l_oe  = r2d_l * atan2 (l_i1, l_i2)
+l_oxe = r2d_l * atan2 (xl * l_i1, l_i2)
+#endif
+#ifdef __GFC_REAL_16__
+q_oe  = r2d_q * atan2 (q_i1, q_i2)
+q_oxe = r2d_q * atan2 (xq * q_i1, q_i2)
+#endif
 
 ! Actual
-f_oa =     atan2d (f_i1, f_i2)
-f_oc =     atan2d (2.679676_4, 1.0_4)
-f_ox = xf* atan2d (f_i1, f_i2)
-d_oa =    datan2d (d_i1, d_i2)
-d_oc =    datan2d (2.679676_8, 1.0_8)
-d_ox = xd* atan2d (d_i1, d_i2)
-
-call cmpf(f_oe,  f_oa, f_tol, "( ) atan2d")
-call cmpf(f_oe,  f_oc, f_tol, "(c) atan2d")
-call cmpf(f_oxe, f_ox, f_tol, "(x) atan2d")
-call cmpd(d_oe,  d_oa, d_tol, "( ) datan2d")
-call cmpd(d_oe,  d_oc, d_tol, "(c) datan2d")
-call cmpd(d_oxe, d_ox, d_tol, "(x) atan2d")
+f_oa = atan2d (f_i1, f_i2)
+f_oc = atan2d (1.0_4, 2.0_4)
+f_ox = atan2d (xf * f_i1, f_i2)
+d_oa = atan2d (d_i1, d_i2)
+d_oc = atan2d (1.0_8, 2.0_8)
+d_ox = atan2d (xd * d_i1, d_i2)
+#ifdef __GFC_REAL_10__
+l_oa = atan2d (l_i1, l_i2)
+l_oc = atan2d (1.0_10, 2.0_10)
+l_ox = atan2d (xl * l_i1, l_i2)
+#endif
+#ifdef __GFC_REAL_16__
+q_oa = atan2d (q_i1, q_i2)
+q_oc = atan2d (1.0_16, 2.0_16)
+q_ox = atan2d (xq * q_i1, q_i2)
+#endif
+
+call cmpf(f_i1, f_oe,  f_oa, f_tol, "( ) fatan2d")
+call cmpf(f_i1, f_oe,  f_oc, f_tol, "(c) fatan2d")
+call cmpf(f_i1, f_oxe, f_ox, f_tol, "(x) fatan2d")
+call cmpd(d_i1, d_oe,  d_oa, d_tol, "( ) datan2d")
+call cmpd(d_i1, d_oe,  d_oc, d_tol, "(c) datan2d")
+call cmpd(d_i1, d_oxe, d_ox, d_tol, "(x) atan2d")
+#ifdef __GFC_REAL_10__
+call cmpl(l_i1, l_oe,  l_oa, l_tol, "( ) latan2d")
+call cmpl(l_i1, l_oe,  l_oc, l_tol, "(c) latan2d")
+call cmpl(l_i1, l_oxe, l_ox, l_tol, "(x) latan2d")
+#endif
+#ifdef __GFC_REAL_16__
+call cmpq(q_i1, q_oe,  q_oa, q_tol, "( ) qatan2d")
+call cmpq(q_i1, q_oe,  q_oc, q_tol, "(c) qatan2d")
+call cmpq(q_i1, q_oxe, q_ox, q_tol, "(x) qatan2d")
+#endif
 
 ! Input
-f_i1 = 1.5874993_4
-d_i1 = 1.5874993_8
+f_i1 = 1.55741_4
+d_i1 = 1.5574077246549_8
+#ifdef __GFC_REAL_10__
+l_i1 = 1.55740772465490229_10
+#endif
+#ifdef __GFC_REAL_16__
+q_i1 = 1.55740772465490229237161656783428_16
+#endif
 
 ! Expected
-f_oe  =    r2d_f*atan (f_i1)
-f_oxe = xf*r2d_f*atan (f_i1)
-d_oe  =    r2d_d*atan (d_i1)
-d_oxe = xd*r2d_d*atan (d_i1)
+f_oe  = r2d_f * atan (f_i1)
+f_oxe = r2d_f * atan (xf * f_i1)
+d_oe  = r2d_d * atan (d_i1)
+d_oxe = r2d_d * atan (xd * d_i1)
+#ifdef __GFC_REAL_10__
+l_oe  = r2d_l * atan (l_i1)
+l_oxe = r2d_l * atan (xl * l_i1)
+#endif
+#ifdef __GFC_REAL_16__
+q_oe  = r2d_q * atan (q_i1)
+q_oxe = r2d_q * atan (xq * q_i1)
+#endif
 
 ! Actual
-f_oa =     atand (f_i1)
-f_oc =     atand (1.5874993_4)
-f_ox = xf* atand (f_i1)
-d_oa =    datand (d_i1)
-d_oc =    datand (1.5874993_8)
-d_ox = xd* atand (d_i1)
-
-call cmpf(f_oe,  f_oa, f_tol, "( ) atand")
-call cmpf(f_oe,  f_oc, f_tol, "(c) atand")
-call cmpf(f_oxe, f_ox, f_tol, "(x) atand")
-call cmpd(d_oe,  d_oa, d_tol, "( ) datand")
-call cmpd(d_oe,  d_oc, d_tol, "(c) datand")
-call cmpd(d_oxe, d_ox, d_tol, "(x) atand")
+f_oa = atand (f_i1)
+f_oc = atand (1.55741_4)
+f_ox = atand (xf * f_i1)
+d_oa = atand (d_i1)
+d_oc = atand (1.5574077246549_8)
+d_ox = atand (xd * d_i1)
+#ifdef __GFC_REAL_10__
+l_oa = atand (l_i1)
+l_oc = atand (1.55740772465490229_10)
+l_ox = atand (xl * l_i1)
+#endif
+#ifdef __GFC_REAL_16__
+q_oa = atand (q_i1)
+q_oc = atand (1.55740772465490229237161656783428_16)
+q_ox = atand (xq * q_i1)
+#endif
+
+call cmpf(f_i1, f_oe,  f_oa, f_tol, "( ) fatand")
+call cmpf(f_i1, f_oe,  f_oc, f_tol, "(c) fatand")
+call cmpf(f_i1, f_oxe, f_ox, f_tol, "(x) fatand")
+call cmpd(d_i1, d_oe,  d_oa, d_tol, "( ) datand")
+call cmpd(d_i1, d_oe,  d_oc, d_tol, "(c) datand")
+call cmpd(d_i1, d_oxe, d_ox, d_tol, "(x) atand")
+#ifdef __GFC_REAL_10__
+call cmpl(l_i1, l_oe,  l_oa, l_tol, "( ) latand")
+call cmpl(l_i1, l_oe,  l_oc, l_tol, "(c) latand")
+call cmpl(l_i1, l_oxe, l_ox, l_tol, "(x) latand")
+#endif
+#ifdef __GFC_REAL_16__
+call cmpq(q_i1, q_oe,  q_oa, q_tol, "( ) qatand")
+call cmpq(q_i1, q_oe,  q_oc, q_tol, "(c) qatand")
+call cmpq(q_i1, q_oxe, q_ox, q_tol, "(x) qatand")
+#endif
 
 ! Input
-f_i1 = 0.6_4
-d_i1 = 0.6_8
+f_i1 = 34.3775_4
+d_i1 = 34.3774677078494_8
+#ifdef __GFC_REAL_10__
+l_i1 = 34.3774677078493909_10
+#endif
+#ifdef __GFC_REAL_16__
+q_i1 = 34.3774677078493908766176900826395_16
+#endif
 
 ! Expected
-f_oe  =    cotan (d2r_f*f_i1)
-f_oxe = xf*cotan (d2r_f*f_i1)
-d_oe  =    cotan (d2r_d*d_i1)
-d_oxe = xd*cotan (d2r_d*d_i1)
+f_oe  = 1.0_4/tan (f_i1)
+f_oxe = 1.0_4/tan (xf * f_i1)
+d_oe  = 1.0_8/tan (d_i1)
+d_oxe = 1.0_8/tan (xd * d_i1)
+#ifdef __GFC_REAL_10__
+l_oe  = 1.0_10/tan (l_i1)
+l_oxe = 1.0_10/tan (xl * l_i1)
+#endif
+#ifdef __GFC_REAL_16__
+q_oe  = 1.0_16/tan (q_i1)
+q_oxe = 1.0_16/tan (xq * q_i1)
+#endif
 
 ! Actual
-f_oa =     cotand (f_i1)
-f_oc =     cotand (0.6_4)
-f_ox = xf* cotand (f_i1)
-d_oa =    dcotand (d_i1)
-d_oc =    dcotand (0.6_8)
-d_ox = xd* cotand (d_i1)
-
-call cmpf(f_oe,  f_oa, f_tol, "( ) cotand")
-call cmpf(f_oe,  f_oc, f_tol, "(c) cotand")
-call cmpf(f_oxe, f_ox, f_tol, "(x) cotand")
-call cmpd(d_oe,  d_oa, d_tol, "( ) dcotand")
-call cmpd(d_oe,  d_oc, d_tol, "(c) dcotand")
-call cmpd(d_oxe, d_ox, d_tol, "(x) cotand")
+f_oa = cotan (f_i1)
+f_oc = cotan (34.3775_4)
+f_ox = cotan (xf * f_i1)
+d_oa = cotan (d_i1)
+d_oc = cotan (34.3774677078494_8)
+d_ox = cotan (xd * d_i1)
+#ifdef __GFC_REAL_10__
+l_oa = cotan (l_i1)
+l_oc = cotan (34.3774677078493909_10)
+l_ox = cotan (xl * l_i1)
+#endif
+#ifdef __GFC_REAL_16__
+q_oa = cotan (q_i1)
+q_oc = cotan (34.3774677078493908766176900826395_16)
+q_ox = cotan (xq * q_i1)
+#endif
+
+call cmpf(f_i1, f_oe,  f_oa, f_tol, "( ) fcotan")
+call cmpf(f_i1, f_oe,  f_oc, f_tol, "(c) fcotan")
+call cmpf(f_i1, f_oxe, f_ox, f_tol, "(x) fcotan")
+call cmpd(d_i1, d_oe,  d_oa, d_tol, "( ) dcotan")
+call cmpd(d_i1, d_oe,  d_oc, d_tol, "(c) dcotan")
+call cmpd(d_i1, d_oxe, d_ox, d_tol, "(x) cotan")
+#ifdef __GFC_REAL_10__
+call cmpl(l_i1, l_oe,  l_oa, l_tol, "( ) lcotan")
+call cmpl(l_i1, l_oe,  l_oc, l_tol, "(c) lcotan")
+call cmpl(l_i1, l_oxe, l_ox, l_tol, "(x) lcotan")
+#endif
+#ifdef __GFC_REAL_16__
+call cmpq(q_i1, q_oe,  q_oa, q_tol, "( ) qcotan")
+call cmpq(q_i1, q_oe,  q_oc, q_tol, "(c) qcotan")
+call cmpq(q_i1, q_oxe, q_ox, q_tol, "(x) qcotan")
+#endif
 
 ! Input
 f_i1 = 0.6_4
 d_i1 = 0.6_8
+#ifdef __GFC_REAL_10__
+l_i1 = 0.6_10
+#endif
+#ifdef __GFC_REAL_16__
+q_i1 = 0.6_16
+#endif
 
 ! Expected
-f_oe  =     1.0_4/tan (f_i1)
-f_oxe = xf* 1.0_4/tan (f_i1)
-d_oe  =    1.0_8/dtan (d_i1)
-d_oxe = xd*1.0_8/dtan (d_i1)
+f_oe  = cotan (d2rf(f_i1))
+f_oxe = cotan (d2rf(xf * f_i1))
+d_oe  = cotan (d2rd(d_i1))
+d_oxe = cotan (d2rd(xd * d_i1))
+#ifdef __GFC_REAL_10__
+l_oe  = cotan (d2rl(l_i1))
+l_oxe = cotan (d2rl(xl * l_i1))
+#endif
+#ifdef __GFC_REAL_16__
+q_oe  = cotan (d2rq(q_i1))
+q_oxe = cotan (d2rq(xq * q_i1))
+#endif
 
 ! Actual
-f_oa =     cotan (f_i1)
-f_oc =     cotan (0.6_4)
-f_ox = xf* cotan (f_i1)
-d_oa =    dcotan (d_i1)
-d_oc =    dcotan (0.6_8)
-d_ox = xd* cotan (d_i1)
-
-call cmpf(f_oe,  f_oa, f_tol, "( ) cotan")
-call cmpf(f_oe,  f_oc, f_tol, "(c) cotan")
-call cmpf(f_oxe, f_ox, f_tol, "(x) cotan")
-call cmpd(d_oe,  d_oa, d_tol, "( ) dcotan")
-call cmpd(d_oe,  d_oc, d_tol, "(c) dcotan")
-call cmpd(d_oxe, d_ox, d_tol, "(x) cotan")
+f_oa = cotand (f_i1)
+f_oc = cotand (0.6_4)
+f_ox = cotand (xf * f_i1)
+d_oa = cotand (d_i1)
+d_oc = cotand (0.6_8)
+d_ox = cotand (xd * d_i1)
+#ifdef __GFC_REAL_10__
+l_oa = cotand (l_i1)
+l_oc = cotand (0.6_10)
+l_ox = cotand (xl * l_i1)
+#endif
+#ifdef __GFC_REAL_16__
+q_oa = cotand (q_i1)
+q_oc = cotand (0.6_16)
+q_ox = cotand (xq * q_i1)
+#endif
+
+call cmpf(f_i1, f_oe,  f_oa, f_tol, "( ) fcotand")
+call cmpf(f_i1, f_oe,  f_oc, f_tol, "(c) fcotand")
+call cmpf(f_i1, f_oxe, f_ox, f_tol, "(x) fcotand")
+call cmpd(d_i1, d_oe,  d_oa, d_tol, "( ) dcotand")
+call cmpd(d_i1, d_oe,  d_oc, d_tol, "(c) dcotand")
+call cmpd(d_i1, d_oxe, d_ox, d_tol, "(x) cotand")
+#ifdef __GFC_REAL_10__
+call cmpl(l_i1, l_oe,  l_oa, l_tol, "( ) lcotand")
+call cmpl(l_i1, l_oe,  l_oc, l_tol, "(c) lcotand")
+call cmpl(l_i1, l_oxe, l_ox, l_tol, "(x) lcotand")
+#endif
+#ifdef __GFC_REAL_16__
+call cmpq(q_i1, q_oe,  q_oa, q_tol, "( ) qcotand")
+call cmpq(q_i1, q_oe,  q_oc, q_tol, "(c) qcotand")
+call cmpq(q_i1, q_oxe, q_ox, q_tol, "(x) qcotand")
+#endif
 
 ! Input
 f_i1 = 60.0_4
 d_i1 = 60.0_8
+#ifdef __GFC_REAL_10__
+l_i1 = 60.0_10
+#endif
+#ifdef __GFC_REAL_16__
+q_i1 = 60.0_16
+#endif
 
 ! Expected
-f_oe  =    tan (d2r_f*f_i1)
-f_oxe = xf*tan (d2r_f*f_i1)
-d_oe  =    tan (d2r_d*d_i1)
-d_oxe = xd*tan (d2r_d*d_i1)
+f_oe  = tan (d2rf(f_i1))
+f_oxe = tan (d2rf(xf * f_i1))
+d_oe  = tan (d2rd(d_i1))
+d_oxe = tan (d2rd(xd * d_i1))
+#ifdef __GFC_REAL_10__
+l_oe  = tan (d2rl(l_i1))
+l_oxe = tan (d2rl(xl * l_i1))
+#endif
+#ifdef __GFC_REAL_16__
+q_oe  = tan (d2rq(q_i1))
+q_oxe = tan (d2rq(xq * q_i1))
+#endif
 
 ! Actual
-f_oa =     tand (f_i1)
-f_oc =     tand (60.0_4)
-f_ox = xf* tand (f_i1)
-d_oa =    dtand (d_i1)
-d_oc =    dtand (60.0_8)
-d_ox = xd* tand (d_i1)
-
-call cmpf(f_oe,  f_oa, f_tol, "( ) tand")
-call cmpf(f_oe,  f_oc, f_tol, "(c) tand")
-call cmpf(f_oxe, f_ox, f_tol, "(x) tand")
-call cmpd(d_oe,  d_oa, d_tol, "( ) dtand")
-call cmpd(d_oe,  d_oc, d_tol, "(c) dtand")
-call cmpd(d_oxe, d_ox, d_tol, "(x) tand")
+f_oa = tand (f_i1)
+f_oc = tand (60.0_4)
+f_ox = tand (xf * f_i1)
+d_oa = tand (d_i1)
+d_oc = tand (60.0_8)
+d_ox = tand (xd * d_i1)
+#ifdef __GFC_REAL_10__
+l_oa = tand (l_i1)
+l_oc = tand (60.0_10)
+l_ox = tand (xl * l_i1)
+#endif
+#ifdef __GFC_REAL_16__
+q_oa = tand (q_i1)
+q_oc = tand (60.0_16)
+q_ox = tand (xq * q_i1)
+#endif
+
+call cmpf(f_i1, f_oe,  f_oa, f_tol, "( ) ftand")
+call cmpf(f_i1, f_oe,  f_oc, f_tol, "(c) ftand")
+call cmpf(f_i1, f_oxe, f_ox, f_tol, "(x) ftand")
+call cmpd(d_i1, d_oe,  d_oa, d_tol, "( ) dtand")
+call cmpd(d_i1, d_oe,  d_oc, d_tol, "(c) dtand")
+call cmpd(d_i1, d_oxe, d_ox, d_tol, "(x) dtand")
+#ifdef __GFC_REAL_10__
+call cmpl(l_i1, l_oe,  l_oa, l_tol, "( ) ltand")
+call cmpl(l_i1, l_oe,  l_oc, l_tol, "(c) ltand")
+call cmpl(l_i1, l_oxe, l_ox, l_tol, "(x) ltand")
+#endif
+#ifdef __GFC_REAL_16__
+call cmpq(q_i1, q_oe,  q_oa, q_tol, "( ) qtand")
+call cmpq(q_i1, q_oe,  q_oc, q_tol, "(c) qtand")
+call cmpq(q_i1, q_oxe, q_ox, q_tol, "(x) qtand")
+#endif
 
 end
diff --git a/gcc/testsuite/gfortran.dg/dec_math_2.f90 b/gcc/testsuite/gfortran.dg/dec_math_2.f90
new file mode 100644
index 0000000..ac49e27
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/dec_math_2.f90
@@ -0,0 +1,14 @@
+! { dg-options "-fdec-math" }
+! { dg-do compile }
+!
+! Ensure extra math intrinsics formerly offered by -fdec-math
+! are still available with -fdec-math.
+!
+
+print *, sind(0.0)
+print *, cosd(0.0)
+print *, tand(0.0)
+print *, cotan(1.0)
+print *, cotand(90.0)
+
+end
diff --git a/gcc/testsuite/gfortran.dg/dec_math_3.f90 b/gcc/testsuite/gfortran.dg/dec_math_3.f90
new file mode 100644
index 0000000..5bf4398
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/dec_math_3.f90
@@ -0,0 +1,8 @@
+! { dg-options "-std=gnu" }
+! { dg-do compile }
+
+! Former ICE when simplifying asind, plus wrong function name in error message
+real, parameter :: d = asind(1.1) ! { dg-error "Argument of ASIND at.*must be between -1 and 1" }
+print *, d
+
+end
diff --git a/gcc/testsuite/gfortran.dg/dec_math_4.f90 b/gcc/testsuite/gfortran.dg/dec_math_4.f90
new file mode 100644
index 0000000..f83210a
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/dec_math_4.f90
@@ -0,0 +1,8 @@
+! { dg-options "-std=gnu" }
+! { dg-do compile }
+
+! Former ICE when simplifying complex cotan
+complex, parameter :: z = cotan((1., 1.))
+print *, z
+
+end
diff --git a/gcc/testsuite/gfortran.dg/dec_math_5.f90 b/gcc/testsuite/gfortran.dg/dec_math_5.f90
new file mode 100644
index 0000000..d761e03
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/dec_math_5.f90
@@ -0,0 +1,228 @@
+! { dg-options "-cpp -std=gnu" }
+! { dg-do run }
+!
+! Test values for degree-valued trigonometric intrinsics.
+!
+
+module dec_math_5
+
+
+  ! Use the highest precision available.
+  ! Note however that if both __GFC_REAL_10__ and __GFC_REAL_16__ are defined,
+  ! the size of real(16) is actually that of REAL(10) (80 bits) in which case
+  ! we should not over-estimate the precision available, or the test will fail.
+#if defined(__GFC_REAL_10__)
+  integer, parameter :: real_kind = 10
+  real(real_kind), parameter :: eps = 5e-11_10
+
+  real(real_kind), parameter :: pi_2 = 1.57079632679489656_10
+  real(real_kind), parameter :: pi = 3.14159265358979312_10
+  real(real_kind), parameter :: tau = 6.28318530717958623_10
+
+#elif defined(__GFC_REAL_16__)
+  integer, parameter :: real_kind = 16
+  real(real_kind), parameter :: eps = 5e-16_16
+
+  real(real_kind), parameter :: pi_2 = 1.5707963267948966192313216916397514_16
+  real(real_kind), parameter :: pi = 3.1415926535897932384626433832795_16
+  real(real_kind), parameter :: tau = 6.28318530717958647692528676655900559_16
+
+#else
+  integer, parameter :: real_kind = 8
+  real(real_kind), parameter :: eps = 5e-10_8
+
+  real(real_kind), parameter :: pi_2 = 1.57079632679490_8
+  real(real_kind), parameter :: pi = 3.14159265358979_8
+  real(real_kind), parameter :: tau = 6.28318530717959_8
+
+#endif
+
+  ! Important angles in canonical form.
+
+  integer, parameter :: nangle = 16
+
+  real(real_kind), dimension(nangle), parameter :: degrees = (/ &
+      0, & !  180 * 0
+     30, & ! 180 * 1/6
+     45, & ! 180 * 1/4
+     60, & ! 180 * 1/3
+     90, & ! 180 * 1/2
+    120, & ! 180 * 2/3
+    135, & ! 180 * 3/4
+    150, & ! 180 * 5/6
+    180, & ! 180
+    210, & ! 180 * 7/6
+    225, & ! 180 * 5/4
+    240, & ! 180 * 4/3
+    270, & ! 180 * 3/2
+    300, & ! 180 * 5/3
+    315, & ! 180 * 7/4
+    330  & ! 180 * 11/6
+  /)
+
+  real(real_kind), dimension(nangle), parameter :: radians = (/ &
+#ifdef __GFC_REAL_10__
+    0.000000000000000000_10, & ! pi * 0
+    0.523598775598298873_10, & ! pi * 1/6
+    0.785398163397448310_10, & ! pi * 1/4
+    1.047197551196597750_10, & ! pi * 1/3
+    1.570796326794896620_10, & ! pi * 1/2
+    2.094395102393195490_10, & ! pi * 2/3
+    2.356194490192344930_10, & ! pi * 3/4
+    2.617993877991494370_10, & ! pi * 5/6
+    3.141592653589793240_10, & ! pi
+    3.665191429188092110_10, & ! pi * 7/6
+    3.926990816987241550_10, & ! pi * 5/4
+    4.188790204786390980_10, & ! pi * 4/3
+    4.712388980384689860_10, & ! pi * 3/2
+    5.235987755982988730_10, & ! pi * 5/3
+    5.497787143782138170_10, & ! pi * 7/4
+    5.759586531581287600_10  & ! pi * 11/6
+
+#elif defined(__GFC_REAL_16__)
+    0.000000000000000000000000000000000_16, & ! pi * 0
+    0.523598775598298873077107230546584_16, & ! pi * 1/6
+    0.785398163397448309615660845819876_16, & ! pi * 1/4
+    1.047197551196597746154214461093170_16, & ! pi * 1/3
+    1.570796326794896619231321691639750_16, & ! pi * 1/2
+    2.094395102393195492308428922186330_16, & ! pi * 2/3
+    2.356194490192344928846982537459630_16, & ! pi * 3/4
+    2.617993877991494365385536152732920_16, & ! pi * 5/6
+    3.141592653589793238462643383279500_16, & ! pi
+    3.665191429188092111539750613826090_16, & ! pi * 7/6
+    3.926990816987241548078304229099380_16, & ! pi * 5/4
+    4.188790204786390984616857844372670_16, & ! pi * 4/3
+    4.712388980384689857693965074919250_16, & ! pi * 3/2
+    5.235987755982988730771072305465840_16, & ! pi * 5/3
+    5.497787143782138167309625920739130_16, & ! pi * 7/4
+    5.759586531581287603848179536012420_16  & ! pi * 11/6
+
+#else
+    0.000000000000000_8, & ! pi * 0
+    0.523598775598299_8, & ! pi * 1/6
+    0.785398163397448_8, & ! pi * 1/4
+    1.047197551196600_8, & ! pi * 1/3
+    1.570796326794900_8, & ! pi * 1/2
+    2.094395102393200_8, & ! pi * 2/3
+    2.356194490192340_8, & ! pi * 3/4
+    2.617993877991490_8, & ! pi * 5/6
+    3.141592653589790_8, & ! pi
+    3.665191429188090_8, & ! pi * 7/6
+    3.926990816987240_8, & ! pi * 5/4
+    4.188790204786390_8, & ! pi * 4/3
+    4.712388980384690_8, & ! pi * 3/2
+    5.235987755982990_8, & ! pi * 5/3
+    5.497787143782140_8, & ! pi * 7/4
+    5.759586531581290_8  & ! pi * 11/6
+#endif
+  /)
+
+  ! sind, cosd, tand, cotand
+
+  ! Ensure precision degrades minimally for large values.
+  integer, parameter :: nphase = 5
+
+  integer, dimension(nphase), parameter :: phases = (/ &
+    0, 1, 5, 100, 10000  &
+  /)
+
+contains
+
+  subroutine compare(strl, xl_in, xl_out, strr, xr_in, xr_out, eps)
+    use ieee_arithmetic
+    implicit none
+    character(*), intent(in) :: strl, strr
+    real(real_kind), intent(in) :: xl_in, xl_out, xr_in, xr_out, eps
+
+    if ((ieee_is_nan(xl_out) .neqv. ieee_is_nan(xr_out)) &
+        .or. (ieee_is_finite(xl_out) .neqv. ieee_is_finite(xr_out)) &
+        .or. (abs(xl_out - xr_out) .gt. eps)) then
+      write (*, 100) strl, "(", xl_in, "): ", xl_out
+      write (*, 100) strr, "(", xr_in, "): ", xr_out
+
+      if ((ieee_is_nan(xl_out) .eqv. ieee_is_nan(xr_out)) &
+          .and. ieee_is_finite(xl_out) .and. ieee_is_finite(xr_out)) then
+        write (*, 300) "|xl - xr| = ", abs(xl_out - xr_out)
+        write (*, 300) "    > eps = ", eps
+      endif
+
+      call abort()
+    endif
+
+#ifdef __GFC_REAL_16__
+    100  format((A8,A,F34.30,A,F34.30,F34.30))
+    200  format((A12,F34.30))
+    !500  format((A8,A,G34.29,A,G34.29,G34.29))
+#elif defined(__GFC_REAL_10__)
+    100  format((A8,A,F21.17,A,F21.17,F21.17))
+    200  format((A12,F21.17))
+    !500  format((A8,A,G21.16,A,G21.16,G21.16))
+#else
+    100  format((A8,A,F18.14,A,F18.14,F18.14))
+    200  format((A12,F18.14))
+    !500  format((A8,A,G18.13,A,G18.13,G18.13))
+#endif
+    300  format((A12,G8.2))
+  endsubroutine
+
+endmodule
+
+use dec_math_5
+use ieee_arithmetic
+implicit none
+
+integer :: phase_index, angle_index, phase
+real(real_kind) :: deg_in, deg_out, deg_out2, rad_in, rad_out
+
+! Try every value in degrees, and make sure they are correct compared to the
+! corresponding radian function.
+
+do phase_index = 1, size(phases)
+  phase = phases(phase_index)
+
+  do angle_index = 1, size(degrees)
+    ! eqv to degrees(angle_index) modulo 360
+    deg_in = degrees(angle_index) + phase * 360
+    rad_in = radians(angle_index) + phase * tau
+
+    ! sind vs. sin
+    deg_out = sind(deg_in)
+    rad_out = sin(rad_in)
+    call compare("sind", deg_in, deg_out, "sin", rad_in, rad_out, eps)
+
+    ! cosd vs. cos
+    deg_out = cosd(deg_in)
+    rad_out = cos(rad_in)
+    call compare("cosd", deg_in, deg_out, "cos", rad_in, rad_out, eps)
+
+    ! tand vs. tan
+    deg_out = tand(deg_in)
+    rad_out = tan(rad_in)
+    if ( ieee_is_finite(deg_out) ) then
+      call compare("tand", deg_in, deg_out, "tan", rad_in, rad_out, eps)
+    endif
+
+    ! cotand vs. cotan
+    deg_out = cotand(deg_in)
+    rad_out = cotan(rad_in)
+
+    ! Skip comparing infinity, because cotan does not return infinity
+    if ( ieee_is_finite(deg_out) ) then
+      call compare("cotand", deg_in, deg_out, "cotan", rad_in, rad_out, eps)
+
+      ! cotand vs. tand
+      deg_out = cotand(deg_in)
+      deg_out2 = -tand(deg_in + 90)
+
+      call compare("cotand", deg_in, deg_out, "-tand+90", deg_in, deg_out2, eps)
+      deg_out2 = 1 / tand(deg_in)
+      call compare("cotand", deg_in, deg_out, "1/tand", deg_in, deg_out2, eps)
+    endif
+
+  enddo
+
+
+enddo
+
+
+end