# Fortran erfinv and erfinv functions

**URL:** <https://fortran-lang.discourse.group/t/fortran-erfinv-and-erfinv-functions/6860>\
**Category:** Uncategorized\
**Created:** [November 26, 2023, 2:59am UTC](https://fortran-lang.discourse.group/t/fortran-erfinv-and-erfinv-functions/6860 "2023-11-26T02:59:53Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\
**Post date:** [November 26, 2023, 2:59am UTC](https://fortran-lang.discourse.group/t/fortran-erfinv-and-erfinv-functions/6860/1 "2023-11-26T02:59:53Z")

</div>

Dear all,

I was checking the methods of using uniform random number to generate gaussian random number, since I try to replace the uniform random number with the quasi-uniform random number (such as Sobol sequence) to generate quasi gaussian random number. It may generate some more ‘gaussian’ distributed gaussian random numbers.

Now, a simple approach to do so, can be using the erfinv and erfcinv functions as simply described in MARLAB page as below,

> **[Inverse error function - MATLAB erfinv](https://www.mathworks.com/help/matlab/ref/erfinv.html#bup2lvj)**
>
> This MATLAB function returns the Inverse Error Function evaluated for each element of x.

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

So they use the inverse of error function erfinc, and the inverse of the complimentary error function erfcinv.

I know Fortran has erf and erfc implemented already, but are there some Fortran implementations of erfinv and erfcinv functions?

Many thanks in advance!

---

<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:** [November 26, 2023, 7:46am UTC](https://fortran-lang.discourse.group/t/fortran-erfinv-and-erfinv-functions/6860/2 "2023-11-26T07:46:43Z")

</div>

The functions sought are in the JPL [Math77](https://netlib.org/math/) library. See [https://netlib.org/math/docpdf/ch02-13.pdf](https://netlib.org/math/docpdf/ch02-13.pdf) .

---

<div class="post-metadata">

**Author:** ![DavidB](https://avatars.discourse-cdn.com/v4/letter/d/76d3ee/32.png) [@DavidB](https://fortran-lang.discourse.group/u/DavidB)\
**Post date:** [November 26, 2023, 8:10am UTC](https://fortran-lang.discourse.group/t/fortran-erfinv-and-erfinv-functions/6860/3 "2023-11-26T08:10:51Z")

</div>

A reasonable recent paper, with references to prior work, is

J. S. Dyer and S. A. Dyer, “Approximations to inverse error functions,” in _IEEE Instrumentation & Measurement Magazine_, vol. 11, no. 5, pp. 32-36, October 2008, doi: 10.1109/MIM.2008.4630740.

Abstract: We have presented four approximations to the inverse cumulative distribution function of a standard normal random variable. The first two were meant purely as a tool to motivate the other approximations, although either could be useful when extreme accuracy is not required and one doesn’t venture too far into the tails of the distribution. There are many other such “simple” approximations to F X -1 , suitable for use on a hand calculator.

---

<div class="post-metadata">

**Author:** ![MarDie](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mardie/32/54_2.png) [@MarDie](https://fortran-lang.discourse.group/u/MarDie)\
**Post date:** [November 26, 2023, 12:41pm UTC](https://fortran-lang.discourse.group/t/fortran-erfinv-and-erfinv-functions/6860/4 "2023-11-26T12:41:48Z")

</div>

there is also [stats\_distribution\_normal – Fortran-lang/stdlib](https://stdlib.fortran-lang.org/page/specs/stdlib_stats_distribution_normal.html)

---

<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:** [November 26, 2023, 12:55pm UTC](https://fortran-lang.discourse.group/t/fortran-erfinv-and-erfinv-functions/6860/5 "2023-11-26T12:55:21Z")

</div>

I ported the `erfinv` function from Julia to Fortran, but lacked the motivation to push it across the finish line for stdlib. You can find it here:

> <https://github.com/fortran-lang/stdlib/issues/179#issuecomment-647178168>
>
> Here is a paper from 2007 that was submitted to the Fortran Committee, but ultim…ately rejected:
> 
> https://wg5-fortran.org/N1651-N1700/N1688.pdf
> 
> The functions there seem to be exactly in the scope of stdlib, so we should include them and we can use this paper as a starting point.

The method in Julia uses the tables from Blair et al. (1976):

> Blair, J. M., Edwards, C. A., & Johnson, J. H. (1976). Rational Chebyshev approximations for the inverse of the error function. _Mathematics of Computation_, _30_(136), 827-830. [AMS :: Mathematics of Computation](https://doi.org/10.1090/S0025-5718-1976-0421040-7)

I’ve also seen an implementation for GPU before:

> Giles, M. (2012). Approximating the erfinv function. In _GPU Computing Gems Jade Edition_ (pp. 109-116). Morgan Kaufmann. [Redirecting](https://doi.org/10.1016/B978-0-12-385963-1.00010-1)

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [November 26, 2023, 1:05pm UTC](https://fortran-lang.discourse.group/t/fortran-erfinv-and-erfinv-functions/6860/6 "2023-11-26T13:05:34Z")

</div>

Alan Miller’s [as241.f90](https://jblevins.org/mirror/amiller/as241.f90) does this. Here it is with a short driver.

```auto
module normdev_mod
implicit none
private
public :: dp, nprob, normdev, normal
integer, parameter :: dp = selected_real_kind(15, 307)
contains
elemental function normal(z) result(prob)
real (dp), intent(in) :: z
real (dp) :: prob
call nprob(z, p=prob)
end function normal
!
elemental subroutine nprob(z, p, q, pdf)
! For a given number (z) of standard deviations from the mean, the
! probabilities to the left (p) and right (q) of z are calculated,
! and the probability density (pdf).

! Programmer - Alan J. Miller
! Latest revision - 7 August 1997
! This version is compatible with the ELF90 subset of Fortran 90.

! REFERENCE: ADAMS,A.G. AREAS UNDER THE NORMAL CURVE,
! ALGORITHM 39, COMPUTER J., VOL. 12, 197-8, 1969.

real(dp), intent(in) :: z
real(dp), intent(out), optional :: p, q, pdf

! Local variables
real(dp), parameter :: a0 = 0.5d0, a1 = 0.398942280444d0, a2 = 0.399903438504d0, &
                    a3 = 5.75885480458d0, a4 = 29.8213557808d0, &
                    a5 = 2.62433121679d0, a6 = 48.6959930692d0, &
                    a7 = 5.92885724438d0, b0 = 0.398942280385d0, &
                    b1 = 3.8052d-8, b2 = 1.00000615302d0, b3 = 3.98064794d-4, &
                    b4 = 1.98615381364d0, b5 = 0.151679116635d0, &
                    b6 = 5.29330324926d0, b7 = 4.8385912808d0, &
                    b8 = 15.1508972451d0, b9 = 0.742380924027d0, &
                    b10 = 30.789933034d0, b11 = 3.99019417011d0, &
                    zero = 0.d0, one = 1.d0
real(kind=dp) :: zabs, y, pp, qq, ppdf
zabs = abs(z)
if (zabs < 12.7d0) then
  y = a0*z*z
  ppdf = exp(-y)*b0
  if (present(pdf)) pdf = ppdf

! Z BETWEEN -1.28 AND +1.28

  if (zabs < 1.28) then
    qq = a0 - zabs*(a1-a2*y/(y+a3-a4/(y+a5+a6/(y+a7))))
    if(z < zero) then
      pp = qq
      qq = one - pp
    else
      pp = one - qq
    end if
    if (present(p)) p = pp
    if (present(q)) q = qq
    return
  end if

! ZABS BETWEEN 1.28 AND 12.7

  qq = ppdf/(zabs-b1+b2/(zabs+b3+b4/(zabs-b5+b6/(zabs+b7-b8/ &
            (zabs+b9+b10/(zabs+b11))))))
  IF(z < zero) THEN
    pp = qq
    qq = one - pp
  ELSE
    pp = one - qq
  END IF
  IF (PRESENT(p)) p = pp
  IF (PRESENT(q)) q = qq
  RETURN

! Z FAR OUT IN TAIL

ELSE
  ppdf = zero
  IF(z < zero) THEN
    pp = zero
    qq = one
  ELSE
    pp = one
    qq = zero
  END IF
  IF (PRESENT(p)) p = pp
  IF (PRESENT(q)) q = qq
  IF (PRESENT(pdf)) pdf = ppdf
  RETURN
END IF
END SUBROUTINE nprob
!
elemental FUNCTION normdev(p) RESULT(fn_val)
REAL (dp), INTENT(IN) :: p
REAL (dp) :: fn_val
! Local variable
INTEGER :: ifault
CALL ppnd16(p, fn_val, ifault)
if (ifault /= 0) fn_val = huge(fn_val)
end function normdev
!
elemental SUBROUTINE ppnd16 (p, normal_dev, ifault)

! ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3

! Produces the normal deviate Z corresponding to a given lower
! tail area of P; Z is accurate to about 1 part in 10**16.

! The hash sums below are the sums of the mantissas of the
! coefficients. They are included for use in checking
! transcription.

! This ELF90-compatible version by Alan Miller - 20 August 1996
! N.B. The original algorithm is as a function; this is a subroutine

REAL (dp), INTENT(IN) :: p
INTEGER, INTENT(OUT) :: ifault
REAL (dp), INTENT(OUT) :: normal_dev

! Local variables

REAL (dp), parameter :: zero = 0.d0, one = 1.d0, half = 0.5d0, &
                    split1 = 0.425d0, split2 = 5.d0, const1 = 0.180625d0, &
                    const2 = 1.6d0

! Coefficients for P close to 0.5

REAL (dp), parameter :: a0 = 3.3871328727963666080D0, &
                    a1 = 1.3314166789178437745D+2, &
                    a2 = 1.9715909503065514427D+3, &
                    a3 = 1.3731693765509461125D+4, &
                    a4 = 4.5921953931549871457D+4, &
                    a5 = 6.7265770927008700853D+4, &
                    a6 = 3.3430575583588128105D+4, &
                    a7 = 2.5090809287301226727D+3, &
                    b1 = 4.2313330701600911252D+1, &
                    b2 = 6.8718700749205790830D+2, &
                    b3 = 5.3941960214247511077D+3, &
                    b4 = 2.1213794301586595867D+4, &
                    b5 = 3.9307895800092710610D+4, &
                    b6 = 2.8729085735721942674D+4, &
                    b7 = 5.2264952788528545610D+3
! HASH SUM AB 55.8831928806149014439

! Coefficients for P not close to 0, 0.5 or 1.

REAL (dp), parameter :: c0 = 1.42343711074968357734D0, &
                    c1 = 4.63033784615654529590D0, &
                    c2 = 5.76949722146069140550D0, &
                    c3 = 3.64784832476320460504D0, &
                    c4 = 1.27045825245236838258D0, &
                    c5 = 2.41780725177450611770D-1, &
                    c6 = 2.27238449892691845833D-2, &
                    c7 = 7.74545014278341407640D-4, &
                    d1 = 2.05319162663775882187D0, &
                    d2 = 1.67638483018380384940D0, &
                    d3 = 6.89767334985100004550D-1, &
                    d4 = 1.48103976427480074590D-1, &
                    d5 = 1.51986665636164571966D-2, &
                    d6 = 5.47593808499534494600D-4, &
                    d7 = 1.05075007164441684324D-9
! HASH SUM CD 49.33206503301610289036

! Coefficients for P near 0 or 1.

REAL (dp), parameter :: e0 = 6.65790464350110377720D0, &
                    e1 = 5.46378491116411436990D0, &
                    e2 = 1.78482653991729133580D0, &
                    e3 = 2.96560571828504891230D-1, &
                    e4 = 2.65321895265761230930D-2, &
                    e5 = 1.24266094738807843860D-3, &
                    e6 = 2.71155556874348757815D-5, &
                    e7 = 2.01033439929228813265D-7, &
                    f1 = 5.99832206555887937690D-1, &
                    f2 = 1.36929880922735805310D-1, &
                    f3 = 1.48753612908506148525D-2, &
                    f4 = 7.86869131145613259100D-4, &
                    f5 = 1.84631831751005468180D-5, &
                    f6 = 1.42151175831644588870D-7, &
                    f7 = 2.04426310338993978564D-15
real(kind=dp) :: q, r
! HASH SUM EF 47.52583317549289671629

ifault = 0
q = p - half
if (abs(q) <= split1) then
  r = const1 - q * q
  normal_dev = q * (((((((a7*r + a6)*r + a5)*r + a4)*r + a3)*r + a2)*r + a1)*r + a0) / &
           (((((((b7*r + b6)*r + b5)*r + b4)*r + b3)*r + b2)*r + b1)*r + one)
  return
else
  if (q < zero) then
    r = p
  else
    r = one - p
  end if
  if (r <= zero) then
    ifault = 1
    normal_dev = zero
    return
  end if
  r = sqrt(-log(r))
  if (r <= split2) then
    r = r - const2
    normal_dev = (((((((c7*r + c6)*r + c5)*r + c4)*r + c3)*r + c2)*r + c1)*r + c0) / &
             (((((((d7*r + d6)*r + d5)*r + d4)*r + d3)*r + d2)*r + d1)*r + one)
  else
    r = r - split2
    normal_dev = (((((((e7*r + e6)*r + e5)*r + e4)*r + e3)*r + e2)*r + e1)*r + e0) / &
             (((((((f7*r + f6)*r + f5)*r + f4)*r + f3)*r + f2)*r + f1)*r + one)
  end if
  if (q < zero) normal_dev = - normal_dev
  return
end if
end subroutine ppnd16
end module normdev_mod
program xnormdev
! 12/27/2021 07:56 PM driver for normdev
use normdev_mod, only: dp, normdev
implicit none
real(kind=dp), parameter :: p(*) = [0.4_dp,0.5_dp,0.6_dp,0.7_dp,0.99_dp, &
   0.999_dp,0.9999_dp,0.99999_dp,0.999999_dp]
print "(/,a5,*(1x,f9.6))","p",p
print "(a5,*(1x,f9.6))","z",normdev(p)
end program xnormdev

```

output with gfortran or ifort:

```auto
    p 0.400000 0.500000 0.600000 0.700000 0.990000 0.999000 0.999900 0.999990 0.999999
    z -0.253347 0.000000 0.253347 0.524401 2.326348 3.090232 3.719016 4.264891 4.753424

```

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [November 28, 2023, 8:40pm UTC](https://fortran-lang.discourse.group/t/fortran-erfinv-and-erfinv-functions/6860/7 "2023-11-28T20:40:02Z")

</div>

> [@CRquantum](#):
>
> I was checking the methods of using uniform random number to generate gaussian random number, since I try to replace the uniform random number with the quasi-uniform random number (such as Sobol sequence) to generate quasi gaussian random number. It may generate some more ‘gaussian’ distributed gaussian random numbers.

According to [Generating low-discrepancy sequences from the normal distribution: Box–Muller or inverse transform? (2011) by Ökten and Göncü](https://www.sciencedirect.com/science/article/pii/S0895717710005935) you can use Box-Muller on quasi-random uniform variates to generate quasi-random normal variates.