PROGRAM dcitest

!   DESCRIPTION:

!   This program tests the performance of DOUBLE-PRECISION codes which
!   calculate the cosine-integral special function, defined as

!       DCOSINT(x) = gamma + log(abs(x)) +
!                  integral from 0 to x of {[cos(t)-1]/t} dt

!   The program compares the value of DCOSINT(x+h) with the Taylor
!   series value about x. The derivatives of DCOSINT at x are
!   generated by a standard recurrence relation.


!   The program assumes the cosine-integral is calculated in a
!   function with the name DCOSINT. Thus users should append their
!   Ci code at the end of this file and edit the small DCOSINT
!   function to include the name of their code.


!   The program uses Cody's MACHAR subroutine to determine certain
!   machine parameters. This can be avoided if values for the following
!   parameters are available: EPS, IT, IBETA, MACHEP, XMIN, XMAX.
!   Values for certain machine-compiler combinations are given in

!       W.J. CODY  Algorithm 665: MACHAR: A subroutine to dynamically
!       determine machine parameters, ACM Trans. Math. Soft. 14 (1988) 303-311.

!   The MACHAR code and a random generator code written by Cody are included
!   in the file SICISUBD.FOR, though with the MACHAR code called DMACHR.

!   To ensure that MACHAR and the argument purification section work as
!   intended, this code should be compiled with any optimisation switched OFF.

!   Users should note that MACHAR tests the machine arithmetic
!   at its very limits. Thus it is possible that warning messages
!   might appear about underflows and similar events. These DO
!   NOT affect the test results. Such warnings can be avoided
!   by explicitly declaring the necessary machine parameters.

!   MODULES USED:

!      This code calls the following modules:

!      (a) DCIDER contained in this file

!      (b) DCOSINT contained in the file DCOSINT.F

!      (c) DMACHR, DREN and DTERMS contained in the file SICISUBD.F


!   AUTHOR: Allan MacLeod
!           Dept. of Mathematics and Statistics
!           University of Paisley
!           Scotland
!            (e-mail: macl_ms0@paisley.ac.uk)

IMPLICIT NONE
INTEGER, PARAMETER    :: dp = SELECTED_REAL_KIND(14, 60)

INTEGER   :: i, ibeta, ipt, it, itest, ndec, nterms, nt1, nt2, numeq,  &
             numla, numsm, tstart(8), tend(8)
REAL (dp) :: cider(0:100), ait, al, albeta, au, cix, ciy, delta, div, erloss, &
             ERR, hint, maxerr, maxpt, relerr, rmserr, sum, temp, x, y, z
REAL      :: dren, time_taken

!   Define constants

REAL (dp), PARAMETER :: alow(5) = (/ 0.25_dp, 2.0_dp, 4.0_dp, 8.0_dp, &
                                     13.5_dp /)
REAL (dp), PARAMETER :: aup(5) = (/ 0.375_dp, 2.5_dp, 4.5_dp, 8.5_dp, &
                                    14.0_dp /)
REAL (dp), PARAMETER :: pt0625 = 0.0625_dp, sixten = 16.0_dp, half = 0.5_dp
REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp, eps = EPSILON(one)
REAL (dp), PARAMETER :: gamma = 0.57721566490153286061_dp
REAL (dp), PARAMETER :: xmin = TINY(one), xmax = HUGE(one)

!   IOUT = output channel number

!   NPTS = number of test points per interval

INTEGER, PARAMETER :: iout = 6, npts = 5000, ntest = 5

INTERFACE
  FUNCTION dcosint(xvalue) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, PARAMETER    :: dp = SELECTED_REAL_KIND(14,60)
    REAL (dp), INTENT(IN) :: xvalue
    REAL (dp)             :: fn_val
  END FUNCTION dcosint
END INTERFACE

!   DETERMINE MACHINE PARAMETERS

ibeta = RADIX( one )
it = DIGITS( one )
ait = DBLE(it)
albeta = LOG(DBLE(ibeta))
ndec = ABS( INT( LOG10( DBLE(ibeta) * eps ) - half ) )

CALL DATE_AND_TIME( VALUES=tstart )

!   START TESTS

DO itest = 1 , ntest
  al = alow(itest)
  au = aup(itest)
  hint = ( au - al ) / DBLE(npts)
  
!   DETERMINE HOW MANY DERIVATIVES ARE NEEDED
  
  CALL dterms(ABS(al), eps, nt1)
  CALL dterms(ABS(au), eps, nt2)
  nterms = nt1
  IF ( nt2 > nt1 ) nterms = nt2
  
!   SET UP STATISTICS COLLECTION
  
  maxerr = zero
  rmserr = zero
  maxpt = zero
  numsm = 0
  numla = 0
  delta = pt0625
  IF ( itest == 2 .OR. itest == 4 ) delta = -pt0625
  DO ipt = 1 , npts
    CALL RANDOM_NUMBER( dren )
    x = al + hint * dren
    
!   PURIFY ARGUMENT
    
    z = sixten * x
    temp = z + x
    x = temp - z
    y = x + delta
    
