#! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'READ.ME' <<'END_OF_FILE' X *************************************************************************** X * All the software contained in this library is protected by copyright. * X * Permission to use, copy, modify, and distribute this software for any * X * purpose without fee is hereby granted, provided that this entire notice * X * is included in all copies of any software which is or includes a copy * X * or modification of this software and in all copies of the supporting * X * documentation for such software. * X *************************************************************************** X * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED * X * WARRANTY. IN NO EVENT, NEITHER THE AUTHORS, NOR THE PUBLISHER, NOR ANY * X * MEMBER OF THE EDITORIAL BOARD OF THE JOURNAL "NUMERICAL ALGORITHMS", * X * NOR ITS EDITOR-IN-CHIEF, BE LIABLE FOR ANY ERROR IN THE SOFTWARE, ANY * X * MISUSE OF IT OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF * X * USING THE SOFTWARE LIES WITH THE PARTY DOING SO. * X *************************************************************************** X * ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE * X * ABOVE STATEMENT. * X *************************************************************************** X X AUTHOR: X X ALLAN J. MACLEOD X UNIVERSITY OF PAISLEY, SCOTLAND X E-MAIL: macl-ms0@paisley.ac.uk X X REFERENCE: X X - RATIONAL APPROXIMATIONS, SOFTWARE AND TEST METHODS FOR SINE AND X COSINE INTEGRALS X NUMERICAL ALGORITHMS, 12 (1996), PP. 259-272 X X SOFTWARE REVISION DATE: X X MARCH 11, 1996 X X SOFTWARE LANGUAGE: X X FORTRAN X END_OF_FILE if test 1776 -ne `wc -c <'READ.ME'`; then echo shar: \"'READ.ME'\" unpacked with wrong size! fi # end of 'READ.ME' fi if test -f 'makefile' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'makefile'\" else echo shar: Extracting \"'makefile'\" \(2956 characters\) sed "s/^X//" >'makefile' <<'END_OF_FILE' X# run either X# make citest X# or X# make sitest X# or X# make dcitest X# or X# make dsitest X# X# Next the single precision or double precision test code may be run. X# X# citest is a test code for a general cosine integral software. X# Here the new single precision function cosint is used. X# X# sitest is a test code for a general sine integral software X# Here the new single precision function sinint is used. X# X# dcitest is a test code for a general cosine integral software. X# Here the new double precision function dcosint is used. X# X# dsitest is a test code for a general sine integral software X# Here the new double precision function dsinint is used. X# X#INSTRUCTIONS FOR IMPLEMENTING TEST CODES X# X#a) For testing the author's functions (Ci and Si), simply use the X# makefile provided e.g. (make citest, make sitest,....) X# X#b) For testing the author's functions without the makefile X# * Combine (for Ci single precision) CITESTSP.F, SICISUBS.F, CINT.F and X# COSINT.F. X# * Compile, link and run. X# and similarly for the other test programs. X# X#c) For testing another function (Ci or Si) with the makefile, X# (for Ci single precision, as an example) X# * Edit CINT.F and change the name COSINT into the name of the new X# Ci function. X# * Edit the MAKEFILE and changes all the occurences of COSINT.F with X# the name of the new file provided, containing the new function Ci to X# be tested. X# * Use the commands make citest. X# and similarly for the other test programs. X# X#d) For testing another function (Ci or Si) without the makefile, X# (for Ci single precision, as an example) X# * Edit CINT.F and change the name COSINT into the name of the new X# Ci function. X# * Combine CITESTSP.F, SICISUBS.F, CINT.F and the new file containing X# the new function Ci to be tested. X# * Compile, link and run. X# and similarly for the others test programs. X# Xcitest: citestsp.o cint.o sicisubs.o cosint.o X f77 -o citest citestsp.o cint.o sicisubs.o cosint.o Xsitest: sitestsp.o sint.o sicisubs.o sinint.o X f77 -o sitest sitestsp.o sint.o sicisubs.o sinint.o Xdcitest: citestdp.o dpcint.o sicisubd.o dcosint.o X f77 -o dcitest citestdp.o dpcint.o sicisubd.o dcosint.o Xdsitest: sitestdp.o dpsint.o sicisubd.o dsinint.o X f77 -o dsitest sitestdp.o dpsint.o sicisubd.o dsinint.o Xcitestsp.o: citestsp.f X f77 -O1 -c citestsp.f Xcint.o: cint.f X f77 -O1 -c cint.f Xcosint.o: cosint.f X f77 -O1 -c cosint.f Xsitestsp.o: sitestsp.f X f77 -O1 -c sitestsp.f Xsint.o: sint.f X f77 -O1 -c sint.f Xsinint.o: sinint.f X f77 -O1 -c sinint.f Xsicisubs.o: sicisubs.f X f77 -O1 -c sicisubs.f Xcitestdp.o: citestdp.f X f77 -O1 -c citestdp.f Xdpcint.o: dpcint.f X f77 -O1 -c dpcint.f Xdcosint.o: dcosint.f X f77 -O1 -c dcosint.f Xsitestdp.o: sitestdp.f X f77 -O1 -c sitestdp.f Xdpsint.o: dpsint.f X f77 -O1 -c dpsint.f Xdsinint.o: dsinint.f X f77 -O1 -c dsinint.f Xsicisubd.o: sicisubd.f X f77 -O1 -c sicisubd.f END_OF_FILE if test 2956 -ne `wc -c <'makefile'`; then echo shar: \"'makefile'\" unpacked with wrong size! fi # end of 'makefile' fi if test -f 'cint.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'cint.f'\" else echo shar: Extracting \"'cint.f'\" \(730 characters\) sed "s/^X//" >'cint.f' <<'END_OF_FILE' X REAL FUNCTION CINT(X) XC XC This function enables the user to tell the test-code XC the particular name of the Ci function being tested. XC In this test code, REAL FUNCTION COSINT will be used. XC To test other Ci function software, the user should XC replace COSINT with the desired actual name of the code. XC Users should note that this format assumes a single XC parameter to the function. XC XC XC INPUT PARAMETER: XC XC X - REAL - The argument to the function XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X REAL COSINT, X X CINT = COSINT(X) X RETURN X END END_OF_FILE if test 730 -ne `wc -c <'cint.f'`; then echo shar: \"'cint.f'\" unpacked with wrong size! fi # end of 'cint.f' fi if test -f 'citestdp.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'citestdp.f'\" else echo shar: Extracting \"'citestdp.f'\" \(8766 characters\) sed "s/^X//" >'citestdp.f' <<'END_OF_FILE' X PROGRAM DCITEST XC XC DESCRIPTION: XC XC This program tests the performance of DOUBLE-PRECISION codes which XC calculate the cosine-integral special function, defined as XC XC DPCINT(x) = gamma + log(abs(x)) + XC integral from 0 to x of {[cos(t)-1]/t} dt XC XC The program compares the value of DPCINT(x+h) with the Taylor XC series value about x. The derivatives of DPCINT at x are XC generated by a standard recurrence relation. XC XC XC The program assumes the cosine-integral is calculated in a XC function with the name DPCINT. Thus users should append their XC Ci code at the end of this file and edit the small DPCINT XC function to include the name of their code. XC XC XC The program uses Cody's MACHAR subroutine to determine certain XC machine parameters. This can be avoided if values for the following XC parameters are available: EPS, IT, IBETA, MACHEP, XMIN, XMAX. XC Values for certain machine-compiler combinations are given in XC XC W.J. CODY Algorithm 665: MACHAR: A subroutine to dynamically XC determine machine parameters, ACM Trans. Math. Soft. 14 (1988) 303-311. XC XC The MACHAR code and a random generator code written by Cody XC are included in the file SICISUBD.FOR, though with the MACHAR XC code called DMACHR. XC XC To ensure that MACHAR and the argument purification section work XC as intended, this code should be compiled with any optimisation XC switched OFF. XC XC Users should note that MACHAR tests the machine arithmetic XC at its very limits. Thus it is possible that warning messages XC might appear about underflows and similar events. These DO XC NOT affect the test results. Such warnings can be avoided XC by explicitly declaring the necessary machine parameters. XC XC MODULES USED: XC XC This code calls the following modules: XC XC (a) DCIDER contained in this file XC XC (b) DPCINT contained in the file DPCINT.F XC XC (c) DMACHR, DREN and DTERMS contained in the file SICISUBD.F XC XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER I,IBETA,IEXP,IOUT,IPT,IRND,IT,ITEST,MACHEP,MAXEXP, X 1 MINEXP,NDEC,NEGEP,NGRD,NPTS,NREN,NTERMS,NTEST, X 2 NT1,NT2,NUMEQ,NUMLA,NUMSM X DOUBLE PRECISION ALOW(5),AUP(5),CIDER(0:100), X 2 AIT,AL,ALBETA,AU,CIX,CIY,DELTA,DIV, X 3 DPCINT,DREN,EPS,EPSNEG,ERLOSS,ERR,GAMMA, X 4 HINT,MAXERR,MAXPT,ONE,PT0625,RELERR,RMSERR, X 5 SIXTEN,SUM,TEMP,X,XMAX,XMIN,Y,Z,ZERO XC XC Define constants XC X DATA ALOW/0.25D0,2.0D0,4.0D0,8.0D0,13.5D0/ X DATA AUP/0.375D0,2.5D0,4.5D0,8.5D0,14.0D0/ X DATA PT0625,SIXTEN/0.0625D0,16.0D0/ X DATA ZERO,ONE/0.0D0,1.0D0/ X DATA GAMMA/0.57721566490153286061D0/ XC XC IOUT = output channel number XC XC NPTS = number of test points per interval XC X DATA IOUT/6/ X DATA NPTS,NTEST/5000,5/ XC XC DETERMINE MACHINE PARAMETERS XC X CALL DMACHR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, X 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) X AIT = DBLE(IT) X ALBETA = LOG(DBLE(IBETA)) X NDEC = ABS(INT(LOG10(DBLE(IBETA))*DBLE(MACHEP))) X NREN = 1234567 XC XC START TESTS XC X DO 300 ITEST = 1 , NTEST X AL = ALOW(ITEST) X AU = AUP(ITEST) X HINT = ( AU - AL ) / DBLE(NPTS) XC XC DETERMINE HOW MANY DERIVATIVES ARE NEEDED XC X CALL DTERMS(ABS(AL),EPS,NT1) X CALL DTERMS(ABS(AU),EPS,NT2) X NTERMS = NT1 X IF ( NT2 .GT. NT1 ) NTERMS = NT2 XC XC SET UP STATISTICS COLLECTION XC X MAXERR = ZERO X RMSERR = ZERO X MAXPT = ZERO X NUMSM = 0 X NUMLA = 0 X DELTA = PT0625 X IF ( ITEST .EQ. 2 .OR. ITEST .EQ. 4 ) DELTA = -PT0625 X DO 200 IPT = 1 , NPTS X X = AL + HINT * DREN(NREN) XC XC PURIFY ARGUMENT XC X Z = SIXTEN * X X TEMP = Z + X X X = TEMP - Z X Y = X + DELTA XC XC CALCULATE CI(Y) AND CI(X+DELTA) XC X CIY = DPCINT(Y) X CIX = DPCINT(X) X CALL DCIDER(X,NTERMS,NDEC,CIDER) X SUM = CIDER(NTERMS) X DIV = DBLE(NTERMS)+ONE X DO 100 I = 0 , NTERMS-1 X SUM = SUM * DELTA / DIV + CIDER(NTERMS-1-I) X DIV = DIV - ONE X 100 CONTINUE X SUM = DELTA * SUM XC XC DETERMINE ERROR AND UPDATE STATS XC X ERR = ( CIY - CIX ) - SUM X IF ( ERR .LT. ZERO ) NUMSM = NUMSM + 1 X IF ( ERR .GT. ZERO ) NUMLA = NUMLA + 1 X RELERR = ABS(ERR/CIY) X RMSERR = RMSERR + RELERR * RELERR X IF ( RELERR .GT. MAXERR ) THEN X MAXERR = RELERR X MAXPT = X X ENDIF X AL = AL + HINT X 200 CONTINUE XC XC OUTPUT RESULTS FOR EACH TEST INTERVAL XC X RMSERR = SQRT(RMSERR/DBLE(NPTS)) X WRITE(IOUT,1000) X WRITE(IOUT,1010)ALOW(ITEST),AUP(ITEST),NPTS X IF ( RMSERR .GT. ZERO ) THEN X ERLOSS = AIT + LOG(RMSERR) / ALBETA X ELSE X ERLOSS = ZERO X ENDIF X IF ( ERLOSS .LT. ZERO ) ERLOSS = ZERO X WRITE(IOUT,1020)RMSERR,ERLOSS X WRITE(IOUT,1030)MAXERR,MAXPT X IF ( MAXERR .GT. ZERO ) THEN X ERLOSS = AIT + LOG(MAXERR) / ALBETA X ELSE X ERLOSS = ZERO X ENDIF X IF ( ERLOSS .LT. ZERO ) ERLOSS = ZERO X WRITE(IOUT,1040)ERLOSS X NUMEQ = NPTS - ( NUMSM + NUMLA ) X WRITE(IOUT,1050)NUMSM,NUMEQ,NUMLA X 300 CONTINUE XC XC PRINT VALUES OF CI(XMAX) = ZERO AND CI(XMIN) = LN(XMIN) XC X WRITE(IOUT,1100) X WRITE(IOUT,1200) X Y = DPCINT(XMAX) X WRITE(IOUT,1210)Y X X = LOG(XMIN) + GAMMA X WRITE(IOUT,1220)X X Y = DPCINT(XMIN) X WRITE(IOUT,1230)Y XC XC FORMAT STATEMENTS XC X 1000 FORMAT(/////10X,'TEST OF CI(X+DELTA) AGAINST TAYLOR SERIES') X 1010 FORMAT(/10X,'TEST CARRIED OUT ON ',F9.4,' TO ',F9.4,' WITH ', X & I5,' ARGUMENTS') X 1020 FORMAT(//10X,'RMS ERROR = ',D15.6,' IS A LOSS OF', X & F8.3,' SIGNIFICANT DIGITS') X 1030 FORMAT(//10X,'MAX ERROR = ',D15.6,' OCCURRED AT X = ',F10.4) X 1040 FORMAT(/10X,'THIS IS A LOSS OF ',F8.3, X & ' SIGNIFICANT DIGITS') X 1050 FORMAT(///10X,I6,' ARGUMENTS GAVE POSITIVE ERROR',/10X, X & I6,' ARGUMENTS GAVE NO ERROR',/10X,I6, X & ' ARGUMENTS GAVE NEGATIVE ERROR') X 1100 FORMAT(/////10X,'SPECIAL TESTS') X 1200 FORMAT(/////10X,'CALCULATE CI(XMAX) - SHOULD RETURN ZERO') X 1210 FORMAT(/10X,'VALUE OF CI(XMAX) = ',D15.5) X 1220 FORMAT(///10X,'CALCULATE CI(XMIN) - SHOULD RETURN ',D15.5) X 1230 FORMAT(/10X,'VALUE OF CI(XMIN) = ',D15.5) X END X X SUBROUTINE DCIDER(X,NMAX,NDIGS,CIDER) XC XC This subroutine calculates the necessary number of XC derivatives of Ci by using the recurrence relation XC code of Gautschi and Klein from Comm. ACM., vol. 13 XC 1970 pp53-54. XC XC INPUT PARAMETERS: XC XC X - DOUBLE PRECISION - Value of argument XC XC NMAX - INTEGER - Max. no of derivatives XC XC NDIGS - INTEGER - No. of decimal digits of accuracy wanted XC XC OUTPUT PARAMETERS: XC XC CIDER - DOUBLE PRECISION ARRAY - The values of the derivatives XC XC XC MODULES CALLED: XC XC This subroutine calls the function DTLNTI which is contained XC in the file SICISUBD.FOR XC XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER MINN,N,NDIGS,NLIM,NMAX,N0,N1 X DOUBLE PRECISION CIDER(0:NMAX),DTLNTI,D1,EXPONE, X 1 ONE,Q,SIGMA(4),S1,TEN,THREE,TWO,X,X1 X DATA ONE,TWO,THREE,TEN/1.0D0,2.0D0,3.0D0,10.0D0/ X X1 = ABS(X) X SIGMA(1) = -SIN(X) X SIGMA(2) = -COS(X) X SIGMA(3) = -SIGMA(1) X SIGMA(4) = -SIGMA(2) X N0 = INT(X1) X CIDER(0) = SIGMA(4) / X X IF ( N0 .LT. NMAX ) THEN X MINN = N0 X ELSE X MINN = NMAX X ENDIF X IF ( X1 .LE. THREE ) THEN X NLIM = NMAX X ELSE X NLIM = MINN X ENDIF X DO 100 N = 1 , NLIM X CIDER(N) = (SIGMA(N-4*((N-1)/4))-DBLE(N)*CIDER(N-1))/X X 100 CONTINUE X IF ( X1 .GT. THREE .AND. N0 .LT. NMAX ) THEN X EXPONE = EXP(ONE) X S1 = EXPONE * X1 X D1 = DBLE(NDIGS) * LOG(TEN) + LOG(TWO) X N1 = INT(S1*DTLNTI(D1/S1)-ONE) X IF ( N1 .LT. NMAX ) N1 = NMAX X Q = ONE / X X DO 200 N = 1 , N1+1 X Q = -DBLE(N) * Q / X X 200 CONTINUE X DO 300 N = N1 , (N0+1) , -1 X Q = (SIGMA(N+1-4*(N/4))-X*Q)/DBLE(N+1) X IF ( N .LE. NMAX ) CIDER(N) = Q X 300 CONTINUE X ENDIF X RETURN X END X END_OF_FILE if test 8766 -ne `wc -c <'citestdp.f'`; then echo shar: \"'citestdp.f'\" unpacked with wrong size! fi # end of 'citestdp.f' fi if test -f 'citestsp.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'citestsp.f'\" else echo shar: Extracting \"'citestsp.f'\" \(8609 characters\) sed "s/^X//" >'citestsp.f' <<'END_OF_FILE' X PROGRAM CITEST XC XC DESCRIPTION: XC XC This program tests the performance of SINGLE-PRECISION codes which XC calculate the cosine-integral special function, defined as XC XC CINT(x) = gamma + log(abs(x)) + XC integral from 0 to x of {[cos(t)-1]/t} dt XC XC The program compares the value of CINT(x+h) with the Taylor XC series value about x. The derivatives of CINT at x are XC generated by a standard recurrence relation. XC XC XC The program assumes the cosine-integral is calculated in a XC function with the name CINT. Thus users should append their XC Ci code at the end of this file and edit the small CINT XC function to include the name of their code. XC XC XC The program uses Cody's MACHAR subroutine to determine certain XC machine parameters. This can be avoided if values for the following XC parameters are available: EPS, IT, IBETA, MACHEP, XMIN, XMAX. XC Values for certain machine-compiler combinations are given in XC XC W.J. CODY Algorithm 665: MACHAR: A subroutine to dynamically XC determine machine parameters, ACM Trans. Math. Soft. 14 (1988) 303-311. XC XC The MACHAR code and a random generator code written by Cody XC are included in this file. XC XC To ensure that MACHAR and the argument purification works, XC this program should be compiled with any optimisation XC switched OFF. XC XC Users should note that MACHAR tests the machine arithmetic XC at its very limits. Thus it is possible that warning messages XC might appear about underflows and similar events. These DO XC NOT affect the test results. Such warnings can be avoided XC by explicitly declaring the necessary machine parameters. XC XC MODULES USED: XC XC This code calls the following modules: XC XC (a) CIDERS contained in this file XC XC (b) CINT contained in the file CINT.F XC XC (c) MACHAR, REN and TERMEQ contained in the file SICISUBS.F XC XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER I,IBETA,IEXP,IOUT,IPT,IRND,IT,ITEST,MACHEP,MAXEXP, X 1 MINEXP,NDEC,NEGEP,NGRD,NPTS,NREN,NTERMS,NTEST, X 2 NT1,NT2,NUMEQ,NUMLA,NUMSM X REAL ALOW(5),AUP(5),CIDER(0:100), X 2 AIT,AL,ALBETA,AU,CINT,CIX,CIY,DELTA,DIV,EPS,EPSNEG, X 3 ERLOSS,ERR,GAMMA,HINT,MAXERR,MAXPT,ONE,PT0625,RELERR, X 4 REN,RMSERR,SIXTEN,SUM,TEMP,X,XMAX,XMIN,Y,Z,ZERO XC XC Define constants XC X DATA ALOW/0.25E0,2.0E0,4.0E0,8.0E0,13.5E0/ X DATA AUP/0.375E0,2.5E0,4.5E0,8.5E0,14.0E0/ X DATA PT0625,SIXTEN/0.0625E0,16.0E0/ X DATA ZERO,ONE/0.0E0,1.0E0/ X DATA GAMMA/0.57721566490153286061E0/ XC XC IOUT = output channel number XC XC NPTS = number of test points per interval XC X DATA IOUT/6/ X DATA NPTS,NTEST/5000,5/ XC XC DETERMINE MACHINE PARAMETERS XC X CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, X 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) X AIT = REAL(IT) X ALBETA = LOG(REAL(IBETA)) X NDEC = ABS(INT(LOG10(REAL(IBETA))*REAL(MACHEP))) X NREN = 1234567 XC XC START TESTS XC X DO 300 ITEST = 1 , NTEST X AL = ALOW(ITEST) X AU = AUP(ITEST) X HINT = ( AU - AL ) / REAL(NPTS) XC XC DETERMINE HOW MANY DERIVATIVES ARE NEEDED XC X CALL TERMEQ(ABS(AL),EPS,NT1) X CALL TERMEQ(ABS(AU),EPS,NT2) X NTERMS = NT1 X IF ( NT2 .GT. NT1 ) NTERMS = NT2 XC XC SET UP STATISTICS COLLECTION XC X MAXERR = ZERO X RMSERR = ZERO X MAXPT = ZERO X NUMSM = 0 X NUMLA = 0 X DELTA = PT0625 X IF ( ITEST .EQ. 2 .OR. ITEST .EQ. 4 ) DELTA = -PT0625 X DO 200 IPT = 1 , NPTS X X = AL + HINT * REN(NREN) XC XC PURIFY ARGUMENT XC X Z = SIXTEN * X X TEMP = Z + X X X = TEMP - Z X Y = X + DELTA XC XC CALCULATE CI(Y) AND CI(X+DELTA) XC X CIY = CINT(Y) X CIX = CINT(X) X CALL CIDERS(X,NTERMS,NDEC,CIDER) X SUM = CIDER(NTERMS) X DIV = REAL(NTERMS)+ONE X DO 100 I = 0 , NTERMS-1 X SUM = SUM * DELTA / DIV + CIDER(NTERMS-1-I) X DIV = DIV - ONE X 100 CONTINUE X SUM = DELTA * SUM XC XC DETERMINE ERROR AND UPDATE STATS XC X ERR = ( CIY - CIX ) - SUM X IF ( ERR .LT. ZERO ) NUMSM = NUMSM + 1 X IF ( ERR .GT. ZERO ) NUMLA = NUMLA + 1 X RELERR = ABS(ERR/CIY) X RMSERR = RMSERR + RELERR * RELERR X IF ( RELERR .GT. MAXERR ) THEN X MAXERR = RELERR X MAXPT = X X ENDIF X AL = AL + HINT X 200 CONTINUE XC XC OUTPUT RESULTS FOR EACH TEST INTERVAL XC X RMSERR = SQRT(RMSERR/REAL(NPTS)) X WRITE(IOUT,1000) X WRITE(IOUT,1010)ALOW(ITEST),AUP(ITEST),NPTS X IF ( RMSERR .GT. ZERO ) THEN X ERLOSS = AIT + LOG(RMSERR) / ALBETA X ELSE X ERLOSS = ZERO X ENDIF X IF ( ERLOSS .LT. ZERO ) ERLOSS = ZERO X WRITE(IOUT,1020)RMSERR,ERLOSS X WRITE(IOUT,1030)MAXERR,MAXPT X IF ( MAXERR .GT. ZERO ) THEN X ERLOSS = AIT + LOG(MAXERR) / ALBETA X ELSE X ERLOSS = ZERO X ENDIF X IF ( ERLOSS .LT. ZERO ) ERLOSS = ZERO X WRITE(IOUT,1040)ERLOSS X NUMEQ = NPTS - ( NUMSM + NUMLA ) X WRITE(IOUT,1050)NUMSM,NUMEQ,NUMLA X 300 CONTINUE XC XC PRINT VALUES OF CI(XMAX) = ZERO AND CI(XMIN) = LN(XMIN) XC X WRITE(IOUT,1100) X WRITE(IOUT,1200) X Y = CINT(XMAX) X WRITE(IOUT,1210)Y X X = LOG(XMIN) + GAMMA X WRITE(IOUT,1220)X X Y = CINT(XMIN) X WRITE(IOUT,1230)Y XC XC FORMAT STATEMENTS XC X 1000 FORMAT(/////10X,'TEST OF CI(X+DELTA) AGAINST TAYLOR SERIES') X 1010 FORMAT(/10X,'TEST CARRIED OUT ON ',F9.4,' TO ',F9.4,' WITH ', X & I5,' ARGUMENTS') X 1020 FORMAT(//10X,'RMS ERROR = ',E15.6,' IS A LOSS OF', X & F8.3,' SIGNIFICANT DIGITS') X 1030 FORMAT(//10X,'MAX ERROR = ',E15.6,' OCCURRED AT X = ',F10.4) X 1040 FORMAT(/10X,'THIS IS A LOSS OF ',F8.3, X & ' SIGNIFICANT DIGITS') X 1050 FORMAT(///10X,I6,' ARGUMENTS GAVE POSITIVE ERROR',/10X, X & I6,' ARGUMENTS GAVE NO ERROR',/10X,I6, X & ' ARGUMENTS GAVE NEGATIVE ERROR') X 1100 FORMAT(/////10X,'SPECIAL TESTS') X 1200 FORMAT(/////10X,'CALCULATE CINT(XMAX) - SHOULD RETURN ZERO') X 1210 FORMAT(/10X,'VALUE OF CINT(XMAX) = ',E15.5) X 1220 FORMAT(///10X,'CALCULATE CINT(XMIN) - SHOULD RETURN ',E15.5) X 1230 FORMAT(/10X,'VALUE OF CINT(XMIN) = ',E15.5) X END X X SUBROUTINE CIDERS(X,NMAX,NDIGS,CIDER) XC XC This subroutine calculates the necessary number of XC derivatives of Ci by using the recurrence relation XC code of Gautschi and Klein from Comm. ACM., vol. 13 XC 1970 pp53-54. XC XC INPUT PARAMETERS: XC XC X - REAL - Value of argument XC XC NMAX - INTEGER - Max. no of derivatives XC XC NDIGS - INTEGER - No. of decimal digits of accuracy wanted XC XC OUTPUT PARAMETERS: XC XC CIDER - REAL ARRAY - The values of the derivatives XC XC XC MODULES CALLED: XC XC This subroutine calls the function TLNTIN which is contained XC in the file SICISUBS.FOR XC XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER MINN,N,NDIGS,NLIM,NMAX,N0,N1 X REAL CIDER(0:NMAX),REAL,D1,EXPONE,ONE,Q,SIGMA(4), X 1 S1,TEN,THREE,TLNTIN,TWO,X,X1 X DATA ONE,TWO,THREE,TEN/1.0E0,2.0E0,3.0E0,10.0E0/ X X1 = ABS(X) X SIGMA(1) = -SIN(X) X SIGMA(2) = -COS(X) X SIGMA(3) = -SIGMA(1) X SIGMA(4) = -SIGMA(2) X N0 = INT(X1) X CIDER(0) = SIGMA(4) / X X IF ( N0 .LT. NMAX ) THEN X MINN = N0 X ELSE X MINN = NMAX X ENDIF X IF ( X1 .LE. THREE ) THEN X NLIM = NMAX X ELSE X NLIM = MINN X ENDIF X DO 100 N = 1 , NLIM X CIDER(N) = (SIGMA(N-4*((N-1)/4))-REAL(N)*CIDER(N-1))/X X 100 CONTINUE X IF ( X1 .GT. THREE .AND. N0 .LT. NMAX ) THEN X EXPONE = EXP(ONE) X S1 = EXPONE * X1 X D1 = REAL(NDIGS) * LOG(TEN) + LOG(TWO) X N1 = INT(S1*TLNTIN(D1/S1)-ONE) X IF ( N1 .LT. NMAX ) N1 = NMAX X Q = ONE / X X DO 200 N = 1 , N1+1 X Q = -REAL(N) * Q / X X 200 CONTINUE X DO 300 N = N1 , (N0+1) , -1 X Q = (SIGMA(N+1-4*(N/4))-X*Q)/REAL(N+1) X IF ( N .LE. NMAX ) CIDER(N) = Q X 300 CONTINUE X ENDIF X RETURN X END X END_OF_FILE if test 8609 -ne `wc -c <'citestsp.f'`; then echo shar: \"'citestsp.f'\" unpacked with wrong size! fi # end of 'citestsp.f' fi if test -f 'cosint.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'cosint.f'\" else echo shar: Extracting \"'cosint.f'\" \(9988 characters\) sed "s/^X//" >'cosint.f' <<'END_OF_FILE' X REAL FUNCTION COSINT(XVALUE) XC XC XC This program calculates the value of the cosine-integral XC defined as XC XC COSINT(x) = Gamma + Ln(x) + Integral (0 to x) [cos(t)-1]/t dt XC XC where Gamma is Euler' constant. XC XC The code uses rational approximations with a maximum XC accuracy of 20sf. XC XC INPUT PARAMETER: XC XC XVALUE - REAL - The argument to the function XC XC XC MACHINE-DEPENDENT PARAMETERS: XC XC XLOW - REAL - The absolute value below which XC COSINT( x ) = gamma + LN(x) , to machine precision. XC The recommended value is SQRT(2*EPSNEG) XC XC XHIGH1 - REAL - The value above which XC COSINT(x) = sin(x)/x - cos(x)/x^2 XC to machine precision. XC The recommended value is SQRT(6/EPSNEG) XC XC XHIGH2 - REAL - The value above which the trig. functions cannot be XC accurately determined. The value of the function is XC COSINT(x) = 0.0 XC The recommended value is pi/EPS. XC XC Values of EPS and EPSNEG for certain machine/compiler XC combinations can be found in the paper XC XC W.J. CODY Algorithm 665: MACHAR: A subroutine to dynamically XC determine machine parameters, ACM Trans. Math. Soft. 14 (1988) 303-311. XC XC The current code gives numerical values for XLOW,XHIGH1,XHIGH2 XC suitable for machines whose arithmetic conforms to the IEEE XC standard. The codes will probably work on other machines XC but might overflow or underflow for certain arguments. XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER I X REAL AFD1(0:7),AFD2(0:7),AFN1(0:7),AFN2(0:7),AGD1(0:8), X 2 AGD2(0:8),AGN1(0:8),AGN2(0:8),AC1D(0:5),AC1N(0:5), X 3 AC2D(0:6),AC2N(0:7),CX,DIF,FIVAL,GAMMA,GIVAL,LOGL1, X 4 LOGL2,LOGVAL,ONE,P(0:2),Q(0:1),ROOT,RT1D, X 5 RT1N,RT1R,RT2D,RT2N,RT2R,SIX,SUM,SUMDEN,SUMNUM,SX, X 6 THREE,TWELVE,X,XHIGH1,XHIGH2,XLOW,XSQ,XVALUE,ZERO XC XC DATA VALUES XC X DATA ZERO,ONE,THREE,SIX,TWELVE/0.0E0,1.0E0,3.0E0,6.0E0,12.0E0/ X DATA GAMMA/0.57721566490153286060E0/ X DATA LOGL1,LOGL2/0.46875E-1,0.25E0/ XC XC MACHINE-DEPENDENT PARAMETERS (SUITABLE FOR IEEE MACHINES) XC X DATA XLOW,XHIGH1,XHIGH2/1.0358E-3,10033.11E0,2.635359E7/ XC XC VALUES FOR COS-INTEGRAL FOR 0 < X <= 3 XC X DATA AC1N/-0.24607411378767540707E0,0.72113492241301534559E-2, X 1 -0.11867127836204767056E-3,0.90542655466969866243E-6, X 2 -0.34322242412444409037E-8,0.51950683460656886834E-11/ X DATA AC1D/1.0E0,0.12670095552700637845E-1, X 1 0.78168450570724148921E-4,0.29959200177005821677E-6, X 2 0.73191677761328838216E-9,0.94351174530907529061E-12/ XC XC VALUES FOR COS-INTEGRAL FOR 3 < X <= 6 XC X DATA AC2N/-0.15684781827145408780E0,0.66253165609605468916E-2, X 1 -0.12822297297864512864E-3,0.12360964097729408891E-5, X 2 -0.66450975112876224532E-8,0.20326936466803159446E-10, X 3 -0.33590883135343844613E-13,0.23686934961435015119E-16/ X DATA AC2D/1.0E0,0.96166044388828741188E-2, X 1 0.45257514591257035006E-4,0.13544922659627723233E-6, X 2 0.27715365686570002081E-9,0.37718676301688932926E-12, X 3 0.27706844497155995398E-15/ XC XC VALUES FOR FI(X) FOR 6 <= X <= 12 XC X DATA AFN1/0.99999999962173909991E0,0.36451060338631902917E3, X 1 0.44218548041288440874E5,0.22467569405961151887E7, X 2 0.49315316723035561922E8,0.43186795279670283193E9, X 3 0.11847992519956804350E10,0.45573267593795103181E9/ X DATA AFD1/1.0E0,0.36651060273229347594E3, X 1 0.44927569814970692777E5,0.23285354882204041700E7, X 2 0.53117852017228262911E8,0.50335310667241870372E9, X 3 0.16575285015623175410E10,0.11746532837038341076E10/ XC XC VALUES OF GI(X) FOR 6 <= X <=12 XC X DATA AGN1/0.99999999920484901956E0,0.51385504875307321394E3, X 1 0.92293483452013810811E5,0.74071341863359841727E7, X 2 0.28142356162841356551E9,0.49280890357734623984E10, X 3 0.35524762685554302472E11,0.79194271662085049376E11, X 4 0.17942522624413898907E11/ X DATA AGD1/1.0E0,0.51985504708814870209E3, X 1 0.95292615508125947321E5,0.79215459679762667578E7, X 2 0.31977567790733781460E9,0.62273134702439012114E10, X 3 0.54570971054996441467E11,0.18241750166645704670E12, X 4 0.15407148148861454434E12/ XC XC VALUES FOR FI(X) FOR X > 12 XC X DATA AFN2/0.19999999999999978257E1,0.22206119380434958727E4, X 1 0.84749007623988236808E6,0.13959267954823943232E9, X 2 0.10197205463267975592E11,0.30229865264524075951E12, X 3 0.27504053804288471142E13,0.21818989704686874983E13/ X DATA AFD2/1.0E0,0.11223059690217167788E4, X 1 0.43685270974851313242E6,0.74654702140658116258E8, X 2 0.58580034751805687471E10,0.20157980379272098841E12, X 3 0.26229141857684496445E13,0.87852907334918467516E13/ XC XC VALUES FOR GI(X) FOR X > 12 XC X DATA AGN2/0.59999999999999993089E1,0.96527746044997139158E4, X 1 0.56077626996568834185E7,0.15022667718927317198E10, X 2 0.19644271064733088465E12,0.12191368281163225043E14, X 3 0.31924389898645609533E15,0.25876053010027485934E16, X 4 0.12754978896268878403E16/ X DATA AGD2/1.0E0,0.16287957674166143196E4, X 1 0.96636303195787870963E6,0.26839734750950667021E9, X 2 0.37388510548029219241E11,0.26028585666152144496E13, X 3 0.85134283716950697226E14,0.11304079361627952930E16, X 4 0.42519841479489798424E16/ XC XC VALUES FOR AN APPROXIMATION TO LN(X/ROOT) XC X DATA P/0.83930008362695945726E1,-0.65306663899493304675E1, X 1 0.569155722227490223E0/ X DATA Q/0.41965004181347972847E1,-0.46641666676862479585E1/ XC XC VALUES OF THE FIRST TWO ROOTS OF THE COSINE-INTEGRAL XC X DATA RT1N,RT1D,RT1R/631.0E0,1024.0E0,0.29454812071623379711E-3/ X DATA RT2N,RT2D,RT2R/3465.0E0,1024.0E0,0.39136005118642639785E-3/ XC XC START COMPUTATION XC X X = XVALUE X IF ( X .LE. ZERO ) THEN X COSINT = ZERO X RETURN X ENDIF X IF ( X .LE. SIX ) THEN XC XC CODE FOR 3 < X < = 6 XC X IF ( X .GT. THREE ) THEN X SUMNUM = ZERO X SUMDEN = ZERO X XSQ = X * X X DO 100 I = 7 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AC2N( I ) X 100 CONTINUE X DO 120 I = 6 , 0 , -1 X SUMDEN = SUMDEN * XSQ + AC2D( I ) X 120 CONTINUE X ROOT = RT2N / RT2D X DIF = ( X - ROOT ) - RT2R X SUM = ROOT + RT2R X IF ( ABS(DIF) .LT. LOGL2 ) THEN X CX = DIF / ( SUM + X ) X XSQ = CX * CX X SX = P(0) + XSQ * ( P(1) + XSQ * P(2) ) X SX = SX / ( Q(0) + XSQ * ( Q(1) + XSQ ) ) X LOGVAL = CX * SX X ELSE X LOGVAL = LOG( X / SUM ) X ENDIF X COSINT = LOGVAL + DIF * ( X + SUM ) * SUMNUM / SUMDEN X ELSE XC XC CODE FOR 0 < X < = 3 XC X IF ( X .GT. XLOW ) THEN X SUMNUM = ZERO X SUMDEN = ZERO X XSQ = X * X X DO 150 I = 5 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AC1N( I ) X SUMDEN = SUMDEN * XSQ + AC1D( I ) X 150 CONTINUE X ROOT = RT1N / RT1D X DIF = ( X - ROOT ) - RT1R X SUM = ROOT + RT1R X IF ( ABS(DIF) .LT. LOGL1 ) THEN X CX = DIF / ( SUM + X ) X XSQ = CX * CX X SX = P(0) + XSQ * ( P(1) + XSQ * P(2) ) X SX = SX / ( Q(0) + XSQ * ( Q(1) + XSQ ) ) X LOGVAL = CX * SX X ELSE X LOGVAL = LOG( X / SUM ) X ENDIF X COSINT = LOGVAL + DIF * ( X + SUM ) * SUMNUM / SUMDEN X ELSE X COSINT = GAMMA + LOG( X ) X ENDIF X ENDIF X ENDIF XC XC CODE FOR 6 < X < = 12 XC X IF ( X .GT. SIX .AND. X .LE. TWELVE ) THEN X SUMNUM = ZERO X SUMDEN = ZERO X XSQ = ONE / ( X * X ) X DO 200 I = 7 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AFN1( I ) X SUMDEN = SUMDEN * XSQ + AFD1( I ) X 200 CONTINUE X FIVAL = SUMNUM / ( X * SUMDEN ) X SUMNUM = ZERO X SUMDEN = ZERO X DO 300 I = 8 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AGN1( I ) X SUMDEN = SUMDEN * XSQ + AGD1( I ) X 300 CONTINUE X GIVAL = XSQ * SUMNUM / SUMDEN X COSINT = FIVAL * SIN( X ) - GIVAL * COS( X ) X ENDIF XC XC CODE FOR X > 12 XC X IF ( X .GT. TWELVE ) THEN X IF ( X .GT. XHIGH2 ) THEN X COSINT = ZERO X ELSE X CX = COS( X ) X SX = SIN( X ) X XSQ = ONE / ( X * X ) X IF ( X .GT. XHIGH1 ) THEN X COSINT = SX / X - CX * XSQ X ELSE X SUMNUM = ZERO X SUMDEN = ZERO X DO 400 I = 7 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AFN2( I ) X SUMDEN = SUMDEN * XSQ + AFD2( I ) X 400 CONTINUE X FIVAL = ( ONE - XSQ * SUMNUM / SUMDEN ) / X X SUMNUM = ZERO X SUMDEN = ZERO X DO 500 I = 8 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AGN2( I ) X SUMDEN = SUMDEN * XSQ + AGD2( I ) X 500 CONTINUE X GIVAL = ( ONE - XSQ * SUMNUM / SUMDEN ) * XSQ X COSINT = SX * FIVAL - CX * GIVAL X ENDIF X ENDIF X ENDIF X RETURN X END X END_OF_FILE if test 9988 -ne `wc -c <'cosint.f'`; then echo shar: \"'cosint.f'\" unpacked with wrong size! fi # end of 'cosint.f' fi if test -f 'dcosint.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dcosint.f'\" else echo shar: Extracting \"'dcosint.f'\" \(10240 characters\) sed "s/^X//" >'dcosint.f' <<'END_OF_FILE' X DOUBLE PRECISION FUNCTION DCOSINT(XVALUE) XC XC XC This program calculates the value of the cosine-integral XC defined as XC XC DCOSINT(x) = Gamma + Ln(x) + Integral (0 to x) [cos(t)-1]/t dt XC XC where Gamma is Euler's constant. XC XC The code uses rational approximations with a XC maximum accuracy of 20sf. XC XC XC INPUT PARAMETER: XC XC XVALUE - DOUBLE PRECISION - The argument to the function XC XC XC MACHINE-DEPENDENT PARAMETERS: XC XC XLOW - DOUBLE PRECISION - The absolute value below which XC DCOSINT( x ) = gamma + LN(x) , XC to machine precision. XC The recommended value is SQRT(2*EPSNEG) XC XC XHIGH1 - DOUBLE PRECISION - The value above which XC DCOSINT(x) = sin(x)/x - cos(x)/x^2 XC to machine precision. XC The recommended value is SQRT(6/EPSNEG) XC XC XHIGH2 - DOUBLE PRECISION - The value above which the trig. functions XC cannot be accurately determined. XC The value of the function is XC DCOSINT(x) = 0.0 XC The recommended value is pi/EPS. XC XC Values of EPS and EPSNEG for certain machine/compiler XC combinations can be found in the paper XC XC W.J. CODY Algorithm 665: MACHAR: A subroutine to dynamically XC determine machine parameters, ACM Trans. Math. Soft. 14 (1988) 303-311. XC XC The current code gives numerical values for XLOW,XHIGH1,XHIGH2 XC suitable for machines whose arithmetic conforms to the IEEE XC standard. The codes will probably work on other machines XC but might overflow or underflow for certain arguments. XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER I X DOUBLE PRECISION AFD1(0:7),AFD2(0:7),AFN1(0:7),AFN2(0:7), X 2 AGD1(0:8),AGD2(0:8),AGN1(0:8),AGN2(0:8),AC1D(0:5),AC1N(0:5), X 3 AC2D(0:6),AC2N(0:7),CX,DIF,FIVAL,GAMMA,GIVAL,LOGL1, X 4 LOGL2,LOGVAL,ONE,P(0:2),Q(0:1),ROOT,RT1D, X 5 RT1N,RT1R,RT2D,RT2N,RT2R,SIX,SUM,SUMDEN,SUMNUM,SX, X 6 THREE,TWELVE,X,XHIGH1,XHIGH2,XLOW,XSQ,XVALUE,ZERO XC XC DATA VALUES XC X DATA ZERO,ONE,THREE,SIX,TWELVE/0.0D0,1.0D0,3.0D0,6.0D0,12.0D0/ X DATA GAMMA/0.57721566490153286060D0/ X DATA LOGL1,LOGL2/0.46875D-1,0.25D0/ XC XC MACHINE-DEPENDENT PARAMETERS (SUITABLE FOR IEEE MACHINES) XC X DATA XLOW,XHIGH1,XHIGH2/1.48996D-8,2.324953D8,1.4148475D16/ XC XC VALUES FOR COS-INTEGRAL FOR 0 < X <= 3 XC X DATA AC1N/-0.24607411378767540707D0,0.72113492241301534559D-2, X 1 -0.11867127836204767056D-3,0.90542655466969866243D-6, X 2 -0.34322242412444409037D-8,0.51950683460656886834D-11/ X DATA AC1D/1.0D0,0.12670095552700637845D-1, X 1 0.78168450570724148921D-4,0.29959200177005821677D-6, X 2 0.73191677761328838216D-9,0.94351174530907529061D-12/ XC XC VALUES FOR COS-INTEGRAL FOR 3 < X <= 6 XC X DATA AC2N/-0.15684781827145408780D0,0.66253165609605468916D-2, X 1 -0.12822297297864512864D-3,0.12360964097729408891D-5, X 2 -0.66450975112876224532D-8,0.20326936466803159446D-10, X 3 -0.33590883135343844613D-13,0.23686934961435015119D-16/ X DATA AC2D/1.0D0,0.96166044388828741188D-2, X 1 0.45257514591257035006D-4,0.13544922659627723233D-6, X 2 0.27715365686570002081D-9,0.37718676301688932926D-12, X 3 0.27706844497155995398D-15/ XC XC VALUES FOR FI(X) FOR 6 <= X <= 12 XC X DATA AFN1/0.99999999962173909991D0,0.36451060338631902917D3, X 1 0.44218548041288440874D5,0.22467569405961151887D7, X 2 0.49315316723035561922D8,0.43186795279670283193D9, X 3 0.11847992519956804350D10,0.45573267593795103181D9/ X DATA AFD1/1.0D0,0.36651060273229347594D3, X 1 0.44927569814970692777D5,0.23285354882204041700D7, X 2 0.53117852017228262911D8,0.50335310667241870372D9, X 3 0.16575285015623175410D10,0.11746532837038341076D10/ XC XC VALUES OF GI(X) FOR 6 <= X <=12 XC X DATA AGN1/0.99999999920484901956D0,0.51385504875307321394D3, X 1 0.92293483452013810811D5,0.74071341863359841727D7, X 2 0.28142356162841356551D9,0.49280890357734623984D10, X 3 0.35524762685554302472D11,0.79194271662085049376D11, X 4 0.17942522624413898907D11/ X DATA AGD1/1.0D0,0.51985504708814870209D3, X 1 0.95292615508125947321D5,0.79215459679762667578D7, X 2 0.31977567790733781460D9,0.62273134702439012114D10, X 3 0.54570971054996441467D11,0.18241750166645704670D12, X 4 0.15407148148861454434D12/ XC XC VALUES FOR FI(X) FOR X > 12 XC X DATA AFN2/0.19999999999999978257D1,0.22206119380434958727D4, X 1 0.84749007623988236808D6,0.13959267954823943232D9, X 2 0.10197205463267975592D11,0.30229865264524075951D12, X 3 0.27504053804288471142D13,0.21818989704686874983D13/ X DATA AFD2/1.0D0,0.11223059690217167788D4, X 1 0.43685270974851313242D6,0.74654702140658116258D8, X 2 0.58580034751805687471D10,0.20157980379272098841D12, X 3 0.26229141857684496445D13,0.87852907334918467516D13/ XC XC VALUES FOR GI(X) FOR X > 12 XC X DATA AGN2/0.59999999999999993089D1,0.96527746044997139158D4, X 1 0.56077626996568834185D7,0.15022667718927317198D10, X 2 0.19644271064733088465D12,0.12191368281163225043D14, X 3 0.31924389898645609533D15,0.25876053010027485934D16, X 4 0.12754978896268878403D16/ X DATA AGD2/1.0D0,0.16287957674166143196D4, X 1 0.96636303195787870963D6,0.26839734750950667021D9, X 2 0.37388510548029219241D11,0.26028585666152144496D13, X 3 0.85134283716950697226D14,0.11304079361627952930D16, X 4 0.42519841479489798424D16/ XC XC VALUES FOR AN APPROXIMATION TO LN(X/ROOT) XC X DATA P/0.83930008362695945726D1,-0.65306663899493304675D1, X 1 0.569155722227490223D0/ X DATA Q/0.41965004181347972847D1,-0.46641666676862479585D1/ XC XC VALUES OF THE FIRST TWO ROOTS OF THE COSINE-INTEGRAL XC X DATA RT1N,RT1D,RT1R/631.0D0,1024.0D0,0.29454812071623379711D-3/ X DATA RT2N,RT2D,RT2R/3465.0D0,1024.0D0,0.39136005118642639785D-3/ XC XC START COMPUTATION XC X X = XVALUE X IF ( X .LE. ZERO ) THEN X DCOSINT = ZERO X RETURN X ENDIF X IF ( X .LE. SIX ) THEN XC XC CODE FOR 3 < X < = 6 XC X IF ( X .GT. THREE ) THEN X SUMNUM = ZERO X SUMDEN = ZERO X XSQ = X * X X DO 100 I = 7 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AC2N( I ) X 100 CONTINUE X DO 120 I = 6 , 0 , -1 X SUMDEN = SUMDEN * XSQ + AC2D( I ) X 120 CONTINUE X ROOT = RT2N / RT2D X DIF = ( X - ROOT ) - RT2R X SUM = ROOT + RT2R X IF ( ABS(DIF) .LT. LOGL2 ) THEN X CX = DIF / ( SUM + X ) X XSQ = CX * CX X SX = P(0) + XSQ * ( P(1) + XSQ * P(2) ) X SX = SX / ( Q(0) + XSQ * ( Q(1) + XSQ ) ) X LOGVAL = CX * SX X ELSE X LOGVAL = LOG( X / SUM ) X ENDIF X DCOSINT = LOGVAL + DIF * ( X + SUM ) * SUMNUM / SUMDEN X ELSE XC XC CODE FOR 0 < X < = 3 XC X IF ( X .GT. XLOW ) THEN X SUMNUM = ZERO X SUMDEN = ZERO X XSQ = X * X X DO 150 I = 5 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AC1N( I ) X SUMDEN = SUMDEN * XSQ + AC1D( I ) X 150 CONTINUE X ROOT = RT1N / RT1D X DIF = ( X - ROOT ) - RT1R X SUM = ROOT + RT1R X IF ( ABS(DIF) .LT. LOGL1 ) THEN X CX = DIF / ( SUM + X ) X XSQ = CX * CX X SX = P(0) + XSQ * ( P(1) + XSQ * P(2) ) X SX = SX / ( Q(0) + XSQ * ( Q(1) + XSQ ) ) X LOGVAL = CX * SX X ELSE X LOGVAL = LOG( X / SUM ) X ENDIF X DCOSINT = LOGVAL + DIF * ( X + SUM ) * SUMNUM / SUMDEN X ELSE X DCOSINT = GAMMA + LOG( X ) X ENDIF X ENDIF X ENDIF XC XC CODE FOR 6 < X < = 12 XC X IF ( X .GT. SIX .AND. X .LE. TWELVE ) THEN X SUMNUM = ZERO X SUMDEN = ZERO X XSQ = ONE / ( X * X ) X DO 200 I = 7 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AFN1( I ) X SUMDEN = SUMDEN * XSQ + AFD1( I ) X 200 CONTINUE X FIVAL = SUMNUM / ( X * SUMDEN ) X SUMNUM = ZERO X SUMDEN = ZERO X DO 300 I = 8 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AGN1( I ) X SUMDEN = SUMDEN * XSQ + AGD1( I ) X 300 CONTINUE X GIVAL = XSQ * SUMNUM / SUMDEN X DCOSINT = FIVAL * SIN( X ) - GIVAL * COS( X ) X ENDIF XC XC CODE FOR X > 12 XC X IF ( X .GT. TWELVE ) THEN X IF ( X .GT. XHIGH2 ) THEN X DCOSINT = ZERO X ELSE X CX = COS( X ) X SX = SIN( X ) X XSQ = ONE / ( X * X ) X IF ( X .GT. XHIGH1 ) THEN X DCOSINT = SX / X - CX * XSQ X ELSE X SUMNUM = ZERO X SUMDEN = ZERO X DO 400 I = 7 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AFN2( I ) X SUMDEN = SUMDEN * XSQ + AFD2( I ) X 400 CONTINUE X FIVAL = ( ONE - XSQ * SUMNUM / SUMDEN ) / X X SUMNUM = ZERO X SUMDEN = ZERO X DO 500 I = 8 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AGN2( I ) X SUMDEN = SUMDEN * XSQ + AGD2( I ) X 500 CONTINUE X GIVAL = ( ONE - XSQ * SUMNUM / SUMDEN ) * XSQ X DCOSINT = SX * FIVAL - CX * GIVAL X ENDIF X ENDIF X ENDIF X RETURN X END END_OF_FILE if test 10240 -ne `wc -c <'dcosint.f'`; then echo shar: \"'dcosint.f'\" unpacked with wrong size! fi # end of 'dcosint.f' fi if test -f 'dpcint.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dpcint.f'\" else echo shar: Extracting \"'dpcint.f'\" \(791 characters\) sed "s/^X//" >'dpcint.f' <<'END_OF_FILE' X DOUBLE PRECISION FUNCTION DPCINT(X) XC XC This function enables the user to tell the test-code XC the particular name of the Ci function being tested. XC In this test code, DOUBLE PRECISION FUNCTION DCOSINT XC will be used. XC To test other Ci function software, the user should XC replace DCOSINT with the desired actual name of the code. XC Users should note that this format assumes a single XC parameter to the function. XC XC XC INPUT PARAMETER: XC XC X - DOUBLE PRECISION - The argument to the function XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X DOUBLE PRECISION DCOSINT, X X DPCINT = DCOSINT(X) X RETURN X END END_OF_FILE if test 791 -ne `wc -c <'dpcint.f'`; then echo shar: \"'dpcint.f'\" unpacked with wrong size! fi # end of 'dpcint.f' fi if test -f 'dpsint.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dpsint.f'\" else echo shar: Extracting \"'dpsint.f'\" \(790 characters\) sed "s/^X//" >'dpsint.f' <<'END_OF_FILE' X DOUBLE PRECISION FUNCTION DPSINT(X) XC XC This function enables the user to tell the test-code XC the particular name of the Si function being tested. XC In this test code, DOUBLE PRECISION FUNCTION DSININT XC will be used. XC To test other Si function software, the user should XC replace DSININT with the desired actual name of the code. XC Users should note that this format assumes a single XC parameter to the function. XC XC XC INPUT PARAMETER: XC XC X - DOUBLE PRECISION - The argument to the function XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X DOUBLE PRECISION DSININT,X X DPSINT = DSININT(X) X RETURN X END END_OF_FILE if test 790 -ne `wc -c <'dpsint.f'`; then echo shar: \"'dpsint.f'\" unpacked with wrong size! fi # end of 'dpsint.f' fi if test -f 'dsinint.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dsinint.f'\" else echo shar: Extracting \"'dsinint.f'\" \(8088 characters\) sed "s/^X//" >'dsinint.f' <<'END_OF_FILE' X DOUBLE PRECISION FUNCTION DSININT(XVALUE) XC XC DEFINITION: XC XC This program calculates the value of the sine-integral XC defined as XC XC DSININT(x) = Integral (0 to x) sin(t)/t dt XC XC The program uses rational approximations with the coefficients XC given to 20sf. accuracy. XC XC XC INPUT PARAMETER: XC XC XVALUE - DOUBLE PRECISION - The argument to the function XC XC XC MACHINE-DEPENDENT PARAMETERS: XC XC XLOW - DOUBLE PRECISION - The absolute value below which XC DSININT( x ) = x, XC to machine precision. XC The recommended value is XC SQRT(18*EPSNEG) XC XC XHIGH1 - DOUBLE PRECISION - The value above which XC DSININT(x) = pi/2 - cos(x)/x -sin(x)/x*x XC to machine precision. XC The recommended value is XC SQRT(6/EPSNEG) XC XC XHIGH2 - DOUBLE PRECISION - The value above which XC DSININT(x) = pi/2 XC to machine precision. XC The recommended value is XC 2 / min(EPS,EPSNEG) XC XC XHIGH3 - DOUBLE PRECISION - The value above which it is not sensible XC to compute COS or SIN. The recommended XC value is pi/EPS XC XC XC Values of EPS and EPSNEG for certain machine/compiler XC combinations can be found in the paper XC XC W.J. CODY Algorithm 665: MACHAR: A subroutine to dynamically XC determine machine parameters, ACM Trans. Math. Soft. 14 (1988) 303-311. XC XC The current code gives numerical values for XLOW,XHIGH1,XHIGH2,XHIGH3, XC suitable for machines whose arithmetic conforms to the IEEE XC standard. The codes will probably work on other machines XC but might overflow or underflow for certain arguments. XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER I,INDSGN X DOUBLE PRECISION AFD1(0:7),AFD2(0:7),AFN1(0:7), X 1 AFN2(0:7),AGD1(0:8),AGD2(0:8),AGN1(0:8), X 2 AGN2(0:8),ASINTD(0:7),ASINTN(0:7), X 3 CX,FIVAL,GIVAL,ONE,PIBY2,SIX,SUMDEN,SUMNUM,SX,TWELVE, X 4 X,XHIGH,XHIGH1,XHIGH2,XHIGH3,XLOW,XSQ,XVALUE,ZERO XC XC DATA VALUES XC X DATA ZERO,ONE,SIX,TWELVE/0.0D0,1.0D0,6.0D0,12.0D0/ X DATA PIBY2/1.5707963267948966192D0/ XC XC MACHINE-DEPENDENT PARAMETERS (SUITABLE FOR IEEE MACHINES) XC X DATA XLOW,XHIGH1/4.47D-8,2.32472D8/ X DATA XHIGH2,XHIGH3/9.0072D15,1.4148475D16/ XC XC VALUES FOR SINE-INTEGRAL FOR 0 <= !X! <= 6 XC X DATA ASINTN/1.0D0,-0.44663998931312457298D-1, X 1 0.11209146443112369449D-2,-0.13276124407928422367D-4, X 2 0.85118014179823463879D-7,-0.29989314303147656479D-9, X 3 0.55401971660186204711D-12,-0.42406353433133212926D-15/ X DATA ASINTD/1.0D0,0.10891556624243098264D-1, X 1 0.59334456769186835896D-4,0.21231112954641805908D-6, X 2 0.54747121846510390750D-9,0.10378561511331814674D-11, X 3 0.13754880327250272679D-14,0.10223981202236205703D-17/ XC XC VALUES FOR FI(X) FOR 6 <= X <= 12 XC X DATA AFN1/0.99999999962173909991D0,0.36451060338631902917D3, X 1 0.44218548041288440874D5,0.22467569405961151887D7, X 2 0.49315316723035561922D8,0.43186795279670283193D9, X 3 0.11847992519956804350D10,0.45573267593795103181D9/ X DATA AFD1/1.0D0,0.36651060273229347594D3, X 1 0.44927569814970692777D5,0.23285354882204041700D7, X 2 0.53117852017228262911D8,0.50335310667241870372D9, X 3 0.16575285015623175410D10,0.11746532837038341076D10/ XC XC VALUES OF GI(X) FOR 6 <= X <=12 XC X DATA AGN1/0.99999999920484901956D0,0.51385504875307321394D3, X 1 0.92293483452013810811D5,0.74071341863359841727D7, X 2 0.28142356162841356551D9,0.49280890357734623984D10, X 3 0.35524762685554302472D11,0.79194271662085049376D11, X 4 0.17942522624413898907D11/ X DATA AGD1/1.0D0,0.51985504708814870209D3, X 1 0.95292615508125947321D5,0.79215459679762667578D7, X 2 0.31977567790733781460D9,0.62273134702439012114D10, X 3 0.54570971054996441467D11,0.18241750166645704670D12, X 4 0.15407148148861454434D12/ XC XC VALUES FOR FI(X) FOR X > 12 XC X DATA AFN2/0.19999999999999978257D1,0.22206119380434958727D4, X 1 0.84749007623988236808D6,0.13959267954823943232D9, X 2 0.10197205463267975592D11,0.30229865264524075951D12, X 3 0.27504053804288471142D13,0.21818989704686874983D13/ X DATA AFD2/1.0D0,0.11223059690217167788D4, X 1 0.43685270974851313242D6,0.74654702140658116258D8, X 2 0.58580034751805687471D10,0.20157980379272098841D12, X 3 0.26229141857684496445D13,0.87852907334918467516D13/ XC XC VALUES FOR GI(X) FOR X > 12 XC X DATA AGN2/0.59999999999999993089D1,0.96527746044997139158D4, X 1 0.56077626996568834185D7,0.15022667718927317198D10, X 2 0.19644271064733088465D12,0.12191368281163225043D14, X 3 0.31924389898645609533D15,0.25876053010027485934D16, X 4 0.12754978896268878403D16/ X DATA AGD2/1.0D0,0.16287957674166143196D4, X 1 0.96636303195787870963D6,0.26839734750950667021D9, X 2 0.37388510548029219241D11,0.26028585666152144496D13, X 3 0.85134283716950697226D14,0.11304079361627952930D16, X 4 0.42519841479489798424D16/ XC XC START COMPUTATION XC X X = XVALUE X INDSGN = 1 X IF ( X .LT. ZERO ) THEN X X = -X X INDSGN = -1 X ENDIF XC XC CODE FOR 0 <= |X| <= 6 XC X IF ( X .LE. SIX ) THEN X IF ( X .LT. XLOW ) THEN X DSININT = X X ELSE X SUMNUM = ZERO X SUMDEN = ZERO X XSQ = X * X X DO 100 I = 7 , 0 , -1 X SUMNUM = SUMNUM * XSQ + ASINTN(I) X SUMDEN = SUMDEN * XSQ + ASINTD(I) X 100 CONTINUE X DSININT = X * SUMNUM / SUMDEN X ENDIF X ENDIF XC XC CODE FOR 6 < |X| <= 12 XC X IF ( X .GT. SIX .AND. X .LE. TWELVE ) THEN X SUMNUM = ZERO X SUMDEN = ZERO X XSQ = ONE / ( X * X ) X DO 200 I = 7 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AFN1(I) X SUMDEN = SUMDEN * XSQ + AFD1(I) X 200 CONTINUE X FIVAL = SUMNUM / ( X * SUMDEN ) X SUMNUM = ZERO X SUMDEN = ZERO X DO 300 I = 8 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AGN1(I) X SUMDEN = SUMDEN * XSQ + AGD1(I) X 300 CONTINUE X GIVAL = XSQ * SUMNUM / SUMDEN X DSININT = PIBY2 - FIVAL * COS(X) - GIVAL * SIN(X) X ENDIF XC XC CODE FOR |X| > 12 XC X IF ( X .GT. TWELVE ) THEN X XHIGH = MIN(XHIGH2,XHIGH3) X IF ( X .GT. XHIGH ) THEN X DSININT = PIBY2 X ELSE X CX = COS(X) X SX = SIN(X) X XSQ = ONE / ( X * X ) X IF ( X .GT. XHIGH1 ) THEN X DSININT = PIBY2 - CX / X - SX * XSQ X ELSE X SUMNUM = ZERO X SUMDEN = ZERO X DO 400 I = 7 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AFN2(I) X SUMDEN = SUMDEN * XSQ + AFD2(I) X 400 CONTINUE X FIVAL = ( ONE - XSQ * SUMNUM / SUMDEN ) / X X SUMNUM = ZERO X SUMDEN = ZERO X DO 500 I = 8 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AGN2(I) X SUMDEN = SUMDEN * XSQ + AGD2(I) X 500 CONTINUE X GIVAL = ( ONE - XSQ * SUMNUM / SUMDEN ) * XSQ X DSININT = PIBY2 - CX * FIVAL - SX * GIVAL X ENDIF X ENDIF X ENDIF X IF ( INDSGN .EQ. -1 ) DSININT = -DSININT X RETURN X END X END_OF_FILE if test 8088 -ne `wc -c <'dsinint.f'`; then echo shar: \"'dsinint.f'\" unpacked with wrong size! fi # end of 'dsinint.f' fi if test -f 'sicisubd.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'sicisubd.f'\" else echo shar: Extracting \"'sicisubd.f'\" \(13731 characters\) sed "s/^X//" >'sicisubd.f' <<'END_OF_FILE' X SUBROUTINE DTERMS(X,EPS,MVAL) XC XC This subroutine determines how many derivatives of Si XC or Ci are needed at the point X, so that the Taylor XC series with this no. of terms has a truncation error XC less than EPS in the relative sense. The no. of terms XC needed is stored in MVAL. XC XC XC INPUT PARAMETERS: XC XC X - DOUBLE PRECISION - The value at which the derivative is wanted. XC XC EPS - DOUBLE PRECISION - The relative error wanted in the Taylor series. XC XC OUTPUT PARAMETERS: XC XC MVAL - INTEGER - The required no. of terms XC XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC High St. XC Paisley XC SCOTLAND XC PA1 2BE XC XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER MVAL X DOUBLE PRECISION CON1,CON2,EPS,FMID,FOUR,LOWER, X 2 MID,ONE,SIXTEN,TOL,TWO,UPPER,X,ZERO X DATA ZERO,TOL,ONE/0.0D0,0.1D0,1.0D0/ X DATA TWO,FOUR,SIXTEN/2.0D0,4.0D0,16.0D0/ X CON1 = LOG ( FOUR * EXP(X) / EPS ) X CON2 = LOG ( SIXTEN * X ) X UPPER = CON1 / CON2 X IF ( UPPER .LE. ONE ) THEN X MVAL = 1 X RETURN X ENDIF X LOWER = ONE X 100 MID = LOWER + ( UPPER - LOWER ) / TWO X IF ( ABS(MID-LOWER) .LT. TOL ) THEN X MVAL = INT(MID) + 1 X GOTO 200 X ENDIF X FMID = LOG(MID) + MID * CON2 - CON1 X IF ( FMID .GE. ZERO ) THEN X UPPER = MID X ELSE X LOWER = MID X ENDIF X GOTO 100 X 200 RETURN X END X X DOUBLE PRECISION FUNCTION DTLNTI(Y) XC XC This function solves the equation XC XC y = t ln(t) XC XC for a given value of y. It is not intended to be very accurate. XC XC XC INPUT PARAMETER: XC XC Y - DOUBLE PRECISION - Given value of y. XC XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC High St. XC Paisley XC SCOTLAND XC PA1 2BE XC XC (e-mail: macl_ms0@paisley.ac.uk) XC X DOUBLE PRECISION F0,FD,ONE,TOL,T0,T1,Y X DATA ONE,TOL/1.0D0,0.05D0/ X T0 = ( Y + Y ) / ( ONE + LOG(Y) ) X 100 F0 = T0 * LOG(T0) - Y X FD = ONE + LOG(T0) X T1 = T0 - F0 / FD X IF ( ABS(T0-T1) .GT. TOL ) THEN X T0 = T1 X GOTO 100 X ENDIF X DTLNTI = T1 X RETURN X END X X SUBROUTINE DMACHR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, X 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) XC---------------------------------------------------------------------- XC This Fortran 77 subroutine is intended to determine the parameters XC of the floating-point arithmetic system specified below. The XC determination of the first three uses an extension of an algorithm XC due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some, XC but not all, of the improvements suggested by M. Gentleman and S. XC Marovich, CACM 17 (1974), pp. 276-277. An earlier version of this XC program was published in the book Software Manual for the XC Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall, XC Englewood Cliffs, NJ, 1980. XC XC XC Parameter values reported are as follows: XC XC IBETA - the radix for the floating-point representation XC IT - the number of base IBETA digits in the floating-point XC significand XC IRND - 0 if floating-point addition chops XC 1 if floating-point addition rounds, but not in the XC IEEE style XC 2 if floating-point addition rounds in the IEEE style XC 3 if floating-point addition chops, and there is XC partial underflow XC 4 if floating-point addition rounds, but not in the XC IEEE style, and there is partial underflow XC 5 if floating-point addition rounds in the IEEE style, XC and there is partial underflow XC NGRD - the number of guard digits for multiplication with XC truncating arithmetic. It is XC 0 if floating-point arithmetic rounds, or if it XC truncates and only IT base IBETA digits XC participate in the post-normalization shift of the XC floating-point significand in multiplication; XC 1 if floating-point arithmetic truncates and more XC than IT base IBETA digits participate in the XC post-normalization shift of the floating-point XC significand in multiplication. XC MACHEP - the largest negative integer such that XC 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that XC MACHEP is bounded below by -(IT+3) XC NEGEPS - the largest negative integer such that XC 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that XC NEGEPS is bounded below by -(IT+3) XC IEXP - the number of bits (decimal places if IBETA = 10) XC reserved for the representation of the exponent XC (including the bias or sign) of a floating-point XC number XC MINEXP - the largest in magnitude negative integer such that XC DBLE(IBETA)**MINEXP is positive and normalized XC MAXEXP - the smallest positive power of BETA that overflows XC EPS - DBLE(IBETA)**MACHEP. XC EPSNEG - DBLE(IBETA)**NEGEPS. XC XMIN - the smallest non-vanishing normalized floating-point XC power of the radix, i.e., XMIN = DBLE(IBETA)**MINEXP XC XMAX - the largest finite floating-point number. In XC particular XMAX = (1.0-EPSNEG)*DBLE(IBETA)**MAXEXP XC Note - on some machines XMAX will be only the XC second, or perhaps third, largest number, being XC too small by 1 or 2 units in the last digit of XC the significand. XC XC Latest modification: May 30, 1989 XC XC Author: W. J. Cody XC Mathematics and Computer Science Division XC Argonne National Laboratory XC Argonne, IL 60439 XC XC---------------------------------------------------------------------- X INTEGER I,IBETA,IEXP,IRND,IT,ITEMP,IZ,J,K,MACHEP,MAXEXP, X 1 MINEXP,MX,NEGEP,NGRD,NXRES X DOUBLE PRECISION A,B,BETA,BETAIN,BETAH,EPS,EPSNEG,ONE,T, X 2 TEMP,TEMPA,TEMP1,TWO,XMAX,XMIN,Y,Z,ZERO XC---------------------------------------------------------------------- X ONE = DBLE(1) X TWO = ONE + ONE X ZERO = ONE - ONE XC---------------------------------------------------------------------- XC Determine IBETA, BETA ala Malcolm. XC---------------------------------------------------------------------- X A = ONE X 10 A = A + A X TEMP = A+ONE X TEMP1 = TEMP-A X IF (TEMP1-ONE .EQ. ZERO) GO TO 10 X B = ONE X 20 B = B + B X TEMP = A+B X ITEMP = INT(TEMP-A) X IF (ITEMP .EQ. 0) GO TO 20 X IBETA = ITEMP X BETA = DBLE(IBETA) XC---------------------------------------------------------------------- XC Determine IT, IRND. XC---------------------------------------------------------------------- X IT = 0 X B = ONE X 100 IT = IT + 1 X B = B * BETA X TEMP = B+ONE X TEMP1 = TEMP-B X IF (TEMP1-ONE .EQ. ZERO) GO TO 100 X IRND = 0 X BETAH = BETA / TWO X TEMP = A+BETAH X IF (TEMP-A .NE. ZERO) IRND = 1 X TEMPA = A + BETA X TEMP = TEMPA+BETAH X IF ((IRND .EQ. 0) .AND. (TEMP-TEMPA .NE. ZERO)) IRND = 2 XC---------------------------------------------------------------------- XC Determine NEGEP, EPSNEG. XC---------------------------------------------------------------------- X NEGEP = IT + 3 X BETAIN = ONE / BETA X A = ONE X DO 200 I = 1, NEGEP X A = A * BETAIN X 200 CONTINUE X B = A X 210 TEMP = ONE-A X IF (TEMP-ONE .NE. ZERO) GO TO 220 X A = A * BETA X NEGEP = NEGEP - 1 X GO TO 210 X 220 NEGEP = -NEGEP X EPSNEG = A XC---------------------------------------------------------------------- XC Determine MACHEP, EPS. XC---------------------------------------------------------------------- X MACHEP = -IT - 3 X A = B X 300 TEMP = ONE+A X IF (TEMP-ONE .NE. ZERO) GO TO 320 X A = A * BETA X MACHEP = MACHEP + 1 X GO TO 300 X 320 EPS = A XC---------------------------------------------------------------------- XC Determine NGRD. XC---------------------------------------------------------------------- X NGRD = 0 X TEMP = ONE+EPS X IF ((IRND .EQ. 0) .AND. (TEMP*ONE-ONE .NE. ZERO)) NGRD = 1 XC---------------------------------------------------------------------- XC Determine IEXP, MINEXP, XMIN. XC XC Loop to determine largest I and K = 2**I such that XC (1/BETA) ** (2**(I)) XC does not underflow. XC Exit from loop is signaled by an underflow. XC---------------------------------------------------------------------- X I = 0 X K = 1 X Z = BETAIN X T = ONE + EPS X NXRES = 0 X 400 Y = Z X Z = Y * Y XC---------------------------------------------------------------------- XC Check for underflow here. XC---------------------------------------------------------------------- X A = Z * ONE X TEMP = Z * T X IF ((A+A .EQ. ZERO) .OR. (ABS(Z) .GE. Y)) GO TO 410 X TEMP1 = TEMP * BETAIN X IF (TEMP1*BETA .EQ. Z) GO TO 410 X I = I + 1 X K = K + K X GO TO 400 X 410 IF (IBETA .EQ. 10) GO TO 420 X IEXP = I + 1 X MX = K + K X GO TO 450 XC---------------------------------------------------------------------- XC This segment is for decimal machines only. XC---------------------------------------------------------------------- X 420 IEXP = 2 X IZ = IBETA X 430 IF (K .LT. IZ) GO TO 440 X IZ = IZ * IBETA X IEXP = IEXP + 1 X GO TO 430 X 440 MX = IZ + IZ - 1 XC---------------------------------------------------------------------- XC Loop to determine MINEXP, XMIN. XC Exit from loop is signaled by an underflow. XC---------------------------------------------------------------------- X 450 XMIN = Y X Y = Y * BETAIN XC---------------------------------------------------------------------- XC Check for underflow here. XC---------------------------------------------------------------------- X A = Y * ONE X TEMP = Y * T X IF (((A+A) .EQ. ZERO) .OR. (ABS(Y) .GE. XMIN)) GO TO 460 X K = K + 1 X TEMP1 = TEMP * BETAIN X IF ((TEMP1*BETA .NE. Y) .OR. (TEMP .EQ. Y)) THEN X GO TO 450 X ELSE X NXRES = 3 X XMIN = Y X END IF X 460 MINEXP = -K XC---------------------------------------------------------------------- XC Determine MAXEXP, XMAX. XC---------------------------------------------------------------------- X IF ((MX .GT. K+K-3) .OR. (IBETA .EQ. 10)) GO TO 500 X MX = MX + MX X IEXP = IEXP + 1 X 500 MAXEXP = MX + MINEXP XC---------------------------------------------------------------------- XC Adjust IRND to reflect partial underflow. XC---------------------------------------------------------------------- X IRND = IRND + NXRES XC---------------------------------------------------------------------- XC Adjust for IEEE-style machines. XC---------------------------------------------------------------------- X IF (IRND .GE. 2) MAXEXP = MAXEXP - 2 XC---------------------------------------------------------------------- XC Adjust for machines with implicit leading bit in binary XC significand, and machines with radix point at extreme XC right of significand. XC---------------------------------------------------------------------- X I = MAXEXP + MINEXP X IF ((IBETA .EQ. 2) .AND. (I .EQ. 0)) MAXEXP = MAXEXP - 1 X IF (I .GT. 20) MAXEXP = MAXEXP - 1 X IF (A .NE. Y) MAXEXP = MAXEXP - 2 X XMAX = ONE - EPSNEG X IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG X XMAX = XMAX / (BETA * BETA * BETA * XMIN) X I = MAXEXP + MINEXP + 3 X IF (I .LE. 0) GO TO 520 X DO 510 J = 1, I X IF (IBETA .EQ. 2) XMAX = XMAX + XMAX X IF (IBETA .NE. 2) XMAX = XMAX * BETA X 510 CONTINUE X 520 RETURN XC---------- Last line of DMACHR ---------- X END X X DOUBLE PRECISION FUNCTION DREN(K) XC--------------------------------------------------------------------- XC Random number generator - based on Algorithm 266 by Pike and XC Hill (modified by Hansson), Communications of the ACM, XC Vol. 8, No. 10, October 1965. XC XC This subprogram is intended for use on computers with XC fixed point wordlength of at least 29 bits. It is XC best if the floating-point significand has at most XC 29 bits. XC XC Latest modification: May 30, 1989 XC XC Author: W. J. Cody XC Mathematics and Computer Science Division XC Argonne National Laboratory XC Argonne, IL 60439 XC XC--------------------------------------------------------------------- X INTEGER IY,J,K X DOUBLE PRECISION C1,C2,C3,ONE X DATA IY/100001/ X DATA ONE,C1,C2,C3/1.0D0,2796203.0D0,1.0D-6,1.0D-12/ XC--------------------------------------------------------------------- XC Statement functions for conversion between integer and float XC--------------------------------------------------------------------- X J = K X IY = IY * 125 X IY = IY - (IY/2796203) * 2796203 X DREN = DBLE(IY) / C1 * (ONE + C2 + C3) X RETURN XC---------- Last card of DREN ---------- X END END_OF_FILE if test 13731 -ne `wc -c <'sicisubd.f'`; then echo shar: \"'sicisubd.f'\" unpacked with wrong size! fi # end of 'sicisubd.f' fi if test -f 'sicisubs.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'sicisubs.f'\" else echo shar: Extracting \"'sicisubs.f'\" \(13627 characters\) sed "s/^X//" >'sicisubs.f' <<'END_OF_FILE' X SUBROUTINE TERMEQ(X,EPS,MVAL) XC XC This subroutine determines how many derivatives of Si XC or Ci are needed at the point X, so that the Taylor XC series with this no. of terms has a truncation error XC less than EPS in the relative sense. The no. of terms XC needed is stored in MVAL. XC XC XC INPUT PARAMETERS: XC XC X - REAL - The value at which the derivative is wanted. XC XC EPS - REAL - The relative error wanted in the Taylor series. XC XC OUTPUT PARAMETERS: XC XC MVAL - INTEGER - The required no. of terms XC XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC High St. XC Paisley XC SCOTLAND XC PA1 2BE XC XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER MVAL X REAL CON1,CON2,EPS,FMID,FOUR,LOWER, X 2 MID,ONE,SIXTEN,TOL,TWO,UPPER,X,ZERO X DATA ZERO,TOL,ONE/0.0E0,0.1E0,1.0E0/ X DATA TWO,FOUR,SIXTEN/2.0E0,4.0E0,16.0E0/ X CON1 = LOG ( FOUR * EXP(X) / EPS ) X CON2 = LOG ( SIXTEN * X ) X UPPER = CON1 / CON2 X IF ( UPPER .LE. ONE ) THEN X MVAL = 1 X RETURN X ENDIF X LOWER = ONE X 100 MID = LOWER + ( UPPER - LOWER ) / TWO X IF ( ABS(MID-LOWER) .LT. TOL ) THEN X MVAL = INT(MID) + 1 X GOTO 200 X ENDIF X FMID = LOG(MID) + MID * CON2 - CON1 X IF ( FMID .GE. ZERO ) THEN X UPPER = MID X ELSE X LOWER = MID X ENDIF X GOTO 100 X 200 RETURN X END X X REAL FUNCTION TLNTIN(Y) XC XC This function solves the equation XC XC y = t ln(t) XC XC for a given value of y. It is not intended to be very accurate. XC XC XC INPUT PARAMETER: XC XC Y - REAL - Given value of y. XC XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC High St. XC Paisley XC SCOTLAND XC PA1 2BE XC XC (e-mail: macl_ms0@paisley.ac.uk) XC X X X REAL F0,FD,ONE,TOL,T0,T1,Y X DATA ONE,TOL/1.0E0,0.05E0/ X T0 = ( Y + Y ) / ( ONE + LOG(Y) ) X 100 F0 = T0 * LOG(T0) - Y X FD = ONE + LOG(T0) X T1 = T0 - F0 / FD X IF ( ABS(T0-T1) .GT. TOL ) THEN X T0 = T1 X GOTO 100 X ENDIF X TLNTIN = T1 X RETURN X END X X SUBROUTINE MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, X 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) XC---------------------------------------------------------------------- XC This Fortran 77 subroutine is intended to determine the parameters XC of the floating-point arithmetic system specified below. The XC determination of the first three uses an extension of an algorithm XC due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some, XC but not all, of the improvements suggested by M. Gentleman and S. XC Marovich, CACM 17 (1974), pp. 276-277. An earlier version of this XC program was published in the book Software Manual for the XC Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall, XC Englewood Cliffs, NJ, 1980. XC XC Parameter values reported are as follows: XC XC IBETA - the radix for the floating-point representation XC IT - the number of base IBETA digits in the floating-point XC significand XC IRND - 0 if floating-point addition chops XC 1 if floating-point addition rounds, but not in the XC IEEE style XC 2 if floating-point addition rounds in the IEEE style XC 3 if floating-point addition chops, and there is XC partial underflow XC 4 if floating-point addition rounds, but not in the XC IEEE style, and there is partial underflow XC 5 if floating-point addition rounds in the IEEE style, XC and there is partial underflow XC NGRD - the number of guard digits for multiplication with XC truncating arithmetic. It is XC 0 if floating-point arithmetic rounds, or if it XC truncates and only IT base IBETA digits XC participate in the post-normalization shift of the XC floating-point significand in multiplication; XC 1 if floating-point arithmetic truncates and more XC than IT base IBETA digits participate in the XC post-normalization shift of the floating-point XC significand in multiplication. XC MACHEP - the largest negative integer such that XC 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that XC MACHEP is bounded below by -(IT+3) XC NEGEPS - the largest negative integer such that XC 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that XC NEGEPS is bounded below by -(IT+3) XC IEXP - the number of bits (decimal places if IBETA = 10) XC reserved for the representation of the exponent XC (including the bias or sign) of a floating-point XC number XC MINEXP - the largest in magnitude negative integer such that XC FLOAT(IBETA)**MINEXP is positive and normalized XC MAXEXP - the smallest positive power of BETA that overflows XC EPS - FLOAT(IBETA)**MACHEP. XC EPSNEG - FLOAT(IBETA)**NEGEPS. XC XMIN - the smallest non-vanishing normalized floating-point XC power of the radix, i.e., XMIN = FLOAT(IBETA)**MINEXP XC XMAX - the largest finite floating-point number. In XC particular XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP XC Note - on some machines XMAX will be only the XC second, or perhaps third, largest number, being XC too small by 1 or 2 units in the last digit of XC the significand. XC XC Latest modification: May 30, 1989 XC XC Author: W. J. Cody XC Mathematics and Computer Science Division XC Argonne National Laboratory XC Argonne, IL 60439 XC XC---------------------------------------------------------------------- X INTEGER I,IBETA,IEXP,IRND,IT,ITEMP,IZ,J,K,MACHEP,MAXEXP, X 1 MINEXP,MX,NEGEP,NGRD,NXRES X REAL A,B,BETA,BETAIN,BETAH,EPS,EPSNEG,ONE,T,TEMP,TEMPA, X 2 TEMP1,TWO,XMAX,XMIN,Y,Z,ZERO XC---------------------------------------------------------------------- X ONE = REAL(1) X TWO = ONE + ONE X ZERO = ONE - ONE XC---------------------------------------------------------------------- XC Determine IBETA, BETA ala Malcolm. XC---------------------------------------------------------------------- X A = ONE X 10 A = A + A X TEMP = A+ONE X TEMP1 = TEMP-A X IF (TEMP1-ONE .EQ. ZERO) GO TO 10 X B = ONE X 20 B = B + B X TEMP = A+B X ITEMP = INT(TEMP-A) X IF (ITEMP .EQ. 0) GO TO 20 X IBETA = ITEMP X BETA = REAL(IBETA) XC---------------------------------------------------------------------- XC Determine IT, IRND. XC---------------------------------------------------------------------- X IT = 0 X B = ONE X 100 IT = IT + 1 X B = B * BETA X TEMP = B+ONE X TEMP1 = TEMP-B X IF (TEMP1-ONE .EQ. ZERO) GO TO 100 X IRND = 0 X BETAH = BETA / TWO X TEMP = A+BETAH X IF (TEMP-A .NE. ZERO) IRND = 1 X TEMPA = A + BETA X TEMP = TEMPA+BETAH X IF ((IRND .EQ. 0) .AND. (TEMP-TEMPA .NE. ZERO)) IRND = 2 XC---------------------------------------------------------------------- XC Determine NEGEP, EPSNEG. XC---------------------------------------------------------------------- X NEGEP = IT + 3 X BETAIN = ONE / BETA X A = ONE X DO 200 I = 1, NEGEP X A = A * BETAIN X 200 CONTINUE X B = A X 210 TEMP = ONE-A X IF (TEMP-ONE .NE. ZERO) GO TO 220 X A = A * BETA X NEGEP = NEGEP - 1 X GO TO 210 X 220 NEGEP = -NEGEP X EPSNEG = A XC---------------------------------------------------------------------- XC Determine MACHEP, EPS. XC---------------------------------------------------------------------- X MACHEP = -IT - 3 X A = B X 300 TEMP = ONE+A X IF (TEMP-ONE .NE. ZERO) GO TO 320 X A = A * BETA X MACHEP = MACHEP + 1 X GO TO 300 X 320 EPS = A XC---------------------------------------------------------------------- XC Determine NGRD. XC---------------------------------------------------------------------- X NGRD = 0 X TEMP = ONE+EPS X IF ((IRND .EQ. 0) .AND. (TEMP*ONE-ONE .NE. ZERO)) NGRD = 1 XC---------------------------------------------------------------------- XC Determine IEXP, MINEXP, XMIN. XC XC Loop to determine largest I and K = 2**I such that XC (1/BETA) ** (2**(I)) XC does not underflow. XC Exit from loop is signaled by an underflow. XC---------------------------------------------------------------------- X I = 0 X K = 1 X Z = BETAIN X T = ONE + EPS X NXRES = 0 X 400 Y = Z X Z = Y * Y XC---------------------------------------------------------------------- XC Check for underflow here. XC---------------------------------------------------------------------- X A = Z * ONE X TEMP = Z * T X IF ((A+A .EQ. ZERO) .OR. (ABS(Z) .GE. Y)) GO TO 410 X TEMP1 = TEMP * BETAIN X IF (TEMP1*BETA .EQ. Z) GO TO 410 X I = I + 1 X K = K + K X GO TO 400 X 410 IF (IBETA .EQ. 10) GO TO 420 X IEXP = I + 1 X MX = K + K X GO TO 450 XC---------------------------------------------------------------------- XC This segment is for decimal machines only. XC---------------------------------------------------------------------- X 420 IEXP = 2 X IZ = IBETA X 430 IF (K .LT. IZ) GO TO 440 X IZ = IZ * IBETA X IEXP = IEXP + 1 X GO TO 430 X 440 MX = IZ + IZ - 1 XC---------------------------------------------------------------------- XC Loop to determine MINEXP, XMIN. XC Exit from loop is signaled by an underflow. XC---------------------------------------------------------------------- X 450 XMIN = Y X Y = Y * BETAIN XC---------------------------------------------------------------------- XC Check for underflow here. XC---------------------------------------------------------------------- X A = Y * ONE X TEMP = Y * T X IF (((A+A) .EQ. ZERO) .OR. (ABS(Y) .GE. XMIN)) GO TO 460 X K = K + 1 X TEMP1 = TEMP * BETAIN X IF ((TEMP1*BETA .NE. Y) .OR. (TEMP .EQ. Y)) THEN X GO TO 450 X ELSE X NXRES = 3 X XMIN = Y X END IF X 460 MINEXP = -K XC---------------------------------------------------------------------- XC Determine MAXEXP, XMAX. XC---------------------------------------------------------------------- X IF ((MX .GT. K+K-3) .OR. (IBETA .EQ. 10)) GO TO 500 X MX = MX + MX X IEXP = IEXP + 1 X 500 MAXEXP = MX + MINEXP XC---------------------------------------------------------------------- XC Adjust IRND to reflect partial underflow. XC---------------------------------------------------------------------- X IRND = IRND + NXRES XC---------------------------------------------------------------------- XC Adjust for IEEE-style machines. XC---------------------------------------------------------------------- X IF (IRND .GE. 2) MAXEXP = MAXEXP - 2 XC---------------------------------------------------------------------- XC Adjust for machines with implicit leading bit in binary XC significand, and machines with radix point at extreme XC right of significand. XC---------------------------------------------------------------------- X I = MAXEXP + MINEXP X IF ((IBETA .EQ. 2) .AND. (I .EQ. 0)) MAXEXP = MAXEXP - 1 X IF (I .GT. 20) MAXEXP = MAXEXP - 1 X IF (A .NE. Y) MAXEXP = MAXEXP - 2 X XMAX = ONE - EPSNEG X IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG X XMAX = XMAX / (BETA * BETA * BETA * XMIN) X I = MAXEXP + MINEXP + 3 X IF (I .LE. 0) GO TO 520 X DO 510 J = 1, I X IF (IBETA .EQ. 2) XMAX = XMAX + XMAX X IF (IBETA .NE. 2) XMAX = XMAX * BETA X 510 CONTINUE X 520 RETURN XC---------- Last line of MACHAR ---------- X END X X REAL FUNCTION REN(K) XC--------------------------------------------------------------------- XC Random number generator - based on Algorithm 266 by Pike and XC Hill (modified by Hansson), Communications of the ACM, XC Vol. 8, No. 10, October 1965. XC XC This subprogram is intended for use on computers with XC fixed point wordlength of at least 29 bits. It is XC best if the floating-point significand has at most XC 29 bits. XC XC Latest modification: May 30, 1989 XC XC Author: W. J. Cody XC Mathematics and Computer Science Division XC Argonne National Laboratory XC Argonne, IL 60439 XC XC--------------------------------------------------------------------- X INTEGER IY,J,K X REAL C1,C2,C3,ONE X DATA IY/100001/ X DATA ONE,C1,C2,C3/1.0E0,2796203.0E0,1.0E-6,1.0E-12/ XC--------------------------------------------------------------------- XC Statement functions for conversion between integer and float XC--------------------------------------------------------------------- X J = K X IY = IY * 125 X IY = IY - (IY/2796203) * 2796203 X REN = REAL(IY) / C1 * (ONE + C2 + C3) X RETURN XC---------- Last card of REN ---------- X END X X X END_OF_FILE if test 13627 -ne `wc -c <'sicisubs.f'`; then echo shar: \"'sicisubs.f'\" unpacked with wrong size! fi # end of 'sicisubs.f' fi if test -f 'sinint.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'sinint.f'\" else echo shar: Extracting \"'sinint.f'\" \(7771 characters\) sed "s/^X//" >'sinint.f' <<'END_OF_FILE' X REAL FUNCTION SININT(XVALUE) XC XC DEFINITION: XC XC This program calculates the value of the sine-integral XC defined as XC XC SININT(x) = Integral (0 to x) sin(t)/t dt XC XC The program uses rational approximations with the coefficients XC given to 20sf. accuracy. XC XC XC INPUT PARAMETER: XC XC XVALUE - REAL - The argument to the function XC XC XC MACHINE-DEPENDENT PARAMETERS: XC XC XLOW - REAL - The absolute value below which XC SININT( x ) = x, to machine precision. XC The recommended value is XC SQRT(18*EPSNEG) XC XC XHIGH1 - REAL - The value above which XC SININT(x) = pi/2 - cos(x)/x -sin(x)/x*x XC to machine precision. XC The recommended value is XC SQRT(6/EPSNEG) XC XC XHIGH2 - REAL - The value above which XC SININT(x) = pi/2 to machine precision. XC The recommended value is XC 2 / min(EPS,EPSNEG) XC XC XHIGH3 - REAL - The value above which it is not sensible XC to compute COS or SIN. The recommended value is XC pi/EPS XC XC XC Values of EPS and EPSNEG for certain machine/compiler XC combinations can be found in the paper XC XC W.J. CODY Algorithm 665: MACHAR: A subroutine to dynamically XC determine machine parameters, ACM Trans. Math. Soft. 14 (1988) 303-311. XC XC The current code gives numerical values for XLOW,XHIGH1,XHIGH2,XHIGH3, XC suitable for machines whose arithmetic conforms to the IEEE XC standard. The codes will probably work on other machines XC but might overflow or underflow for certain arguments. XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER I,INDSGN X REAL AFD1(0:7),AFD2(0:7),AFN1(0:7),AFN2(0:7),AGD1(0:8), X 2 AGD2(0:8),AGN1(0:8),AGN2(0:8),ASINTD(0:7),ASINTN(0:7), X 3 CX,FIVAL,GIVAL,ONE,PIBY2,SIX,SUMDEN,SUMNUM,SX,TWELVE, X 4 X,XHIGH,XHIGH1,XHIGH2,XHIGH3,XLOW,XSQ,XVALUE,ZERO XC XC DATA VALUES XC X DATA ZERO,ONE,SIX,TWELVE/0.0E0,1.0E0,6.0E0,12.0E0/ X DATA PIBY2/1.5707963267948966192E0/ XC XC MACHINE-DEPENDENT PARAMETERS (SUITABLE FOR IEEE MACHINES) XC X DATA XLOW,XHIGH1/1.0358E-3,10033.11E0/ X DATA XHIGH2,XHIGH3/1.67722E7,2.635359E7/ XC XC VALUES FOR SINE-INTEGRAL FOR 0 <= |X| <= 6 XC X DATA ASINTN/1.0E0,-0.44663998931312457298E-1, X 1 0.11209146443112369449E-2,-0.13276124407928422367E-4, X 2 0.85118014179823463879E-7,-0.29989314303147656479E-9, X 3 0.55401971660186204711E-12,-0.42406353433133212926E-15/ X DATA ASINTD/1.0E0,0.10891556624243098264E-1, X 1 0.59334456769186835896E-4,0.21231112954641805908E-6, X 2 0.54747121846510390750E-9,0.10378561511331814674E-11, X 3 0.13754880327250272679E-14,0.10223981202236205703E-17/ XC XC VALUES FOR FI(X) FOR 6 <= X <= 12 XC X DATA AFN1/0.99999999962173909991E0,0.36451060338631902917E3, X 1 0.44218548041288440874E5,0.22467569405961151887E7, X 2 0.49315316723035561922E8,0.43186795279670283193E9, X 3 0.11847992519956804350E10,0.45573267593795103181E9/ X DATA AFD1/1.0E0,0.36651060273229347594E3, X 1 0.44927569814970692777E5,0.23285354882204041700E7, X 2 0.53117852017228262911E8,0.50335310667241870372E9, X 3 0.16575285015623175410E10,0.11746532837038341076E10/ XC XC VALUES OF GI(X) FOR 6 <= X <=12 XC X DATA AGN1/0.99999999920484901956E0,0.51385504875307321394E3, X 1 0.92293483452013810811E5,0.74071341863359841727E7, X 2 0.28142356162841356551E9,0.49280890357734623984E10, X 3 0.35524762685554302472E11,0.79194271662085049376E11, X 4 0.17942522624413898907E11/ X DATA AGD1/1.0E0,0.51985504708814870209E3, X 1 0.95292615508125947321E5,0.79215459679762667578E7, X 2 0.31977567790733781460E9,0.62273134702439012114E10, X 3 0.54570971054996441467E11,0.18241750166645704670E12, X 4 0.15407148148861454434E12/ XC XC VALUES FOR FI(X) FOR X > 12 XC X DATA AFN2/0.19999999999999978257E1,0.22206119380434958727E4, X 1 0.84749007623988236808E6,0.13959267954823943232E9, X 2 0.10197205463267975592E11,0.30229865264524075951E12, X 3 0.27504053804288471142E13,0.21818989704686874983E13/ X DATA AFD2/1.0E0,0.11223059690217167788E4, X 1 0.43685270974851313242E6,0.74654702140658116258E8, X 2 0.58580034751805687471E10,0.20157980379272098841E12, X 3 0.26229141857684496445E13,0.87852907334918467516E13/ XC XC VALUES FOR GI(X) FOR X > 12 XC X DATA AGN2/0.59999999999999993089E1,0.96527746044997139158E4, X 1 0.56077626996568834185E7,0.15022667718927317198E10, X 2 0.19644271064733088465E12,0.12191368281163225043E14, X 3 0.31924389898645609533E15,0.25876053010027485934E16, X 4 0.12754978896268878403E16/ X DATA AGD2/1.0E0,0.16287957674166143196E4, X 1 0.96636303195787870963E6,0.26839734750950667021E9, X 2 0.37388510548029219241E11,0.26028585666152144496E13, X 3 0.85134283716950697226E14,0.11304079361627952930E16, X 4 0.42519841479489798424E16/ XC XC START COMPUTATION XC X X = XVALUE X INDSGN = 1 X IF ( X .LT. ZERO ) THEN X X = -X X INDSGN = -1 X ENDIF XC XC CODE FOR 0 <= |X| <= 6 XC X IF ( X .LE. SIX ) THEN X IF ( X .LT. XLOW ) THEN X SININT = X X ELSE X SUMNUM = ZERO X SUMDEN = ZERO X XSQ = X * X X DO 100 I = 7 , 0 , -1 X SUMNUM = SUMNUM * XSQ + ASINTN(I) X SUMDEN = SUMDEN * XSQ + ASINTD(I) X 100 CONTINUE X SININT = X * SUMNUM / SUMDEN X ENDIF X ENDIF XC XC CODE FOR 6 < |X| <= 12 XC X IF ( X .GT. SIX .AND. X .LE. TWELVE ) THEN X SUMNUM = ZERO X SUMDEN = ZERO X XSQ = ONE / ( X * X ) X DO 200 I = 7 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AFN1(I) X SUMDEN = SUMDEN * XSQ + AFD1(I) X 200 CONTINUE X FIVAL = SUMNUM / ( X * SUMDEN ) X SUMNUM = ZERO X SUMDEN = ZERO X DO 300 I = 8 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AGN1(I) X SUMDEN = SUMDEN * XSQ + AGD1(I) X 300 CONTINUE X GIVAL = XSQ * SUMNUM / SUMDEN X SININT = PIBY2 - FIVAL * COS(X) - GIVAL * SIN(X) X ENDIF XC XC CODES FOR |X| > 12 XC X IF ( X .GT. TWELVE ) THEN X XHIGH = MIN(XHIGH2,XHIGH3) X IF ( X .GT. XHIGH ) THEN X SININT = PIBY2 X ELSE X CX = COS(X) X SX = SIN(X) X XSQ = ONE / ( X * X ) X IF ( X .GT. XHIGH1 ) THEN X SININT = PIBY2 - CX / X - SX * XSQ X ELSE X SUMNUM = ZERO X SUMDEN = ZERO X DO 400 I = 7 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AFN2(I) X SUMDEN = SUMDEN * XSQ + AFD2(I) X 400 CONTINUE X FIVAL = ( ONE - XSQ * SUMNUM / SUMDEN ) / X X SUMNUM = ZERO X SUMDEN = ZERO X DO 500 I = 8 , 0 , -1 X SUMNUM = SUMNUM * XSQ + AGN2(I) X SUMDEN = SUMDEN * XSQ + AGD2(I) X 500 CONTINUE X GIVAL = ( ONE - XSQ * SUMNUM / SUMDEN ) * XSQ X SININT = PIBY2 - CX * FIVAL - SX * GIVAL X ENDIF X ENDIF X ENDIF X IF ( INDSGN .EQ. -1 ) SININT = -SININT X RETURN X END X END_OF_FILE if test 7771 -ne `wc -c <'sinint.f'`; then echo shar: \"'sinint.f'\" unpacked with wrong size! fi # end of 'sinint.f' fi if test -f 'sint.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'sint.f'\" else echo shar: Extracting \"'sint.f'\" \(731 characters\) sed "s/^X//" >'sint.f' <<'END_OF_FILE' X REAL FUNCTION SINT(X) XC XC This function enables the user to tell the test-code XC the particular name of the Si function being tested. XC In this test code, REAL FUNCTION SININT will be used. XC To test other Si function software, the user should XC replace SININT with the desired actual name of the code. XC Users should note that this format assumes a single XC parameter to the function. XC XC XC INPUT PARAMETER: XC XC X - REAL - The argument to the function XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X REAL SININT, X X SINT = SININT(X) X RETURN X END END_OF_FILE if test 731 -ne `wc -c <'sint.f'`; then echo shar: \"'sint.f'\" unpacked with wrong size! fi # end of 'sint.f' fi if test -f 'sitestdp.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'sitestdp.f'\" else echo shar: Extracting \"'sitestdp.f'\" \(8896 characters\) sed "s/^X//" >'sitestdp.f' <<'END_OF_FILE' X PROGRAM DSITEST XC XC XC This program tests the performance of DOUBLE-PRECISION codes which XC calculate the sine-integral special function, defined as XC XC DPSINT(x) = integral from 0 to x of {sin(t)/t} dt XC XC The program compares the value of DPSINT(x+h) with the Taylor XC series value about x. The derivatives of DPSINT at x are XC generated by a standard recurrence relation. XC XC The program assumes the sine-integral is calculated in a XC function with the name DPSINT. Thus users should append their XC Si code at the end of this file and edit the small DPSINT XC function to include the name of their code. XC XC The program uses Cody's MACHAR subroutine to determine certain XC machine parameters. This can be avoided if values for the following XC parameters are available: EPS, IT, IBETA, MACHEP, XMAX, XMIN. XC Values for certain machine-compiler combinations are given in XC XC W.J. CODY Algorithm 665: MACHAR: A subroutine to dynamically XC determine machine parameters, ACM Trans. Math. Soft. 14 (1988) 303-311. XC XC The MACHAR code and a random generator code written by Cody XC are included in this file, though the name has been changed XC to DMACHR. XC XC To ensure that DMACHR and the argument purification section XC work as intended, this code should be compiled with any XC optimisation switched OFF. XC XC Users should note that MACHAR tests the machine arithmetic XC at its very limits. Thus it is possible that warning messages XC might appear about underflows and similar events. These DO XC NOT affect the test results. Such warnings can be avoided XC by explicitly declaring the necessary machine parameters. XC XC MODULES USED: XC XC This code calls the following modules: XC XC (a) DSIDER contained in this file XC XC (b) DPSINT contained in the file DPSINT.F XC XC (c) DMACHR, DREN and DTERMS contained in the file SICISUBD.F XC XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER I,IBETA,IEXP,IOUT,IPT,IRND,IT,ITEST,MACHEP,MAXEXP, X 1 MINEXP,NDEC,NEGEP,NGRD,NPTS,NREN,NTERMS,NTEST, X 2 NT1,NT2,NUMEQ,NUMLA,NUMSM X DOUBLE PRECISION ALOW(5),AUP(5),SDER(0:100), X 2 AIT,AL,ALBETA,AU,DELTA,DIV,DPSINT,DREN, X 3 EPS,EPSNEG,ERR,ERLOSS,HINT,MAXERR,MAXPT,ONE, X 4 PT0625,RELERR,RMSERR,SIX,SIXTEN,SIY,SUM,TEMP, X 5 X,XMAX,XMIN,Y,Z,ZERO X DATA ALOW/0.5D0,-2.5D0,5.0D0,-11.0D0,14.0D0/ X DATA AUP/1.0D0,-2.0D0,6.0D0,-10.0D0,15.0D0/ X DATA PT0625,SIXTEN/0.0625D0,16.0D0/ X DATA ZERO,ONE/0.0D0,1.0D0/ XC XC IOUT = STANDARD OUTPUT STREAM NUMBER XC XC NPTS = NO. OF TEST ARGUMENTS PER INTERVAL XC X DATA IOUT/6/ X DATA NPTS,NTEST/5000,5/ XC XC DETERMINE MACHINE PARAMETERS XC X CALL DMACHR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, X 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) X AIT = DBLE(IT) X ALBETA = LOG(DBLE(IBETA)) X NDEC = ABS(INT(LOG10(DBLE(IBETA))*DBLE(MACHEP))) X NREN = 1234567 XC XC START TESTS XC X DO 300 ITEST = 1 , NTEST X AL = ALOW(ITEST) X AU = AUP(ITEST) X HINT = ( AU - AL ) / DBLE(NPTS) XC XC DETERMINE HOW MANY DERIVATIVES ARE NEEDED XC X CALL DTERMS(ABS(AL),EPS,NT1) X CALL DTERMS(ABS(AU),EPS,NT2) X NTERMS = NT1 X IF ( NT2 .GT. NT1 ) NTERMS = NT2 XC XC SET UP STATISTICS COLLECTION XC X MAXERR = ZERO X RMSERR = ZERO X MAXPT = ZERO X NUMSM = 0 X NUMLA = 0 X DELTA = PT0625 X IF ( ITEST .EQ. 2 .OR. ITEST .EQ. 3 ) DELTA = -PT0625 X DO 200 IPT = 1 , NPTS X X = AL + HINT * DREN(NREN) XC XC PURIFY ARGUMENT SO THAT ( X + DELTA ) IS EXACT XC X Z = SIXTEN * X X TEMP = Z + X X X = TEMP - Z X Y = X + DELTA XC XC CALCULATE SI(Y) AND SI(X+DELTA) XC X SIY = DPSINT(Y) X SIX = DPSINT(X) X CALL DSIDER(X,NTERMS,NDEC,SDER) X SUM = SDER(NTERMS) X DIV = DBLE(NTERMS) + ONE X DO 100 I = 0 , NTERMS-1 X SUM = SUM * DELTA / DIV + SDER(NTERMS-1-I) X DIV = DIV - ONE X 100 CONTINUE X SUM = DELTA * SUM XC XC DETERMINE ERROR AND UPDATE STATS XC X ERR = ( SIY - SIX ) - SUM X IF ( ERR .LT. ZERO ) NUMSM = NUMSM + 1 X IF ( ERR .GT. ZERO ) NUMLA = NUMLA + 1 X RELERR = ABS(ERR/SIY) X RMSERR = RMSERR + RELERR * RELERR X IF ( RELERR .GT. MAXERR ) THEN X MAXERR = RELERR X MAXPT = X X ENDIF X AL = AL + HINT X 200 CONTINUE XC XC OUTPUT RESULTS FOR EACH TEST INTERVAL XC X RMSERR = SQRT(RMSERR/DBLE(NPTS)) X WRITE(IOUT,1000) X WRITE(IOUT,1010)ALOW(ITEST),AUP(ITEST),NPTS X ERLOSS = AIT + LOG(RMSERR) / ALBETA X IF(ERLOSS.LT.ZERO)ERLOSS = ZERO X WRITE(IOUT,1020)RMSERR,ERLOSS X WRITE(IOUT,1030)MAXERR,MAXPT X ERLOSS = AIT + LOG(MAXERR) / ALBETA X IF(ERLOSS.LT.ZERO)ERLOSS = ZERO X WRITE(IOUT,1040)ERLOSS X NUMEQ = NPTS - ( NUMSM + NUMLA ) X WRITE(IOUT,1050)NUMSM,NUMEQ,NUMLA X 300 CONTINUE XC XC TEST SI(X) + SI(-X) = 0.0 XC X WRITE(IOUT,1100) X WRITE(IOUT,1110) X DO 400 I = 1 , 20 X X = DREN(NREN) + DBLE(I-1) X Y = DPSINT(X) + DPSINT(-X) X WRITE(IOUT,1120)X,Y X 400 CONTINUE XC XC PRINT VALUES OF SI(XMAX) = PI/2 AND SI(XMIN) = XMIN XC X WRITE(IOUT,1200) X Y = DPSINT(XMAX) X WRITE(IOUT,1210)Y X WRITE(IOUT,1220)XMIN X Y = DPSINT(XMIN) X WRITE(IOUT,1230)Y XC XC FORMAT STATEMENTS XC X 1000 FORMAT(/////10X,'TEST OF SI(X+DELTA) AGAINST TAYLOR SERIES') X 1010 FORMAT(/10X,'TEST CARRIED OUT ON ',F9.4,' TO ',F9.4,' WITH ', X & I5,' ARGUMENTS') X 1020 FORMAT(//10X,'RMS ERROR = ',D15.6,' IS A LOSS OF', X & F8.3,' SIGNIFICANT DIGITS') X 1030 FORMAT(//10X,'MAX ERROR = ',D15.6,' OCCURRED AT X = ',F10.4) X 1040 FORMAT(/10X,'THIS IS A LOSS OF ',F8.3, X & ' SIGNIFICANT DIGITS') X 1050 FORMAT(///10X,I6,' ARGUMENTS GAVE POSITIVE ERROR',/10X, X & I6,' ARGUMENTS GAVE NO ERROR',/10X,I6, X & ' ARGUMENTS GAVE NEGATIVE ERROR') X 1100 FORMAT(/////10X,'SPECIAL TESTS') X 1110 FORMAT(///10X,'TEST OF SI(X) + SI(-X) = 0.0',//10X, X & 'VALUE OF X',10X,'VALUE OF SI(X)+SI(-X)') X 1120 FORMAT(/10X,F10.4,15X,D16.5) X 1200 FORMAT(/////10X,'CALCULATE SI(XMAX) - SHOULD RETURN PI/2') X 1210 FORMAT(/10X,'VALUE OF SI(XMAX) = ',D15.5) X 1220 FORMAT(///10X,'CALCULATE SI(XMIN) - SHOULD RETURN XMIN =',D15.5) X 1230 FORMAT(/10X,'VALUE OF SI(XMIN) = ',D15.5) X END X X SUBROUTINE DSIDER(X,NMAX,NDIGS,SDER) XC XC This subroutine calculates the necessary number of XC derivatives of Si by using the recurrence relation XC code of Gautschi and Klein from Comm. ACM., vol. 13 XC 1970 pp53-54. XC XC INPUT PARAMETERS: XC XC X - DOUBLE PRECISION - Value of argument XC XC NMAX - INTEGER - Max. no of derivatives XC XC NDIGS - INTEGER - No. of decimal digits of accuracy wanted XC XC OUTPUT PARAMETERS: XC XC SDER - DOUBLE PRECISION ARRAY - The values of the derivatives XC XC XC MODULES CALLED: XC XC This subroutine calls the function DTLNTI which is contained XC in the file SICISUBD.FOR XC XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER N,NDIGS,NLIM,NMAX,NU,N0 X DOUBLE PRECISION DTLNTI,D1,EXPONE,ONE,SDER(0:NMAX), X 1 SIGMA(4),S1,TEN,TWO,X,X1,ZERO X DATA ZERO,ONE,TWO,TEN/0.0D0,1.0D0,2.0D0,10.0D0/ XC XC SET UP INITIAL PARAMETERS FOR RECURRENCE XC X X1 = ABS(X) X SIGMA(1) = COS(X) X SIGMA(2) = -SIN(X) X SIGMA(3) = -SIGMA(1) X SIGMA(4) = -SIGMA(2) X N0 = INT(X1) X IF ( X .NE. ZERO ) THEN X SDER(0) = SIGMA(4) / X X ELSE X SDER(0) = ONE X ENDIF X IF ( N0 .LE. NMAX ) THEN X NLIM = N0 X ELSE X NLIM = NMAX X ENDIF XC XC USE FORWARD RECURRENCE UP TO A CERTAIN VALUE XC X DO 100 N = 1 , NLIM X SDER(N) = ( SIGMA(N-4*((N-1)/4)) - DBLE(N) * SDER(N-1)) / X X 100 CONTINUE XC XC USE BACKWARD RECURRENCE FOR REMAINDER XC X IF ( N0 .LT. NMAX ) THEN X EXPONE = EXP(ONE) X S1 = ZERO X D1 = DBLE(NDIGS) * LOG(TEN) + LOG(TWO) X IF ( NMAX .LE. INT(EXPONE*X1) ) THEN X NU = 1 + INT ( EXPONE * X1 * DTLNTI(D1/(EXPONE*X1)) ) X ELSE X NU = 1 + INT ( DBLE(NMAX) * DTLNTI(D1/DBLE(NMAX)) ) X ENDIF X DO 200 N = NU , N0+2 , -1 X S1 = ( SIGMA(N-4*((N-1)/4)) - X * S1 ) / DBLE(N) X IF ( N .LE. (NMAX+1) ) SDER(N-1) = S1 X 200 CONTINUE X ENDIF X RETURN X END X END_OF_FILE if test 8896 -ne `wc -c <'sitestdp.f'`; then echo shar: \"'sitestdp.f'\" unpacked with wrong size! fi # end of 'sitestdp.f' fi if test -f 'sitestsp.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'sitestsp.f'\" else echo shar: Extracting \"'sitestsp.f'\" \(8758 characters\) sed "s/^X//" >'sitestsp.f' <<'END_OF_FILE' X PROGRAM SITEST XC XC XC This program tests the performance of SINGLE-PRECISION codes which XC calculate the sine-integral special function, defined as XC XC SINT(x) = integral from 0 to x of {sin(t)/t} dt XC XC The program compares the value of SINT(x+h) with the Taylor XC series value about x. The derivatives of SINT at x are XC generated by a standard recurrence relation. XC XC The program assumes the sine-integral is calculated in a XC function with the name SINT. Thus users should append their XC Si code at the end of this file and edit the small SINT XC function to include the name of their code. XC XC The program uses Cody's MACHAR subroutine to determine certain XC machine parameters. This can be avoided if values for the following XC parameters are available: EPS, IT, IBETA, MACHEP, XMAX, XMIN. XC Values for certain machine-compiler combinations are given in XC XC W.J. CODY Algorithm 665: MACHAR: A subroutine to dynamically XC determine machine parameters, ACM Trans. Math. Soft. 14 (1988) 303-311. XC XC The MACHAR code and a random generator code written by Cody XC are included in this file. XC XC To ensure that MACHAR and the argument purification work as XC expected this code should be compiled with any optimisation XC swithced OFF. XC XC Users should note that MACHAR tests the machine arithmetic XC at its very limits. Thus it is possible that warning messages XC might appear about underflows and similar events. These DO XC NOT affect the test results. Such warnings can be avoided XC by explicitly declaring the necessary machine parameters. XC XC MODULES USED: XC XC This code calls the following modules: XC XC (a) SIDERS contained in this file XC XC (b) SINT contained in the file SINT.F XC XC (c) MACHAR, REN and TERMEQ contained in the file SICISUBS.F XC XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER I,IBETA,IEXP,IOUT,IPT,IRND,IT,ITEST,MACHEP,MAXEXP, X 1 MINEXP,NDEC,NEGEP,NGRD,NPTS,NREN,NTERMS,NTEST, X 2 NT1,NT2,NUMEQ,NUMLA,NUMSM X REAL ALOW(5),AUP(5),SDER(0:100), X 2 AIT,AL,ALBETA,AU,DELTA,DIV,EPS,EPSNEG,ERR,ERLOSS, X 3 HINT,MAXERR,MAXPT,ONE,PT0625,RELERR,REN,RMSERR,SINT, X 4 SIX,SIXTEN,SIY,SUM,TEMP,X,XMAX,XMIN,Y,Z,ZERO X DATA ALOW/0.5E0,-2.5E0,5.0E0,-11.0E0,14.0E0/ X DATA AUP/1.0E0,-2.0E0,6.0E0,-10.0E0,15.0E0/ X DATA PT0625,SIXTEN/0.0625E0,16.0E0/ X DATA ZERO,ONE/0.0E0,1.0E0/ XC XC IOUT = STANDARD OUTPUT STREAM NUMBER XC XC NPTS = NO. OF TEST ARGUMENTS PER INTERVAL XC X DATA IOUT/6/ X DATA NPTS,NTEST/5000,5/ XC XC DETERMINE MACHINE PARAMETERS XC X CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, X 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) X AIT = REAL(IT) X ALBETA = LOG(REAL(IBETA)) X NDEC = ABS(INT(LOG10(REAL(IBETA))*REAL(MACHEP))) X NREN = 1234567 XC XC START TESTS XC X DO 300 ITEST = 1 , NTEST X AL = ALOW(ITEST) X AU = AUP(ITEST) X HINT = ( AU - AL ) / REAL(NPTS) XC XC DETERMINE HOW MANY DERIVATIVES ARE NEEDED XC X CALL TERMEQ(ABS(AL),EPS,NT1) X CALL TERMEQ(ABS(AU),EPS,NT2) X NTERMS = NT1 X IF ( NT2 .GT. NT1 ) NTERMS = NT2 XC XC SET UP STATISTICS COLLECTION XC X MAXERR = ZERO X RMSERR = ZERO X MAXPT = ZERO X NUMSM = 0 X NUMLA = 0 X DELTA = PT0625 X IF ( ITEST .EQ. 2 .OR. ITEST .EQ. 3 ) DELTA = -PT0625 X DO 200 IPT = 1 , NPTS X X = AL + HINT * REN(NREN) XC XC PURIFY ARGUMENT SO THAT ( X + DELTA ) IS EXACT XC X Z = SIXTEN * X X TEMP = Z + X X X = TEMP - Z X Y = X + DELTA XC XC CALCULATE SI(Y) AND SI(X+DELTA) XC X SIY = SINT(Y) X SIX = SINT(X) X CALL SIDERS(X,NTERMS,NDEC,SDER) X SUM = SDER(NTERMS) X DIV = REAL(NTERMS) + ONE X DO 100 I = 0 , NTERMS-1 X SUM = SUM * DELTA / DIV + SDER(NTERMS-1-I) X DIV = DIV - ONE X 100 CONTINUE X SUM = DELTA * SUM XC XC DETERMINE ERROR AND UPDATE STATS XC X ERR = ( SIY - SIX ) - SUM X IF ( ERR .LT. ZERO ) NUMSM = NUMSM + 1 X IF ( ERR .GT. ZERO ) NUMLA = NUMLA + 1 X RELERR = ABS(ERR/SIY) X RMSERR = RMSERR + RELERR * RELERR X IF ( RELERR .GT. MAXERR ) THEN X MAXERR = RELERR X MAXPT = X X ENDIF X AL = AL + HINT X 200 CONTINUE XC XC OUTPUT RESULTS FOR EACH TEST INTERVAL XC X RMSERR = SQRT(RMSERR/REAL(NPTS)) X WRITE(IOUT,1000) X WRITE(IOUT,1010)ALOW(ITEST),AUP(ITEST),NPTS X ERLOSS = AIT + LOG(RMSERR) / ALBETA X IF(ERLOSS.LT.ZERO)ERLOSS = ZERO X WRITE(IOUT,1020)RMSERR,ERLOSS X WRITE(IOUT,1030)MAXERR,MAXPT X ERLOSS = AIT + LOG(MAXERR) / ALBETA X IF(ERLOSS.LT.ZERO)ERLOSS = ZERO X WRITE(IOUT,1040)ERLOSS X NUMEQ = NPTS - ( NUMSM + NUMLA ) X WRITE(IOUT,1050)NUMSM,NUMEQ,NUMLA X 300 CONTINUE XC XC TEST SI(X) + SI(-X) = 0.0 XC X WRITE(IOUT,1100) X WRITE(IOUT,1110) X DO 400 I = 1 , 20 X X = REN(NREN) + REAL(I-1) X Y = SINT(X) + SINT(-X) X WRITE(IOUT,1120)X,Y X 400 CONTINUE XC XC PRINT VALUES OF SI(XMAX) = PI/2 AND SI(XMIN) = XMIN XC X WRITE(IOUT,1200) X Y = SINT(XMAX) X WRITE(IOUT,1210)Y X WRITE(IOUT,1220)XMIN X Y = SINT(XMIN) X WRITE(IOUT,1230)Y XC XC FORMAT STATEMENTS XC X 1000 FORMAT(/////10X,'TEST OF SI(X+DELTA) AGAINST TAYLOR SERIES') X 1010 FORMAT(/10X,'TEST CARRIED OUT ON ',F9.4,' TO ',F9.4,' WITH ', X & I5,' ARGUMENTS') X 1020 FORMAT(//10X,'RMS ERROR = ',E15.6,' IS A LOSS OF', X & F8.3,' SIGNIFICANT DIGITS') X 1030 FORMAT(//10X,'MAX ERROR = ',E15.6,' OCCURRED AT X = ',F10.4) X 1040 FORMAT(/10X,'THIS IS A LOSS OF ',F8.3, X & ' SIGNIFICANT DIGITS') X 1050 FORMAT(///10X,I6,' ARGUMENTS GAVE POSITIVE ERROR',/10X, X & I6,' ARGUMENTS GAVE NO ERROR',/10X,I6, X & ' ARGUMENTS GAVE NEGATIVE ERROR') X 1100 FORMAT(/////10X,'SPECIAL TESTS') X 1110 FORMAT(///10X,'TEST OF SI(X) + SI(-X) = 0.0',//10X, X & 'VALUE OF X',10X,'VALUE OF SI(X)+SI(-X)') X 1120 FORMAT(/10X,F10.4,15X,E16.5) X 1200 FORMAT(/////10X,'CALCULATE SINT(XMAX) - SHOULD RETURN PI/2') X 1210 FORMAT(/10X,'VALUE OF SINT(XMAX) = ',E15.5) X 1220 FORMAT(///10X,'CALCULATE SINT(XMIN) - SHOULD RETURN XMIN =',E15.5) X 1230 FORMAT(/10X,'VALUE OF SINT(XMIN) = ',E15.5) X END X X SUBROUTINE SIDERS(X,NMAX,NDIGS,SDER) XC XC This subroutine calculates the necessary number of XC derivatives of Si by using the recurrence relation XC code of Gautschi and Klein from Comm. ACM., vol. 13 XC 1970 pp53-54. XC XC INPUT PARAMETERS: XC XC X - REAL - Value of argument XC XC NMAX - INTEGER - Max. no of derivatives XC XC NDIGS - INTEGER - No. of decimal digits of accuracy wanted XC XC OUTPUT PARAMETERS: XC XC SDER - REAL ARRAY - The values of the derivatives XC XC XC MODULES CALLED: XC XC This subroutine calls the function TLNTIN which is contained XC in the file SICISUBS.FOR XC XC XC AUTHOR: Allan MacLeod XC Dept. of Mathematics and Statistics XC University of Paisley XC Scotland XC (e-mail: macl_ms0@paisley.ac.uk) XC X INTEGER N,NDIGS,NLIM,NMAX,NU,N0 X REAL D1,EXPONE,ONE,SDER(0:NMAX),SIGMA(4), X 1 S1,TEN,TLNTIN,TWO,X,X1,ZERO X DATA ZERO,ONE,TWO,TEN/0.0E0,1.0E0,2.0E0,10.0E0/ XC XC SET UP INITIAL PARAMETERS FOR RECURRENCE XC X X1 = ABS(X) X SIGMA(1) = COS(X) X SIGMA(2) = -SIN(X) X SIGMA(3) = -SIGMA(1) X SIGMA(4) = -SIGMA(2) X N0 = INT(X1) X IF ( X .NE. ZERO ) THEN X SDER(0) = SIGMA(4) / X X ELSE X SDER(0) = ONE X ENDIF X IF ( N0 .LE. NMAX ) THEN X NLIM = N0 X ELSE X NLIM = NMAX X ENDIF XC XC USE FORWARD RECURRENCE UP TO A CERTAIN VALUE XC X DO 100 N = 1 , NLIM X SDER(N) = ( SIGMA(N-4*((N-1)/4)) - REAL(N) * SDER(N-1)) / X X 100 CONTINUE XC XC USE BACKWARD RECURRENCE FOR REMAINDER XC X IF ( N0 .LT. NMAX ) THEN X EXPONE = EXP(ONE) X S1 = ZERO X D1 = REAL(NDIGS) * LOG(TEN) + LOG(TWO) X IF ( NMAX .LE. INT(EXPONE*X1) ) THEN X NU = 1 + INT ( EXPONE * X1 * TLNTIN(D1/(EXPONE*X1)) ) X ELSE X NU = 1 + INT ( REAL(NMAX) * TLNTIN(D1/REAL(NMAX)) ) X ENDIF X DO 200 N = NU , N0+2 , -1 X S1 = ( SIGMA(N-4*((N-1)/4)) - X * S1 ) / REAL(N) X IF ( N .LE. (NMAX+1) ) SDER(N-1) = S1 X 200 CONTINUE X ENDIF X RETURN X END X END_OF_FILE if test 8758 -ne `wc -c <'sitestsp.f'`; then echo shar: \"'sitestsp.f'\" unpacked with wrong size! fi # end of 'sitestsp.f' fi echo shar: End of shell archive. exit 0