# Accurate computation of exp(-x\*\*2)/erfc(x)

**URL:** <https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360>\
**Category:** Help\
**Created:** [December 2, 2021, 12:20pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360 "2021-12-02T12:20:30Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [December 2, 2021, 12:20pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/1 "2021-12-02T12:20:30Z")

</div>

I’m calculating some integrals that include the function

f(x) = \frac{e^{-x^2}}{\mathrm{erfc}(x)}

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

```fortran
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:

\mathrm{erfcx(x)} = e^{x^2} \mathrm{erfc} (x)

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

f(x) = \exp \left(-x^2 - \underbrace{\ln(\mathrm{erfc}(x))}\_{\ln(\mathrm{erfcx}(x)) - x^2}\right)

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

```fortran
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:

f(x) = \frac{e^{x^2}}{e^{x^2}} \frac{e^{-x^2}}{\mathrm{erfc}(x)} = \frac{1}{\mathrm{erfcx}(x)}

which I’ll call method B, in Fortran code

```fortran
! 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**
>
> ```fortran
> 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
> 
> ```

---

<div class="post-metadata">

**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\
**Post date:** [December 2, 2021, 2:55pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/2 "2021-12-02T14:55:51Z")

</div>

You did not state which integral you want to evaluate, but it may be that the well-known asymptotic behavior of the complementary error function may be of use.

The expression _√π x_ exp(_x_^2 ) erfc(_x_) approaches 1 as _x_ approaches infinity, and there is an asymptotic series in terms of (1/2_x_^2) for this expression, see for example [https://young.physics.ucsc.edu/116A/erfc.pdf](https://young.physics.ucsc.edu/116A/erfc.pdf) .

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [December 2, 2021, 3:21pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/3 "2021-12-02T15:21:26Z")

</div>

Well the integral are slightly more complex. The problem comes from applying the [method of weighted residuals](http://faculty.washington.edu/finlayso/MWR-AReview.pdf) to a distribution function W(x) = \sum\_{i=0}^N C\_i \phi\_i(x) expressed on the basis of Laguerre functions \phi\_n(x) = e^{-x} L\_n(x) where L\_n is the n-th Laguerre polynomial. The weighting functions are \psi\_j(x) = x^{j-1} e^{-\lambda x}. Next, defining the inner product as

\langle f, g \rangle = \int\_0^\infty f(x) \cdot g(x) dx

I am looking for the matrix

\langle \psi\_j, x M(x) \phi\_n(x)\rangle, \\ j = 1,...,N-1,\; n = 0,...,N

where

M(x) = \frac{2}{\sqrt\pi} \frac{e^{-(x - x\_c)^2}}{\mathrm{erfc}(x - x\_c)}

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [December 2, 2021, 3:22pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/4 "2021-12-02T15:22:07Z")

</div>

I think you can also compute quadrature weights and points exactly for your function \exp(-x^2)\over\mathrm{erfc}(x), you can use sympy or any other arbitrary precision package to get exact double precision points and weights. Then it should work nicely.

---

<div class="post-metadata">

**Author:** ![septc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/septc/32/77_2.png) [@septc](https://fortran-lang.discourse.group/u/septc)\
**Post date:** [December 2, 2021, 3:35pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/5 "2021-12-02T15:35:07Z")

</div>

Isn’t Method-A less accurate because of the cancellation of large values (`-x**2` and `logerfc(x)`) for large `x`? According to the page below, `erfcx(x)` behaves like `1/x` for large x, so Method-B may not suffer from such cancellation? (If the built-in routine is accurate enough…)

[https://www.intel.com/content/www/us/en/develop/documentation/fortran-compiler-oneapi-dev-guide-and-reference/top/language-reference/a-to-z-reference/e-to-f/erfc-scaled.html](https://www.intel.com/content/www/us/en/develop/documentation/fortran-compiler-oneapi-dev-guide-and-reference/top/language-reference/a-to-z-reference/e-to-f/erfc-scaled.html)

> The complementary error function is asymptotic to ![note graphic a](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/d/d6ca0c4c0d6097b2123462b9ebb1a54213fb787a.gif). As such it underflows for ![note graphicb](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/8/81fdd2515850d3db0f6cf3af3400151084434566.gif) when using IEEE single precision arithmetic. The exponentially-scaled complementary error function is asymptotic to ![note graphic c](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/0/01b0fb68fe92aac17fd67eb83a3ee9ce7d736f05.gif). As such it does not underflow until ![note graphic d](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/3/380ef3c27e26b26b107271338b0e587a8c80a02a.gif).

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [December 2, 2021, 3:44pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/6 "2021-12-02T15:44:40Z")

</div>

This would involve finding the recurrence relation of orthogonal polynomials for the custom weight function? Then using the one of the procedures from the [OPQ](https://www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html) package (or similar tools in SymPy) to find the custom quadrature rule. Sounds doable but complicated.

The computations we are trying to reproduce where originally done in 1969 on an [IBM 7044](https://www.computerhistory.org/collections/catalog/102646610), however the numerical values are not given so we have no means of error-checking. The authors said they verified their Gauss-Laguerre results with a Newton-Cotes rule, the two different integration methods matching to three digits. The accuracy might not have been so crucial at the time, given that the model contains some experimental parameters, and the constraints on computer time (~20 hrs). Still I’d like to be sure we’ve done the integrations as accurately as possible with modern tools.

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [December 2, 2021, 4:14pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/8 "2021-12-02T16:14:09Z")

</div>

The authors state they used the 32-node Gauss-Laguerre quadrature rule, with the abscissas taken from the book “Approximate calculation of integrals” by V. I. Krylov, published in 1962 (translation by Arthur Stroud). These are the same values provided above, meaning the domain was [0, 111.75…]

The computations were performed at the Computer Center of [IIT Kanpur](https://en.wikipedia.org/wiki/IIT_Kanpur) in India. The paper is given here: [On the solution of statistical models of cell populations - ScienceDirect](https://www.sciencedirect.com/science/article/abs/pii/0025556471900502) (however it’s missing many details of the calculation procedures).

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [December 3, 2021, 2:46pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/11 "2021-12-03T14:46:41Z")

</div>

> [@ivanpribec](#):
>
> This would involve finding the recurrence relation of orthogonal polynomials for the custom weight function? Then using the one of the procedures from the [OPQ](https://www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html) package (or similar tools in SymPy) to find the custom quadrature rule. Sounds doable but complicated.

Yes. I think this task is actually very common, so having an easy to use package (say in SymPy) would be very helpful (I created an issue for this: [Gaussian quadrature rule for arbitrary weight function and interval · Issue #22594 · sympy/sympy · GitHub](https://github.com/sympy/sympy/issues/22594)). On the Fortran side, it would be very simple to use the new points and weights.

I don’t see the functionality here: [https://github.com/sympy/sympy/tree/master/sympy/integrals](https://github.com/sympy/sympy/tree/master/sympy/integrals), but we should add it.

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://avatars.discourse-cdn.com/v4/letter/p/958977/32.png) [@PeterSimon](https://fortran-lang.discourse.group/u/PeterSimon)\
**Post date:** [December 3, 2021, 5:30pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/12 "2021-12-03T17:30:34Z")

</div>

There is an easy to use function for generating arbitrary precision Gaussian quadrature nodes and weights from an arbitrary weight function in the Julia [QuadGK](https://github.com/JuliaMath/QuadGK.jl) package. See [this section](https://github.com/JuliaMath/QuadGK.jl#gaussian-quadrature-and-arbitrary-weight-functions) of the README.md file and the [documentation](https://github.com/JuliaMath/QuadGK.jl/blob/9f9dc68a90e42dc99351f564dbf2c9dc5f10db19/src/weightedgauss.jl#L42) for the `gauss` method. It does not require one to first find the recurrence relation of the associated orthogonal polynomials.

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [December 3, 2021, 5:53pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/13 "2021-12-03T17:53:15Z")

</div>

> [@PeterSimon](#):
>
> There is an easy to use function for generating arbitrary precision Gaussian quadrature nodes and weights from an arbitrary weight function in the Julia [QuadGK](https://github.com/JuliaMath/QuadGK.jl) package.

Indeed. Unfortunately it seems numerical, so it works for smaller intervals, but it’s hard to get it to converge to get high accuracy:

```auto
julia> using SpecialFunctions

julia> using QuadGK

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, 0, 20, rtol=1e-9)
([0.38768537755077603, 1.7568858761982542, 3.803903391026576, 6.336780964014735, 9.152644231849653, 12.025290131154804, 14.722945139994678, 17.027550936537043, 18.752707347179335, 19.759208293813685], [1.385352625729292, 6.137392669224711, 16.158630107464138, 30.815101135175087, 47.04694093046928, 60.38252780711042, 66.31501257350119, 61.74935694664822, 46.08036828553858, 21.583335269112947])

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, 0, 30, rtol=1e-9)
ERROR: DomainError with 15.0:
integrand produced NaN in the interval (0, 30)
Stacktrace:
  [1] evalrule(f::var"#3#4", a::Int64, b::Int64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, nrm::typeof(LinearAlgebra.norm))
    @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/evalrule.jl:41
  [2] #2
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:10 [inlined]
  [3] ntuple
    @ ./ntuple.jl:48 [inlined]
  [4] do_quadgk(f::var"#3#4", s::Tuple{Int64, Int64}, n::Int64, atol::Nothing, rtol::Float64, maxevals::Int64, nrm::typeof(LinearAlgebra.norm))
    @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:10
  [5] (::QuadGK.var"#28#29"{Nothing, Float64, Int64, Int64, typeof(LinearAlgebra.norm)})(f::Function, s::Tuple{Int64, Int64}, #unused#::Function)
    @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:179
  [6] handle_infinities
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:113 [inlined]
  [7] #quadgk#27
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:177 [inlined]
  [8] _gauss(W::var"#3#4", N::Int64, a::Int64, b::Int64, rtol::Float64, quad::Function)
    @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/weightedgauss.jl:74
  [9] gauss(W::Function, N::Int64, a::Int64, b::Int64; rtol::Float64, quad::Function)
    @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/weightedgauss.jl:62
 [10] top-level scope
    @ REPL[4]:1

```

I think this has to be done symbolically (say in SymPy) to obtain accurate answers.

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://avatars.discourse-cdn.com/v4/letter/p/958977/32.png) [@PeterSimon](https://fortran-lang.discourse.group/u/PeterSimon)\
**Post date:** [December 3, 2021, 6:19pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/14 "2021-12-03T18:19:19Z")

</div>

You can get quite a bit further using extended precision float types:

```auto
julia> using SpecialFunctions, QuadGK, Quadmath

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, Float128(0), Float128(20), rtol=1e-9)
(Float128[3.87685377550762212241773092779845364e-01, 1.75688587619823879175737923069572830e+00, 3.80390339102656684792616744050199912e+00, 6.33678096401472421945915155325746325e+00, 9.15264423184965727190966650300169598e+00, 1.20252901311547957434572683103166754e+01, 1.47229451399946729266320462093104682e+01, 1.70275509365370377100260435382257308e+01, 1.87527073471793350183947750993046031e+01, 1.97592082938136925783941627417989462e+01], Float128[1.38535262572932963828648215216764593e+00, 6.13739266922469909867619904385367344e+00, 1.61586301074638846605178484544102100e+01, 3.08151011351749197195145122889382423e+01, 4.70469409304694732526895358042952104e+01, 6.03825278071104610553729784281794292e+01, 6.63150125735013793052707456707147675e+01, 6.17493569466472859903879566530466625e+01, 4.60803682855379145652097342145010420e+01, 2.15833352691123531123418668325872807e+01])

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, Float128(0), Float128(30), rtol=1e-9)
(Float128[6.34659407224465463142196297858512046e-01, 2.71171768735281275985578726914630361e+00, 5.77172336757198750829571201010836730e+00, 9.55841360601610044207575638130383600e+00, 1.37700767174041187617258255999648973e+01, 1.80676417080122640594540138116366018e+01, 2.21038634233279372174925678952818276e+01, 2.55522054067187704666617804877047030e+01, 2.81336073698157885360344273804487407e+01, 2.96396882894937045382468107761158695e+01], Float128[2.67466117744153237486666477119578476e+00, 1.33106909436555609478674658058197336e+01, 3.60073195706674527496696891058388764e+01, 6.90461357829318821070728103163382012e+01, 1.05537609676403258061350223554796575e+02, 1.35475607326787583431931249161960996e+02, 1.48777295551138076209416609062028994e+02, 1.38521426665211170169908953165660962e+02, 1.03363635159087066450004869773698807e+02, 4.84118206425137524135593371145666744e+01])

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, Float128(0), Float128(40), rtol=1e-9)
(Float128[9.00780115683006498306807582292482960e-01, 3.67737225112924308072550023859929769e+00, 7.74713912436374194731955235756185976e+00, 1.27862534586589845322831579891730792e+01, 1.83924496689391570427540879394873073e+01, 2.41136476243153953855474693368597716e+01, 2.94872120320537425920291275042218204e+01, 3.40782336562836041685666433805137078e+01, 3.75150851457707441115157831550554809e+01, 3.95202799518310128320584312641343963e+01], Float128[4.42137875691497530987515277003128497e+00, 2.34095476609499634285705096610784881e+01, 6.38949552209232300492328233141610626e+01, 1.22642636137265081544108285968894329e+02, 1.87438369155112030182459335689592933e+02, 2.40553300756315163130237215575937092e+02, 2.64120951575474784096709157734308366e+02, 2.45880067727169222441715410450082335e+02, 1.83457612796487614807187590772431716e+02, 8.59209671882741275499940977858120291e+01])

