C ALGORITHM 453, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM C VOL. 16, NO. 8, August, 1973, PP.486--487. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Fortran/ # Fortran/Sp/ # Fortran/Sp/Drivers/ # Fortran/Sp/Drivers/Makefile # Fortran/Sp/Drivers/driver.f # Fortran/Sp/Drivers/res # Fortran/Sp/Src/ # Fortran/Sp/Src/src.f # This archive created: Thu Dec 15 13:28:05 2005 export PATH; PATH=/bin:$PATH if test ! -d 'Fortran' then mkdir 'Fortran' fi cd 'Fortran' if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << "SHAR_EOF" > 'Makefile' all: Res src.o: src.f $(F77) $(F77OPTS) -c src.f driver.o: driver.f $(F77) $(F77OPTS) -c driver.f DRIVERS= driver RESULTS= Res Objs1= driver.o src.o driver: $(Objs1) $(F77) $(F77OPTS) -o driver $(Objs1) $(SRCLIBS) Res: driver ./driver >Res diffres:Res res echo "Differences in results from driver" $(DIFF) Res res clean: rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS) SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << "SHAR_EOF" > 'driver.f' program main c*********************************************************************** c cc TOMS453_PRB tests BROMIN, ACM TOMS algorithm 453. c implicit none write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOMS453_PRB' write ( *, '(a)' ) ' Test TOMS algorithm 453, which computes' write ( *, '(a)' ) ' the complex abscissas and weights for ' write ( *, '(a)' ) ' Gaussian quadrature formulas for ' write ( *, '(a)' ) ' Bromwich''s integral.' call test01 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOMS453_PRB' write ( *, '(a)' ) ' Normal end of execution.' stop end subroutine test01 c*********************************************************************** c cc TEST01 tests BROMIN, ACM TOMS algorithm 453. c implicit none integer nhalf_max integer n_num integer s_num parameter ( nhalf_max = 12 ) parameter ( n_num = 3 ) parameter ( s_num = 4 ) double precision eps integer i integer ier integer j integer k integer n integer n_half integer n_vec(n_num) double precision s double precision s_vec(s_num) double precision tol double precision wi(nhalf_max) double precision wr(nhalf_max) double precision xi(nhalf_max) double precision xr(nhalf_max) save n_vec save s_vec data n_vec / 6, 9, 12 / data s_vec / 0.0D+00, 0.1D+00, 1.0D+00, 4.0D+00 / tol = 0.1D-8 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST01' write ( *, '(a)' ) ' Determine abscissas and weights for' write ( *, '(a)' ) ' a variety of values of S and N.' do i = 1, n_num n = n_vec(i) n_half = ( n + 1 ) / 2 do j = 1, s_num s = s_vec(j) write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' N = ', n write ( *, '(a,g14.6)' ) ' S = ', s call bromin ( n, s, tol, xr, xi, wr, wi, eps, ier ) if ( 0 .lt. ier ) then write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'BROMIN returned IER = ', ier else if ( ier .eq. -1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Note that the requested accuracy' write ( *, '(a)' ) ' was not achieved.' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' XR XI WR WI' write ( *, '(a)' ) ' ' do k = 1, n_half write ( *, '(4(2x,g14.6))' ) xr(k), xi(k), wr(k), wi(k) end do end if end do end do return end CS REAL FUNCTION GAMMA(X) DOUBLE PRECISION FUNCTION DGAMMA(X) C---------------------------------------------------------------------- C C This routine calculates the GAMMA function for a real argument X. C Computation is based on an algorithm outlined in reference 1. C The program uses rational functions that approximate the GAMMA C function to at least 20 significant decimal digits. Coefficients C for the approximation over the interval (1,2) are unpublished. C Those for the approximation for X .GE. 12 are from reference 2. C The accuracy achieved depends on the arithmetic system, the C compiler, the intrinsic functions, and proper selection of the C machine-dependent constants. C C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants C C beta - radix for the floating-point representation C maxexp - the smallest positive power of beta that overflows C XBIG - the largest argument for which GAMMA(X) is representable C in the machine, i.e., the solution to the equation C GAMMA(XBIG) = beta**maxexp C XINF - the largest machine representable floating-point number; C approximately beta**maxexp C EPS - the smallest positive floating-point number such that C 1.0+EPS .GT. 1.0 C XMININ - the smallest positive floating-point number such that C 1/XMININ is machine representable C C Approximate values for some important machines are: C C beta maxexp XBIG C C CRAY-1 (S.P.) 2 8191 966.961 C Cyber 180/855 C under NOS (S.P.) 2 1070 177.803 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 128 35.040 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 1024 171.624 C IBM 3033 (D.P.) 16 63 57.574 C VAX D-Format (D.P.) 2 127 34.844 C VAX G-Format (D.P.) 2 1023 171.489 C C XINF EPS XMININ C C CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 C Cyber 180/855 C under NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 C IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76 C VAX D-Format (D.P.) 1.70D+38 1.39D-17 5.88D-39 C VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.12D-308 C C******************************************************************* C******************************************************************* C C Error returns C C The program returns the value XINF for singularities or C when overflow would occur. The computation is believed C to be free of underflow and overflow. C C C Intrinsic functions required are: C C INT, DBLE, EXP, LOG, REAL, SIN C C C References: "An Overview of Software Development for Special C Functions", W. J. Cody, Lecture Notes in Mathematics, C 506, Numerical Analysis Dundee, 1975, G. A. Watson C (ed.), Springer Verlag, Berlin, 1976. C C Computer Approximations, Hart, Et. Al., Wiley and C sons, New York, 1968. C C Latest modification: October 12, 1989 C C Authors: W. J. Cody and L. Stoltz C Applied Mathematics Division C Argonne National Laboratory C Argonne, IL 60439 C C---------------------------------------------------------------------- INTEGER I,N LOGICAL PARITY CS REAL DOUBLE PRECISION 1 C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, 2 TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO DIMENSION C(7),P(8),Q(8) C---------------------------------------------------------------------- C Mathematical constants C---------------------------------------------------------------------- CS DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/, CS 1 SQRTPI/0.9189385332046727417803297E0/, CS 2 PI/3.1415926535897932384626434E0/ DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/, 1 SQRTPI/0.9189385332046727417803297D0/, 2 PI/3.1415926535897932384626434D0/ C---------------------------------------------------------------------- C Machine dependent parameters C---------------------------------------------------------------------- CS DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/, CS 1 XINF/3.4E38/ DATA XBIG,XMININ,EPS/34.844D0,5.88D-39,1.39D-17/, 1 XINF/1.70D+38/ C---------------------------------------------------------------------- C Numerator and denominator coefficients for rational minimax C approximation over (1,2). C---------------------------------------------------------------------- CS DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, CS 1 -3.79804256470945635097577E+2,6.29331155312818442661052E+2, CS 2 8.66966202790413211295064E+2,-3.14512729688483675254357E+4, CS 3 -3.61444134186911729807069E+4,6.64561438202405440627855E+4/ CS DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, CS 1 -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, CS 2 2.25381184209801510330112E+4,4.75584627752788110767815E+3, CS 3 -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/ DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1, 1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2, 2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4, 3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/ DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2, 1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3, 2 2.25381184209801510330112D+4,4.75584627752788110767815D+3, 3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/ C---------------------------------------------------------------------- C Coefficients for minimax approximation over (12, INF). C---------------------------------------------------------------------- CS DATA C/-1.910444077728E-03,8.4171387781295E-04, CS 1 -5.952379913043012E-04,7.93650793500350248E-04, CS 2 -2.777777777777681622553E-03,8.333333333333333331554247E-02, CS 3 5.7083835261E-03/ DATA C/-1.910444077728D-03,8.4171387781295D-04, 1 -5.952379913043012D-04,7.93650793500350248D-04, 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02, 3 5.7083835261D-03/ C---------------------------------------------------------------------- C Statement functions for conversion between integer and float C---------------------------------------------------------------------- CS CONV(I) = REAL(I) CONV(I) = DBLE(I) PARITY = .FALSE. FACT = ONE N = 0 Y = X IF (Y .LE. ZERO) THEN C---------------------------------------------------------------------- C Argument is negative C---------------------------------------------------------------------- Y = -X Y1 = AINT(Y) RES = Y - Y1 IF (RES .NE. ZERO) THEN IF (Y1 .NE. AINT(Y1*HALF)*TWO) PARITY = .TRUE. FACT = -PI / SIN(PI*RES) Y = Y + ONE ELSE RES = XINF GO TO 900 END IF END IF C---------------------------------------------------------------------- C Argument is positive C---------------------------------------------------------------------- IF (Y .LT. EPS) THEN C---------------------------------------------------------------------- C Argument .LT. EPS C---------------------------------------------------------------------- IF (Y .GE. XMININ) THEN RES = ONE / Y ELSE RES = XINF GO TO 900 END IF ELSE IF (Y .LT. TWELVE) THEN Y1 = Y IF (Y .LT. ONE) THEN C---------------------------------------------------------------------- C 0.0 .LT. argument .LT. 1.0 C---------------------------------------------------------------------- Z = Y Y = Y + ONE ELSE C---------------------------------------------------------------------- C 1.0 .LT. argument .LT. 12.0, reduce argument if necessary C---------------------------------------------------------------------- N = INT(Y) - 1 Y = Y - CONV(N) Z = Y - ONE END IF C---------------------------------------------------------------------- C Evaluate approximation for 1.0 .LT. argument .LT. 2.0 C---------------------------------------------------------------------- XNUM = ZERO XDEN = ONE DO 260 I = 1, 8 XNUM = (XNUM + P(I)) * Z XDEN = XDEN * Z + Q(I) 260 CONTINUE RES = XNUM / XDEN + ONE IF (Y1 .LT. Y) THEN C---------------------------------------------------------------------- C Adjust result for case 0.0 .LT. argument .LT. 1.0 C---------------------------------------------------------------------- RES = RES / Y1 ELSE IF (Y1 .GT. Y) THEN C---------------------------------------------------------------------- C Adjust result for case 2.0 .LT. argument .LT. 12.0 C---------------------------------------------------------------------- DO 290 I = 1, N RES = RES * Y Y = Y + ONE 290 CONTINUE END IF ELSE C---------------------------------------------------------------------- C Evaluate for argument .GE. 12.0, C---------------------------------------------------------------------- IF (Y .LE. XBIG) THEN YSQ = Y * Y SUM = C(7) DO 350 I = 1, 6 SUM = SUM / YSQ + C(I) 350 CONTINUE SUM = SUM/Y - Y + SQRTPI SUM = SUM + (Y-HALF)*LOG(Y) RES = EXP(SUM) ELSE RES = XINF GO TO 900 END IF END IF C---------------------------------------------------------------------- C Final adjustments and return C---------------------------------------------------------------------- IF (PARITY) RES = -RES IF (FACT .NE. ONE) RES = FACT / RES CS900 GAMMA = RES 900 DGAMMA = RES RETURN C ---------- Last line of GAMMA ---------- END SHAR_EOF fi # end of overwriting check if test -f 'res' then echo shar: will not over-write existing file "'res'" else cat << "SHAR_EOF" > 'res' TOMS453_PRB Test TOMS algorithm 453, which computes the complex abscissas and weights for Gaussian quadrature formulas for Bromwich's integral. TEST01 Determine abscissas and weights for a variety of values of S and N. N = 6 S = 0.00000 XR XI WR WI 6.48253 1.49931 -27.8412 315.065 5.46926 4.51993 37.1570 -124.634 3.04821 7.65179 -9.31577 13.1079 N = 6 S = 0.100000 XR XI WR WI 6.58346 1.51198 -22.4305 290.584 5.56941 4.55940 30.2273 -114.988 3.14713 7.72402 -7.74421 12.1904 N = 6 S = 1.00000 XR XI WR WI 7.49064 1.62150 1.67441 122.163 6.47051 4.90012 -0.490540 -46.7814 4.03885 8.34560 -0.683871 4.88323 N = 6 S = 4.00000 XR XI WR WI 10.5053 1.94212 0.405800 1.70332 9.47263 5.89402 -0.357778 -0.458555 7.02212 10.1435 0.353115E-01 0.271604E-01 N = 9 S = 0.00000 XR XI WR WI 10.5793 0.00000 19740.0 -0.00000 10.2467 3.15936 -14158.4 -3267.65 9.20449 6.34804 4931.37 2600.30 7.28383 9.61625 -650.747 -710.452 3.97536 13.1066 7.79286 48.4906 N = 9 S = 0.100000 XR XI WR WI 10.6802 0.00000 17320.4 -0.00000 10.3474 3.17594 -12442.5 -2686.32 9.30475 6.38204 4360.66 2145.19 7.38342 9.66975 -587.140 -590.179 4.07435 13.1844 8.82134 40.9010 N = 9 S = 1.00000 XR XI WR WI 11.5874 0.00000 4906.08 -0.00000 11.2533 3.32134 -3542.61 -328.959 10.2069 6.68014 1272.27 283.871 8.28004 10.1384 -187.912 -88.0612 4.96613 13.8647 5.70963 7.30219 N = 9 S = 4.00000 XR XI WR WI 14.6045 0.00000 28.9272 -0.00000 14.2673 3.76533 -19.7412 4.62697 13.2119 7.58890 6.10477 -2.73269 11.2719 11.5626 -0.771534 0.429792 7.94671 15.9216 0.276524E-01 -0.931263E-02 N = 12 S = 0.00000 XR XI WR WI 14.4930 1.61733 -125948. 0.100930E+07 13.9836 4.86081 234619. -584473. 12.9259 8.13362 -144785. 182670. 11.2246 11.4674 40404.5 -24658.7 8.67056 14.9254 -4375.78 283.201 4.70228 18.6891 84.9833 77.5368 N = 12 S = 0.100000 XR XI WR WI 14.5938 1.62344 -102398. 857076. 14.0843 4.87931 190985. -498186. 13.0262 8.16504 -118182. 157306. 11.3245 11.5128 33151.5 -21898.5 8.76993 14.9867 -3628.70 400.168 4.80135 18.7702 73.2458 59.9454 N = 12 S = 1.00000 XR XI WR WI 15.5004 1.67741 -13676.4 184762. 14.9895 5.04267 25882.9 -109815. 13.9287 8.44250 -16513.1 36822.0 12.2232 11.9134 4876.46 -6045.66 9.66460 15.5270 -584.577 318.896 5.69358 19.4846 15.2259 3.43448 N = 12 S = 4.00000 XR XI WR WI 18.5175 1.84584 29.0190 540.588 18.0029 5.55230 -47.3026 -313.027 16.9351 9.30721 21.7632 101.905 15.2200 13.1600 -3.45609 -17.2573 12.6508 17.2045 0.498979E-01 1.24831 8.67366 21.6952 0.988037E-02 -0.199247E-01 TOMS453_PRB Normal end of execution. SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << "SHAR_EOF" > 'src.f' SUBROUTINE BROMIN ( N, S, TOL, XR, XI, WR, WI, EPS, IER ) DOUBLE PRECISION AK, AN, ARG, CI, CR, D, D1, D2, E, EPS, & FAC, FACTI, FACTR, PI, PR, QI, QR, RI, RR, S, T1, T2, & TOL, U, V, WI, WR, XI, XR, YI, YR, Z DOUBLE PRECISION DGAMMA INTEGER IER, J, K, L, N, N1, NUM, NUP, IGNAL DIMENSION XR(N), XI(N), WR(N), WI(N) C THIS SUBROUTINE CALCULATES ABSCISSAS AND WEIGHTS OF THE C GAUSSIAN QUADRATURE FORMULA OF ORDER N FOR THE BROMWICH C INTEGRAL. ONLY THE ABSCISSAS OF THE FIRST QUADRANT OF C THE COMPLEX PLANE, THE REAL ABSCISSA (IF N IS ODD) AND C THE CORRESPONDING WEIGHTS ARE CALCULATED. THE OTHER C ABSCISSAS AND WEIGHTS ARE COMPLEX CONJUGATES. C INPUT PARAMETERS C N, ORDER OF THE QUADRATURE FORMULA. C N MUST BE GREATER THAN 2. C TOL, REQUESTED RELATIVE ACCURACY OF THE ABSCISSAS. C S, PARAMETERS OF THE WEIGHT FUNCTION. C OUTPUT PARAMETERS C XR AND XI CONTAIN THE REAL AND IMAGINARY PARTS OF C THE ABSCISSAS. IF N IS ODD, THE REAL ABSCISSA C IS XR(1). C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS OF C THE CORRESPONDING WEIGHTS. C EPS IS A CRUDE ESTIMATION OF THE OBTAINED RELATIVE C ACCURACY OF THE ABSCISSAS. C IER IS AN ERROR CODE. C IF IER = 0, THE COMPUTATION IS CARRIED OUT TO C THE REQUESTED ACCURACY. C IF IER.GT.0, THE IER-TH ABSCISSA IS NOT FOUND. C IF IER = -1, THE COMPUTATIONS ARE CARRIED OUT, C BUT THE REQUESTED ACCURACY IS NOT C ACHIEVED. C IF IER = -2, N IS LESS THAN 3. C FUNCTION PROGRAMS REQUIRED C FUNCTION GAMMA(X), WHICH EVALUATES THE GAMMA C FUNCTION FOR POSITIVE X. IER = -2 IF ( N .LT. 3 ) RETURN N1 = ( N + 1 ) / 2 L = N - 1 AN = N IER = 0 EPS = TOL ARG = 0.034D0 * ( 30.D0 + AN + AN ) / ( AN - 1.D0 ) FACTR = DCOS ( ARG ) FACTI = DSIN ( ARG ) FAC = 1.D0 AK = 0.D0 DO 10 K = 1, L AK = AK + 1.D0 FAC = -FAC * AK 10 CONTINUE FAC = FAC * ( AN+AN+S-2.D0)**2 / ( AN * DGAMMA ( AN+S-1.D0 ) ) C CALCULATION OF AN APPROXIMATION OF THE FIRST ABSCISSA. YR = 1.333D0 * AN + S - 1.5D0 YI = 0.0D0 IF ( N - N1 - N1 ) 30, 20, 20 20 YI = YI + 1.6D0 + 0.07D0 * S C START MAIN LOOP 30 DO 140 K = 1, N1 E = TOL IGNAL = 0 NUM = 0 NUP = 0 C NEWTON-RAPHSON METHOD. D = YR * YR + YI * YI YR = YR / D YI = -YI / D GO TO 50 40 IGNAL = 1 50 QR = S * YR - 1.D0 QI = S * YI PR = (S+1.D0) * ( (S+2.D0) * ( YR*YR-YI*YI ) - 2.D0*YR ) + 1.D0 PI = 2.D0 * ( S + 1.D0 ) * YI * ( ( S + 2.D0 ) * YR - 1.D0 ) Z = 2.D0 DO 60 J = 3, N RR = QR RI = QI QR = PR QI = PI Z = Z + 1.D0 U = Z + S - 2.D0 V = U + Z D = ( V * YR + ( 2.D0 - S ) / ( V - 2.D0 ) ) / U D1 = ( Z - 1.D0 ) * V / ( U * ( V - 2.D0 ) ) D2 = V * YI / U PR = ( V - 1.D0 ) * ( QR * D - QI * D2 ) + D1 * RR PI = ( V - 1.D0 ) * ( QI * D + QR * D2 ) + D1 * RI 60 CONTINUE IF ( IGNAL .EQ. 1 ) GO TO 100 D = ( YR * YR + YI * YI ) * V D1 = ( ( PR + QR ) * YR + ( PI + QI ) * YI ) / D + PR D2 = ( ( PI + QI ) * YR - ( PR + QR ) * YI ) / D + PI D = ( D1 * D1 + D2 * D2 ) * AN T1 = PR * YR - PI * YI T2 = PI * YR + PR * YI CR = ( T1 * D1 + T2 * D2 ) / D CI = ( T2 * D1 - T1 * D2 ) / D YR = YR - CR YI = YI - CI NUM = NUM + 1 C TEST OF CONVERGENCE OF ITERATION PROCESS. IF ( CR*CR + CI*CI - E*E * ( YR*YR + YI*YI ) ) 40, 40, 70 C TEST OF NUMBER OF ITERATION STEPS. 70 IF ( NUM - 10 ) 50, 50, 80 80 E = E * 10.D0 IER = -1 NUP = NUP + 1 IF ( NUP - 5 ) 50, 50, 90 90 IER = K RETURN C CALCULATION OF WEIGHTS. 100 IF ( EPS .GE. E ) GO TO 110 EPS = E 110 D = ( QR * QR + QI * QI )**2 D1 = YR * QR + YI * QI D2 = YI * QR - YR * QI WR(K) = FAC * ( D1 * D1 - D2 * D2 ) / D WI(K) = 2.D0 * FAC * D2 * D1 / D D = YR * YR + YI * YI XR(K) = YR / D XI(K) = -YI / D IF ( K + 1 - N1 ) 130, 120, 150 120 FACTR = DCOS ( 1.5D0 * ARG ) FACTI = DSIN ( 1.5D0 * ARG ) C CALCULATION OF AN APPROXIMATION OF THE (K+1)-TH ABSCISSA. 130 YR = ( XR(K) + 0.67D0*AN ) * FACTR - XI(K) * FACTI - 0.67D0*AN YI = ( XR(K) + 0.67D0*AN ) * FACTI + XI(K) * FACTR 140 CONTINUE 150 RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. # End of shell archive exit 0