!   CALCULATE CI(Y) AND CI(X+DELTA)
    
    ciy = dcosint(y)
    cix = dcosint(x)
    CALL dcider(x, nterms, ndec, cider)
    sum = cider(nterms)
    div = DBLE(nterms) + one
    DO i = 0 , nterms-1
      sum = sum * delta / div + cider(nterms-1-i)
      div = div - one
    END DO
    sum = delta * sum
    
!   DETERMINE ERROR AND UPDATE STATS
    
    ERR = ( ciy - cix ) - sum
    IF ( ERR < zero ) numsm = numsm + 1
    IF ( ERR > zero ) numla = numla + 1
    relerr = ABS(ERR/ciy)
    rmserr = rmserr + relerr * relerr
    IF ( relerr > maxerr ) THEN
      maxerr = relerr
      maxpt = x
    END IF
    al = al + hint
  END DO
  
!   OUTPUT RESULTS FOR EACH TEST INTERVAL
  
  rmserr = SQRT(rmserr/DBLE(npts))
  WRITE(iout, 1000)
  WRITE(iout, 1010) alow(itest), aup(itest), npts
  IF ( rmserr > zero ) THEN
    erloss = ait + LOG(rmserr) / albeta
  ELSE
    erloss = zero
  END IF
  IF ( erloss < zero ) erloss = zero
  WRITE(iout, 1020) rmserr, erloss
  WRITE(iout, 1030) maxerr, maxpt
  IF ( maxerr > zero ) THEN
    erloss = ait + LOG(maxerr) / albeta
  ELSE
    erloss = zero
  END IF
  IF ( erloss < zero ) erloss = zero
  WRITE(iout, 1040) erloss
  numeq = npts - ( numsm + numla )
  WRITE(iout, 1050) numsm, numeq, numla
END DO

!   PRINT VALUES OF CI(XMAX) = ZERO AND CI(XMIN) = LN(XMIN)

WRITE(iout, 1100)
WRITE(iout, 1200)
y = dcosint(xmax)
WRITE(iout, 1210) y
x = LOG(xmin) + gamma
WRITE(iout, 1220) x
y = dcosint(xmin)
WRITE(iout, 1230) y

CALL DATE_AND_TIME( VALUES=tend )
time_taken = 3600.*(tend(5) - tstart(5)) + 60.*(tend(6) - tstart(6)) +  &
             (tend(7) - tstart(7)) + 0.001*(tend(8) - tstart(8))
WRITE(iout, 1240) time_taken

STOP

!   FORMAT STATEMENTS

