#include "fpreproc/library_names.h"
#include "fpreproc/f77_dcomplx.h"
#include "fpreproc/fp_const.h"
*DECK DCOT
      REAL*8 FUNCTION DCOT (X)
C***BEGIN PROLOGUE  DCOT
C***PURPOSE  Compute the cotangent.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C4A
C***TYPE      REAL*8 (COT-S, DCOT-D, CCOT-C)
C***KEYWORDS  COTANGENT, ELEMENTARY FUNCTIONS, FNLIB, TRIGONOMETRIC
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DCOT(X) calculates the double precision trigonometric cotangent
C for double precision argument X.  X is in units of radians.
C
C Series for COT        on the interval  0.          to  6.25000E-02
C                                        with weighted error   5.52E-34
C                                         log weighted error  33.26
C                               significant figures required  32.34
C                                    decimal places required  33.85
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, DCSEVL, INITDS, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   920618  Removed space from variable names.  (RWC, WRB)
C***END PROLOGUE  DCOT
C============
C idecl:  explicitize implicit INTEGER declarations:
      IMPLICIT NONE
      INTEGER nterms,initds,ifn
C============
      REAL*8 X, COTCS(15), AINTY, AINTY2, PI2REC, SQEPS,
     1  XMAX, XMIN, XSML, Y, YREM, PRODBG, DCSEVL, D1MACH
      LOGICAL FIRST
      SAVE COTCS, PI2REC, NTERMS, XMAX, XSML, XMIN, SQEPS, FIRST
      DATA COTCS(  1) / +CONST(.240259160982956302509553617744970,+0)/
      DATA COTCS(  2) / -CONST(.165330316015002278454746025255758,-1)/
      DATA COTCS(  3) / -CONST(.429983919317240189356476228239895,-4)/
      DATA COTCS(  4) / -CONST(.159283223327541046023490851122445,-6)/
      DATA COTCS(  5) / -CONST(.619109313512934872588620579343187,-9)/
      DATA COTCS(  6) / -CONST(.243019741507264604331702590579575,-11)/
      DATA COTCS(  7) / -CONST(.956093675880008098427062083100000,-14)/
      DATA COTCS(  8) / -CONST(.376353798194580580416291539706666,-16)/
      DATA COTCS(  9) / -CONST(.148166574646746578852176794666666,-18)/
      DATA COTCS( 10) / -CONST(.583335658903666579477984000000000,-21)/
      DATA COTCS( 11) / -CONST(.229662646964645773928533333333333,-23)/
      DATA COTCS( 12) / -CONST(.904197057307483326719999999999999,-26)/
      DATA COTCS( 13) / -CONST(.355988551920600064000000000000000,-28)/
      DATA COTCS( 14) / -CONST(.140155139824298666666666666666666,-30)/
      DATA COTCS( 15) / -CONST(.551800436872533333333333333333333,-33)/
      DATA PI2REC / CONST(.011619772367581343075535053490057,0)/
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DCOT
      IF (FIRST) THEN
         NTERMS = INITDS (COTCS, 15, CONST(0.1,0)*AREAL(D1MACH(3)) )
         XMAX = CONST(1.0,0)/D1MACH(4)
         XSML = SQRT(CONST(3.0,0)*D1MACH(3))
         XMIN = EXP (MAX(LOG(D1MACH(1)),
     >      -LOG(D1MACH(2))) + CONST(0.01,0))
         SQEPS = SQRT(D1MACH(4))
      ENDIF
      FIRST = .FALSE.
C
      Y = ABS(X)
      IF (Y .LT. XMIN) CALL XERMSG ('SLATEC', 'DCOT',
     +   'ABS(X) IS ZERO OR SO SMALL DCOT OVERFLOWS', 2, 2)
      IF (Y .GT. XMAX) CALL XERMSG ('SLATEC', 'DCOT',
     +   'NO PRECISION BECAUSE ABS(X) IS TOO BIG', 3, 2)
C
C CAREFULLY COMPUTE Y * (2/PI) = (AINT(Y) + REM(Y)) * (.625 + PI2REC)
C = AINT(.625*Y) + REM(.625*Y) + Y*PI2REC  =  AINT(.625*Y) + Z
C = AINT(.625*Y) + AINT(Z) + REM(Z)
C
      AINTY = AINT (Y)
      YREM = Y - AINTY
      PRODBG = CONST(0.625,0)*AINTY
      AINTY = AINT (PRODBG)
      Y = (PRODBG-AINTY) + CONST(0.625,0)*YREM + PI2REC*Y
      AINTY2 = AINT (Y)
      AINTY = AINTY + AINTY2
      Y = Y - AINTY2
C
      IFN = MOD (AINTY, CONST(2.0,0))
      IF (IFN.EQ.1) Y = CONST(1.0,0)- Y
C
      IF (ABS(X) .GT. CONST(0.5,0).AND. Y .LT. ABS(X)*SQEPS) CALL XERMSG
     +   ('SLATEC', 'DCOT',
     +   'ANSWER LT HALF PRECISION, ABS(X) TOO BIG OR X NEAR N*PI ' //
     +   '(N.NE.0)', 1, 1)
C
      IF (Y.GT.CONST(0.25,0)) GO TO 20
      DCOT = CONST(1.0,0)/X
      IF (Y.GT.XSML) DCOT = (CONST(0.5,0)+
     1   DCSEVL (CONST(32.0,0)*Y*Y-CONST(1.,0), COTCS,
     1   NTERMS)) / Y
      GO TO 40
C
 20   IF (Y.GT.CONST(0.5,0)) GO TO 30
      DCOT = (CONST(0.5,0)+ DCSEVL (CONST(8.,0)*Y*Y-CONST(1.,0),
     >   COTCS, NTERMS))/(CONST(0.5,0)*Y)
      DCOT = (DCOT*DCOT-CONST(1.,0))*CONST(0.5,0)/DCOT
      GO TO 40
C
 30   DCOT = (CONST(0.5,0)+ DCSEVL (CONST(2.,0)*Y*Y-CONST(1.,0),
     >   COTCS, NTERMS))/(CONST(.25,0)*Y)
      DCOT = (DCOT*DCOT-CONST(1.,0))*CONST(0.5,0)/DCOT
      DCOT = (DCOT*DCOT-CONST(1.,0))*CONST(0.5,0)/DCOT
C
 40   IF (X.NE.CONST(0.,0)) DCOT = SIGN (DCOT, X)
      IF (IFN.EQ.1) DCOT = -DCOT
C
      RETURN
      END
! 22Jun2000 fgtok -s r8_precision.sub "r8con.csh conversion"
! 22Jun2000 fgtok
! 22Jun2000 fgtok
! 22Jun2000 fgtok
