# The recursive adaptive quadrature in MS Fortran-77

## 1 Introduction

As it was shown in  the application of recursion makes it possible to create compact, explicit and effective integration programs. In the mentioned papers the C++ version of such a routine is presented. However, it is historically formed that a large number of science and engineering Fortran-77 codes have been accumulated by now in the form of applied libraries and packages. That is one of the reasons why Fortran-77 is still quite popular in the applied programming. From this standpoint it seems to be very useful to use such an effective recursive integration algorithm in Fortran-77. There exist at least two possibilities to realize it. The first one is described in  where the interface for calling mentioned C++ recursive integration function from MS Fortran-77 is presented. The second possibility consist in constructing the recursive subroutine by means of Fortran-77 only. This is the particular subject of the paper where the possibility and benefits of recursion strategy in MS Fortran-77 is discussed.

## 2 Recursion in MS Fortran-77

The direct transformation of the mentioned C++ code is not possible mainly due to the formal inhibition of the recursion in Fortran-77. However, Microsoft extensions of Fortran-77 (e.g. MS Fortran V.5.0, Fortran Power Station) allow to make indirect recursive calls. It means that subprogram can call itself through intermediate subprogram. If anybody doubts he can immediately try:

|  |
| --- |
| call rec(1.0) |
| end |
| subroutine rec(hh) |
| integer i/0/ |
| i = i + 1 |
| h = 0.5\*hh |
| write(\*,\*) i, h |
| if (i.lt.3) call mediator(h) |
| write(\*,\*) i, h |
| end |
| subroutine mediator(h) |
| call rec(h) |
| end |

and get the following results:

|  |  |
| --- | --- |
| 1 | 5.000000E-01 |
| 2 | 2.500000E-01 |
| 3 | 1.250000E-01 |
| 3 | 1.250000E-01 |
| 3 | 1.250000E-01 |
| 3 | 1.250000E-01 |

But this is not a true recursion because no mechanism is supplied for restoring the values of the internal variables of the subroutine after its returning from recursion. The last requirement can be fulfilled by the forced storing of the internal variables into the program stack. The AUTOMATIC description of variables provides such possibility in MS Fortran-77. Taking this into account, the above example can be rewritten:

|  |
| --- |
| call rec(1.0) |
| end |
| subroutine rec(hh) |
| integer i/0/ |
| automatic h, i |
| i = i + 1 |
| h = 0.5\*hh |
| write(\*,\*) i, h |
| if (i.lt.3) call mediator(h) |
| write(\*,\*) i, h |
| end |
| subroutine mediator(h) |
| call rec(h) |
| end |

that yields:

|  |  |
| --- | --- |
| 1 | 5.000000E-01 |
| 2 | 2.500000E-01 |
| 3 | 1.250000E-01 |
| 3 | 1.250000E-01 |
| 3 | 2.500000E-01 |
| 3 | 5.000000E-01 |

Here the values of h are restored after each returning from recursion because it is saved in the stack before the recursive call. Note, that although the i variable is described as AUTOMATIC nonetheless its value is not saved.

## 3 Recursive adaptive quadrature algorithm

The described possibilities allow to employ effective recursion strategy for creating adaptive quadrature subroutine in MS Fortran-77.

The presented algorithm consists of two independent parts: adaptive subroutine and quadrature formula. The adaptive subroutine uses recursive algorithm to implement standard bisection method (see fig.1). For reaching desired relative accuracy $`\epsilon `$ of the integration the integral estimation $`I_{whole}`$ over \[$`a_i`$,$`b_i`$\] subinterval on the i-th step of bisection is compared with the sum of $`I_{left}`$ and $`I_{right}`$ integral values that are evaluated over left and right halves of the considered subinterval. The comparison rule was chosen in the form:

$$I_{left}+I_{right}I_{whole}\epsilon I,$$
(1)   
where $`I`$ denotes the integral sum over whole integration interval \[a,b\]. The value of $`I`$ is accumulated and adjusted on each step of bisection.

Should (1) be not fulfilled the adaptive procedure is called recursively for both (left and right) subintervals. Evaluation of the integral sums on each step of bisection is performed by means of quadrature formula. There are no restrictions on the type of quadratures used during integration. This makes the code to be very flexible and applicable to a wide range of integration problems.

The form (1) of the chosen comparison rule does not pretend on effectiveness rather on simplicity and generality. Really it seems to be very common and does not depend on the integrand as well as quadrature type. At the same time the use of (1) in some cases can result in overestimation of the calculated integral that consequently leads to more integrand function calls. One certainly can get some gains using, for instance, definite quadratures with different number or/and equidistant points or Gauss-Kronrod quadrature  etc. The comparison rule in the later cases becomes more effective but complex, intricate and sometimes less common. Whatever the case, the choice of comparison rule as well as the problems connected with it lie outside the subject of the publication.

