import mpmath


def mode(p, b, loc=0, scale=1):
    """
    Mode of the generalized inverse Gaussian distribution.

    The mode is:

                    p - 1 + sqrt((p - 1)**2 + b**2)
        loc + scale -------------------------------
                                  b
    """
    p = mpmath.mpf(p)
    b = mpmath.mpf(b)
    loc = mpmath.mpf(loc)
    scale = mpmath.mpf(scale)

    return loc + scale*(p - 1 + mpmath.sqrt((p - 1)**2 + b**2))/b