[SciPy-user] exponential integral behaviour near -20

[David Joyner]

Thanks very much Robert!


---- Original message ----
>Date: Tue, 16 Dec 2008 00:20:44 -0600
>From: "Robert Kern" <robert.kern at gmail.com>  
>Subject: Re: [SciPy-user] exponential integral behaviour near -20  
>To: "SciPy Users List" <scipy-user at scipy.org>
>
>On Mon, Dec 15, 2008 at 23:41, David Cournapeau
><david at 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:
>
>          IMPLICIT DOUBLE PRECISION (A,D-H,O-Y)
>
>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
>SciPy-user at scipy.org
>http://projects.scipy.org/mailman/listinfo/scipy-user