import mpmath


def logpdf(x, p, b, loc=0, scale=1):
    """
    Log of the PDF of the generalized inverse Gaussian distribution.

    The PDF for x > loc is:

        z**(p - 1) * exp(-b*(z + 1/z)/2))
        ---------------------------------
                 scale * K_p(b)

    where s is the scale, z = (x - loc)/s, and K_p(b) is the modified Bessel
    function of the second kind.  For x <= loc, the PDF is zero.
    """
    x = mpmath.mpf(x)
    p = mpmath.mpf(p)
    b = mpmath.mpf(b)
    loc = mpmath.mpf(loc)
    scale = mpmath.mpf(scale)

    if x <= loc:
        return -mpmath.mp.inf
    z = (x - loc)/scale
    return ((p - 1)*mpmath.log(z)
            - b*(z + 1/z)/2
            - mpmath.log(2*mpmath.besselk(p, b))
            - mpmath.log(scale))