C	shcrest

C	Sven Ivansson 1990




	COMPLEX*16 FUNCTION SHCREST (Z)
C	to compute  shcrest = sinh(z) / z
C	accurately, even for small z

	IMPLICIT NONE

	COMPLEX*16 Z

	INTEGER N, I
	REAL*8 EPS
	COMPLEX*16 Z2, SREL(20), CSLASK1,CSLASK2

	DATA EPS /1.0D-15/  ! desired relative accuracy (approximately), 
C	strict analysis possible with MacLaurin remainder term

	Z2 = Z*Z

	IF ( DMAX1(DABS(DREAL(Z)),DABS(DIMAG(Z))) .GT. 0.5D0 ) THEN
	  CSLASK1 = CDEXP(Z)
	  CSLASK2 = 1/CSLASK1
	  SHCREST = (CSLASK1-CSLASK2) / (2*Z)
	ELSE
C	  shcrest = 1 + z**2/6 + z**4/120 + ...
	  N = 1
	  I = 2
	  SREL(1) = Z2/6
	  DO WHILE ( DMAX1(DABS(DREAL(SREL(N))),DABS(DIMAG(SREL(N)))) 
     $	  .GT. EPS )
	    N = N + 1
	    I = I + 2
	    SREL(N) = SREL(N-1) * Z2 / (I*(I+1))  !  z**(2*n) / ((2*n+1)!)
	  END DO
	  SHCREST = (0.0D0,0.0D0)
	  DO I = N, 1, -1
	    SHCREST = SHCREST + SREL(I)
	  END DO
	  SHCREST = (1.0D0,0.0D0) + SHCREST
	END IF

	RETURN
	END