import mpmath


def pdf(x, dfn, dfd):
    """
    Probability density function of the F distribution.

    `dfn` and `dfd` are the numerator and denominator degrees of freedom, resp.
    """
    if x <= 0:
        return mpmath.mp.zero
    with mpmath.mp.extradps(5):
        x = mpmath.mp.mpf(x)

        dfn = mpmath.mp.mpf(dfn)
        dfd = mpmath.mp.mpf(dfd)
        r = dfn / dfd
        hdfn = dfn / 2
        hdfd = dfd / 2
        p = (r**hdfn
             * x**(hdfn - 1)
             * (1 + r*x)**(-(hdfn + hdfd))
             / mpmath.beta(hdfn, hdfd))
        return p