C ALGORITHM 614 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 2, C JUN., 1984, P. 152-160. SUBROUTINE INTHP(A, B, D, F, M, P, EPS, INF, QUADR) INT 10 C C THIS SUBROUTINE COMPUTES INTEGRAL OF FUNCTIONS WHICH C MAY HAVE SINGULARITIES AT ONE OR BOTH END-POINTS OF AN C INTERVAL (A,B), SEE [1, 2]. IT CONTAINS FOUR DIFFERENT C QUADRATURE ROUTINES: ONE OVER A FINITE INTERVAL (A,B), C TWO OVER (A,+INFINITY), AND ONE OVER (-INFINITY,+INFINITY). C OF THE TWO FORMULAS OVER (A,+INFINITY), THE FIRST (INF=2 C BELOW) IS MORE SUITED TO NON-OSCILLATORY INTEGRANDS, WHILE C THE SECOND (INF=3) IS MORE SUITED TO OSCILLATORY INTEGRANDS. C THE USER SUPPLIES THE INTEGRAND FUNCTION, HE SPECIFIES THE C INTERVAL, AS WELL AS THE RELATIVE ERROR TO WHICH THE INTE- C GRAL IS TO BE EVALUATED. C THE FORMULAS ARE OPTIMAL IN CERTAIN HARDY SPACES H(P,DD), C SEE [1, 2]. HERE DD IS AN OPEN DOMAIN IN THE COMPLEX PLANE, C A AND B BELONG TO THE BOUNDARY OF DD AND H(P,DD), P.GT.1, IS C THE SET OF ALL ANALYTIC FUNCTONS IN DD WHOSE P-TH NORM DEFI- C NED AS IN [2] IS FINITE. C IF THE USER IS UNABLE TO SPECIFY THE PARAMETERS P AND D C OF THE SPACE H(P,DD) TO WHICH HIS INTEGRAND BELONGS, THE C ALGORITHM TERMINATES ACCORDING TO A HEURISTIC CRITERION, SEE C [2] AND COMMENTS TO EPS. C IF THE USER CAN SPECIFY THE PARAMETERS P AND D OF THE C SPACE H(P,DD) TO WHICH HIS INTEGRAND BELONGS, THE ALGORITHM C TERMINATES WITH AN ANSWER HAVING A GUARANTEED ACCURACY ( DE- C TEMINISTIC CRITERION, SEE [1, 2] AND COMMENTS TO EPS). C C C C C INPUT PARAMETERS C C C A = LOWER LIMIT OF INTEGRATION (SEE COMMENTS TO INF). C C B = UPPER LIMIT OF INTEGRATION (SEE COMMENTS TO INF). C C D = A PARAMETER OF THE CLASS H(P,DD) (SEE COMMENTS TO C INF). C C USER SETS D: C C HEURISTIC TERMINATION C = ANY REAL NUMBER C C DETERMINISTIC TERMINATION C = A NUMBER IN THE RANGE 0.LT.D.LE.PI/2. C C F = A NAME OF AN EXTERNAL INTEGRAND FUNCTION TO BE C SUPPLIED BY THE USER. F(X) COMPUTES THE VALUE OF C A FUNCTION F AT A POINT X. THE STATEMENT C ...EXTERNAL F... MUST APPEAR IN THE MAIN PROGRAM. C C M = MAXIMAL NUMBER OF FUNCTION EVALUATIONS ALLOWED IN C THE COMPUTATIONS, M.GE.3.( ALTERED ON EXIT ). C C P = 0, 1, .GT.1 A PARAMETER OF THE CLASS H(P,DD). C C USER SETS P: C = 0 - HEURISTIC TERMINATION. C = 1 - DETERMINISTIC TERMINATION WITH THE INFINITY C NORM. C .GT.1 -DETERMINISTIC TERMINATION WITH THE P-TH NORM. C C EPS = A REAL NUMBER - THE RELATIVE ERROR BOUND - SEE C REMARKS BELOW. ( ALTERED ON EXIT ). C C C INF = 1, 2, 3, 4 - INFORMATION PARAMETER. ( ALTERED ON EXIT ). C C = 1 SIGNIFIES AN INFINITE INTERVAL (A,B)=REAL LINE, C A AND B ANY NUMBERS. C DETERMINISTIC TERMINATION - C DD=STRIP(Z:ABS(IM(Z)).LT.D). C C = 2 SIGNIFIES A SEMI-INFINITE INTERVAL (A, +INFINITY) C USER SUPPLIES A, B ANY NUMBER. C QUADRATURE SUITED TO NON-OSCILLATORY INTEGRANDS. C DETERMINISTIC TERMINATION - C DD=SECTOR(Z:ABS(ARG(Z-A)).LT.D). C C = 3 SIGNIFIES A SEMI INFINITE INTERVAL (A,+INFINITY) C USER SUPPLIES A, B ANY NUMBER. C QUADRATURE SUITED TO OSCILLATORY INTEGRANDS. C DETERMINISTIC TERMINATION - C DD=REGION(Z:ABS(ARG(SINH(Z-A))).LT.D). C C = 4 SIGNIFIES A FINITE INTERVAL (A,B). C USER SUPPLIES A AND B. C DETERMINISTIC TERMINATION - C DD=LENS REGION(Z:ABS(ARG((Z-A)/(B-Z))).LT.D). C C C C C OUTPUT PARAMETERS C C C M = THE NUMBER OF FUNCTION EVALUATIONS USED IN THE C QUADRATURE. C C EPS = THE RELATIVE ERROR BOUND (SEE REMARKS BELOW). C C DETERMINISTIC TERMINATION C C = THE RELATIVE ERROR REXA BOUND, I.E., C REXA(F,M(OUTPUT)) .LE. EPS. C C HEURISTIC TERMINATION C C = MAX(EPS(INPUT),MACHEP). C C INF = 0, 1 - DETERMINISTIC TERMINATION C C = 0 COMPUTED QUADRATURE QCOM(F,M(EPS)), SEE REMARKS C BELOW. C C = 1 COMPUTED QUADRATURE QCOM(F,M1), SEE REMARKS C BELOW. C C INF = 2, 3, 4 - HEURISTIC TERMINATION. C C = 2 INTEGRATION COMPLETED WITH EPS=MAX(EPS(INPUT), C MACHEP). WE CAN EXPECT THE RELATIVE ERROR C REXA TO BE OF THE ORDER OF EPS (FOR SOME P.GE.1). C C = 3 INTEGRATION NOT COMPLETED. ATTEMPT TO EXCEED THE C MAXIMAL ALLOWED NUMBER OF FUNCTION EVALUATIONS M. C TRUNCATION CONDITIONS (SEE [2]) SATISFIED. QUADR C SET TO BE EQUAL TO THE LAST TRAPEZOIDAL APPRO- C XIMATION. IT IS LIKELY THAT QUADR APPROXIMATES THE C INTEGRAL IF M IS LARGE. C C = 4 INTEGRATION NOT COMPLETED. ATTEMPT TO EXCEED THE C MAXIMAL ALLOWED NUMBER OF FUNCTION EVALUATIONS M. C TRUNCATION CONDITIONS (SEE [2]) NOT SATISFIED. C QUADR SET TO BE EQUAL TO THE COMPUTED TRAPEZOIDAL C APPROXIMATION. IT IS UNLIKELY THAT QUADR APPROXIMATES C THE INTEGRAL. C C INF = 10, 11, 12, 13 - INCORRECT INPUT C C = 10 M.LT.3. C C = 11 P DOES NOT SATISFY P=0, P=1 OR P.GT.1 OR IN THE C CASE OF DETERMINISTIC TERMINATION D DOES NOT C SATISFY 0.LT.D.LE.PI/2. C C = 12 A.GE.B IN CASE OF A FINITE INTERVAL. C C = 13 INF NOT EQUAL TO 1, 2, 3, OR 4. C C C QUADR = THE COMPUTED VALUE OF QUADRATURE. C C C C C REMARKS: C C LET QEXA(F,M) ( QCOM(F,M) ) BE THE EXACT (COMPUTED) C VALUE OF THE QUADRATURE WITH M FUNCTION EVALUATIONS, C AND LET REXA(F,M) ( RCOM(F,M) ) BE THE RELATIVE ERROR C OF QEXA (QCOM) ,I.E., C REXA(F,M)=ABS(INTEGRAL(F)-QEXA(F,M))/NORM(F), C RCOM(F,M)=ABS(INTEGRAL(F)-QCOM(F,M))/NORM(F), C WITH THE NOTATION 0/0=0. C DUE TO THE ROUNDOFF ONE CANNOT EXPECT THE ERROR C RCOM TO BE LESS THAN THE RELATIVE MACHINE PRECISION C MACHEP. THEREFORE THE INPUT VALUE OF EPS IS CHANGED C ACCORDING TO THE FORMULA C EPS=MAX(EPS,MACHEP). C C DETERMINISTIC TERMINATON CASE C C THE NUMBER OF FUNCTON EVALUATIONS M(EPS) IS COMPUTED C SO THAT THE ERROR REXA IS NO GREATER THAN EPS,I.E., C C (*) REXA(F,M(EPS)) .LE. EPS . C C IF M(EPS).LE.M THEN THE QUADRATURE QCOM(F,M(EPS)) IS COM- C PUTED. OTHERWISE, WHICH MEANS THAT EPS IS TOO SMALL WITH C RESPECT TO M, THE QUADRATURE QCOM(F,M1) IS COMPUTED, WHERE C M1=2*INT((M-1)/2)+1. IN THIS CASE EPS IS CHANGED TO THE C SMALLEST NUMBER FOR WHICH THE ESTIMATE (*) HOLDS WITH C M(EPS)=M1 FUNCTION EVALUATIONS. C C HEURISTIC TERMINATION CASE C C WE CAN EXPECT THE RELATIVE ERROR REXA TO BE OF THE C ORDER OF EPS, SEE [2]. IF EPS IS TOO SMALL WITH RESPECT C TO M THEN THE QUADRATURE QCOM(F,M) IS COMPUTED. C C ROUNDOFF ERRORS C C IN BOTH DETERMINISTIC AND HEURISTIC CASES THE ROUND- C OFF ERROR C ROFF=ABS(QEXA(F,M)-QCOM(F,M)) C CAN BE ESTIMATED BY C C (**) ROFF .LE. 3*C1*R*MACHEP, C C WHERE R=QCOM(ABS(F),M)+(1+2*C2)/3*SUM(W(I),I=1,2,...M) C AND C1 IS OF THE ORDER OF UNITY, C1=1/(1-3*MACHEP), W(I) C ARE THE WEIGHTS OF THE QUADRATURE, SEE [2], AND C2 IS C A CONSTANT ESTIMATING THE ACCURACY OF COMPUTING FUNCTION C VALUES, I.E., C ABS(EXACT(F(X))-COMPUTED(F(X))).LE.C2*MACHEP. C IF THE INTEGRAND VALUES ARE COMPUTED INACCURATELY, I.E., C C2 IS LARGE, THEN THE ESTIMATE (**) IS LARGE AND ONE CAN C EXPECT THE ACTUAL ERROR ROFF TO BE LARGE. NUMERICAL TESTS C INDICATE THAT THIS HAPPENS ESPECIALLY WHEN THE INTEGRAND C IS EVALUATED INACCURATELY NEAR A SINGULARITY. THE WAYS OF C CIRCUMVENTING SUCH PITFALLS ARE EXPLAINED IN [2]. C C REFERENCES: C C [1] SIKORSKI,K., OPTIMAL QUADRATURE ALGORITHMS IN HP C SPACES, NUM. MATH., 39, 405-410 (1982). C [2] SIKORSKI,K., STENGER,F., OPTIMAL QUADRATURES IN C HP SPACES, ACM TOMS. C C C C C INTEGER I, I1, INF, K, L, L1, M, M1, M2, N, N1 REAL A, ALFA, B, BA, C, C0, COR, D, E1, EPS, EPS3, EXPH, EXPH0, * F, H, H0 REAL H1, P, PI, QUADR, S, S1, SR, SQ2, SUM, SUM1, SUM2, U, T, V, * V0, V1 REAL V2, W, W1, W2, W3, W4 LOGICAL INF1, INF2 C PI = 4.*ATAN(1.0) C C CHECK THE INPUT DATA C IF (INF.NE.1 .AND. INF.NE.2 .AND. INF.NE.3 .AND. INF.NE.4) GO TO * 300 IF (M.LT.3) GO TO 270 IF (P.LT.1. .AND. P.NE.0.) GO TO 280 IF (P.GE.1. .AND. (D.LE.0. .OR. D.GT.PI/2.)) GO TO 280 IF (INF.EQ.4 .AND. A.GE.B) GO TO 290 C SQ2 = SQRT(2.0) I1 = INF - 2 BA = B - A N1 = 0 C C COMPUTE THE RELATIVE MACHINE PRECISION AND CHECK C THE VALUE OF EPS. CAUTION...THIS LOOP MAY NOT WORK ON A C MACHINE THAT HAS AN ACCURATED ARITHMETIC PROCESS COMPARED C TO THE STORAGE PRECISION. THE VALUE OF U MAY NEED TO BE C SIMPLY DEFINED AS THE RELATIVE ACCURACY OF STORAGE PRECISION. C U = 1. 10 U = U/10. T = 1. + U IF (1..NE.T) GO TO 10 U = U*10. IF (EPS.LT.U) EPS = U C IF (P.EQ.0.) GO TO 40 C C SET UP DATA FOR THE DETERMINISTIC TERMINATION C IF (P.EQ.1.) ALFA = 1. IF (P.GT.1.) ALFA = (P-1.)/P C = 2.*PI/(1.-1./EXP(PI*SQRT(ALFA))) + 4.**ALFA/ALFA W = ALOG(C/EPS) W1 = 1./(PI*PI*ALFA)*W*W N = INT(W1) IF (W1.GT.FLOAT(N)) N = N + 1 IF (W1.EQ.0.) N = 1 N1 = 2*N + 1 SR = SQRT(ALFA*FLOAT(N)) IF (N1.LE.M) GO TO 20 C C EPS TOO SMALL WITH RESPECT TO M. COMPUTE THE NEW EPS C GUARANTEED BY THE VALUE OF M. C N1 = 1 N = INT(FLOAT((M-1)/2)) SR = SQRT(ALFA*FLOAT(N)) M = 2*N + 1 EPS = C/EXP(PI*SR) GO TO 30 C 20 M = N1 N1 = 0 30 H = 2.*D/SR SUM2 = 0. L1 = N K = N INF1 = .FALSE. INF2 = .FALSE. H0 = H GO TO 50 C C SET UP DATA FOR THE HEURISTIC TERMINATION C 40 H = 1. H0 = 1. EPS3 = EPS/3. SR = SQRT(EPS) V1 = EPS*10. V2 = V1 M1 = M - 1 N = INT(FLOAT(M1/2)) M2 = N L1 = 0 INF1 = .TRUE. INF2 = .FALSE. C C INITIALIZE THE QUADRATURE C 50 I = 0 IF (INF.EQ.1) SUM = F(0.) IF (INF.EQ.2) SUM = F(A+1.) IF (INF.EQ.3) SUM = F(A+ALOG(1.+SQ2))/SQ2 IF (INF.EQ.4) SUM = F((A+B)/2.)/4.*BA C C COMPUTE WEIGHTS, NODES AND FUNCTION VALUES C 60 EXPH = EXP(H) EXPH0 = EXP(H0) H1 = H0 E1 = EXPH0 U = 0. COR = 0. C 70 IF (I1) 80, 90, 100 C 80 V = F(H1) H1 = H1 + H GO TO 150 C 90 V = E1*F(A+E1) E1 = E1*EXPH GO TO 150 C 100 IF (INF.EQ.4) GO TO 140 W1 = SQRT(E1+1./E1) W2 = SQRT(E1) IF (E1.LT.0.1) GO TO 110 S = ALOG(E1+W1*W2) GO TO 130 110 W3 = E1 W4 = E1*E1 C0 = 1. S = E1 S1 = E1 T = 0. 120 C0 = -C0*(0.5+T)*(2.*T+1.)/(2.*T+3.)/(T+1.) T = T + 1. W3 = W3*W4 S = S + C0*W3 IF (S.EQ.S1) GO TO 130 S1 = S GO TO 120 130 V = W2/W1*F(A+S) E1 = E1*EXPH GO TO 150 C 140 W1 = E1 + 1. V = E1/W1/W1*F((A+B*E1)/W1)*BA E1 = E1*EXPH C C SUMMATION ALGORITHM C 150 I = I + 1 SUM1 = U + V IF (ABS(U).LT.ABS(V)) GO TO 160 COR = V - (SUM1-U) + COR GO TO 170 160 COR = U - (SUM1-V) + COR 170 U = SUM1 IF (I.LT.L1) GO TO 70 C C SWITCH TO CHECK TRUNCATION CONDITION ( HEURISTIC C TERMINATION) C IF (INF1) GO TO 190 C C SWITCH TO COMPUTE THE MIDORDINATE APPROXIMATION C ( HEURISTIC TERMINATION ) OR TO STOP ( DETERMINIS- C TIC TERMINATION) C IF (INF2) GO TO 210 C C SET UP PARAMETERS TO CONTINUE SUMMATION C L1 = K 180 INF2 = .TRUE. I = 0. EXPH = 1./EXPH H0 = -H0 E1 = 1./EXPH0 H1 = H0 H = -H GO TO 70 C C TRUNCATION CONDITION C 190 V0 = V1 V1 = V2 V2 = ABS(V) IF (V0+V1+V2.LE.EPS3) GO TO 200 IF (I.LT.M2) GO TO 70 N1 = 5 200 IF (INF2) K = I IF (.NOT.INF2) L = I V1 = 10.*EPS V2 = V1 M2 = M1 - L IF (.NOT.INF2) GO TO 180 C C N1=5 - TRUNCATION CONDITION NOT SATISFIED C IF (N1.EQ.5) GO TO 260 C C TRUNCATION CONDITION SATISFIED, SUM2=TRAPEZOIDAL C APPROXIMATION C SUM2 = SUM1 + COR + SUM M2 = 2*(K+L) C C CHECK THE NUMBER OF FUNCTION EVALUATIONS C IF (M2.GT.M1) GO TO 240 C C INITIALIZE ITERATION C INF1 = .FALSE. INF2 = .FALSE. L1 = L I = 0 H = -H H0 = H/2. GO TO 60 C C P.GE.1 = DETERMINISTIC TERMINATION C 210 IF (P.GE.1.) GO TO 220 C C COMPUTE THE MIDORDINATE APPROXIMATION SUM1 C H = -H SUM1 = (SUM1+COR)*H W1 = (SUM1+SUM2)/2. C C TERMINATION CONDITION C IF (ABS(SUM1-SUM2).LE.SR) GO TO 230 C C SET UP DATA FOR THE NEXT ITERATION C M2 = 2*M2 IF (M2.GT.M1) GO TO 250 I = 0 K = 2*K L = 2*L L1 = L H = H/2. H0 = H/2. SUM2 = W1 INF2 = .FALSE. GO TO 60 C C FINAL RESULTS C 220 QUADR = -H*(SUM1+COR+SUM) INF = N1 RETURN C 230 QUADR = W1 INF = 2 M = M2 + 1 RETURN C 240 QUADR = SUM2 INF = 3 M = K + L + 1 RETURN C 250 QUADR = W1 INF = 3 M = M2/2 + 1 RETURN C 260 QUADR = U + COR + SUM INF = 4 M = K + L + 1 RETURN C 270 INF = 10 RETURN C 280 INF = 11 RETURN C 290 INF = 12 RETURN C 300 INF = 13 RETURN C END C PROGRAM DRIVER MAN 10 C DRIVING PROGRAM FOR INTHP MAN 20 C MAN 30 C TEST OF 8 INTEGRANDS REPRESENTING CLASSES OF PROBLEMS MAN 40 C THAT INTHP SOLVES. EACH INTEGRAND FUNCTION INCLUDES MAN 50 C COMMENTS ON THE CLASS OF PROBLEMS IT BELONGS TO. MAN 60 C THE VALUES OF PARAMETERS A,B,D,M,P,EPS,INF ARE MAN 70 C READ FROM THE DATA FILES FOR33.DAT AND FOR34.DAT. MAN 80 C MAN 90 INTEGER I, J, M, INF, NTIN, NTOUT MAN 100 REAL A, B, D, P, EPS, QUADR MAN 110 EXTERNAL F1, F2, F3, F4, F5, F6, F7, F8 MAN 120 C MAN 130 NTIN = 32 MAN 140 NTOUT = 6 MAN 150 C NTIN, NTOUT SPECIFY INPUT OUTPUT DEVICES MAN 160 C MAN 170 C MAIN LOOP WITH SWITCHES TO CALL INTHP MAN 180 C MAN 190 DO 110 J=1,2 MAN 200 C MAN 210 C FIRST EXECUTION (J=1) WITH DETERMINISTIC TERMINATION - MAN 220 C NTIN=33. MAN 230 C SECOND EXECUTION (J=2) WITH HEURISTIC TERMINATION - MAN 240 C NTIN=34. MAN 250 C MAN 260 NTIN = NTIN + 1 MAN 270 C MAN 280 DO 100 I=1,8 MAN 290 C MAN 300 WRITE (NTOUT,99999) I MAN 310 C MAN 320 C INPUT PARAMETERS A,B,D,M,P,EPS,INF MAN 330 C MAN 340 READ (NTIN,99998) A, B, D, M, P, EPS, INF MAN 350 WRITE (NTOUT,99995) MAN 360 WRITE (NTOUT,99994) A, B, D, M MAN 370 WRITE (NTOUT,99993) P, EPS, INF MAN 380 C MAN 390 C SWITCHES TO CALL SUBROUTINE INTHP MAN 400 C MAN 410 IF (I.EQ.1) GO TO 10 MAN 420 IF (I.EQ.2) GO TO 20 MAN 430 IF (I.EQ.3) GO TO 30 MAN 440 IF (I.EQ.4) GO TO 40 MAN 450 IF (I.EQ.5) GO TO 50 MAN 460 IF (I.EQ.6) GO TO 60 MAN 470 IF (I.EQ.7) GO TO 70 MAN 480 IF (I.EQ.8) GO TO 80 MAN 490 10 CALL INTHP(A, B, D, F1, M, P, EPS, INF, QUADR) MAN 500 GO TO 90 MAN 510 20 CALL INTHP(A, B, D, F2, M, P, EPS, INF, QUADR) MAN 520 GO TO 90 MAN 530 30 CALL INTHP(A, B, D, F3, M, P, EPS, INF, QUADR) MAN 540 GO TO 90 MAN 550 40 CALL INTHP(A, B, D, F4, M, P, EPS, INF, QUADR) MAN 560 GO TO 90 MAN 570 50 CALL INTHP(A, B, D, F5, M, P, EPS, INF, QUADR) MAN 580 GO TO 90 MAN 590 60 CALL INTHP(A, B, D, F6, M, P, EPS, INF, QUADR) MAN 600 GO TO 90 MAN 610 70 CALL INTHP(A, B, D, F7, M, P, EPS, INF, QUADR) MAN 620 GO TO 90 MAN 630 80 CALL INTHP(A, B, D, F8, M, P, EPS, INF, QUADR) MAN 640 C MAN 650 C FINAL RESULTS MAN 660 C MAN 670 90 WRITE (NTOUT,99992) MAN 680 WRITE (NTOUT,99997) INF, QUADR, EPS MAN 690 WRITE (NTOUT,99996) M MAN 700 C MAN 710 100 CONTINUE MAN 720 110 CONTINUE MAN 730 C MAN 740 99999 FORMAT (//////12H TEST OF THE, I3, 12H - INTEGRAND) MAN 750 99998 FORMAT (2E6.0, F6.4, I7, 2E6.0, I1) MAN 760 99997 FORMAT (5H INF=, I3, 7H QUADR=, E16.8, 5H EPS=, E16.8) MAN 770 99996 FORMAT (5H M=, I7) MAN 780 99995 FORMAT (/28H INPUT VALUES OF PARAMETERS ) MAN 790 99994 FORMAT (/3H A=, E6.0/3H B=, E6.0/3H D=, F6.4/3H M=, I7) MAN 800 99993 FORMAT (/3H P=, E7.1/5H EPS=, E7.1/5H INF=, I1/) MAN 810 99992 FORMAT (/32H OUTPUT VALUES OF PARAMETERS ) MAN 820 C MAN 830 STOP MAN 840 END MAN 850 FUNCTION F1(X) F1 10 REAL X C C F1(X)=EXP(-X*X) C C F1 BELONGS TO THE SPACE H(P,DD) WITH P=+INFINITY C AND DD=STRIP(Z:ABS(IM(Z)).LT.PI/4). F1 DECREASES C EXPOTENTIALLY AS ABS(X) GOES TO INFINITY. C F1 = EXP(-X*X) RETURN END FUNCTION F2(X) F2 10 REAL X C C F2(X)=2./(SQRT(X)*(X+1)) C C F2 BELONGS TO THE SPACE H(P,DD) WITH P=2 AND C DD=LENS REGION(Z:ABS(ARG(Z/(1-Z)).LT.PI/2). C F2 HAS ALGEBRAIC TYPE SINGULARITY AT X=0. C F2 = 2./SQRT(X)/(X+1.) RETURN END FUNCTION F3(X) F3 10 REAL X C C F3(X)=1/(SQRT(X)*EXP(X)) C C F3 BELONGS TO THE SPACE H(P,DD) WITH P=2 AND C DD=REGION(Z:ABS(ARG(SINH(Z))).LT.PI/2). F3 HAS C ALGEBRAIC TYPE SINGULARITY AT X=0 AND DECREASES C EXPOTENTIALLY AS X GOES TO +INFINITY. C F3 = 1./SQRT(X)/EXP(X) RETURN END FUNCTION F4(X) F4 10 REAL X C C F4(X)=1/(SQRT(X)*(X+1)) C C F4 BELONGS TO THE SPACE H(P,DD) WITH P=2 AND C DD=SECTOR(Z:ABS(ARG(Z)).LT.PI/2). F4 HAS ALGEBRAIC C TYPE SINGULARITY AT X=0 AND DECREASES WITH ALGEBRAIC C RATE AS X GOES TO +INFINITY. C F4 = 1./SQRT(X)/(X+1.) RETURN END FUNCTION F5(X) F5 10 REAL F, S, T, X, Z, S1, Z1, SIG C C F5(X)=LOG(1-SIN(X)/X)/EXP(X) C C F5 BELONGS TO THE SPACE H(P,DD) WITH P=100 AND C DD=REGION(Z:ABS(ARG(SINH(Z))).LT.PI/2). F5 IS AN C OSCILLATORY INTEGRAND, HAS LOGARITHMIC TYPE SINGULARITY C AT X=0 AND DECREASES EXPOTENTIALLY AS X GOES TO +INFINITY. C IF (X.LE.0.1) GO TO 10 S = 1. - SIN(X)/X GO TO 30 10 SIG = 1. Z1 = X*X Z = Z1 S = 0. S1 = 0. T = 3. F = 1. 20 F = F*T*(T-1.) S = S + SIG*Z/F IF (S.EQ.S1) GO TO 30 S1 = S Z = Z*Z1 SIG = -SIG T = T + 2. GO TO 20 30 F5 = ALOG(S)/EXP(X) RETURN END FUNCTION F6(X) F6 10 REAL X C C F6(X)=2/SQRT(3-2*X-X*X) C C F6 BELONGS TO THE SPACE H(P,DD) WITH P=2 AND C DD=LENS REGION(Z:ABS(ARG((Z+1)/(1-Z))).LT.PI/2). F6 HAS C ALGEBRAIC TYPE SINGULARITY AT X=1. C THIS IS AN EXAMPLE FOR LOSS OF SIGNIFICANT FIGURES DUE C TO INACCURATE EVALUATION OF AN INTEGRAND NEAR THE SINGU- C LARITY-SEE [2] AND INTEGRAND F7. C F6 = 2./SQRT(3.-2.*X-X*X) RETURN END FUNCTION F7(X) F7 10 REAL X C C F7(X)=2/SQRT(4*X-X*X) C C F7 BELONGS TO THE SPACE H(P,DD) WITH P=2 AND C DD=LENS REGION(Z:ABS(ARG(Z/(2-Z)).LT.PI/2). F7 HAS C ALGEBRAIC TYPE SINGULARITY AT X=0. C TO OVERCOME NUMERICAL PROBLEMS IN COMPUTING F6 NEAR C THE SINGULARITY WE REPLACE X IN F6 BY 1-T AND GET F7. C NOW WE CAN ACHIEVE ALL SIGNIFICANT FIGURES-SEE [2]. C F7 = 2./SQRT(4.*X-X*X) RETURN END FUNCTION F8(X) F8 10 REAL X, Y C C F8(X)=1/(EXP(X/2)+EXP(-X/2)) C C F8 BELONGS TO THE SPACE H(P,DD) WITH P=2 AND C DD=STRIP(Z:ABS(IM(Z)).LT.PI/2). F8 DECREASES EXPOTENTIALLY C AS ABS(X) GOES TO INFINITY. THE RATE OF DECREASE IS SLOWER C THAN THAT OF F1. C Y = ABS(X/2.) Y = EXP(Y) F8 = 1./(Y+1./Y) RETURN END 0.E0000.E0000.785400100001.E0001.E-101 0.E0001.E0001.570600100002.E0001.E-104 0.E0000.E0001.570600100002.E0001.E-103 0.E0000.E0001.570600100002.E0001.E-102 0.E0000.E0001.570600100001.E0021.E-103 -1.E001.E0001.570600100002.E0005.E-034 0.E0002.E0001.570600100002.E0001.E-104 0.E0000.E0001.570600100002.E0001.E-101 TEST OF THE 1 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B=0. D= .7854 M= 10000 P= .1E+01 EPS= .1E-09 INF=1 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .17724539E+01 EPS= .10000000E-09 M= 133 TEST OF THE 2 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B= E+01 D=1.5706 M= 10000 P= .2E+01 EPS= .1E-09 INF=4 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .31415927E+01 EPS= .10000000E-09 M= 265 TEST OF THE 3 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B=0. D=1.5706 M= 10000 P= .2E+01 EPS= .1E-09 INF=3 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .17724539E+01 EPS= .10000000E-09 M= 265 TEST OF THE 4 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B=0. D=1.5706 M= 10000 P= .2E+01 EPS= .1E-09 INF=2 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .31415927E+01 EPS= .10000000E-09 M= 265 TEST OF THE 5 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B=0. D=1.5706 M= 10000 P= .1E+03 EPS= .1E-09 INF=3 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= -.30456893E+01 EPS= .10000000E-09 M= 133 TEST OF THE 6 - INTEGRAND INPUT VALUES OF PARAMETERS A= -E+01 B= E+01 D=1.5706 M= 10000 P= .2E+01 EPS= .5E-02 INF=4 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .31409151E+01 EPS= .50000000E-02 M= 27 TEST OF THE 7 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B= E+01 D=1.5706 M= 10000 P= .2E+01 EPS= .1E-09 INF=4 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .31415927E+01 EPS= .10000000E-09 M= 265 TEST OF THE 8 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B=0. D=1.5706 M= 10000 P= .2E+01 EPS= .1E-09 INF=1 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .31415927E+01 EPS= .10000000E-09 M= 265 TEST OF THE 1 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B=0. D= .7854 M= 10000 P= .1E+01 EPS= .1E-09 INF=1 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .17724539E+01 EPS= .10000000E-09 M= 133 TEST OF THE 2 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B= E+01 D=1.5706 M= 10000 P= .2E+01 EPS= .1E-09 INF=4 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .31415927E+01 EPS= .10000000E-09 M= 265 TEST OF THE 3 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B=0. D=1.5706 M= 10000 P= .2E+01 EPS= .1E-09 INF=3 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .17724539E+01 EPS= .10000000E-09 M= 265 TEST OF THE 4 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B=0. D=1.5706 M= 10000 P= .2E+01 EPS= .1E-09 INF=2 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .31415927E+01 EPS= .10000000E-09 M= 265 TEST OF THE 5 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B=0. D=1.5706 M= 10000 P= .1E+03 EPS= .1E-09 INF=3 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= -.30456893E+01 EPS= .10000000E-09 M= 133 TEST OF THE 6 - INTEGRAND INPUT VALUES OF PARAMETERS A= -E+01 B= E+01 D=1.5706 M= 10000 P= .2E+01 EPS= .5E-02 INF=4 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .31409151E+01 EPS= .50000000E-02 M= 27 TEST OF THE 7 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B= E+01 D=1.5706 M= 10000 P= .2E+01 EPS= .1E-09 INF=4 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .31415927E+01 EPS= .10000000E-09 M= 265 TEST OF THE 8 - INTEGRAND INPUT VALUES OF PARAMETERS A=0. B=0. D=1.5706 M= 10000 P= .2E+01 EPS= .1E-09 INF=1 OUTPUT VALUES OF PARAMETERS INF= 0 QUADR= .31415927E+01 EPS= .10000000E-09 M= 265