Does Silverfrost have the FORTRAN ERFC function?

ERFC - complementary error function - could not be found in the FTN95 library.

I get the build error: This expression cannot receive a value, in the call to ERFC, with argument 1 ('X') which has been declared as INTENT(INOUT)

Can someone verify whether this function is avaiable in the Silverfrost provided libraries?

Thanks, Alex.

I am not aware of the function being available in FTN95.

You can use the approximation provided in Wikipedia http://en.wikipedia.org/wiki/Error_function This gives about 4 figures of accuracy.

Excel also has the function, and you could use this to generate a look-up table for log(erfc(x)) to provide improved accuracy, say using quadratic or cubic interpolation. I've done this for the normal distribution.

John

My approach using a look-up table. ! test erfc by dumping values into a spreadsheet !
real8 random, erfc, x, y integer4 i external random, erfc ! do i = 1,10000 x = random () * 3.0 y = erfc (x) write (,) i,x,y end do end

      real*8 function erfc (x)
!
!  calculate the value of erfc for x  using a quadratic lookup table for x in the range 0:6
!  use table of ln(erfc(x)) as this is closer to quadratic variation
!
      real*8 x
!
      real*8   quadratic_function, y 
      external quadratic_function   !  (x, table, x0, dx, n)
      real*8   x0, dx, log_erfc_table(-1:120)
      integer*4 n
!
      data n, x0, dx / 120, 0, 0.05 /
!
!  look up table for ln(erfc(x)) in steps of 0.05
      data log_erfc_table /          &   !  ln(erfc.precise(x)) from excel_2010
     0.054840375,               0,   &   !     -0.05,  0
    -0.058023235,    -0.119304974,   &   !       0.1
    -0.183917996,    -0.251932234,   &   !       0.2
    -0.323414803,    -0.398430051,   &   !       0.3
     -0.47703961,    -0.559302458,   &   !       0.4
    -0.645274999,     -0.73501113,   &   !       0.5
    -0.828562327,    -0.925977727,   &   !       0.6
    -1.027304216,    -1.132586514,   &   !       0.7
    -1.241867261,    -1.355187107,   &   !       0.8
    -1.472584795,    -1.594097247,   &   !       0.9
    -1.719759644,     -1.84960551,   &   !         1
    -1.983666787,    -2.121973911,   &   !       1.1
    -2.264555887,    -2.411440355,   &   !       1.2
    -2.562653662,    -2.718220922,   &   !       1.3
    -2.878166081,    -3.042511975,   &   !       1.4
    -3.211280384,     -3.38449209,   &   !       1.5
    -3.562166921,    -3.744323809,   &   !       1.6
    -3.930980826,    -4.122155234,   &   !       1.7
    -4.317863525,     -4.51812146,   &   !       1.8
    -4.722944103,    -4.932345863,   &   !       1.9
    -5.146340521,    -5.364941265,   &   !         2
    -5.588160718,    -5.816010969,   &   !       2.1
    -6.048503596,    -6.285649693,   &   !       2.2
    -6.527459896,    -6.773944403,   &   !       2.3
    -7.025112994,    -7.280975056,   &   !       2.4
    -7.541539599,    -7.806815273,   &   !       2.5
    -8.076810388,    -8.351532931,   &   !       2.6
    -8.630990575,    -8.915190702,   &   !       2.7
     -9.20414041,    -9.497846531,   &   !       2.8
    -9.796315639,    -10.09955407,   &   !       2.9
    -10.40756791,    -10.72036304,   &   !         3
    -11.03794513,    -11.36031963,   &   !       3.1
    -11.68749181,    -12.01946675,   &   !       3.2
    -12.35624936,    -12.69784435,   &   !       3.3
    -13.04425632,    -13.39548966,   &   !       3.4
    -13.75154865,     -14.1124374,   &   !       3.5
    -14.47815991,    -14.84872003,   &   !       3.6
    -15.22412149,    -15.60436789,   &   !       3.7
    -15.98946274,    -16.37940941,   &   !       3.8
    -16.77421118,    -17.17387122,   &   !       3.9
     -17.5783926,    -17.98777831,   &   !         4
    -18.40203124,    -18.82115418,   &   !       4.1
    -19.24514986,    -19.67402092,   &   !       4.2
     -20.1077699,     -20.5463993,   &   !       4.3
    -20.98991154,    -21.43830894,   &   !       4.4
    -21.89159379,     -22.3497683,   &   !       4.5
    -22.81283463,    -23.28079485,   &   !       4.6
      -23.753651,    -24.23140506,   &   !       4.7
    -24.71405894,    -25.20161453,   &   !       4.8
    -25.69407363,    -26.19143803,   &   !       4.9
    -26.69370944,    -27.20088955,   &   !         5

    -25.69407363,    -26.19143803,   &   !       4.9
    -26.69370944,    -27.20088955,   &   !         5
    -27.71297998,    -28.22998235,   &   !       5.1
    -28.75189819,    -29.27872902,   &   !       5.2
    -29.81047631,    -30.34714152,   &   !       5.3
    -30.88872603,    -31.43523121,   &   !       5.4
     -31.9866584,    -32.54300891,   &   !       5.5
    -33.10428399,     -33.6704849,   &   !       5.6
    -34.24161284,    -34.81766899,   &   !       5.7
    -35.39865451,    -35.98457053,   &   !       5.8
    -36.57541814,    -37.17119842,   &   !       5.9
    -37.77191242,    -38.37756117  /     !         6
