I have a fortran code from scipy that looks like this:
erf.f
DOUBLE PRECISION FUNCTION ERF(x)
C-----------------------------------------------------------------------
C EVALUATION OF THE REAL ERROR FUNCTION
C-----------------------------------------------------------------------
C .. Scalar Arguments ..
DOUBLE PRECISION x
C ..
C .. Local Scalars ..
DOUBLE PRECISION ax,bot,c,t,top,x2
C ..
C .. Local Arrays ..
DOUBLE PRECISION a(5),b(3),p(8),q(8),r(5),s(4)
C ..
C .. Intrinsic Functions ..
INTRINSIC abs,exp,sign
C ..
C .. Data statements ..
C-------------------------
C-------------------------
C-------------------------
C-------------------------
DATA c/.564189583547756D0/
DATA a(1)/.771058495001320D-04/,a(2)/-.133733772997339D-02/,
+ a(3)/.323076579225834D-01/,a(4)/.479137145607681D-01/,
+ a(5)/.128379167095513D+00/
DATA b(1)/.301048631703895D-02/,b(2)/.538971687740286D-01/,
+ b(3)/.375795757275549D+00/
DATA p(1)/-1.36864857382717D-07/,p(2)/5.64195517478974D-01/,
+ p(3)/7.21175825088309D+00/,p(4)/4.31622272220567D+01/,
+ p(5)/1.52989285046940D+02/,p(6)/3.39320816734344D+02/,
+ p(7)/4.51918953711873D+02/,p(8)/3.00459261020162D+02/
DATA q(1)/1.00000000000000D+00/,q(2)/1.27827273196294D+01/,
+ q(3)/7.70001529352295D+01/,q(4)/2.77585444743988D+02/,
+ q(5)/6.38980264465631D+02/,q(6)/9.31354094850610D+02/,
+ q(7)/7.90950925327898D+02/,q(8)/3.00459260956983D+02/
DATA r(1)/2.10144126479064D+00/,r(2)/2.62370141675169D+01/,
+ r(3)/2.13688200555087D+01/,r(4)/4.65807828718470D+00/,
+ r(5)/2.82094791773523D-01/
DATA s(1)/9.41537750555460D+01/,s(2)/1.87114811799590D+02/,
+ s(3)/9.90191814623914D+01/,s(4)/1.80124575948747D+01/
C ..
C .. Executable Statements ..
C-------------------------
ax = abs(x)
IF (ax.GT.0.5D0) GO TO 10
t = x*x
top = ((((a(1)*t+a(2))*t+a(3))*t+a(4))*t+a(5)) + 1.0D0
bot = ((b(1)*t+b(2))*t+b(3))*t + 1.0D0
erf = x* (top/bot)
RETURN
C
10 IF (ax.GT.4.0D0) GO TO 20
top = ((((((p(1)*ax+p(2))*ax+p(3))*ax+p(4))*ax+p(5))*ax+p(6))*ax+
+ p(7))*ax + p(8)
bot = ((((((q(1)*ax+q(2))*ax+q(3))*ax+q(4))*ax+q(5))*ax+q(6))*ax+
+ q(7))*ax + q(8)
erf = 0.5D0 + (0.5D0-exp(-x*x)*top/bot)
IF (x.LT.0.0D0) erf2 = -erf2
RETURN
C
20 IF (ax.GE.5.8D0) GO TO 30
x2 = x*x
t = 1.0D0/x2
top = (((r(1)*t+r(2))*t+r(3))*t+r(4))*t + r(5)
bot = (((s(1)*t+s(2))*t+s(3))*t+s(4))*t + 1.0D0
erf = (c-top/ (x2*bot))/ax
erf = 0.5D0 + (0.5D0-exp(-x2)*erf)
IF (x.LT.0.0D0) erf = -erf
RETURN
C
30 erf = sign(1.0D0,x)
RETURN
END
I'm making a module in python and I want this function to work with numpy arrays too, like scipy does. The only way I found that make this work is creating a subroutine above the code which takes an array and every element is passed to the erf function, and then compile with f2py.
erfmod.f
subroutine erf(a,n)
implicit none
integer :: n,i
real*8 :: a(n)
Cf2py intent(in,out,copy) :: a
cf2py integer intent(hide),depend(a) :: n=len(a)
do i=1, n
a(i) = erf2(a(i))
end do
contains
DOUBLE PRECISION FUNCTION erf2(x)
C-----------------------------------------------------------------------
C EVALUATION OF THE REAL ERROR FUNCTION
C-----------------------------------------------------------------------
C .. Scalar Arguments ..
DOUBLE PRECISION x
C ..
C .. Local Scalars ..
DOUBLE PRECISION ax,bot,c,t,top,x2
C ..
C .. Local Arrays ..
DOUBLE PRECISION a(5),b(3),p(8),q(8),r(5),s(4)
C ..
C .. Intrinsic Functions ..
INTRINSIC abs,exp,sign
C ..
C .. Data statements ..
C-------------------------
C-------------------------
C-------------------------
C-------------------------
DATA c/.564189583547756D0/
DATA a(1)/.771058495001320D-04/,a(2)/-.133733772997339D-02/,
+ a(3)/.323076579225834D-01/,a(4)/.479137145607681D-01/,
+ a(5)/.128379167095513D+00/
DATA b(1)/.301048631703895D-02/,b(2)/.538971687740286D-01/,
+ b(3)/.375795757275549D+00/
DATA p(1)/-1.36864857382717D-07/,p(2)/5.64195517478974D-01/,
+ p(3)/7.21175825088309D+00/,p(4)/4.31622272220567D+01/,
+ p(5)/1.52989285046940D+02/,p(6)/3.39320816734344D+02/,
+ p(7)/4.51918953711873D+02/,p(8)/3.00459261020162D+02/
DATA q(1)/1.00000000000000D+00/,q(2)/1.27827273196294D+01/,
+ q(3)/7.70001529352295D+01/,q(4)/2.77585444743988D+02/,
+ q(5)/6.38980264465631D+02/,q(6)/9.31354094850610D+02/,
+ q(7)/7.90950925327898D+02/,q(8)/3.00459260956983D+02/
DATA r(1)/2.10144126479064D+00/,r(2)/2.62370141675169D+01/,
+ r(3)/2.13688200555087D+01/,r(4)/4.65807828718470D+00/,
+ r(5)/2.82094791773523D-01/
DATA s(1)/9.41537750555460D+01/,s(2)/1.87114811799590D+02/,
+ s(3)/9.90191814623914D+01/,s(4)/1.80124575948747D+01/
C ..
C .. Executable Statements ..
C-------------------------
ax = abs(x)
IF (ax.GT.0.5D0) GO TO 10
t = x*x
top = ((((a(1)*t+a(2))*t+a(3))*t+a(4))*t+a(5)) + 1.0D0
bot = ((b(1)*t+b(2))*t+b(3))*t + 1.0D0
erf2 = x* (top/bot)
RETURN
C
10 IF (ax.GT.4.0D0) GO TO 20
top = ((((((p(1)*ax+p(2))*ax+p(3))*ax+p(4))*ax+p(5))*ax+p(6))*ax+
+ p(7))*ax + p(8)
bot = ((((((q(1)*ax+q(2))*ax+q(3))*ax+q(4))*ax+q(5))*ax+q(6))*ax+
+ q(7))*ax + q(8)
erf2 = 0.5D0 + (0.5D0-exp(-x*x)*top/bot)
IF (x.LT.0.0D0) erf2 = -erf2
RETURN
C
20 IF (ax.GE.5.8D0) GO TO 30
x2 = x*x
t = 1.0D0/x2
top = (((r(1)*t+r(2))*t+r(3))*t+r(4))*t + r(5)
bot = (((s(1)*t+s(2))*t+s(3))*t+s(4))*t + 1.0D0
erf2 = (c-top/ (x2*bot))/ax
erf2 = 0.5D0 + (0.5D0-exp(-x2)*erf2)
IF (x.LT.0.0D0) erf2 = -erf2
RETURN
C
30 erf2 = sign(1.0D0,x)
RETURN
end function erf2
end subroutine
after compiling with f2py
import module
print(module.erfmod(0.5))
print(module.erfmod(np.array([0.5])))
>>> array(0.52049988)
>>> array(0.52049988)
shoud look like this:
import module
print(module.erfmod(0.5))
print(module.erfmod(np.array([0.5])))
>>> 0.5204998778130465
>>> array(0.52049988)
But now I lost the precision when I'm passing a float, the result is an array with less digits. Scipy somehow manage to return a float when I'm passing a float, and an array when I pass a numpy array (the second case). How can I get the same result?
There are no precision issues here. You did not lost any precision. Python just decided to print fewer digits for you in the text representation array(0.52049988), that's all. The values are identical.
In [1]: import numpy as np
In [2]: np.array(0.5204998778130465)
Out[2]: array(0.52049988)
In [3]: a = np.array(0.5204998778130465)
In [4]: a
Out[4]: array(0.52049988)
In [5]: a.item()
Out[5]: 0.5204998778130465