/************************************************************************/
/* */
/* ALGORITHM 353 - FILON QUADRATURE */
/* */
/* Stephen M. Chase and Lloyd D. Fosdick 7 July 1967. */
/* With modifications suggested by Bo Einarsson, */
/* Research Institute of National Defence, Tumba, Sweden, */
/* 8 December 1969. */
/* */
/************************************************************************/
/* Translated to Fortran 77 by R. A. Vowels. 27 June 2006. */
/* Translated to Fortran 90 by R. A. Vowels, 13 May 2007. */
/* Translated to PL/I by R. A. Vowels, 13 May 2007. */
/* Fortran 77 version works 27 June 2006. */
/* PL/I version works 2 June 2007. */
/* FSER1 evaluates the integrals */
/* 1 1 */
/* C = S F(X) cos (M pi X) dX, S = S F(X) sin (M pi X) dX */
/* 0 0 */
/* using the Filon quadrature algorithm. */
/* The user may request an evaluation of C only, S only, or both C and S. */
/* With the function F(x) = x**2*(1-x), and M = 64, MAX = 20, LC = 1, */
/* LS = 1 and EPS = 1e-10: */
/* Exact solutions are : */
/* C = -1/(64*pi)**2 = -2.473661710e-5 */
/* S = -6/(64*pi)**2 = -7.381790413e-7 */
/* Actual results obtained were: */
/* C = -2.473661709349E-05 */
/* S = -7.381790410980E-07 */
/* LC = 9 LS = 9 */
/* With EPS changed to 1E-20, the actual results were: */
/* C = -2.473661710018E-05 */
/* S= -7.381790412808E-07 */
/* LC = 20 LS = 20 */
(SUBSCRIPTRANGE, FIXEDOVERFLOW, SIZE):
FSER1: PROCEDURE (F, EPS, MAX, M, C, S, LC, LS) OPTIONS (REORDER);
/* F = name of a function to evaluate F(X), supplied by the user. */
/* EPS = error control for algorithm. */
/* MAX = maximum number of halvings of the step size that are allowed. */
/* M = parameter in formula. */
/* C = value of the cosine integral computed by this subroutine. */
/* S = value of the sine integral computed by this subroutine. */
/* LC = 1 if the cosine integral is to be computed (value returned via */
/* C). The value of step size used to compute C is returned in LC. */
/* (The step size h = 2**(-LC). */
/* = 0 if the cosine integral is not to be computed. */
/* LS = 1 if the sine integral is to be computed (value returned via S). */
/* The value of step size used to compute S is returned in LS. */
/* (The step size h = 2**(-LS). */
/* LS = 0 if the sine integral is not to be computed. */
DECLARE F ENTRY (FLOAT (18)) RETURNS (FLOAT (18));
DECLARE ( EPS, C, S ) FLOAT (18);
DECLARE ( MAX, M, LC, LS ) FIXED BINARY (31);
DECLARE ( PI, XM, F0, F1, SUMSIN, SUMCOS, B1, TMAX, H, T, TP, ODSIN, ODCOS,
XI, THA, FRAC, COSIP, TEMP1, TEMP2, TSQ, A, B2, B3, B4, B5, B, G, S1, S1SQ,
C1, P, PVT1, PVT2, T1, T2 ) FLOAT (18);
DECLARE ( I, N, NSTOP, NST, LN, IN, M_ON_2, IN_ON_2 ) FIXED BINARY (31);
DECLARE ( LLC, LLS, MSWITCH ) FIXED BINARY (31);
PI = 3.14159265358979324E0;
XM = M;
M_ON_2 = M/2;
F1 = 1 - 2 * (M - M_ON_2 * 2);
F0 = F(0.00000000000000000E0);
F1 = F(1.00000000000000000E0) * F1;
SUMCOS = (F1 * F0)/2;
SUMSIN = 0;
B1 = 2.00000000000000000E0/3.00000000000000000E0;
/* TMAX is the switchover point in the angle T. */
/* Our analysis indicates that TMAX = 1/6 is the best for a machine */
/* having a 44 bit floating-point mantissa. */
TMAX = 0.16600000000000000E0;
/* N is the number of the iteration. Note that we start at the */
/* fourth iteration step. */
/* Actually, the first evaluation of an integral is at N = 5, and */
/* therefore, the first comparison of values is at N = 6. */
N = LOG(2*XM) / 0.69300000000000000E0;
PUT SKIP LIST ( N);
/* Both TMAX and N may be changed if the machine for which this */
/* routine is intended has greater or less accuracy than a 44-bit mantissa. */
/* If N is changed, then corresponding changes must be made */
/* in the assignments of H and NSTOP. */
H = 1/2.00000000000000000E0**N;
NSTOP = ISLL(1, N) - 1; /* A simple way to do NSTOP = 2**N-1; */
T = H * XM;
TP = T * PI;
NST = 1;
MSWITCH = 67;
/* LLC and LLS are used by the routine in computed-go-to statements. */
/* As soon as LLS and LLC have been defined, we can use LS and LC */
/* as return values (see description of dummy arguments). */
IF LS <= 0 THEN
LLS = 2;
ELSE
DO;
LLS = 1;
LS = MAX;
END;
IF LC <= 0 THEN
LLC = 2;
ELSE
DO;
LLC = 1;
LC = MAX;
END;
LN = 1;
/* All of the above is executed once per call. */
/* Now the iteration begins. */
L10:
ODCOS, ODSIN = 0;
/* Begin summation for ODCOS and ODSIN. */
DO I = 1 TO NSTOP BY NST;
XI = I;
THA = XI * T;
/* THA*PI is the angle used in this I-th term. */
/* CIR(I*T*PI) is calculated here using the identity */
/* CIR (integer multiple of PI + fractional mult of PI ) */
/* = COS(integer*PI) * CIR(FRAC*PI) */
/* = (+ or -) * CIR(FRAC*PI). */
FRAC = THA;
IN = THA;
THA = IN;
FRAC = (FRAC - THA) * PI;
/* THA is a floating-point integer, FRAC is the fractional part * PI. */
IN_ON_2 = IN/2;
COSIP = 1 - 2*(IN - 2*IN_ON_2);
TEMP1 = COSIP * F(XI*H);
/* TEMP1 = COS(integer part) * F(I*H). */
IF LLS = 2 THEN GO TO L55;
ODSIN = TEMP1 * SIN (FRAC) + ODSIN;
L55: IF LLC = 1 THEN
ODCOS = TEMP1 * COS(FRAC) + ODCOS;
END;
IF MSWITCH = 67 THEN GO TO L67; ELSE GO TO L70;
L67:
NST = 2;
/* Now have made up for the first 4 iteration steps, so reset these */
/* two numbers to look like the general case. */
NSTOP = ISLL(1, N); /* A simple way to get NSTOP = 2**N; */
MSWITCH = 70;
GO TO L92;
L70:
TSQ = TP * TP;
IF T <= TMAX THEN
DO;
/* Power series for small T. */
/* (The closed form is used with larger values of T; see the ELSE */
/* below.) */
/* The power series are (with 'TN' = TP ** N) */
/* A = (2/45)*T3 - (2/315)*T5 + (2/4725)*T7 */
/* B = (2/3) + (2/15)*T2 - (4/105)*T4 + (2/567)*T6 - (4/22275)*T8 */
/* G = (4/3) - (2/15)*T2 + (1/210)*T4 - (1/11340)*T6 */
A = TP * TSQ * (1 - TSQ * (1 - TSQ / 15) / 7.00000000000000000E0)
/ 22.5000000000000000E0;
B2 = B1 * TSQ / 5;
B3 = B2 * TSQ * 2 / 7.00000000000000000E0;
B4 = B3 * TSQ / 10.8000000000000000E0;
B5 = B4 * TSQ * 14.0000000000000000E0 / 275.000000000000000E0;
B = B1 + B2 - B3 + B4 - B5;
G = 2*B1 - B2 + B3/8 - B4/40 + 5*B5/896.000000000000000E0;
GO TO L80;
END;
ELSE
DO;
/* Closed form of the coefficients, where again 'TN' means TP**N. */
/* A = 1/TP + COS(TP)*SIN(TP)/T2 - 2*(SIN(TP))**2/TP */
/* B = 2*((1 + (COS(TP))**2/T2 - 2*SIN(TP)*COS(TP)/T3) */
/* G = 4*(SIN(TP)/T3 - COS(TP)/T2) */
IN = T;
IN_ON_2 = IN/2;
TEMP1 = 1 - 2 * (IN - 2*IN_ON_2);
TEMP2 = IN;
/* TEMP1 is COS (integer part of TP), TEMP2 is the fractional part of TP. */
TEMP2 = (T - TEMP2) * PI;
S1 = TEMP1 * SIN (TEMP2);
C1 = TEMP1 * COS (TEMP2);
P = S1 * C1;
S1SQ = S1 * S1;
A = (((-2*S1SQ/TP) + P)/ TP + 1) / TP;
B = 2 * ((-2*P/TP) + 2 - S1SQ) / TSQ;
G = 4 * (S1 / TP - C1) / TSQ;
END;
L80:
IF LLS = 2 THEN GO TO L85;
/* Have calculated the coefficients, now ready for the integration formulas. */
T2 = H * (A * (F0 - F1) + B * SUMSIN + G * ODSIN);
/* ENDT1 is a subroutine that checks for the convergence of the */
/* previous value to within EPS2, where */
/* EPS2 = (1.0 + ABSF(present value))*EPS */
/* EPS is supplied by the user. */
CALL ENDT1 (PVT2, T2, EPS, S, LLS, LN);
IF LLS = 1 THEN GO TO L85;
L84:
LS = N;
L85:
IF LLC = 2 THEN GO TO L90;
/* This is the cosine integral. */
T1 = H * (B * SUMCOS + G * ODCOS);
CALL ENDT1 (PVT1, T1, EPS, C, LLC, LN);
IF LLC = 1 THEN GO TO L90;
LC = N;
L90:
LN = 2;
IF LLC + LLS > 3 THEN GO TO L100;
L92:
N = N + 1;
IF N > MAX THEN GO TO L100;
H = H / 2;
T = T / 2;
TP = TP / 2;
NSTOP = 2 * NSTOP;
SUMSIN = SUMSIN + ODSIN;
SUMCOS = SUMCOS + ODCOS;
GO TO L10;
L100:
S = T2;
C = T1;
RETURN;
END FSER1;
ENDT1: PROCEDURE (PREVQT, QUANT, EPS, VALUE, L1, L2) OPTIONS (REORDER);
DECLARE ( PREVQT, QUANT, EPS, VALUE ) FLOAT (18);
DECLARE ( L1, L2 ) FIXED BINARY (31);
DECLARE ( REPS ) FLOAT (18);
IF L2 = 1 THEN GO TO L29;
REPS = EPS * (1.00000000000000000E0 + ABS (QUANT));
IF ABS(PREVQT - QUANT) <= REPS THEN
DO;
VALUE = QUANT;
L1 = 2;
RETURN;
END;
L29:
PREVQT = QUANT;
END ENDT1;