!
      y = quadratic_function (x, log_erfc_table, x0, dx, n)
!
      erfc = exp (y)
!
      end function erfc

      real*8 function quadratic_function (x, table, x0, dx, n)
!
!  Interpolate the value of x from a table of values using 3-point quadratic interpolation
!
!   x           is the value to be estimated
!   table(-1:n) is a table of values
!   x0          is the function value at table(0)
!   dx          is the spacing of x values for table(i) = f(x0+i*dx)
!
      integer*4 n                ! size of table
      real*8    x                ! value to be interpolated
      real*8    table(-1:n)      ! table of values
      real*8    x0               ! x for table(0)
      real*8    dx               ! spacing of values in table
!
      real*8    fy, y, u     ! , f1,f2,f3, a1,a2,a3
      integer*4 i, sgn
!
      real*8, parameter :: one   = 1
      real*8, parameter :: two   = 2
!
!    Test for -ve value
      if (x < x0) then
         sgn = -1
         y   = (x0-x)/dx
      else
         sgn = 0
         y   = (x-x0)/dx
      end if
!
!  require i+1 <= n ; i < n ; u in range -0.5 to 0.5
      i  = nint (y)                    ! table reference
      u  = (y - dble(i))               ! transform to -1:1
!
!    Shape function formulation
!      f1 = u+one
!      f2 = u/two
!      f3 = u-one
!
!      a1 =     f2*f3                   ! y = -1
!      a2 = -f1*   f3                   ! y =  0
!      a3 =  f1*f2                      ! y =  1
!
!      fy = a1*table(i-1) + a2*table(i) + a3*table(i+1)
!
!    Expanded formulation
      fy = table(i) + u/two * ( (u-one)*table(i-1) - two*u*table(i) + (u+one)*table(i+1) )
       if (sgn /= 0) fy = -fy
!
      quadratic_function = fy
!
      end function quadratic_function
 

I'd be interested if others agree with my approach.

John

John,

