[SciPy-user] exponential integral behaviour near -20

Robert Kern robert.kern@gmail....
Tue Dec 16 00:20:44 CST 2008


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:

          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


More information about the SciPy-user mailing list