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?
CodePudding user response：
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