Your approach is reasonable! (But please add D0 to the constants in the data statement 😃 I may copy your idea for a different distribution.

I usually use a 'standard' approach which involves using rational polynomial approximations for different ranges of the value x. Details are on the web if people want to find it (search AS241 for Fortan code for the normal distribution - you have to do a trivial transformation on the normal distribution to get erfc).

I believe that the error function is a new intrinsic function in Fortran 2008, but I suspect we will have to wait awhile before that becomes easily available. I've been using Cody's algorithm for a while, which you can find here:

http://www.netlib.org/toms/715

I also used a 4 point cubic_function interpolation, again based on FE shape function formulation and also reduced the spacing to 0.025 in the range 0-1 and get 9 figures of accuracy, in comparison to the values reported in excel 2010 function : erfc.precise.

Some points of interest.

1. The 4-point cubic, based on shape functions gives a much more accurate interpolation than compared to bessel's 4 point interpolation, although this is probably only the case for interpolating from a table of very accurate values. If the table entries are approximate (experimental results), then higher order cubics do not work well.
2. The references to other polynomial approaches and also Excel's functions do not give an indication of accuracy. I wonder how accurate they are. I find that interpolating from tables and getting 9 figure accuracy with cubic is much better than 3 to 4 figure accuracy with linear interpolation.
3. the use of log to transform the table to be more suited for quadratic/cubic interpolation works well for this approach. For generating a standard normal distribution from a random number generator, I use a tan transformation.
4. the requirement of including d0, when specifying a value like -38.3775611732 in a data statement for a real8 is a bit of a puzzle to me. The compiler should retain the precision, rather than truncating to a real4 constant.

Interesting problem and I hope my approach is of interest.

John

If I may, I would add some comments on two of your points.

(2) Methods that are published do usually include an error analysis, so the accuracy is know. I have the paper for AS241 and the accuracy is about 16 significant figures for the double precision version (so very high). Excel doesn't tell you how accurate it is, but I know from experimentation it gives the similar accuracy to AS241 (may 1 or 2 digits less).

(4) You are allowed to omit the kind parameter from constants. However, if it is omitted the compiler will treat it as default real kind (usually, this is what we think of as single precision), and then cast or convert the value to the kind of the variable. If the variable has a higher precision that that of default real some truncation and loss of resolution will occur. It's always best to put the precision in (using D0, or _kind numbers).

D

Well it appears that my look up table approach needs some improvement to approach the precision of the Toms polynomials.

I have taken my cubic interpolation and converted it to a quintic (order 5, 6 point). I have based this on the shape functions used for higher order isoparametric shell elements in the 70's and 80's. It provides a clean way of calculating the exact polynomials of required order. (As I have said before, exact high order polyniomials perform very badly if there are errors in the look-up table.)

The 6-point polyniomial now produces about 12 figures of accuracy. I'll have to settle for that ! For those interested, the revised program (without the lookup table that has been taken from excel 2010; post limits!!) is listed below.

!  test erfc by dumping values into a spreadsheet
!  this used a dual table - did not work, 
!  now works for 1/2 spacing in 0:1 !!
!
      real*8    random, erfc, x, y
      integer*4 i
      external  random, erfc
!
      do i = 1,10000
         x = random () * 3.0d0
         y = erfc (x)
         write (*,*) i,x,y
      end do
      end
 
      real*8 function erfc (x)
!
!  calculate the value of erfc for x  using a quadratic lookup table for x in the range 0:8
!  use table of ln(erfc(x)) as this is closer to polynomial variation than erfc(x)
!  use a finer table in range 0-1 as this is where max interpolation error occurs
!
      real*8 x
!
      real*8   quintic_function, y
      external quintic_function   !  (x, table, x0, dx, n)
      real*8   x0, dx,  log_erfc_table_50(-2:160),  & ! table for range 1:8
                   dx2, log_erfc_table_25(-2:44)      ! table for range 0:1
!
      data x0, dx, dx2 / 0, 0.05d0, 0.025d0 /
!
!  look up table for ln(erfc(x)) in steps of 0.05
      data log_erfc_table_50 / 0.10657640058652d0,  &  !  -0.1 ln(erfc.precise(x)) from excel_2010
    0.05484037495972d0,  0.0d0,                &  !   -0.05
   -0.05802323476898d0, -0.11930497373740d0,   &  !    0.05
...
  -64.26333403703940d0,-65.05708478041260d0,   &  !    7.85
  -65.85579727212700d0,-66.65947197080500d0  /    !    7.95
!
!  look up table for ln(erfc(x)) in steps of 0.025
      data log_erfc_table_25 / 0.05484037495972d0, &  !  -0.05  ln(erfc.precise(x)) from excel_2010
    0.02781320511117d0,  0.0d0,                &  !  -0.025
   -0.02860896488336d0, -0.05802323476898d0,   &  !   0.025
...
   -1.91610728531629d0, -1.98366678682610d0,   &  !   1.025
   -2.05228777258070d0, -2.12197391117516d0  /    !   1.075
!
      if ( x > 1 ) then
         y = quintic_function (x, log_erfc_table_50, x0, dx, 160)
      else 
         y = quintic_function (x, log_erfc_table_25, x0, dx2, 44) 
      end if
!
      erfc = exp (y)   !  using log table provides better convergence
!
      end function erfc

The 6-point interpolation routine.

      real*8 function quintic_function (x, table, x0, dx, n)
!
!  Interpolate the value of x from a table of values using 6-point quintic interpolation
!   formulation based on shape functions from isoparametric shell element 
!
!   x           is for value to be estimated  { in range 0 to dx*(n-2) }
!   table(-2:n) is a table of values
!   x0          is the x value at table(0)
!   dx          is the spacing of x values for table(i) = f(x0+i*dx)
!
      integer*4 n                ! size of table
      real*8    x                ! value to be interpolated
      real*8    table(-2:n)      ! table of values
      real*8    x0               ! x for table(0)
      real*8    dx               ! spacing of values in table
!
      real*8    fy, y, u, f0,f1,f2,f3,f4,f5, w0,w1,w2,w3,w4,w5, a0,a1,a2,a3,a4,a5
      integer*4 i, sgn
!
      real*8, parameter :: one   = 1
      real*8, parameter :: two   = 2
      real*8, parameter :: three = 3
!
!    Test for -ve value
      if (x < x0) then
         sgn = -1
         y   = (x0-x)/dx
      else
         sgn = 0
         y   = (x-x0)/dx
      end if
!
!  require i+3 <= n ; i+2 < n ; i < n-2 ; u in range 0-1
      i  = int (y)                     ! table reference
      u  = (y - dble(i))               ! transform to 0:1
!
!    Shape function formulation
      f0 = u+two                       ! u = -2
      f1 = u+one                       ! u = -1
      f2 = u                           ! u =  0
      f3 = u-one                       ! u =  1
      f4 = u-two                       ! u =  2
      f5 = u-three                     ! u =  3
!
      w0 = -120   !     1*2*3*4*5      ! u = -2
      w1 =   24   !   1  *1*2*3*4      ! u = -1
      w2 =  -12   !  -2*1  *1*2*3      ! u =  0
      w3 =   12   !   3*2*1  *1*2      ! u =  1
      w4 =  -24   !  -4*3*2*1  *1      ! u =  2
      w5 =  120   !   5*4*3*2*1        ! u =  3
!
      a0 =     f1*f2*f3*f4*f5 / w0     ! u = -2
      a1 =  f0*   f2*f3*f4*f5 / w1     ! u = -1
      a2 =  f0*f1*   f3*f4*f5 / w2     ! u =  0
      a3 =  f0*f1*f2   *f4*f5 / w3     ! u =  1
      a4 =  f0*f1*f2*f3   *f5 / w4     ! u =  2
      a5 =  f0*f1*f2*f3*f4    / w5     ! u =  3
!
      fy =  a0*table(i-2)    &
          + a1*table(i-1)    &
          + a2*table(i)      &
          + a3*table(i+1)    &
          + a4*table(i+2)    &
          + a5*table(i+3)
       if (sgn /= 0) fy = -fy   !   ( assume F(-x) = -f(x) )
!
      quintic_function = fy
!
      end function quintic_function

Well I think I've beaten the precision problem when interpolating from a table of values. The only requirement is the table must be for a smooth function. I've taken the shape function approach and written it for a general order polynomial. The only limits are: the table must have at least p+1 values, and the array F must be big enough. (it could be removed) It does require a table of values for the required function. I think it is a neat solution where a function is not available in FTN95.

The following is a test example for SIN, providing the test rprogram and the POLY_function

!     Last change:  JDC   1 Aug 2011   10:34 am
!  test general polynomial for interpolating a table
!  polynomial calculation based on shape function
!  used in iosparametric shell element formulation
!
!  use SIN as easy to replicate
!
      real*8 poly_function, random, x, f, e, s, emax
      integer*4 p, n, i
      external poly_function, random
!
      real*8   table(-2:50), x0, dx      ! table for range 1:8
!
      dx = 0.1d0
      x0 = 0
      n  = 50
      do i = -2,n
         x = dble(i) * dx
         table(i) = sin(x)
      end do
!
      write (*,*) ' Test interpolation of SIN using higher order polynomial'
      do p = 1,10
         s = 0
         emax = 0
         do i = 0,100000
            x = random () * 3.14
            f = poly_function (x, table, x0, dx, n, p)
            e = abs ( f - sin(x))
            s = s + e
            if (e > emax) emax = e
         end do
         write (*,*) 'error = ', s/dble(i),emax, p
      end do
!
      end

The general polynomial function

     real*8 function poly_function (x, table, x0, dx, n, p) 
! 
!  Interpolate the value of x from a table of values using multi-point polynomial interpolation 
!   formulation based on shape functions from isoparametric shell element 
! 
!   x           is for value to be estimated  { in range 0 to dx*(n-high) } 
!   table(-2:n) is a table of values
!   x0          is the x value at table(0) 
!   dx          is the spacing of x values for table(i) = f(x0+i*dx) 
!   p           is power of polynomial approximation
!
      integer*4 p                ! polynomial order (1-8)
      integer*4 n                ! size of table 
      real*8    x                ! value to be interpolated 
      real*8    table(-2:n)      ! table of values
      real*8    x0               ! x for table(0) 
      real*8    dx               ! spacing of values in table 
! 
      real*10   fy, y, u, s    ! , f(-6:6)   ! will go to order 12
      integer*4 low,high, i, sgn, j, k  ! , w 
!
!  P  Low:high values   low = -(p/2);  high = (p+1)/2
!  1    0  1  int 
!  2   -1  1  nint
!  3   -1  2  int
!  4   -2  2  nint
!  5   -2  3  int
!  6   -3  3  nint
!  7   -3  4  int
!  ...
!
      low  = -(p/2)
      high = (p+1)/2
!
!   Test for -ve value 
      if (x < x0) then 
         sgn = -1 
         y   = (x0-x)/dx 
      else 
         sgn = 0 
         y   = (x-x0)/dx 
      end if 
! 
!   require i+high <= n ; i+high < n+1 ; i < n+1-high ; u in range 0:1 or -.5:.5
      if (low+high==0) then
         i = nint (y)                    ! table reference about 0
      else
         i = int (y)                     ! table reference in 0:1
      end if
       if (i < -2-low) i = -2-low        ! test lower bound of available table
       if (i > n-high) i = n-high        ! test upper bound of available table 
!
!   Shape function formulation 
      u  = (y - dble(i))                 ! transform wrt i
!      do k = low,high
!         f(k) = u-dble(k)               ! f removed for general polynomial
!      end do
!
!   Calculate Weighting (w) and Coeffecient (s)
      fy = 0          ! value of function
      do k = low,high
!         w = 1        ! weighting product so shape function fk = 1 at u=k
         s = 1        ! table coefficient product
         do j = low,high
            if (k == j) cycle
!            w = w * (k-j)                ! w removed as overflows for p > 12
!            s = s * f(j)                    ! f(j) removed for general polynomial
            s = s * (u-j) / (k-j)
         end do
!         s  = s/w
         fy = fy + table(i+k) * s
      end do
      if (sgn /= 0) fy = -fy   !   ( assume F(-x) = -f(x) )
!
      poly_function = fy 
! 
      end function poly_function

I find this problem very interesting. I have presented a way of using “P” points from a look-up table to fit a polynomial of order “P-1” and then interpolating the value of the function. This works very well if the function is “smooth” and the table of values is available to good accuracy. If this accuracy is not available, an alternative could be to use a least-squares approach and fit a polynomial of say “P-2”. The main part of the calculation is inverting the matrix used for least squares approximation. This matrix is independent of the function values, being only a function of P: the number of points (x). This method does not require a look-up table of equally spaced points. If the table is of equally spaced points, the matrix does not change and can be retained for the same P. The least squares approach also gives an error estimate, which can be good to qualify the accuracy of the function estimate. Then again, much simpler if it is an intrinsic function !

John