PROGRAM dsitest


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

!       DPSINT(x) = integral from 0 to x of {sin(t)/t} dt

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

!   The program assumes the sine-integral is calculated in a
!   function with the name DPSINT. Thus users should append their
!   Si code at the end of this file and edit the small DPSINT
!   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, XMAX, XMIN.
!   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 this file, though the name has been changed
!   to DMACHR.

!   To ensure that DMACHR 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) DSIDER contained in this file

!      (b) DPSINT contained in the file DPSINT.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 dsinint(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 dsinint
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 == 3 ) delta = -pt0625
  DO ipt = 1 , npts
    x = al + hint * dren(nren)
    
!   PURIFY ARGUMENT SO THAT ( X + DELTA ) IS EXACT
    
    z = sixten * x
    temp = z + x
    x = temp - z
    y = x + delta
    
!   CALCULATE SI(Y) AND SI(X+DELTA)
    
    siy = dsinint(y)
    six = dsinint(x)
    CALL dsider(x, nterms, ndec, sder)
    sum = sder(nterms)
    div = DBLE(nterms) + one
    DO i = 0 , nterms-1
      sum = sum * delta / div + sder(nterms-1-i)
      div = div - one
    END DO
    sum = delta * sum
    
!   DETERMINE ERROR AND UPDATE STATS
    
    ERR = ( siy - six ) - sum
    IF ( ERR < zero ) numsm = numsm + 1
    IF ( ERR > zero ) numla = numla + 1
    relerr = ABS(ERR/siy)
    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
  erloss = ait + LOG(rmserr) / albeta
  IF(erloss < zero) erloss = zero
  WRITE(iout, 1020) rmserr, erloss
  WRITE(iout, 1030) maxerr, maxpt
  erloss = ait + LOG(maxerr) / albeta
  IF(erloss < zero) erloss = zero
  WRITE(iout, 1040) erloss
  numeq = npts - ( numsm + numla )
  WRITE(iout, 1050) numsm, numeq, numla
END DO

!   TEST SI(X) + SI(-X)  =  0.0

WRITE(iout, 1100)
WRITE(iout, 1110)
DO i = 1 , 20
  x = dren(nren) + DBLE(i-1)
  y = dsinint(x) + dsinint(-x)
  WRITE(iout, 1120) x, y
END DO

!   PRINT VALUES OF SI(XMAX) = PI/2 AND SI(XMIN) = XMIN

WRITE(iout, 1200)
y = dsinint(xmax)
WRITE(iout, 1210) y
WRITE(iout, 1220) xmin
y = dsinint(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 SI(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')
1110 FORMAT(///t11, 'TEST OF SI(X) + SI(-X)  =  0.0', //  &
            t11, 'VALUE OF X', t31, 'VALUE OF SI(X)+SI(-X)')
1120 FORMAT(/t11, f10.4, t36, d16.5)
1200 FORMAT(/////t11, 'CALCULATE SI(XMAX) - SHOULD RETURN PI/2')
1210 FORMAT(/t11, 'VALUE OF SI(XMAX)  = ', e15.5)
1220 FORMAT(///t11, 'CALCULATE SI(XMIN) - SHOULD RETURN XMIN =', e15.5)
1230 FORMAT(/t11, 'VALUE OF SI(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

REAL (dp) :: con1, con2, fmid, four, lower,  &
    mid, one, sixten, tol, two, upper, zero
DATA zero, tol, one/0.0D0, 0.1D0, 1.0D0/
DATA two, four, sixten/2.0D0, 4.0D0, 16.0D0/

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 dsider(x, nmax, ndigs, sder)

!   This subroutine calculates the necessary number of derivatives
!   of Si 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:

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


!   MODULES CALLED:

!      This subroutine calls the function DTLNTI which is contained
!      in the file SICISUBD.FOR


!   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 OUT)                  :: ndigs
REAL (dp), INTENT(OUT)            :: sder(0:nmax)
INTEGER :: n, nlim, nu, n0
REAL (dp) :: dtlnti, d1, expone, one,  sigma(4), s1, ten, two, x1, zero
DATA zero, one, two, ten/0.0D0, 1.0D0, 2.0D0, 10.0D0/

!   SET UP INITIAL PARAMETERS FOR RECURRENCE

x1 = ABS(x)
sigma(1) = COS(x)
sigma(2) = -SIN(x)
sigma(3) = -sigma(1)
sigma(4) = -sigma(2)
n0 = INT(x1)
IF ( x /= zero ) THEN
  sder(0) = sigma(4) / x
ELSE
  sder(0) = one
END IF
IF ( n0 <= nmax ) THEN
  nlim = n0
ELSE
  nlim = nmax
END IF

!  USE FORWARD RECURRENCE UP TO A CERTAIN VALUE

DO n = 1 , nlim
  sder(n) = ( sigma(n-4*((n-1)/4)) - DBLE(n) * sder(n-1)) / x
END DO

!   USE BACKWARD RECURRENCE FOR REMAINDER

IF ( n0 < nmax ) THEN
  expone = EXP(one)
  s1 = zero
  d1 = DBLE(ndigs) * LOG(ten) + LOG(two)
  IF ( nmax <= INT(expone*x1) ) THEN
    nu = 1 + INT ( expone * x1 * dtlnti(d1/(expone*x1)) )
  ELSE
    nu = 1 + INT ( DBLE(nmax) * dtlnti(d1/DBLE(nmax)) )
  END IF
  DO n = nu , n0+2 , -1
    s1 = ( sigma(n-4*((n-1)/4)) - x * s1 ) / DBLE(n)
    IF ( n <= (nmax+1) ) sder(n-1) = s1
  END DO
END IF
RETURN
END SUBROUTINE dsider



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.


!   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

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

RETURN
END FUNCTION dtlnti

END PROGRAM dsitest