[SciPy-user] exponential integral behaviour near -20
David Cournapeau
david@ar.media.kyoto-u.ac...
Mon Dec 15 23:41:36 CST 2008
Robert Kern wrote:
>
> Hmm. I dunno. Perhaps gfortran's runtime changed a branch cut for
> CDLOG() between 4.2 and 4.3. That would explain the +/- pi*I in the
> output.
Actually, on my machine, the first goes through the code path with
CDLOG, and the other one with the code path using CDEXP. I used a simple
program fortran program:
C
**********************************
SUBROUTINE
E1Z(Z,CE1)
C
C
====================================================
C Purpose: Compute complex exponential integral
E1(z)
C Input : z --- Argument of
E1(z)
C Output: CE1 ---
E1(z)
C
====================================================
C
IMPLICIT COMPLEX*16
(C,Z)
IMPLICIT DOUBLE PRECISION
(D-H,O-Y)
PI=3.141592653589793D0
EL=0.5772156649015328D0
X=DBLE(Z)
A0=CDABS(Z)
IF (A0.EQ.0.0D0)
THEN
CE1=(1.0D+300,0.0D0)
ELSE IF (A0.LE.10.0.OR.X.LT.0.0.AND.A0.LT.20.0)
THEN
CE1=(1.0D0,0.0D0)
CR=(1.0D0,0.0D0)
DO 10
K=1,150
CR=-CR*K*Z/(K+1.0D0)**2
CE1=CE1+CR
IF (CDABS(CR).LE.CDABS(CE1)*1.0D-15) GO TO
15
10
CONTINUE
15
CE1=-EL-CDLOG(Z)+Z*CE1
ELSE
CT0=(0.0D0,0.0D0)
DO 20
K=120,1,-1
CT0=K/(1.0D0+K/(Z+CT0))
20
CONTINUE
CT=1.0D0/(Z+CT0)
CE1=CDEXP(-Z)*CT
IF (X.LE.0.0.AND.DIMAG(Z).EQ.0.0)
CE1=CE1-PI*(0.0D0,1.0D0)
ENDIF
RETURN
END
program
main
complex*16 a,
b
complex*16 ra,
rb
a = (-19.9999990D0,
0D0)
b = (-19.9999991D0,
0D0)
call E1Z(a,
ra)
call E1Z(b,
rb)
print *, DIMAG(ra), DIMAG(rb)
end
In that case, I have the same values for ra and rb complex parts (I have
both gfortran 4.3 and 4.2, and tried both, no difference, even with
optimization flags -O3 -funroll-loops as used by distutils). I noticed
that both values use different code paths, though (the one using CDLOG
for a, the one using CDEXP for b).
I don't know why there is a discrepancy between this fortran program and
scipy - can the C/F transition cause this ?
David
More information about the SciPy-user
mailing list