import math


def mmd_formula(D, k):
    """ Calculate Myrheim-Meyer f_D formula
    
    Parameters
    ----------
    
    D : Minkowski spacetime dimension (space + time)
    k : Length of chains being counted
    
    Returns
    -------
    
    S_k / N^k : Expected number of k-chains per N^k vertices in the interval
    For k=2 this is equal to the ordering fraction"""
    if D < 1:
        raise ValueError('D value %s less than 1' % D)
    if k==2:
        top = math.gamma(D + 1.) * math.gamma(D / 2.)
        bottom = 4. * math.gamma(D * 1.5)
    else:
        top = math.gamma(D / 2.) * math.gamma(D) * math.gamma(D+1.)**(k-1.)
        bottom = 2**(k-1) * k * math.gamma(0.5 * k * D) * math.gamma(0.5 * D * (k + 1))
        
    return top / bottom