
C	hyprest & hyprest2

C	Sven Ivansson 1990




	SUBROUTINE HYPREST (Z, CHYPR,SHYPR)
C	to compute  chypr = (cosh(z)-1) / (z*z/2)
C	            shypr = (sinh(z)-z) / (z*z/2)
C	accurately, even for small z

	IMPLICIT NONE

	COMPLEX*32 Z
	COMPLEX*32 CHYPR,SHYPR

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

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

	Z2 = Z*Z

	IF ( MAX(ABS(REAL(Z)),ABS(AIMAG(Z))) .GT. 0.5Q0 ) THEN
	  CSLASK1 = EXP(Z)
	  CSLASK2 = 1/CSLASK1
	  CHYPR = (CSLASK1+CSLASK2-2.0Q0) / Z2
	  SHYPR = (CSLASK1-CSLASK2-2.0Q0*Z) / Z2
	ELSE
C	  chypr =           1 + z**2/12 + z**4/360 + ...
C	  shypr = (z/3) * ( 1 + z**2/20 + z**4/840 + ... )
	  N = 1
	  I = 4
	  CREL(1) = Z2/12
	  SREL(1) = Z2/20
	  DO WHILE ( MAX(ABS(REAL(CREL(N))),ABS(AIMAG(CREL(N)))) .
     $	  GT. EPS )
	    N = N + 1
	    I = I + 2
	    CREL(N) = CREL(N-1) * Z2 / ((I-1)*I)  !  2 * z**(2*n) / ((2*n+2)!)
	    SREL(N) = SREL(N-1) * Z2 / (I*(I+1))  !  6 * z**(2*n) / ((2*n+3)!)
	  END DO
	  CHYPR = (0.0Q0,0.0Q0)
	  SHYPR = (0.0Q0,0.0Q0)
	  DO I = N, 1, -1
	    CHYPR = CHYPR + CREL(I)
	    SHYPR = SHYPR + SREL(I)
	  END DO
	  CHYPR =           (1.0Q0,0.0Q0) + CHYPR
	  SHYPR = (Z/3) * ( (1.0Q0,0.0Q0) + SHYPR )
	END IF

	RETURN
	END






	SUBROUTINE HYPREST2 (Z, BHYPR,CHYPR,SHYPR)
C	to compute  bhypr = 4 * ( cosh(z) - sinh(z)/z ) / (z*z)
C	            chypr = 4 * ( sinh(2z)/(2z) - (3+cosh(2z))/4 ) / (z*z)
C	                    - bhypr
C	            shypr = 2 * ( cosh(2z) - ((2z)/4)*sinh(2z) - 1 - 
C	                                ((2z)*(2z))/4 ) / (z*z*z)  
C	accurately, even for small z

	IMPLICIT NONE

	COMPLEX*32 Z
	COMPLEX*32 BHYPR,CHYPR,SHYPR

	INTEGER N, I
	REAL*16 EPS, SLASK
	COMPLEX*32 Z2, ZPZ,ZPZ2, BREL(20), CREL(20), SREL(20), 
     $	CSLASK1,CSLASK2, CSLASK, QCMPLX

	DATA EPS /1.0Q-30/  ! desired relative accuracy (approximately)

	Z2 = Z*Z
	ZPZ = 2 * Z
	ZPZ2 = ZPZ*ZPZ

	IF ( MAX(ABS(REAL(ZPZ)),ABS(AIMAG(ZPZ))) .GT. 0.5Q0 ) THEN 
	  CSLASK1 = EXP(Z)
	  CSLASK2 = 1/CSLASK1
	  BHYPR = (2/Z2) * ( (CSLASK1+CSLASK2) - (CSLASK1-CSLASK2)/Z )
	  CSLASK = CSLASK1*CSLASK1
	  CSLASK2 = CSLASK2*CSLASK2
	  CSLASK1 = CSLASK + CSLASK2
	  CSLASK2 = CSLASK - CSLASK2
	  CHYPR = ( 2*CSLASK2/ZPZ - (3+0.5Q0*CSLASK1) ) / Z2 - BHYPR
	  SHYPR = ( 2*(CSLASK1-2)/ZPZ - 0.5Q0*CSLASK2 - ZPZ ) / Z2
	ELSE
C	  bhypr =          4/3 + (2/15)*(z**2)    + ...
C	  chypr + bhypr =  2/3 - (1/30)*((2z)**2) + ...
C	  shypr =      -(4/45)*(z**3) * ( 1 + (1/28)*((2z)**2) +  ... )
	  N = 1
	  I = 5
	  BREL(1) = 2*Z2/15
	  CREL(1) = BREL(1)
	  SREL(1) = ZPZ2/28
	  DO WHILE ( MAX(ABS(REAL(BREL(N))),ABS(AIMAG(BREL(N)))) .
     $	  GT. EPS )
	    N = N + 1
	    I = I + 2
	    SLASK = (I-1)*I
	    BREL(N) = BREL(N-1) * Z2 * ( (N+1)/(N*SLASK) )
	    CREL(N) = CREL(N-1) * ZPZ2 * ( (N-0.5Q0)/((N-1.5Q0)*SLASK) )
	    SLASK = (I+2)*(I+3)
	    SREL(N) = SREL(N-1) * ZPZ2 * ( (N+1)/(N*SLASK) )
	  END DO
	  BHYPR = (0.0Q0,0.0Q0)
	  CHYPR = (0.0Q0,0.0Q0)
	  SHYPR = (0.0Q0,0.0Q0)
	  DO I = N, 1, -1
	    BHYPR = BHYPR + BREL(I)
	    CHYPR = CHYPR + CREL(I)
	    SHYPR = SHYPR + SREL(I)
	  END DO
	  CSLASK = CHYPR + BHYPR
	  BHYPR =           QCMPLX(4/3.0Q0,0.0Q0) + BHYPR
	  CHYPR =          -QCMPLX(2/3.0Q0,0.0Q0) - CSLASK
	  SHYPR = -(4/45.0Q0) * Z*Z2 * ( (1.0Q0,0.0Q0) + SHYPR )
	END IF

	RETURN
	END
