search

haskell / math-functions

Special mathematical functions
http://hackage.haskell.org/package/math-functions
BSD 2-Clause "Simplified" License
40 stars 28 forks source link

Bad asympotic behaviour for incomplete gamma #40

Closed Shimuuar closed 8 years ago

Shimuuar commented 8 years ago 
>>> incompleteGamma 1000 3000
-0.0
Shimuuar commented 8 years ago

Problem lies with Gauss-Legendre integration

    -- For large values of `p' we use 18-point Gauss-Legendre
    -- integration.
    approx
      | ans > 0   = 1 - ans
      | otherwise = -ans
      where
        -- Set upper limit for integration
        xu | x > p1    =         (p1 + 11.5*sqrtP1) `max` (x + 6*sqrtP1)
           | otherwise = max 0 $ (p1 -  7.5*sqrtP1) `min` (x - 5*sqrtP1)
        s = U.sum $ U.zipWith go coefY coefW
        go y w = let t = x + (xu - x)*y
                 in w * exp( -(t-p1) + p1*(log t - lnP1) )
        ans = s * (xu - x) * exp( p1 * (lnP1 - 1) - logGamma p)
        --
        p1     = p - 1
        lnP1   = log  p1
        sqrtP1 = sqrt p1