New flag -fdec-math for COTAN and degree trig intrinsics.

2016-10-11  Fritz Reese  <fritzoreese@gmail.com>

New flag -fdec-math for COTAN and degree trig intrinsics.

	gcc/fortran/
	* lang.opt: New flag -fdec-math.
	* options.c (set_dec_flags): Enable with -fdec.
	* invoke.texi, gfortran.texi, intrinsic.texi: Update documentation.
	* intrinsics.c (add_functions, do_simplify): New intrinsics
	with -fdec-math.
	* gfortran.h (gfc_isym_id): New isym GFC_ISYM_COTAN.
	* gfortran.h (gfc_resolve_atan2d, gfc_resolve_cotan,
	gfc_resolve_trigd, gfc_resolve_atrigd): New prototypes.
	* iresolve.c (resolve_trig_call, get_degrees, get_radians,
	is_trig_resolved, gfc_resolve_cotan, gfc_resolve_trigd,
	gfc_resolve_atrigd, gfc_resolve_atan2d): New functions.
	* intrinsics.h (gfc_simplify_atan2d, gfc_simplify_atrigd,
	gfc_simplify_cotan, gfc_simplify_trigd): New prototypes.
	* simplify.c (simplify_trig_call, degrees_f, radians_f,
	gfc_simplify_cotan, gfc_simplify_trigd, gfc_simplify_atrigd,
	gfc_simplify_atan2d): New functions.

	gcc/testsuite/gfortran.dg/
	* dec_math.f90: New testsuite.

From-SVN: r240989

13 files changed, 1463 insertions, 5 deletions

diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog
index 0c54c6b..907a8ef 100644
--- a/gcc/fortran/ChangeLog
+++ b/gcc/fortran/ChangeLog
@@ -1,3 +1,22 @@
+2016-10-11  Fritz Reese  <fritzoreese@gmail.com>
+
+	* lang.opt: New flag -fdec-math.
+	* options.c (set_dec_flags): Enable with -fdec.
+	* invoke.texi, gfortran.texi, intrinsic.texi: Update documentation.
+	* intrinsics.c (add_functions, do_simplify): New intrinsics
+	with -fdec-math.
+	* gfortran.h (gfc_isym_id): New isym GFC_ISYM_COTAN.
+	* gfortran.h (gfc_resolve_atan2d, gfc_resolve_cotan,
+	gfc_resolve_trigd, gfc_resolve_atrigd): New prototypes.
+	* iresolve.c (resolve_trig_call, get_degrees, get_radians,
+	is_trig_resolved, gfc_resolve_cotan, gfc_resolve_trigd,
+	gfc_resolve_atrigd, gfc_resolve_atan2d): New functions.
+	* intrinsics.h (gfc_simplify_atan2d, gfc_simplify_atrigd,
+	gfc_simplify_cotan, gfc_simplify_trigd): New prototypes.
+	* simplify.c (simplify_trig_call, degrees_f, radians_f,
+	gfc_simplify_cotan, gfc_simplify_trigd, gfc_simplify_atrigd,
+	gfc_simplify_atan2d): New functions.
+
 2016-10-10  Thomas Koenig  <tkoenig@gcc.gnu.org>
 
 	PR fortran/77915
diff --git a/gcc/fortran/gfortran.h b/gcc/fortran/gfortran.h
index 2cac42b..33de0ffb 100644
--- a/gcc/fortran/gfortran.h
+++ b/gcc/fortran/gfortran.h
@@ -390,6 +390,7 @@ enum gfc_isym_id
   GFC_ISYM_CONVERSION,
   GFC_ISYM_COS,
   GFC_ISYM_COSH,
+  GFC_ISYM_COTAN,
   GFC_ISYM_COUNT,
   GFC_ISYM_CPU_TIME,
   GFC_ISYM_CSHIFT,
diff --git a/gcc/fortran/gfortran.texi b/gcc/fortran/gfortran.texi
index 797730c..301c286 100644
--- a/gcc/fortran/gfortran.texi
+++ b/gcc/fortran/gfortran.texi
@@ -1463,6 +1463,7 @@ without warning.
 * UNION and MAP::
 * Type variants for integer intrinsics::
 * AUTOMATIC and STATIC attributes::
+* Extended math intrinsics::
 @end menu
 
 @node Old-style kind specifications
@@ -2472,6 +2473,42 @@ subroutine f
 endsubroutine
 @end example
 
+@node Extended math intrinsics
+@subsection Extended math intrinsics
+@cindex intrinsics, math
+@cindex intrinsics, trigonometric functions
+
+GNU Fortran supports an extended list of mathematical intrinsics with the
+compile flag @option{-fdec-math} for compatability with legacy code.
+These intrinsics are described fully in @ref{Intrinsic Procedures} where it is
+noted that they are extensions and should be avoided whenever possible.
+
+Specifically, @option{-fdec-math} enables the @ref{COTAN} intrinsic, and
+trigonometric intrinsics which accept or produce values in degrees instead of
+radians.  Here is a summary of the new intrinsics:
+
+@multitable @columnfractions .5 .5
+@headitem Radians @tab Degrees
+@item @code{@ref{ACOS}}   @tab @code{@ref{ACOSD}}*
+@item @code{@ref{ASIN}}   @tab @code{@ref{ASIND}}*
+@item @code{@ref{ATAN}}   @tab @code{@ref{ATAND}}*
+@item @code{@ref{ATAN2}}  @tab @code{@ref{ATAN2D}}*
+@item @code{@ref{COS}}    @tab @code{@ref{COSD}}*
+@item @code{@ref{COTAN}}* @tab @code{@ref{COTAND}}*
+@item @code{@ref{SIN}}    @tab @code{@ref{SIND}}*
+@item @code{@ref{TAN}}    @tab @code{@ref{TAND}}*
+@end multitable
+
+* Enabled with @option{-fdec-math}.
+
+For advanced users, it may be important to know the implementation of these
+functions. They are simply wrappers around the standard radian functions, which
+have more accurate builtin versions. These functions convert their arguments
+(or results) to degrees (or radians) by taking the value modulus 360 (or 2*pi)
+and then multiplying it by a constant radian-to-degree (or degree-to-radian)
+factor, as appropriate. The factor is computed at compile-time as 180/pi (or
+pi/180).
+
 
 @node Extensions not implemented in GNU Fortran
 @section Extensions not implemented in GNU Fortran
