[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