[SciPy-user] exponential integral behaviour near -20

David Joyner wdj@usna....
Tue Dec 16 05:51:59 CST 2008

Thanks very much Robert!

---- Original message ----
>Date: Tue, 16 Dec 2008 00:20:44 -0600
>From: "Robert Kern" <robert.kern@gmail.com>  
>Subject: Re: [SciPy-user] exponential integral behaviour near -20  
>To: "SciPy Users List" <scipy-user@scipy.org>
>On Mon, Dec 15, 2008 at 23:41, David Cournapeau
><david@ar.media.kyoto-u.ac.jp> wrote:
>> 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.
>Ah, I think found it using this clue. It's a bug in SPECFUN. The
>"IMPLICIT DOUBLE PRECISION" statement is missing "A" so A0 is REAL
>rather than DOUBLE. Fixing that makes both of them go through the same
>code path. Can you change the line to this:
>in your specfun.f file, and rebuild scipy?
>> 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 ?
>I don't think so.
>Robert Kern
>"I have come to believe that the whole world is an enigma, a harmless
>enigma that is made terrible by our own mad attempt to interpret it as
>though it had an underlying truth."
>  -- Umberto Eco
>SciPy-user mailing list

More information about the SciPy-user mailing list