Conference turris::fortran

1315.0. "erf(x), erfc(x) problem?" by NNTPD::"[email protected]" (ofer) Sun May 25 1997 05:53

Hi,

Angry customer called me about it. Could we have a problem with our erf(x) and
erfc(x) functions? See details and small program below:

Env: 
DU 4.0B, Fortran V4.0, RTL 369

Checks I did: 
1) Checking simple relations - erf(0)=0 ; erf(inf)=1; erf(-x)=-erf(x); That's
worked fine
2) For Small values (x<1.0) , the erf(x) is almost linear (i.e. erf(x) ~
1.13x),
This can be seen from the following relation: erf(x)=
2/sqrt(pi)(x-x**3/x+x**5/10-x**7/42 -...+...)
The results looks wired:
  
         X                      erf(x)                approx. erf(x)

  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
  1.000000000000000E-002   1.00000000000000       1.128341555584961E-002
  2.000000000000000E-002   1.00000000000000       2.256457469184491E-002
  3.000000000000000E-002  -1.00000000000000       3.384122234173529E-002
  4.000000000000000E-002   1.00000000000000       4.511110614512333E-002
  5.000000000000000E-002 -1.792076933179376E-023  5.637197779700635E-002
  6.000000000000000E-002  -1.00000000000000       6.762159439325575E-002
  7.000000000000001E-002  2.202104070867331E-020  7.885771977068003E-002
  8.000000000000000E-002   1.00000000000000       9.007812584031763E-002
  9.000000000000000E-002   1.00000000000000       0.101280593912606     
  0.100000000000000      -1.792076775407195E-023  0.112462916013069     
  0.110000000000000       -1.00000000000000       0.123622896187181     
  0.120000000000000       -1.00000000000000       0.134758351793029     
  0.130000000000000       1.029217597047488E-032  0.145867114780451     
  0.140000000000000       2.202103909308618E-020  0.156947032955268     
  0.150000000000000       4.707964507133511E-008  0.167995971227273     
  0.160000000000000        1.00000000000000       0.179011812840615     
  0.170000000000000        1.00000000000000       0.189992460585226     
  0.180000000000000        1.00000000000000       0.200935837987940     
  0.190000000000000      -8.421070342186311E-036  0.211839890481944     
  0.200000000000000      -1.792077090951557E-023  0.222702586553210     
  0.210000000000000      -3.809458651304887E-011  0.233521918862560     
  0.220000000000000       -1.00000000000000       0.244295905341992     
  0.230000000000000       -1.00000000000000       0.255022590263937     
  0.240000000000000       -1.00000000000000       0.26570004528207
 

The program I used:

	program bbb
	implicit real*8 (a-h,o-z)
	pi=3.1415926535898

c Check of some known relations

	x1=0.0
	x2=10000000000000000000.0
	x3=0.003264
        	 write(6,*) x1,erf(x1)
        	 write(6,*) x2,erf(x2)
        	 write(6,*) x3,erf(-x3), -erf(x3)

c More  checks Using:
c erf(x)= 2/sqrt(pi)(x - x^3/3*1! +x^5/5*2! - x^7/7*3! - ..)

	do x=0.0,1.0, 0.01
		approxerf=2.0/sqrt(pi)* (x- x**3/3 + x**5/10 - x**7/42)
        	 write(6,*) x,erf(x),approxerf
	enddo


      	 stop
      end
[Posted by WWW Notes gateway]

Small mistake in the formula provided. It should be:

erf(X) = 2/sqrt(pi)(X - X**3/3 + X**5/10 - X**7/42 +.....)
[Posted by WWW Notes gateway]

>>Angry customer called me about it. Could we have a problem with our erf(x) and
>>erfc(x) functions?

The problem is with the customer...  The "implicit real*8 (a-h,o-z)" ends up
declaring "erf" as a function that takes a double precision arg and returns
a double precision result.  The external ERF is single precision so ERF was
being passed a bogus bit pattern and returning a bogus result.  Changing "ERF"
to "DERF" in the test program and printing out "DERF(X)-APPROXERF", you will
see that the result of DERF is pretty close to the approximate for numbers close
to 0.0 and the approxmiation gets worse as the arg gets closer to 1.  Same
rsuls if you make everything REAL*4.

/Stan

PS:  Notice that this mismatch of REAL*4 vs *8 args and results between caller
     and callee can go almost unnoticed on VAX/VMS when using the default
     F_ and D_float.  UNIX uses IEEE S_ and T_float which has different formats
     in their fisrt 32 bits {like F_ and G_float on VMS} so you see differences
     when the types don't match.

Ya.. should have noticed it myself. Anyway thanks.

[Posted by WWW Notes gateway]