I’m calculating some integrals that include the function

The full integral will be performed using Gauss-Laguerre quadrature, with the following set of abscissa’s:

```
real(dp), parameter :: x(32) = [&
0.0444893658332670184_dp, &
0.234526109519618537_dp, &
0.576884629301886426_dp, &
1.07244875381781763_dp, &
1.72240877644464544_dp, &
2.52833670642579488_dp, &
3.49221327302199449_dp, &
4.61645676974976739_dp, &
5.90395850417424395_dp, &
7.35812673318624111_dp, &
8.98294092421259610_dp, &
10.7830186325399721_dp, &
12.7636979867427251_dp, &
14.9311397555225573_dp, &
17.2924543367153148_dp, &
19.8558609403360547_dp, &
22.6308890131967745_dp, &
25.6286360224592478_dp, &
28.8621018163234747_dp, &
32.3466291539647370_dp, &
36.1004948057519738_dp, &
40.1457197715394415_dp, &
44.5092079957549380_dp, &
49.2243949873086392_dp, &
54.3337213333969073_dp, &
59.8925091621340182_dp, &
65.9753772879350528_dp, &
72.6876280906627086_dp, &
80.1874469779135231_dp, &
88.7353404178923987_dp, &
98.8295428682839726_dp, &
111.751398097937695_dp]
```

Unfortunately, for larger x I start getting NaN’s. This is due to an arithmetic underflow in the calculation of `erfc`

, which can be caught at compile time if the array of abscissas is given the parameter attribute.

Later I discovered the Fortran intrinsic `erfc_scaled`

sometimes abbreviated as \mathrm{erfcx} which is equal to:

Using some logarithms, it’s possible to compute f(x) as the expression

I’ll call this method A; in terms of Fortran code this becomes:

```
elemental function logerfc(x) result(y)
real(dp), intent(in) :: x
real(dp) :: y
y = log(erfc_scaled(x)) - x**2
end function
! Method A
elemental function f(x) result(y)
real(dp), intent(in) :: x
real(dp) :: y
y = exp(-x**2 - logerfc(x))
end function
```

Another alternative is to compute:

which I’ll call method B, in Fortran code

```
! Method B
elemental function f(x) result(y)
real(dp), intent(in) :: x
real(dp) :: y
y = 1.0_dp/erfc_scaled(x)
end function
```

I’ve tested this and found that A and B agree very well at smaller values, however at the larger quadrature nodes, the differences are of order 1.0E-11. (One caveat, actually I’m computing f(x - x_c), where x_c \approx 7)

Which computation method should I prefer, A or B?

Should I perform a comparison of A and B with a quadruple precision calculation of the original function and take the method with the lower error?

Edit: I am using gfortran v10.3, presumably linked with glibc v2.27.

Edit 2: Full code is provided in the details box below

## Details

```
module lag32
implicit none
private
public :: x
integer, parameter :: dp = kind(1.0d0)
real(dp), parameter :: x(32) = [&
0.0444893658332670184_dp, &
0.234526109519618537_dp, &
0.576884629301886426_dp, &
1.07244875381781763_dp, &
1.72240877644464544_dp, &
2.52833670642579488_dp, &
3.49221327302199449_dp, &
4.61645676974976739_dp, &
5.90395850417424395_dp, &
7.35812673318624111_dp, &
8.98294092421259610_dp, &
10.7830186325399721_dp, &
12.7636979867427251_dp, &
14.9311397555225573_dp, &
17.2924543367153148_dp, &
19.8558609403360547_dp, &
22.6308890131967745_dp, &
25.6286360224592478_dp, &
28.8621018163234747_dp, &
32.3466291539647370_dp, &
36.1004948057519738_dp, &
40.1457197715394415_dp, &
44.5092079957549380_dp, &
49.2243949873086392_dp, &
54.3337213333969073_dp, &
59.8925091621340182_dp, &
65.9753772879350528_dp, &
72.6876280906627086_dp, &
80.1874469779135231_dp, &
88.7353404178923987_dp, &
98.8295428682839726_dp, &
111.751398097937695_dp]
end module
program main
use lag32, only: x
implicit none
integer, parameter :: dp = kind(1.0d0)
real(dp), parameter :: mc = 3.0e-12_dp
real(dp), parameter :: epsilon = 3*sqrt(2.0_dp)*1.0E-13_dp
real(dp), parameter :: xc = mc/epsilon
real(dp), allocatable :: ya(:), yb(:)
integer :: i
! Original
! y = exp(-(x - xc)**2)/erfc(x - xc)
! Method A
ya = exp(-(x - xc)**2 - logerfc(x - xc))
! Method B
yb = 1.0_dp/erfc_scaled(x - xc)
do i = 1, size(x)
print *, x(i), ya(i) ,ya(i) - yb(i)
end do
contains
elemental function logerfc(x) result(y)
real(dp), intent(in) :: x
real(dp) :: y
y = log(erfc_scaled(x)) - x**2
end function
end program
```