heather@MATH.UCLA.EDU (06/01/88)
Recently I submitted the question: >Please forgive me if this question has been asked before. > >One of the staff members here at UCLA is porting a statical analysis >application from an IBM mainframe to his departmental VAX. The erf() >function, which is supplied with the FORTRAN compiler on the mainframe >(VS FORTRAN), is not present under VMS FORTRAN. We believe that this >function is part of the FORTRAN 77 standard and are therefore unsure why >it is not supplied. Be that as it may, DEC Tech Support has claimed no >knowledge of this function. Here are the many responses I received. Thanks to you all! Heather Burris, UCLA ************************************************************************ From: Bjorn Ellertsson <bje@MATH.UCLA.EDU> Do "man erf" on a unix system to find out what it is. The code is given in "Numerical Recipies" available at the bookstore. ************************************************************************ From: heisterb%uxe.cso.uiuc.edu@UXC.CSO.UIUC.EDU(Dave Heisterberg) The error function is NOT part of the FORTRAN 77 standard. I recommend you look at _Numerical Recipes_ by Press, Flannery, Teukolsky, and Vetterling, publisher Cambridge University Press. It contains standard comforming routines for evaluation of erf () and its complement. There are many other numerical computing books that should also have useable source code. The method used above involves incomplete gamma functions, but if you're not familiar with them you can simply type in the code and it will work - about 100 lines. Good luck. Oh, may I also recommend that your department purchase a copy of the FORTRAN 77 standard from ANSI? There is so much misinformation on the net concerning subjects near and dear to my heart as it is! ************************************************************************ Also for "Numercial Recipes" From: quintus!ok@SUN.COM(Richard A. O'Keefe) (The answer is found in section 6.2, pp 163-164.) From: Joe Canter <canter@LL-VLSI.ARPA> From: chpf127@EMX.UTEXAS.EDU(Jello Biafra) ************************************************************************ From: quintus!ok@SUN.COM(Richard A. O'Keefe) An obvious place to start looking is CALGO. (The collected algorithms from the ACM.) The Computer Journal (the journal of the BCS) has an algorithms series. JCAM [J. Comp & Appl Math] has some good stuff. Applied Statistics has a series of Fortran-coded algorithms. ************************************************************************ From: jim@DANDELION.CI.COM(Jim Hurt) Approximations to eight place accuracy are not difficult to code. For such formulae, look in your copy of Abramowitz and Stegun. What, you don't have a copy of Abramowitz and Stegun? Then: Abramowitz and Stegun, Eds., "Handbook of Mathematical Functions - AMS 55", Dover, pages 931-932, equations 26.2.1, 26.2.5 and 26.2.17. Fortran code fragment: Calculate Normal cumulative probability: PARAMETER (P * 0.2316419) PARAMETER (B1 * 0.319381530, B2 * -0.356563782, B3 * 1.781477937, B4 * -1.821255978, B5 * 1.330274429) C If X is too small or too big, then the probability is VERY close C to 0.0 or 1.0. Not doing this round-off would cause an C arithmetic overflow when calculating the exp (-(X*X)/2) C on the next statement. IF (X .LE. -6.0) THEN PROB * 0.0 ELSE IF (X .GE. 6.0) THEN PROB * 1.0 ELSE C X is within the normal range - ZX * exp (-(X*X)/2.0) / sqrt (kTwoPi) T * 1.0 / (1.0 + (P * abs(X))) C Cumulative probability for variable <* - abs (X) PROB * (ZX*(B1+(B2+(B3+(B4+B5*T)*T)*T)*T)*T) C Adjust if x > 0 IF (X .GT. 0.0) PROB * 1.0 - PROB ENDIF This is the Cumulative Normal Probability function, which is closely related to the erf and erfc functions. ************************************************************************ Also for Abramowitz and Stegun: From: shenkin@CUBSUN.BIO.COLUMBIA.EDU(Peter Shenkin) ************************************************************************ For IMSL or NAG (commercial packages): From: jerry@VIOLET.BERKELEY.EDU( Jerry Berkman ) From: dgh@SUN.COM(David Hough) From: chpf127@EMX.UTEXAS.EDU(Jello Biafra) ************************************************************************ From: chpf127@EMX.UTEXAS.EDU(Jello Biafra) I'm assuming the function erf() [erf(x)? with x a real/double or complex number] returns the error function -- 2 x erf(x) * ----------- integral [ exp ( -t**2 ) ] dt sqrt ( pi ) 0 for a given value of x. In any case, it's NOT standard fortran, but is available in fortran from a variety of sources. Some to consider are: [Editor's note: ideas already mentioned have been deleted] Argonne National Laboratory's software library 'Netlib'. This is a free software service -- you can get a wide variety of (mostly numerical, mostly written in fortran) software over the net. Some major packages are available this way (Linpack, Eispack, Odepack, Minpack, etc.) but they strongly discourage getting whole libraries over the network. Since you're only after one routine, this seems like the best bet -- I know several special function libraries which are included in all the stuff they have. To find out what they do currently have, send a message to netlib@anl-mcs.arpa of the form: send index ************************************************************************ Also for the netlib server: From: perry@VU-VLSI.VILLANOVA.EDU(Rick Perry) "send erf from fn core" From: Bill Gropp <gropp-bill@YALE.ARPA> ************************************************************************ Misc. comments ************************************************************************ From: dillon@UHCCUX.UHCC.HAWAII.EDU(Ian Dillon) We too have encountered this problem on a u-VAX machine. The only advice I can give you is to get a copy of the erf () yourself. I am in the process of doing this right now, and if you would like a copy sent to you, I'll send it through e-mail. ************************************************************************ From: shenkin@CUBSUN.BIO.COLUMBIA.EDU(Peter Shenkin) As others will probably be quick to tell you, of course this is not part of the FORTRAN 77 standard. Only sin, cos, exp, log, etc. are. The STAT pakg you mention might have it; otherwise, don't you have the IMSL library? (I don't, actually, but it *must* have it -- at least I believe so as certainly as you believe(d?) that it *must* be part of the standard :-( ) ************************************************************************ From: ucla-an!ondine!steve@EE.UCLA.EDU(Steve Jenkins) I have a large collection of portable Fortran numerical routines. It includes single and double precision ERF routines. I'll send them to you Tuesday. ************************************************************************