forums.silverfrost.com

[alex21]

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.

[JohnCampbell]

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

[JohnCampbell]

My approach using a look-up table.

[JohnCampbell]

[davidb]

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).
_________________
Programmer in: Fortran 77/95/2003/2008, C, C++ (& OpenMP), java, Python, Perl

[simon]

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

[JohnCampbell]

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 real*8 is a bit of a puzzle to me. The compiler should retain the precision, rather than truncating to a real*4 constant.

Interesting problem and I hope my approach is of interest.

John

[davidb]

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
_________________
Programmer in: Fortran 77/95/2003/2008, C, C++ (& OpenMP), java, Python, Perl

[JohnCampbell]

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.

[JohnCampbell]

The 6-point interpolation routine.

[JohnCampbell]

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

[JohnCampbell]

The general polynomial function

[JohnCampbell]

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