```

I guess it’ll run out of gas at some point, but perhaps this is far enough for the OP’s application?

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [December 3, 2021, 6:44pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/15 "2021-12-03T18:44:25Z")

</div>

Thanks. I tried it but I am getting:

```auto
julia> using SpecialFunctions, QuadGK, Quadmath

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, Float128(0), Float128(20), rtol=1e-9)
ERROR: MethodError: no method matching eigen!(::LinearAlgebra.SymTridiagonal{Float128, Vector{Float128}})
Closest candidates are:
  eigen!(::LinearAlgebra.SymTridiagonal{var"#s814", V} where {var"#s814"<:Union{Float32, Float64}, V<:AbstractVector{var"#s814"}}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:280
  eigen!(::LinearAlgebra.SymTridiagonal{var"#s814", V} where {var"#s814"<:Union{Float32, Float64}, V<:AbstractVector{var"#s814"}}, ::UnitRange) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:283
  eigen!(::LinearAlgebra.SymTridiagonal{var"#s814", V} where {var"#s814"<:Union{Float32, Float64}, V<:AbstractVector{var"#s814"}}, ::Real, ::Real) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:288
  ...
Stacktrace:
 [1] eigen(A::LinearAlgebra.SymTridiagonal{Float128, Vector{Float128}})
   @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:281
 [2] _gauss(W::var"#5#6", N::Int64, a::Float128, b::Float128, rtol::Float64, quad::Function)
   @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/weightedgauss.jl:98
 [3] gauss(W::Function, N::Int64, a::Float128, b::Float128; rtol::Float64, quad::Function)
   @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/weightedgauss.jl:62
 [4] top-level scope
   @ REPL[8]:1

```

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://avatars.discourse-cdn.com/v4/letter/p/958977/32.png) [@PeterSimon](https://fortran-lang.discourse.group/u/PeterSimon)\
**Post date:** [December 3, 2021, 9:49pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/16 "2021-12-03T21:49:29Z")

</div>

I am able to reproduce your error when I start from a clean Julia session. Prior to posting my previous result, I had attempted to use the DoubleFloats package instead of the Quadmath package. The DoubleFloats package worked for an upper limit of 20 for the range of integration but failed when I tried 30 as the upper limit. So I then switched to the higher precision Quadmath package, but in the same Julia session. If I start with a clean Julia session I find that the following works (I’ve added some checks on the results versus standard Float64 for the case with an upper limit of 10):

```auto
julia> using SpecialFunctions, QuadGK, DoubleFloats, Quadmath

julia> x, w = gauss(x -> exp(-x^2) / erfc(x), 10, 0, 10, rtol=1e-9)
([0.17002966561186994, 0.8275983996020136, 1.8519610263611523, 3.1268538804877104, 4.544187900979246, 5.989495947115991, 7.346369502028654, 8.5053477390109, 9.372837776505476, 9.878927439535989], [0.5113140688691087, 1.798695873902142, 4.288803672361198, 7.917855524736181, 11.95622529170446, 15.281890838276784, 16.75469046153688, 15.589435499209033, 11.629641986675479, 5.446374524138117])

julia> x2, w2 = gauss(x -> exp(-x^2) / erfc(x), 10, Double64(0), Double64(10), rtol=1e-9)
(Double64[0.17002966561186245, 0.8275983996020075, 1.851961026361147, 3.1268538804877046, 4.544187900979243, 5.9894959471159845, 7.346369502028653, 8.505347739010896, 9.372837776505476, 9.87892743953599], Double64[0.5113140688691216, 1.7986958739021315, 4.288803672361115, 7.91785552473613, 11.956225291704522, 15.28189083827682, 16.754690461536985, 15.589435499208923, 11.629641986675392, 5.446374524138095])

julia> x3, w3 = gauss(x -> exp(-x^2) / erfc(x), 10, Float128(0), Float128(10), rtol=1e-9)
(Float128[1.70029665611862433659224649784303822e-01, 8.27598399602007522716310792972464948e-01, 1.85196102636114703415280335809190867e+00, 3.12685388048770480613305321518244925e+00, 4.54418790097924340654885464452625230e+00, 5.98949594711598487888835436986946672e+00, 7.34636950202865233803120255124583270e+00, 8.50534773901089557533579797067995003e+00, 9.37283777650547568931649112718171529e+00, 9.87892743953598985096491870518316858e+00], Float128[5.11314068869121514855563410584036841e-01, 1.79869587390213138052062743401453301e+00, 4.28880367236111457651850781821224047e+00, 7.91785552473613082702253498409252333e+00, 1.19562252917045218072471809363012327e+01, 1.52818908382768184272135307269923489e+01, 1.67546904615369832026084526020777775e+01, 1.55894354992089224522041996753377941e+01, 1.16296419866753912481246448811522357e+01, 5.44637452413809532657056757168147170e+00])

julia> all(x .≈ x2 .≈ x3)
true

julia> all(w .≈ w2 .≈ w3)
true

```

Clearly, it should not be necessary to use the `DoubleFloat` type first, and the results are therefore suspect, despite the checks above. I will raise this as a question on the Julia discourse group. Meantime, for whatever it’s worth, I’m providing the following result for the OP:

```auto
julia> x32_100, w32_100 = gauss(x -> exp(-x^2) / erfc(x), 32, Float128(0), Float128(100), rtol=1e-9);

julia> for (x,w) in zip(x32_100, w32_100)
         println("(", x, ", ", w, ")")
       end
(1.91963135636113087723944127177350547e-01, 5.91881243864281397404682379530372879e-01)
(9.47784705001635673735635624115440130e-01, 2.28103406732697400070025195716587568e+00)
(2.18462250050000638860517247490728286e+00, 6.15270038294108072957879090736266827e+00)
(3.85626859427070854587320930367100973e+00, 1.32715351122281943804061115615437698e+01)
(5.94303187308467266601038365516999644e+00, 2.44376376900740335278372351151039321e+01)
(8.42638901227391390038722405982762608e+00, 4.02246061438534077839506391402296262e+01)
(1.12844428360141741498899057330258462e+01, 6.09839427552017486920159068027314587e+01)
(1.44916867548752248131178061450934793e+01, 8.68305492503164311358005709379575465e+01)
(1.80192958552878944913476484862152326e+01, 1.17631095079496323732167665455715792e+02)
(2.18354546184332115314525436246271114e+01, 1.53001126546825735713393779220942402e+02)
(2.59056840759581005940415884388064161e+01, 1.92311969911848824229307360626933545e+02)
(3.01931752794688898020601038483737474e+01, 2.34707308606169716260678027779734561e+02)
(3.46591343113510252791447051883222276e+01, 2.79128854746171706162643569323787704e+02)
(3.92631402019628602635706496887434664e+01, 3.24350200145202995678699215690886389e+02)
(4.39635148406255416897901927267930996e+01, 3.69017632568504248917732540166495400e+02)
(4.87177025562547582957593255134124877e+01, 4.11696339913465915607232700011629950e+02)
(5.34826585903904232532942371599805010e+01, 4.50921293491815542882676339208858114e+02)
(5.82152359116053165875036252841666312e+01, 4.85248084612860171671828307608346226e+02)
(6.28725812916780241616703248427819922e+01, 5.13306732970190961236736744924939985e+02)
(6.74125209660953831834501594401474889e+01, 5.33851694739406618351103770667634646e+02)
(7.17939436152740935644911222492100765e+01, 5.45809077392199778016709139664051254e+02)
(7.59771728543853796379063808248813219e+01, 5.48318448338041974701335519348192086e+02)
(7.99243267019881804634226842838910505e+01, 5.40767807257881291911081577349390172e+02)
(8.35996607583697538931447220700354058e+01, 5.22820421115705296761046480758239117e+02)
(8.69698920177915542133292792344160127e+01, 4.94432508255697671603529107031386430e+02)
(9.00045004009175853681883379968963752e+01, 4.55861098417666542493676776663789586e+02)
(9.26760052938985683299002166380371211e+01, 4.07661727779874996996623675559100551e+02)
(9.49602147232202962186959757488623784e+01, 3.50675997908408904323976953687314454e+02)
(9.68364453107542932707453209942753953e+01, 2.86009489026088228757267541583554843e+02)
(9.82877127420278085522421230085364442e+01, 2.15001274170254064551428620628945615e+02)
(9.93009023803250404838548836044167475e+01, 1.39190361543864742545054280921176578e+02)
(9.98670662667736393395987132357698739e+01, 6.03633386886598401030275140720523847e+01)

```

Edit: The usage of `Quadmath` types in a call to `gauss` following a call using `DoubleFloats` types is valid, as explained in [this](https://discourse.julialang.org/t/unexpected-behavior-of-quadgk-with-extended-precision-floats/72684?u=petersimon) post.

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [December 3, 2021, 10:18pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/17 "2021-12-03T22:18:20Z")

</div>

@ivanpribec would you mind trying the points and weights that @PeterSimon provided?

I tested the Julia code and it works for me.

Why are the weights decaying so slowly? I would expect the weights for x=100 to go to zero much faster.

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://avatars.discourse-cdn.com/v4/letter/p/958977/32.png) [@PeterSimon](https://fortran-lang.discourse.group/u/PeterSimon)\
**Post date:** [December 3, 2021, 11:06pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/18 "2021-12-03T23:06:15Z")

</div>

> Why are the weights decaying so slowly? I would expect the weights for x=100 to go to zero much faster.

Looking at this a bit more, the weight function `exp(-x^2) / erfc(x)` grows with increasing x, which is probably why the weights aren’t decaying. I don’t really know much about Gaussian integration theory but I would be surprised if such a weight function is even valid. I like the OP’s original thought of using a Gauss-Laguerre scheme that uses the decaying exponential as the weight function. Some of the other ideas presented in this thread can allow him to evaluate the remaining part of the integrand at large x values.

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [December 4, 2021, 1:05am UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/19 "2021-12-04T01:05:28Z")

</div>

Yeah, right after I posted my comment, I noticed that `exp(-x^2) / erfc(x) ~ O(x) for x->oo`, so that explains it. I would choose a weight function that decays by incorporating into this some part of the function being integrated, such as `exp(-x)` or something. Then the weights should look good.

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [December 4, 2021, 10:09am UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/20 "2021-12-04T10:09:59Z")

</div>

@PeterSimon, @certik, thanks for the demonstration of the custom weight calculation. This is one of the big advantages of Julia, that you can leverage custom types across the whole library ecosystem.

As you’ve both noticed, the function f(x) (or M(x)) grows linearly as x approaches infinity. It’s supposed to model the probability of cell division. Above the critical mass x\_c the function slopes upward.

![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/1/1e1fcef8390957d88fae44cf916bbba8fe329f5b.png)

However the model cell distribution behaves like exp(-x) for large x values (large cells are very rare, because they divide), so the weight in the integral actually looks something like this (depending on the trial function and weighting functions used in the method of weighted residuals):

![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/2/2ab979f5b5e2a5af24816fcabc9030c14ee366f1.png)

My current plan was to first compare the results of Gauss-Laguerre quadrature with some adaptive methods available in QUADPACK, and MATLAB, to determine the accuracy. If not sufficient then I know custom quadrature rules are an option.

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [December 4, 2021, 10:14am UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/21 "2021-12-04T10:14:59Z")

</div>

Thanks for the function using the continued fraction. It appears that in the region x \< 20, the product `exp(-x**2) * erfc(x)` is also slightly more accurate than `erf_scaled`.

I’ve tested the routine and also found the absolute errors of order 10\*\*(-14), in close agreement with the MPFR results. However if I evaluate at the quadrature nodes, instead of the whole numbers (5,10,15,…) the accuracy appears to drop back to 10\*\*(-12) on par with Method A and B. I guess this is just a property of floating-point math?

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [December 4, 2021, 4:58pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/22 "2021-12-04T16:58:13Z")

</div>

Is the effective weight function then something like `exp(-x^2) / erfc(x) * exp(-x)`? That will behave nicely. If so, then you can use the well behaved points and weights for this weight function, and then when you integrate your actual function `f(x)` simply integrate `f(x)*exp(x)`, because we include `exp(-x)` in the new weight function, so that things cancel.

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [December 4, 2021, 5:33pm UTC](https://fortran-lang.discourse.group/t/accurate-computation-of-exp-x-2-erfc-x/2360/23 "2021-12-04T17:33:14Z")

</div>

Yes, precisely.

So after shuffling the terms around, we get something like:

\int\_0^\infty \underbrace{e^{-(\lambda + 1)x} \frac{e^{-(x-x\_c)^2}}{\mathrm{erfc}(x-x\_c)}}\_{\omega(x)\text{ - effective weight}} x^{i + k - 1} dx

where the i + k - 1 is a non-negative integer constant. There’s an additional \lambda factor in the exponent, which provides some control over the conditioning of the moment system.

The factor x^k originates from the terms of the Laguerre polynomial expansion. In the original work they decided that since the Laguerre polynomials are oscillating functions it’s better to integrate the monomial parts individually and then recombine them by multiplying with the coefficients of the closed-form Laguerre polynomial,

L\_n(x) = \sum\_{k=0}^n {n \choose k}\frac{(-1)^k}{k!} x^k.