# [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
```