1000 FORMAT(/////t11, 'TEST OF CI(X+DELTA) AGAINST TAYLOR SERIES')
1010 FORMAT(/t11, 'TEST CARRIED OUT ON ', f9.4,' TO ', f9.4, ' WITH ',  &
            i5, ' ARGUMENTS')
1020 FORMAT(//t11, 'RMS ERROR  =  ', e15.6,' IS A LOSS OF',  &
            f8.3, ' SIGNIFICANT DIGITS')
1030 FORMAT(//t11, 'MAX ERROR  =  ', e15.6, ' OCCURRED AT X  =  ', f10.4)
1040 FORMAT(/t11, 'THIS IS A LOSS OF ', f8.3, ' SIGNIFICANT DIGITS')
1050 FORMAT(///t11, i6, ' ARGUMENTS GAVE POSITIVE ERROR', /  &
            t11, i6, ' ARGUMENTS GAVE NO ERROR', /  &
            t11, i6, ' ARGUMENTS GAVE NEGATIVE ERROR')
1100 FORMAT(/////t11, 'SPECIAL TESTS')
1200 FORMAT(/////t11, 'CALCULATE CI(XMAX) - SHOULD RETURN ZERO')
1210 FORMAT(/t11, 'VALUE OF CI(XMAX)  = ', e15.5)
1220 FORMAT(///t11, 'CALCULATE CI(XMIN) - SHOULD RETURN ', e15.5)
1230 FORMAT(/t11, 'VALUE OF CI(XMIN)  = ', e15.5)
1240 FORMAT(/t11, 'Time taken = ', f8.2, 'sec.')


CONTAINS

SUBROUTINE dterms(x, eps, mval)

!   This subroutine determines how many derivatives of Si or Ci
!   are needed at the point X, so that the Taylor series with this no.
!   of terms has a truncation error less than EPS in the
!   relative sense.  The no. of terms needed is stored in MVAL.


!   INPUT PARAMETERS:

!      X - DOUBLE PRECISION - The value at which the derivative is wanted.

!      EPS - DOUBLE PRECISION - The relative error wanted in the Taylor series.

!   OUTPUT PARAMETERS:

!      MVAL - INTEGER - The required no. of terms


!   AUTHOR:  Allan MacLeod
!            Dept. of Mathematics and Statistics
!            University of Paisley
!            High St.
!            Paisley
!            SCOTLAND
!            PA1 2BE

!            (e-mail: macl_ms0@paisley.ac.uk)


REAL (dp), INTENT(IN) :: x
REAL (dp), INTENT(IN) :: eps
INTEGER, INTENT(OUT)  :: mval

! Local variables
REAL (dp)            :: con1, con2, fmid, lower, mid, upper
REAL (dp), PARAMETER :: zero = 0.0_dp, tol = 0.1_dp, one = 1.0_dp
REAL (dp), PARAMETER :: two = 2.0_dp, four = 4.0_dp, sixten = 16.0_dp

con1 = LOG ( four * EXP(x) / eps )
con2 = LOG ( sixten * x )
upper = con1 / con2
IF ( upper <= one ) THEN
  mval = 1
  RETURN
END IF
lower = one
100 mid = ( upper + lower ) / two
IF ( ABS(mid-lower) < tol ) THEN
  mval = INT(mid) + 1
  GO TO 200
END IF
fmid = LOG(mid) + mid * con2 - con1
IF ( fmid >= zero ) THEN
  upper = mid
ELSE
  lower = mid
END IF
GO TO 100

200 RETURN
END SUBROUTINE dterms



SUBROUTINE dcider(x, nmax, ndigs, cider)

!   This subroutine calculates the necessary number of derivatives of
!   Ci by using the recurrence relation code of Gautschi and Klein
!   from Comm. ACM., vol. 13 1970 pp53-54.

!   INPUT PARAMETERS:

!      X - DOUBLE PRECISION - Value of argument

!      NMAX - INTEGER - Max. no of derivatives

!      NDIGS - INTEGER - No. of decimal digits of accuracy wanted

!   OUTPUT PARAMETERS:

!      CIDER - DOUBLE PRECISION ARRAY - The values of the derivatives


!   MODULES CALLED:

!      This subroutine calls the function DTLNTI


!   AUTHOR: Allan MacLeod
!           Dept. of Mathematics and Statistics
!           University of Paisley
!           Scotland
!            (e-mail: macl_ms0@paisley.ac.uk)

REAL (dp), INTENT(IN)  :: x
INTEGER, INTENT(IN)    :: nmax
INTEGER, INTENT(IN)    :: ndigs
REAL (dp), INTENT(OUT) :: cider(0:)

! Local variables
INTEGER              :: minn, n, nlim, n0, n1
REAL (dp)            :: d1, expone, q, sigma(4), s1, x1
REAL (dp), PARAMETER :: one = 1.0_dp, two = 2.0_dp, three = 3.0_dp,  &
                        ten = 10.0_dp

x1 = ABS(x)
sigma(1) = -SIN(x)
sigma(2) = -COS(x)
sigma(3) = -sigma(1)
sigma(4) = -sigma(2)
n0 = INT(x1)
cider(0) = sigma(4) / x
IF ( n0 < nmax ) THEN
  minn = n0
ELSE
  minn = nmax
END IF
IF ( x1 <= three ) THEN
  nlim = nmax
ELSE
  nlim = minn
END IF
DO n = 1 , nlim
  cider(n) = (sigma(n-4*((n-1)/4)) - DBLE(n)*cider(n-1))/x
END DO
IF ( x1 > three .AND. n0 < nmax ) THEN
  expone = EXP(one)
  s1 = expone * x1
  d1 = DBLE(ndigs) * LOG(ten) + LOG(two)
  n1 = INT(s1*dtlnti(d1/s1) - one)
  IF ( n1 < nmax ) n1 = nmax
  q = one / x
  DO n = 1 , n1+1
    q = -DBLE(n) * q / x
  END DO
  DO n = n1 , n0+1 , -1
    q = (sigma(n+1-4*(n/4)) - x*q) / DBLE(n+1)
    IF ( n <= nmax ) cider(n) = q
  END DO
END IF
RETURN
END SUBROUTINE dcider



FUNCTION dtlnti(y) RESULT(fn_val)

!   This function solves the equation

!          y  =  t ln(t)

!   for a given value of y. It is not intended to be very accurate.
!   N.B. Halley's method has been substituted for Newton's method.

!   INPUT PARAMETER:

!      Y - DOUBLE PRECISION - Given value of y.


!   AUTHOR:  Allan MacLeod
!            Dept. of Mathematics and Statistics
!            University of Paisley
!            High St.
!            Paisley
!            SCOTLAND
!            PA1 2BE

!            (e-mail: macl_ms0@paisley.ac.uk)


REAL (dp), INTENT(IN) :: y
REAL (dp)             :: fn_val

! Local variables
REAL (dp)            :: f0, fd, t0, t1
REAL (dp), PARAMETER :: one = 1.0_dp, tol = 0.05_dp, two = 2.0_dp

t0 = ( y + y ) / ( one + LOG(y) )
DO
  f0 = t0 * LOG(t0) - y
  fd = one + LOG(t0)
  t1 = t0 - f0 / (fd - f0 / (two*t0*fd) )
  IF ( ABS(t0-t1) > tol ) THEN
    t0 = t1
  ELSE
    EXIT
  END IF
END DO
fn_val = t1

RETURN
END FUNCTION dtlnti

END PROGRAM dcitest