diff --git a/gcc/fortran/intrinsic.c b/gcc/fortran/intrinsic.c
index cad54b8..fdc11d8 100644
--- a/gcc/fortran/intrinsic.c
+++ b/gcc/fortran/intrinsic.c
@@ -3139,6 +3139,117 @@ 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);
+
+      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 ("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);
+
+      make_generic ("asind", GFC_ISYM_ASIN, 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 ("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);
+
+      make_generic ("atand", GFC_ISYM_ATAN, 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_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);
+
+      make_generic ("atan2d", GFC_ISYM_ATAN2, 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_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);
+
+      make_generic ("cosd", GFC_ISYM_COS, 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 ("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);
+
+      make_generic ("cotan", GFC_ISYM_COTAN, 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 ("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);
+
+      make_generic ("cotand", GFC_ISYM_COTAN, GFC_STD_GNU);
+
+      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 ("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 ("sind", GFC_ISYM_SIN, GFC_STD_GNU);
+
+      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 ("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 ("tand", GFC_ISYM_TAN, GFC_STD_GNU);
+    }
+
   /* The following function is internally used for coarray libray functions.
      "make_from_module" makes it inaccessible for external users.  */
   add_sym_1 (GFC_PREFIX ("caf_get"), GFC_ISYM_CAF_GET, CLASS_IMPURE, ACTUAL_NO,
@@ -4227,6 +4338,15 @@ 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 f228976..8bba6e0 100644
--- a/gcc/fortran/intrinsic.h
+++ b/gcc/fortran/intrinsic.h
@@ -238,6 +238,7 @@ 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_dint (gfc_expr *);
 gfc_expr *gfc_simplify_anint (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_dnint (gfc_expr *);
@@ -248,6 +249,7 @@ gfc_expr *gfc_simplify_asinh (gfc_expr *);
 gfc_expr *gfc_simplify_atan (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 *);
 gfc_expr *gfc_simplify_bessel_j0 (gfc_expr *);
 gfc_expr *gfc_simplify_bessel_j1 (gfc_expr *);
 gfc_expr *gfc_simplify_bessel_jn (gfc_expr *, gfc_expr *);
@@ -271,6 +273,7 @@ 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_cosh (gfc_expr *);
+gfc_expr *gfc_simplify_cotan (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 *);
@@ -401,6 +404,7 @@ 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 *);
@@ -434,6 +438,7 @@ void gfc_resolve_asinh (gfc_expr *, gfc_expr *);
 void gfc_resolve_atan (gfc_expr *, gfc_expr *);
 void gfc_resolve_atanh (gfc_expr *, gfc_expr *);
 void gfc_resolve_atan2 (gfc_expr *, gfc_expr *, gfc_expr *);
+void gfc_resolve_atan2d (gfc_expr *, gfc_expr *, gfc_expr *);
 void gfc_resolve_atomic_def (gfc_code *);
 void gfc_resolve_atomic_ref (gfc_code *);
 void gfc_resolve_besn (gfc_expr *, gfc_expr *, gfc_expr *);
@@ -452,6 +457,7 @@ 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 *);
@@ -582,6 +588,8 @@ void gfc_resolve_time (gfc_expr *);
 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_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/intrinsic.texi b/gcc/fortran/intrinsic.texi
index 8cca9b1..16e1d5c 100644
--- a/gcc/fortran/intrinsic.texi
+++ b/gcc/fortran/intrinsic.texi
@@ -23,6 +23,9 @@ Some basic guidelines for editing this document:
 @end ignore
 
 @tex
+\gdef\acosd{\mathop{\rm acosd}\nolimits}
+\gdef\asind{\mathop{\rm asind}\nolimits}
+\gdef\atand{\mathop{\rm atand}\nolimits}
 \gdef\acos{\mathop{\rm acos}\nolimits}
 \gdef\asin{\mathop{\rm asin}\nolimits}
 \gdef\atan{\mathop{\rm atan}\nolimits}
@@ -43,6 +46,7 @@ Some basic guidelines for editing this document:
 * @code{ACCESS}:        ACCESS,    Checks file access modes
 * @code{ACHAR}:         ACHAR,     Character in @acronym{ASCII} collating sequence
 * @code{ACOS}:          ACOS,      Arccosine function
+* @code{ACOSD}:         ACOSD,     Arccosine function, degrees
 * @code{ACOSH}:         ACOSH,     Inverse hyperbolic cosine function
 * @code{ADJUSTL}:       ADJUSTL,   Left adjust a string
 * @code{ADJUSTR}:       ADJUSTR,   Right adjust a string
@@ -55,10 +59,13 @@ Some basic guidelines for editing this document:
 * @code{ANINT}:         ANINT,     Nearest whole number
 * @code{ANY}:           ANY,       Determine if any values are true
 * @code{ASIN}:          ASIN,      Arcsine function
+* @code{ASIND}:         ASIND,     Arcsine function, degrees
 * @code{ASINH}:         ASINH,     Inverse hyperbolic sine function
 * @code{ASSOCIATED}:    ASSOCIATED, Status of a pointer or pointer/target pair
 * @code{ATAN}:          ATAN,      Arctangent function
+* @code{ATAND}:         ATAND,     Arctangent function, degrees
 * @code{ATAN2}:         ATAN2,     Arctangent function
+* @code{ATAN2D}:        ATAN2D,    Arctangent function, degrees
 * @code{ATANH}:         ATANH,     Inverse hyperbolic tangent function
 * @code{ATOMIC_ADD}:    ATOMIC_ADD, Atomic ADD operation
 * @code{ATOMIC_AND}:    ATOMIC_AND, Atomic bitwise AND operation
@@ -106,7 +113,10 @@ Some basic guidelines for editing this document:
 * @code{COMPLEX}:       COMPLEX,   Complex conversion function
 * @code{CONJG}:         CONJG,     Complex conjugate function
 * @code{COS}:           COS,       Cosine function
+* @code{COSD}:          COSD,      Cosine function, degrees
 * @code{COSH}:          COSH,      Hyperbolic cosine function
+* @code{COTAN}:         COTAN,     Cotangent function
+* @code{COTAND}:        COTAND,    Cotangent function, degrees
 * @code{COUNT}:         COUNT,     Count occurrences of TRUE in an array
 * @code{CPU_TIME}:      CPU_TIME,  CPU time subroutine
 * @code{CSHIFT}:        CSHIFT,    Circular shift elements of an array
@@ -277,6 +287,7 @@ Some basic guidelines for editing this document:
 * @code{SIGN}:          SIGN,      Sign copying function
 * @code{SIGNAL}:        SIGNAL,    Signal handling subroutine (or function)
 * @code{SIN}:           SIN,       Sine function
+* @code{SIND}:          SIND,      Sine function, degrees
 * @code{SINH}:          SINH,      Hyperbolic sine function
 * @code{SIZE}:          SIZE,      Function to determine the size of an array
 * @code{SIZEOF}:        SIZEOF,    Determine the size in bytes of an expression
@@ -292,6 +303,7 @@ Some basic guidelines for editing this document:
 * @code{SYSTEM}:        SYSTEM,    Execute a shell command
 * @code{SYSTEM_CLOCK}:  SYSTEM_CLOCK, Time function
 * @code{TAN}:           TAN,       Tangent function
+* @code{TAND}:          TAND,      Tangent function, degrees
 * @code{TANH}:          TANH,      Hyperbolic tangent function
 * @code{THIS_IMAGE}:    THIS_IMAGE, Cosubscript index of this image
 * @code{TIME}:          TIME,      Time function
@@ -619,6 +631,65 @@ end program test_acos
 
 @item @emph{See also}:
 Inverse function: @ref{COS}
+Degrees function: @ref{ACOSD}
+
+@end table
+
+
+
+@node ACOSD
+@section @code{ACOSD} --- Arccosine function, degrees
+@fnindex ACOSD
+@fnindex DACOSD
+@cindex trigonometric function, cosine, inverse, degrees
+@cindex cosine, inverse, degrees
+
+@table @asis
+@item @emph{Description}:
+@code{ACOSD(X)} computes the arccosine of @var{X} in degrees (inverse of
+@code{COSD(X)}).
+
+This function is for compatibility only and should be avoided in favor of
+standard constructs wherever possible.
+
+@item @emph{Standard}:
+GNU Extension, enabled with @option{-fdec-math}
+
+@item @emph{Class}:
+Elemental function
+
+@item @emph{Syntax}:
+@code{RESULT = ACOSD(X)}
+
+@item @emph{Arguments}:
+@multitable @columnfractions .15 .70
+@item @var{X} @tab The type shall either be @code{REAL} with a magnitude that is
+less than or equal to one - or the type shall be @code{COMPLEX}.
+@end multitable
+
+@item @emph{Return value}:
+The return value is of the same type and kind as @var{X}.
+The real part of the result is in degrees and lies in the range
+@math{0 \leq \Re \acos(x) \leq 180}.
+
+@item @emph{Example}:
+@smallexample
+program test_acosd
+  real(8) :: x = 0.866_8
+  x = acosd(x)
+end program test_acosd
+@end smallexample
+
+@item @emph{Specific names}:
+@multitable @columnfractions .20 .20 .20 .25
+@item Name            @tab Argument         @tab Return type     @tab Standard
+@item @code{ACOSD(X)}  @tab @code{REAL(4) X} @tab @code{REAL(4)}  @tab GNU Extension
+@item @code{DACOSD(X)} @tab @code{REAL(8) X} @tab @code{REAL(8)}  @tab GNU Extension
+@end multitable
+
+@item @emph{See also}:
+Inverse function: @ref{COSD}
+Radians function: @ref{ACOS}
 
 @end table
 
@@ -1269,6 +1340,65 @@ end program test_asin
 
 @item @emph{See also}:
 Inverse function: @ref{SIN}
+Degrees function: @ref{ASIND}
+
+@end table
+
+
+
+@node ASIND
+@section @code{ASIND} --- Arcsine function, degrees
+@fnindex ASIND
+@fnindex DASIND
+@cindex trigonometric function, sine, inverse, degrees
+@cindex sine, inverse, degrees
+
+@table @asis
+@item @emph{Description}:
+@code{ASIND(X)} computes the arcsine of its @var{X} in degrees (inverse of
+@code{SIND(X)}).
+
+This function is for compatibility only and should be avoided in favor of
+standard constructs wherever possible.
+
+@item @emph{Standard}:
+GNU Extension, enabled with @option{-fdec-math}.
+
+@item @emph{Class}:
+Elemental function
+
+@item @emph{Syntax}:
+@code{RESULT = ASIND(X)}
+
+@item @emph{Arguments}:
+@multitable @columnfractions .15 .70
+@item @var{X} @tab The type shall be either @code{REAL} and a magnitude that is
+less than or equal to one - or be @code{COMPLEX}.
+@end multitable
+
+@item @emph{Return value}:
+The return value is of the same type and kind as @var{X}.
+The real part of the result is in degrees and lies in the range
+@math{-90 \leq \Re \asin(x) \leq 90}.
+
+@item @emph{Example}:
+@smallexample
+program test_asind
+  real(8) :: x = 0.866_8
+  x = asind(x)
+end program test_asind
+@end smallexample
+
+@item @emph{Specific names}:
+@multitable @columnfractions .20 .20 .20 .25
+@item Name            @tab Argument          @tab Return type       @tab Standard
+@item @code{ASIND(X)}  @tab @code{REAL(4) X}  @tab @code{REAL(4)}    @tab GNU Extension
+@item @code{DASIND(X)} @tab @code{REAL(8) X}  @tab @code{REAL(8)}    @tab GNU Extension
+@end multitable
+
+@item @emph{See also}:
+Inverse function: @ref{SIND}
+Radians function: @ref{ASIN}
 
 @end table
 
@@ -1458,6 +1588,71 @@ end program test_atan
 
 @item @emph{See also}:
 Inverse function: @ref{TAN}
+Degrees function: @ref{ATAND}
+
+@end table
+
+
+
+@node ATAND
+@section @code{ATAND} --- Arctangent function, degrees
+@fnindex ATAND
+@fnindex DATAND
+@cindex trigonometric function, tangent, inverse, degrees
+@cindex tangent, inverse, degrees
+
+@table @asis
+@item @emph{Description}:
+@code{ATAND(X)} computes the arctangent of @var{X} in degrees (inverse of
+@ref{TAND}).
+
+This function is for compatibility only and should be avoided in favor of
+standard constructs wherever possible.
+
+@item @emph{Standard}:
+GNU Extension, enabled with @option{-fdec-math}.
+
+@item @emph{Class}:
+Elemental function
+
+@item @emph{Syntax}:
+@multitable @columnfractions .80
+@item @code{RESULT = ATAND(X)}
+@item @code{RESULT = ATAND(Y, X)}
+@end multitable
+
+@item @emph{Arguments}:
+@multitable @columnfractions .15 .70
+@item @var{X} @tab The type shall be @code{REAL} or @code{COMPLEX};
+if @var{Y} is present, @var{X} shall be REAL.
+@item @var{Y} shall be of the same type and kind as @var{X}.
+@end multitable
+
+@item @emph{Return value}:
+The return value is of the same type and kind as @var{X}.
+If @var{Y} is present, the result is identical to @code{ATAND2(Y,X)}.
+Otherwise, it is the arcus tangent of @var{X}, where the real part of
+the result is in degrees and lies in the range
+@math{-90 \leq \Re \atand(x) \leq 90}.
+
+@item @emph{Example}:
+@smallexample
+program test_atand
+  real(8) :: x = 2.866_8
+  x = atand(x)
+end program test_atand
+@end smallexample
+
+@item @emph{Specific names}:
+@multitable @columnfractions .20 .20 .20 .25
+@item Name            @tab Argument          @tab Return type       @tab Standard
+@item @code{ATAND(X)}  @tab @code{REAL(4) X}  @tab @code{REAL(4)}    @tab GNU Extension
+@item @code{DATAND(X)} @tab @code{REAL(8) X}  @tab @code{REAL(8)}    @tab GNU Extension
+@end multitable
+
+@item @emph{See also}:
+Inverse function: @ref{TAND}
+Radians function: @ref{ATAN}
 
 @end table
 
@@ -1473,7 +1668,7 @@ Inverse function: @ref{TAN}
 @table @asis
 @item @emph{Description}:
 @code{ATAN2(Y, X)} computes the principal value of the argument
-function of the complex number @math{X + i Y}. This function can
+function of the complex number @math{X + i Y}.  This function can
 be used to transform from Cartesian into polar coordinates and
 allows to determine the angle in the correct quadrant.
 
@@ -1518,6 +1713,78 @@ end program test_atan2
 @item @code{ATAN2(X, Y)}  @tab @code{REAL(4) X, Y} @tab @code{REAL(4)} @tab Fortran 77 and later
 @item @code{DATAN2(X, Y)} @tab @code{REAL(8) X, Y} @tab @code{REAL(8)} @tab Fortran 77 and later
 @end multitable
+
+@item @emph{See also}:
+Alias: @ref{ATAN}
+Degrees function: @ref{ATAN2D}
+
+@end table
+
+
+
+@node ATAN2D
+@section @code{ATAN2D} --- Arctangent function, degrees
+@fnindex ATAN2D
+@fnindex DATAN2D
+@cindex trigonometric function, tangent, inverse, degrees
+@cindex tangent, inverse, degrees
+
+@table @asis
+@item @emph{Description}:
+@code{ATAN2D(Y, X)} computes the principal value of the argument
+function of the complex number @math{X + i Y} in degrees.  This function can
+be used to transform from Cartesian into polar coordinates and
+allows to determine the angle in the correct quadrant.
+
+This function is for compatibility only and should be avoided in favor of
+standard constructs wherever possible.
+
+@item @emph{Standard}:
+GNU Extension, enabled with @option{-fdec-math}.
+
+@item @emph{Class}:
+Elemental function
+
+@item @emph{Syntax}:
+@code{RESULT = ATAN2D(Y, X)}
+
+@item @emph{Arguments}:
+@multitable @columnfractions .15 .70
+@item @var{Y} @tab The type shall be @code{REAL}.
+@item @var{X} @tab The type and kind type parameter shall be the same as @var{Y}.
+If @var{Y} is zero, then @var{X} must be nonzero.
+@end multitable
+
+@item @emph{Return value}:
+The return value has the same type and kind type parameter as @var{Y}. It
+is the principal value of the complex number @math{X + i Y}.  If @var{X}
+is nonzero, then it lies in the range @math{-180 \le \atan (x) \leq 180}.
+The sign is positive if @var{Y} is positive.  If @var{Y} is zero, then
+the return value is zero if @var{X} is strictly positive, @math{180} if
+@var{X} is negative and @var{Y} is positive zero (or the processor does
+not handle signed zeros), and @math{-180} if @var{X} is negative and
+@var{Y} is negative zero.  Finally, if @var{X} is zero, then the
+magnitude of the result is @math{90}.
+
+@item @emph{Example}:
+@smallexample
+program test_atan2d
+  real(4) :: x = 1.e0_4, y = 0.5e0_4
+  x = atan2d(y,x)
+end program test_atan2d
+@end smallexample
+
+@item @emph{Specific names}:
+@multitable @columnfractions .20 .20 .20 .25
+@item Name                @tab Argument            @tab Return type    @tab Standard
+@item @code{ATAN2D(X, Y)}  @tab @code{REAL(4) X, Y} @tab @code{REAL(4)} @tab GNU Extension
+@item @code{DATAN2D(X, Y)} @tab @code{REAL(8) X, Y} @tab @code{REAL(8)} @tab GNU Extension
+@end multitable
+
+@item @emph{See also}:
+Alias: @ref{ATAND}
+Radians function: @ref{ATAN2}
+
 @end table
 
 
@@ -3895,6 +4162,70 @@ end program test_cos
 
 @item @emph{See also}:
 Inverse function: @ref{ACOS}
+Degrees function: @ref{COSD}
+
+@end table
+
+
+
+@node COSD
+@section @code{COSD} --- Cosine function, degrees
+@fnindex COSD
+@fnindex DCOSD
+@fnindex CCOSD
+@fnindex ZCOSD
+@fnindex CDCOSD
+@cindex trigonometric function, cosine, degrees
+@cindex cosine, degrees
+
+@table @asis
+@item @emph{Description}:
+@code{COSD(X)} computes the cosine of @var{X} in degrees.
+
+This function is for compatibility only and should be avoided in favor of
+standard constructs wherever possible.
+
+@item @emph{Standard}:
+GNU Extension, enabled with @option{-fdec-math}.
+
+@item @emph{Class}:
+Elemental function
+
+@item @emph{Syntax}:
+@code{RESULT = COSD(X)}
+
+@item @emph{Arguments}:
+@multitable @columnfractions .15 .70
+@item @var{X} @tab The type shall be @code{REAL} or
+@code{COMPLEX}.
+@end multitable
+
+@item @emph{Return value}:
+The return value is of the same type and kind as @var{X}. The real part
+of the result is in degrees.  If @var{X} is of the type @code{REAL},
+the return value lies in the range @math{ -1 \leq \cosd (x) \leq 1}.
+
+@item @emph{Example}:
+@smallexample
+program test_cosd
+  real :: x = 0.0
+  x = cosd(x)
+end program test_cosd
+@end smallexample
+
+@item @emph{Specific names}:
+@multitable @columnfractions .20 .20 .20 .25
+@item Name            @tab Argument            @tab Return type       @tab Standard
+@item @code{COSD(X)}   @tab @code{REAL(4) X}    @tab @code{REAL(4)}    @tab GNU Extension
+@item @code{DCOSD(X)}  @tab @code{REAL(8) X}    @tab @code{REAL(8)}    @tab GNU Extension
+@item @code{CCOSD(X)}  @tab @code{COMPLEX(4) X} @tab @code{COMPLEX(4)} @tab GNU Extension
+@item @code{ZCOSD(X)}  @tab @code{COMPLEX(8) X} @tab @code{COMPLEX(8)} @tab GNU extension
+@item @code{CDCOSD(X)} @tab @code{COMPLEX(8) X} @tab @code{COMPLEX(8)} @tab GNU extension
+@end multitable
+
+@item @emph{See also}:
+Inverse function: @ref{ACOSD}
+Radians function: @ref{COS}
 
 @end table
 
@@ -3954,6 +4285,115 @@ Inverse function: @ref{ACOSH}
 
 
 
+@node COTAN
+@section @code{COTAN} --- Cotangent function
+@fnindex COTAN
+@fnindex DCOTAN
+@cindex trigonometric function, cotangent
+@cindex cotangent
+
+@table @asis
+@item @emph{Description}:
+@code{COTAN(X)} computes the cotangent of @var{X}. Equivalent to @code{COS(x)}
+divided by @code{SIN(x)}, or @code{1 / TAN(x)}.
+
+This function is for compatibility only and should be avoided in favor of
+standard constructs wherever possible.
+
+@item @emph{Standard}:
+GNU Extension, enabled with @option{-fdec-math}.
+
+@item @emph{Class}:
+Elemental function
+
+@item @emph{Syntax}:
+@code{RESULT = COTAN(X)}
+
+@item @emph{Arguments}:
+@multitable @columnfractions .15 .70
+@item @var{X} @tab The type shall be @code{REAL} or @code{COMPLEX}.
+@end multitable
+
+@item @emph{Return value}:
+The return value has same type and kind as @var{X}, and its value is in radians.
+
+@item @emph{Example}:
+@smallexample
+program test_cotan
+  real(8) :: x = 0.165_8
+  x = cotan(x)
+end program test_cotan
+@end smallexample
+
+@item @emph{Specific names}:
+@multitable @columnfractions .20 .20 .20 .25
+@item Name            @tab Argument          @tab Return type     @tab Standard
+@item @code{COTAN(X)}   @tab @code{REAL(4) X}  @tab @code{REAL(4)}  @tab GNU Extension
+@item @code{DCOTAN(X)}  @tab @code{REAL(8) X}  @tab @code{REAL(8)}  @tab GNU Extension
+@end multitable
+
+@item @emph{See also}:
+Converse function: @ref{TAN}
+Degrees function: @ref{COTAND}
+@end table
+
+
+
+@node COTAND
+@section @code{COTAND} --- Cotangent function, degrees
+@fnindex COTAND
+@fnindex DCOTAND
+@cindex trigonometric function, cotangent, degrees
+@cindex cotangent, degrees
+
+@table @asis
+@item @emph{Description}:
+@code{COTAND(X)} computes the cotangent of @var{X} in degrees.  Equivalent to
+@code{COSD(x)} divided by @code{SIND(x)}, or @code{1 / TAND(x)}.
+
+@item @emph{Standard}:
+GNU Extension, enabled with @option{-fdec-math}.
+
+This function is for compatibility only and should be avoided in favor of
+standard constructs wherever possible.
+
+@item @emph{Class}:
+Elemental function
+
+@item @emph{Syntax}:
+@code{RESULT = COTAND(X)}
+
+@item @emph{Arguments}:
+@multitable @columnfractions .15 .70
+@item @var{X} @tab The type shall be @code{REAL} or @code{COMPLEX}.
+@end multitable
+
+@item @emph{Return value}:
+The return value has same type and kind as @var{X}, and its value is in degrees.
+
+@item @emph{Example}:
+@smallexample
+program test_cotand
+  real(8) :: x = 0.165_8
+  x = cotand(x)
+end program test_cotand
+@end smallexample
+
+@item @emph{Specific names}:
+@multitable @columnfractions .20 .20 .20 .25
+@item Name            @tab Argument          @tab Return type     @tab Standard
+@item @code{COTAND(X)}   @tab @code{REAL(4) X}  @tab @code{REAL(4)}  @tab GNU Extension
+@item @code{DCOTAND(X)}  @tab @code{REAL(8) X}  @tab @code{REAL(8)}  @tab GNU Extension
+@end multitable
+
+@item @emph{See also}:
+Converse function: @ref{TAND}
+Radians function: @ref{COTAN}
+
+@end table
+
+
+
 @node COUNT
 @section @code{COUNT} --- Count function
 @fnindex COUNT
@@ -12390,7 +12830,69 @@ end program test_sin
 @end multitable
 
 @item @emph{See also}:
-@ref{ASIN}
+Inverse function: @ref{ASIN}
+Degrees function: @ref{SIND}
+@end table
+
+
+
+@node SIND
+@section @code{SIND} --- Sine function, degrees
+@fnindex SIND
+@fnindex DSIND
+@fnindex CSIND
+@fnindex ZSIND
+@fnindex CDSIND
+@cindex trigonometric function, sine, degrees
+@cindex sine, degrees
+
+@table @asis
+@item @emph{Description}:
+@code{SIND(X)} computes the sine of @var{X} in degrees.
+
+This function is for compatibility only and should be avoided in favor of
+standard constructs wherever possible.
+
+@item @emph{Standard}:
+GNU Extension, enabled with @option{-fdec-math}.
+
+@item @emph{Class}:
+Elemental function
+
+@item @emph{Syntax}:
+@code{RESULT = SIND(X)}
+
+@item @emph{Arguments}:
+@multitable @columnfractions .15 .70
+@item @var{X} @tab The type shall be @code{REAL} or
+@code{COMPLEX}.
+@end multitable
+
+@item @emph{Return value}:
+The return value has same type and kind as @var{X}, and its value is in degrees.
+
+@item @emph{Example}:
+@smallexample
+program test_sind
+  real :: x = 0.0
+  x = sind(x)
+end program test_sind
+@end smallexample
+
+@item @emph{Specific names}:
+@multitable @columnfractions .20 .20 .20 .25
+@item Name            @tab Argument             @tab Return type       @tab Standard
+@item @code{SIND(X)}   @tab @code{REAL(4) X}     @tab @code{REAL(4)}    @tab GNU Extension
+@item @code{DSIND(X)}  @tab @code{REAL(8) X}     @tab @code{REAL(8)}    @tab GNU Extension
+@item @code{CSIND(X)}  @tab @code{COMPLEX(4) X}  @tab @code{COMPLEX(4)} @tab GNU Extension
+@item @code{ZSIND(X)}  @tab @code{COMPLEX(8) X}  @tab @code{COMPLEX(8)} @tab GNU Extension
+@item @code{CDSIND(X)} @tab @code{COMPLEX(8) X}  @tab @code{COMPLEX(8)} @tab GNU Extension
+@end multitable
+
+@item @emph{See also}:
+Inverse function: @ref{ASIND}
+Radians function: @ref{SIN}
+
 @end table
 
 
@@ -13151,7 +13653,7 @@ Elemental function
 @end multitable
 
 @item @emph{Return value}:
-The return value has same type and kind as @var{X}.
+The return value has same type and kind as @var{X}, and its value is in radians.
 
 @item @emph{Example}:
 @smallexample
@@ -13169,7 +13671,61 @@ end program test_tan
 @end multitable
 
 @item @emph{See also}:
-@ref{ATAN}
+Inverse function: @ref{ATAN}
+Degrees function: @ref{TAND}
+@end table
+
+
+
+@node TAND
+@section @code{TAND} --- Tangent function, degrees
+@fnindex TAND
+@fnindex DTAND
+@cindex trigonometric function, tangent, degrees
+@cindex tangent, degrees
+
+@table @asis
+@item @emph{Description}:
+@code{TAND(X)} computes the tangent of @var{X} in degrees.
+
+This function is for compatibility only and should be avoided in favor of
+standard constructs wherever possible.
+
+@item @emph{Standard}:
+GNU Extension, enabled with @option{-fdec-math}.
+
+@item @emph{Class}:
+Elemental function
+
+@item @emph{Syntax}:
+@code{RESULT = TAND(X)}
+
+@item @emph{Arguments}:
+@multitable @columnfractions .15 .70
+@item @var{X} @tab The type shall be @code{REAL} or @code{COMPLEX}.
+@end multitable
+
+@item @emph{Return value}:
+The return value has same type and kind as @var{X}, and its value is in degrees.
+
+@item @emph{Example}:
+@smallexample
+program test_tand
+  real(8) :: x = 0.165_8
+  x = tand(x)
+end program test_tand
+@end smallexample
+
+@item @emph{Specific names}:
+@multitable @columnfractions .20 .20 .20 .25
+@item Name            @tab Argument          @tab Return type     @tab Standard
+@item @code{TAND(X)}   @tab @code{REAL(4) X}  @tab @code{REAL(4)}  @tab GNU Extension
+@item @code{DTAND(X)}  @tab @code{REAL(8) X}  @tab @code{REAL(8)}  @tab GNU Extension
+@end multitable
+
+@item @emph{See also}:
+Inverse function: @ref{ATAND}
+Radians function: @ref{TAN}
 @end table
 
 
diff --git a/gcc/fortran/invoke.texi b/gcc/fortran/invoke.texi
index 268d155..655ee6f 100644
--- a/gcc/fortran/invoke.texi
+++ b/gcc/fortran/invoke.texi
@@ -116,7 +116,7 @@ by type.  Explanations are in the following sections.
 @xref{Fortran Dialect Options,,Options controlling Fortran dialect}.
 @gccoptlist{-fall-intrinsics -fbackslash -fcray-pointer -fd-lines-as-code @gol
 -fd-lines-as-comments @gol
--fdec -fdec-structure -fdec-intrinsic-ints -fdec-static @gol
+-fdec -fdec-structure -fdec-intrinsic-ints -fdec-static -fdec-math @gol
 -fdefault-double-8 -fdefault-integer-8 @gol
 -fdefault-real-8 -fdollar-ok -ffixed-line-length-@var{n} @gol
 -ffixed-line-length-none -ffree-form -ffree-line-length-@var{n} @gol
@@ -255,6 +255,11 @@ instead where possible.
 Enable B/I/J/K kind variants of existing integer functions (e.g. BIAND, IIAND,
 JIAND, etc...). For a complete list of intrinsics see the full documentation.
 
+@item -fdec-math
+@opindex @code{fdec-math}
+Enable legacy math intrinsics such as COTAN and degree-valued trigonometric
+functions (e.g. TAND, ATAND, etc...) for compatability with older code.
+
 @item -fdec-static
 @opindex @code{fdec-static}
 Enable DEC-style STATIC and AUTOMATIC attributes to explicitly specify
diff --git a/gcc/fortran/iresolve.c b/gcc/fortran/iresolve.c
index ecea1c3..f4f81b2 100644
--- a/gcc/fortran/iresolve.c
+++ b/gcc/fortran/iresolve.c
@@ -673,6 +673,86 @@ 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
+    && 0 == strncmp ("__", f->value.function.name, 2);
+}
+
+/* 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)
 {
@@ -2578,6 +2658,159 @@ 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_d (factor->value.real, factor->value.real, 180.0, 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;
+
+  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_d (factor->value.real, 180.0, 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:
+	break;
+    }
+
+  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);
+}
+
+
 void
 gfc_resolve_image_index (gfc_expr *f, gfc_expr *array ATTRIBUTE_UNUSED,
 			 gfc_expr *sub ATTRIBUTE_UNUSED)
diff --git a/gcc/fortran/lang.opt b/gcc/fortran/lang.opt
index ef421d3..b563e09 100644
--- a/gcc/fortran/lang.opt
+++ b/gcc/fortran/lang.opt
@@ -428,6 +428,10 @@ fdec-intrinsic-ints
 Fortran Var(flag_dec_intrinsic_ints)
 Enable kind-specific variants of integer intrinsic functions.
 
+fdec-math
+Fortran Var(flag_dec_math)
+Enable legacy math intrinsics for compatibility.
+
 fdec-structure
 Fortran
 Enable support for DEC STRUCTURE/RECORD.
diff --git a/gcc/fortran/options.c b/gcc/fortran/options.c
index 5881a88..93403f7c 100644
--- a/gcc/fortran/options.c
+++ b/gcc/fortran/options.c
@@ -55,6 +55,7 @@ set_dec_flags (int value)
     gfc_option.flag_dec_structure  = value;
     flag_dec_intrinsic_ints = value;
     flag_dec_static = value;
+    flag_dec_math = value;
 }
 
 
diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c
index ad547a1..bf60f74 100644
--- a/gcc/fortran/simplify.c
+++ b/gcc/fortran/simplify.c
@@ -1706,6 +1706,152 @@ 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:
+	 break;
+  }
+
+  gfc_internal_error ("in simplify_trig_call(): Bad intrinsic");
+  return NULL;
+}
+
+/* Convert a floating-point number from radians to degrees.  */
+
+static void
+degrees_f (mpfr_t x, mp_rnd_t rnd_mode)
+{
+    mpfr_t tmp;
+    mpfr_init (tmp);
+
+    /* Set x = x % 2pi to avoid offsets with large angles.  */
+    mpfr_const_pi (tmp, rnd_mode);
+    mpfr_mul_ui (tmp, tmp, 2, rnd_mode);
+    mpfr_fmod (tmp, x, tmp, rnd_mode);
+
+    /* Set x = x * 180.  */
+    mpfr_mul_d (x, x, 180.0, rnd_mode);
+
+    /* Set x = x / pi.  */
+    mpfr_const_pi (tmp, rnd_mode);
+    mpfr_div (x, x, tmp, rnd_mode);
+
+    mpfr_clear (tmp);
+}
+
+/* Convert a floating-point number from degrees to radians.  */
+
+static void
+radians_f (mpfr_t x, mp_rnd_t rnd_mode)
+{
+    mpfr_t tmp;
+    mpfr_init (tmp);
+
+    /* Set x = x % 360 to avoid offsets with large angles.  */
+    mpfr_fmod_d (tmp, x, 360.0, rnd_mode);
+
+    /* Set x = x * pi.  */
+    mpfr_const_pi (tmp, rnd_mode);
+    mpfr_mul (x, x, tmp, rnd_mode);
+
+    /* Set x = x / 180.  */
+    mpfr_div_d (x, x, 180.0, rnd_mode);
+
+    mpfr_clear (tmp);
+}
+
+
+/* Convert argument to radians before calling a trig function.  */
+
+gfc_expr *
+gfc_simplify_trigd (gfc_expr *icall)
+{
+  gfc_expr *arg;
+
+  arg = icall->value.function.actual->expr;
+
+  if (arg->ts.type != BT_REAL)
+    gfc_internal_error ("in gfc_simplify_trigd(): Bad type");
+
+  if (arg->expr_type == EXPR_CONSTANT)
+    /* Convert constant to radians before passing off to simplifier.  */
+    radians_f (arg->value.real, GFC_RND_MODE);
+
+  /* Let the usual simplifier take over - we just simplified the arg.  */
+  return simplify_trig_call (icall);
+}
+
+/* Convert result of an inverse trig function to degrees.  */
+
+gfc_expr *
+gfc_simplify_atrigd (gfc_expr *icall)
+{
+  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 (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;
+  }
+
+  /* Let gfc_resolve_atrigd take care of the non-constant case.  */
+  return NULL;
+}
+
+/* Convert the result of atan2 to degrees.  */
+
+gfc_expr *
+gfc_simplify_atan2d (gfc_expr *y, 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 && y->expr_type == EXPR_CONSTANT)
+    {
+      result = gfc_simplify_atan2 (y, x);
+      if (result != NULL)
+	{
+	  degrees_f (result->value.real, GFC_RND_MODE);
+	  return result;
+	}
+    }
+
+  /* Let gfc_resolve_atan2d take care of the non-constant case.  */
+  return NULL;
+}
 
 gfc_expr *
 gfc_simplify_cos (gfc_expr *x)
@@ -6244,6 +6390,41 @@ gfc_simplify_sum (gfc_expr *array, gfc_expr *dim, gfc_expr *mask)
 
 
 gfc_expr *
+gfc_simplify_cotan (gfc_expr *x)
+{
+  gfc_expr *result;
+  mpc_t swp, *val;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+
+  switch (x->ts.type)
+    {
+      case BT_REAL:
+	mpfr_cot (result->value.real, x->value.real, GFC_RND_MODE);
+	break;
+
+      case BT_COMPLEX:
+	/* 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_div (*val, swp, *val, GFC_MPC_RND_MODE);
+	mpc_clear (swp);
+	break;
+
+      default:
+	gcc_unreachable ();
+    }
+
+  return range_check (result, "COTAN");
+}
+
+
+gfc_expr *
 gfc_simplify_tan (gfc_expr *x)
 {
   gfc_expr *result;
diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog
index 1b6044c..e1ed8e5 100644
--- a/gcc/testsuite/ChangeLog
+++ b/gcc/testsuite/ChangeLog
@@ -1,3 +1,7 @@
+2016-10-11  Fritz Reese  <fritzoreese@gmail.com>
+
+	* gfortran.dg/dec_math.f90: New testsuite.
+
 2016-10-11  Senthil Kumar Selvaraj  <senthil_kumar.selvaraj@atmel.com>
 
 	* gcc.dg/tree-ssa/pr59597.c: Typedef  __INT32_TYPE__ to i32.
diff --git a/gcc/testsuite/gfortran.dg/dec_math.f90 b/gcc/testsuite/gfortran.dg/dec_math.f90
new file mode 100644
index 0000000..857a261
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/dec_math.f90
@@ -0,0 +1,289 @@
+! { dg-options "-fdec-math" }
+! { dg-do run }
+!
+! Test extra math intrinsics offered by -fdec-math.
+!
+
+  subroutine cmpf(f1, f2, tolerance, str)
+    implicit none
+    real(4), intent(in) :: f1, f2, tolerance
+    character(len=*), intent(in) :: str
+    if ( abs(f2 - f1) .gt. tolerance ) then
+      write (*, '(A,F12.6,F12.6)') str, f1, f2
+      call abort()
+    endif
+  endsubroutine
+
+  subroutine cmpd(d1, d2, tolerance, str)
+    implicit none
+    real(8), intent(in) :: d1, d2, tolerance
+    character(len=*), intent(in) :: str
+    if ( dabs(d2 - d1) .gt. tolerance ) then
+      write (*, '(A,F12.6,F12.6)') str, d1, d2
+      call abort()
+    endif
+  endsubroutine
+
+implicit none
+
+  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
+
+! inputs
+real(4) :: f_i1, f_i2
+real(4), volatile :: xf
+real(8) :: d_i1, d_i2
+real(8), volatile :: xd
+
+! expected outputs from (oe) default (oxe) expression
+real(4) :: f_oe, f_oxe
+real(8) :: d_oe, d_oxe
+
+! 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
+
+! 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
+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)
+
+! 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")
+
+! Input
+f_i1 = 60.0_4
+d_i1 = 60.0_8
+
+! 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)
+
+! 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)
+
+! 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")
+
+! Input
+f_i1 = 60.0_4
+d_i1 = 60.0_8
+
+! 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)
+
+! 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")
+
+! Input
+f_i1 = 2.679676_4
+f_i2 = 1.0_4
+d_i1 = 2.679676_8
+d_i2 = 1.0_8
+
+! 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)
+
+! 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")
+
+! Input
+f_i1 = 1.5874993_4
+d_i1 = 1.5874993_8
+
+! 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)
+
+! 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")
+
+! Input
+f_i1 = 0.6_4
+d_i1 = 0.6_8
+
+! 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)
+
+! 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")
+
+! Input
+f_i1 = 0.6_4
+d_i1 = 0.6_8
+
+! 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)
+
+! 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")
+
+! Input
+f_i1 = 60.0_4
+d_i1 = 60.0_8
+
+! 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)
+
+! 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")
+
+end