Let us note some advantages which the application of the recursive call ideology to numerical integration can reveal:

* Very simple and evident algorithm that could result in extremely short as well as easy for further modifications and possible enhancements adaptive code.
* Because of the indicated shortness the adaptive procedure’s own running time has to be diminutive. That could result in its better performance compared to the known programs especially in the cases when the integrand function calculations are not time consuming.
* There is no need to store the integrand function values and pay attention on their correct usage. Besides no longer the control of subinterval bounds is in need. Indicated features permit utmost reduction of the efforts that one has to pay while creating the adaptive code.
* Nothing but program’s stack size sets the restriction on the number of subintervals when the recursive procedure is used (see next section). At the same time for the existing programs the crush level of the primary interval is strictly limited by the dimensions of the static arrays.

## 4 Program realization

Fortran-77 version of adaptive subroutine practically coincides with the corresponding C++ QUADREC (Quadrature used Adaptively and Recursively) function :

```
SUBROUTINE Quadrec(Fun,Left,Right,Estimation)
real*8 Fun, Left, Right, Estimation
real*8 Eps, Result, Section, SumLeft, SumRight, QuadRule
integer*4 RecMax, RecCur, RawInt
common /IP/ Eps, Result, RecMax, RawInt, RecCur
automatic SumLeft, SumRight, Section
external Fun
RecCur = RecCur+1
if (RecCur.le.RecMax) then
  Section = 0.5d0*(Left+Right)
  SumLeft = QuadRule(Fun,Left,Section)
  SumRight = QuadRule(Fun,Section,Right)
  Result = Result+SumLeft+SumRight-Estimation
  if (dabs(SumLeft+SumRight-Estimation).gt.Eps*dabs(Result)) then
    call Mediator(Fun,Left,Section,SumLeft)
    call Mediator(Fun,Section,Right,SumRight)
  end if
else
  RawInt = RawInt+1
end if
RecCur = RecCur-1
return
end

```

Note that subroutine contains only eleven executable statements. The integrand function name, left and right bounds of the integration interval as well as the initial estimation of the integral value are the formal arguments of the subroutine. The IP common block contains the following variables: desired relative accuracy (Eps), the result of the integration (Result), maximum and current levels of recursion (RecMax, RecCur) as well as raw (not processed during integration) subintervals (RawInt).

The Section variable is used for storing a value of midpoint of the current subinterval. The integral sums over its left and right halves are estimated by QuadRule external function and stored in LeftSum and RightSum variables. The last three variables are declared as AUTOMATIC allowing to preserve their values from changing and use them after returning from recursion.

Execution of the subroutine begins with increasing of recursion level counter. If its value does not exceed RecMax, the integral sums over left and right halves of the current subinterval are evaluated, the integration result is updated and accuracy of the integration is checked. If achieved accuracy is not sufficient then subprogram calls itself (with the help of mediator subroutine) over left and right halves of subinterval. In the case the accuracy condition is satisfied the recursion counter decreases and subprogram returns to the previous level of recursion. The number of raw subintervals is increased when desired accuracy is not reached and RecCur is equal to RecMax.

The mediator subroutine has only one executable statement:

```
      subroutine Mediator(Fun,Left,Right,Estimation)
      real*8 Estimation, Fun, Left, Right
      external Fun
      call Quadrec(Fun,Left,Right,Estimation)
      return
      end

```

The main part of the integration program can be look like:

```
      common   /IP/ Eps, Result, RecMax, RawInt, RecCur
      common   /XW/ X, W, N
      integer  RecMax, RawInt, RecCur
      real*8   X(100), W(100), Left/0.0d0/, Right/1.0d0/
      real*8   Result, Eps, QuadRule
      real*8   Integrand
      external Integrand
      Eps = 1.0d-14
      N = 10
      RecCur = 0
      RecMax = 10
      call gauleg(-1.0d0,1.0d0,X,W,N)
      Result = QuadRule(Integrand,Left,Right)
      call Quadrec(Integrand,Left,Right,Result)
      write(*,*) ’ Result = ’,Result,’  RawInt = ’,RawInt
      end

```

The common block XW contains Gaussian abscissas and weights which are calculated with the help of gauleg subroutine for a given number of points N. The text of subroutine, reproduced from , is presented in Appendix A.

The text of QuadRule function is presented below:

```
      real*8 function QuadRule(Integrand,Left,Right)
      common /XW/ X, W, N
      real*8 X(100),W(100),IntSum,Abscissa,Left,Right,Integrand
      IntSum = 0.0d0
      do 1 i = 1, N
      Abscissa = 0.5d0*(Right+Left+(Right-Left)*X(i))
1     IntSum = IntSum + W(i)*Integrand(Abscissa)
      QuadRule = 0.5d0*IntSum*(Right-Left)
      return
      end

```

It is important to note that the number of recursive calls is limited by the size the program stack. This fact obviously sets the limit on the reachable number of the primary integration interval bisections and consequently restricts the integration accuracy. Note that stack size of the program can be enlarged by using /Fxxxx option of MS Fortran-77 compiler.

## 5 Numerical tests

The program testing was performed on four different integrals. In each case the exact values can be found analytically. That made it possible to control the desired and reached accuracy of the integration. Besides the same integrals were obtained with the help of well-known adaptive program QUANC8 reproduced from . It allowed to compare the number of integrand function calls and the number of raw intervals for both programs.

The presented comparison has merely the aim to show that the use of recursion allows to construct very short and simple adaptive quadrature code that is not inferior to such a sophisticated program as QUANC8. Meanwhile the direct comparison of these programs seems to be incorrect because of a number of reasons.

The Newton-Cotes equidistant quadrature formula which is used in QUANC8 allows to make reuse of integrand function values calculated in the previous steps of bisection. That is the reason why QUANC8 has to have higher performance in integrand function calls compared to adaptive programs that use quadratures with non-equidistant points. Since QUADREC is not oriented on the use of definite but any quadrature formula it can be specified as a program of the later type.

At the same time QUANC8 gives bad results for functions with unlimitedly growing derivative and does not work at all for functions that go to infinity at the either of the integration interval endpoints. There are none of the indicated restrictions in QUADREC. Furthermore the opportunity of choosing of quadrature type makes it to be a very flexible tool for integration. Here QUADREC gives a chance to choose quadrature which is the most appropriate to the task (see Section 6).

For integrals in sections 5.1 and 5.2 the optimal numbers of quadrature points were found and used for integration. The 24-point quadrature was applied for integration in sections 5.3 and 5.4.

### 5.1 Sharp peaks at a smooth background

Let us start with the calculation of the integral cited in :

$$\underset{0}{\overset{1}{}}\left(\frac{1}{(xa_1)^2+a_2}+\frac{1}{(xb_1)^2+b_2}c_0\right)𝑑x$$
(2)   
The integrand is the sum of two Lorenz type peaks and a constant background. At the beginning, values of $`a_1`$, $`a_2`$, $`b_1`$, $`b_2`$ and $`c_0`$ parameters were chosen to be the same as in the cited work. Then test was conducted at decreasing values of $`a_2`$ and $`b_2`$, which determine width of the peaks, while both programs satisfied desired accuracy and did not signal about raw subintervals. The results of the test when $`a_2=b_2=10^8`$ are presented in Table 1. Note that only the optimal values are given for QUADREC program. The corresponding optimal numbers of Gaussian quadrature points are indicated.

As it follows from the given data the number of integrand function calls are compared for the both programs in the wide range of desired accuracy. Meanwhile it is interesting to point the attention on the fact that the use of optimal quadratures can give definite profits in the reached accuracy and number of integrand function calls even when the simple comparison rule (1) and no reuse of the integrand values are applied.

At further decreasing of $`a_2`$ and $`b_2`$ parameters QUANC8 informed about raw intervals and did not satisfy desired accuracy. At the same time QUADREC gave correct results for the integral down to parameter values $`a_2=b_2=10^{19}`$. This is mainly due to the differences between static and dynamic memory allocation ideologies which are used in QUANC8 and QUADREC respectively.

### 5.2 Integration at a given absolute accuracy

The next test concerns the integration of the function:

$$f(x)=\frac{(1\alpha x)\mathrm{exp}(\alpha x)}{x^2\mathrm{exp}(2\alpha x)+b^2}$$
(3)   
over whole positive real axes. It is easy to show that its exact value is equal to zero. For reducing the integration interval to  the evident substitution x=t/(1-t) was used. As far as absolute accuracy of the integration was required the fourteenth line in the listed above QUADREC function text was changed to:

`if (dabs(SumLeft+SumRight-Estimation).gt.Eps) then`

The results of the test are presented in Table 2.

### 5.3 Improper integral

The next integrand function:

$$f(x)=x^{\frac{1}{n}1}$$
(4)   
becomes nonintegrable in the  interval when n goes to infinity. Besides function 4 goes to infinity at the low integration limit. That is the reason why QUANC8 can not be directly applied to the problem. To have still the opportunity of comparison, the integration of indicated function was performed from $`10^{10}`$ to 1. The results of the test for number n up to and including 20 and desired relative accuracy of $`10^{14}`$ are listed in Table 3. The second and fifth columns give the number of intervals that were not processed during the integration by both routines. The number of integrand function calls and values of reached relative accuracy are also presented in the table.

Such a low performance of QUANC8, as it follows from the given data, can be explained by the type of comparison rule which it exploits. Really it is not applicable to functions like 4 that have unlimitedly growing derivative in some point(s) of integration interval.

The results of integration over whole interval  by QUADREC are shown in Table 4. The last column contains the maximum recursion level needed to achieve the desired accuracy. As far as it turned out to be a rather big the program stack size was correspondingly increased.

### 5.4 Evidence of program adaptation

Finally, integration of:

$$f(x)=sin(Mx)$$
(5)   
over \[0,2$`\pi `$\] interval is appropriate for demonstrating the ability of an integration routine for adaptation. The exact integral value evidently is equal to zero for any integer M.

During the test a number of large simple integers were assigned to M and absolute integration accuracy of $`10^{10}`$ was required. The output of the test is presented in Table 5. The accommodation of the routine evidently follows from the data. Namely, for M’s going from $`10^5`$ to $`1.210^6`$ the maximum recursion level (integrand function call number) changes from 15 ($`1.610^6`$) to 18 ($`1.210^7`$).

For all chosen values of M the desired accuracy was fulfilled. Furthermore program succeeded down to the accuracy of $`10^{12}`$. Note that the standard stack size was used during the test and higher performance of the QUADREC certainly could be reached had the stack size been enlarged.

## 6 The optimal quadrature

Here we want to take note of the possibility to minimize number of integrand calls by choosing the quadrature with optimal number of abscissas. The fig.2 demonstrates such a possibility.

Particularly, circles present the number of integrand calls versus the number of Gaussian abscissas used for the integration of 2 over  and desired relative accuracy of $`10^{10}`$ and $`a_1=0.3`$, $`a_2=10^9`$, $`b_1=0.9`$, $`b_2=4.010^9`$, $`c_0=6`$. From the presented data one can see that the optimal number of Gaussian abscissas ranges approximately from 7 to 17. The use of more Gaussian abscissas leads to the linear growth of the number of function calls, because it does not result in essential reduction of the integration interval fragmentation. From the other hand the use of less number of Gaussian abscissas results in the significant growth of the number of function calls due to the extremely high fragmentation required. Dependence of the maximum recursion level upon the number of abscissas, presented by triangles, confirms this consideration.

Note that despite significant differences of the testing integrals (2,3) the optimal number of Gaussian abscissas turned out to be in the above limits.

## 7 Conclusion

Thus, the indirect recursion combined with AUTOMATIC variable description allow to employ true recursion mechanism in MS Fortran-77. In particular, the recursion strategy was applied to create effective adaptive quadrature code. Despite the simplicity and extremely small program body it showed good results on rather complex testing integrals. The created subroutine is very flexible and applicable to a wide range of integration problems. In particular, it was applied for constructing effective Hilbert transformation program. The last one was used to restore frequency dependence of refraction coefficient in analysis of optical properties of complex organic compounds. The subroutine can be easily incorporated into existing Fortran programs. Note that the coding trick, described in the paper, is very convenient for constructing multidimensional adaptive quadrature programs.

## 8 Acknowledgments

We express our thanks to Dr.V.K.Basenko for stimulating and useful discussion of the problem.

Appendix

## A Guass-Legendr weights and quadratures

```
      subroutine gauleg(x1,x2,x,w,n)
      integer n
      real*8 x1,x2,x(n),w(n)
      real*8  eps
      parameter (eps=3.d-14)
      integer i,j,m
      real*8  p1,p2,p3,pp,xl,xm,z,z1
      m=(n+1)/2
      xm=0.5d0*(x2+x1)
      xl=0.5d0*(x2-x1)
      do 12 i=1,m
        z=cos(3.141592654d0*(i-.25d0)/(n+.5d0))
1       continue
          p1=1.d0
          p2=0.d0
          do 11 j=1,n
            p3=p2
            p2=p1
            p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j
11        continue
          pp=n*(z*p1-p2)/(z*z-1.d0)
          z1=z
          z=z1-p1/pp
        if(abs(z-z1).gt.eps)goto 1
        x(i)=xm-xl*z
        x(n+1-i)=xm+xl*z
        w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
        w(n+1-i)=w(i)
12    continue
      return
      end

```