import decimal


def RichardsonExtrapolation(fh, fhn, n, k):
    """Compute the Richardson extrapolation based on two approximations of order k
    where the finite difference parameter h is used in fh and h/n in fhn.
    Inputs:
    fh:  Approximation using h
    fhn: Approximation using h/n
    n:   divisor of h
    k:   original order of approximation
    
    Returns:
    Richardson estimate of order k+1"""
    n = decimal.Decimal(n)
    k = decimal.Decimal(k)
    numerator = decimal.Decimal(n**k * decimal.Decimal(fhn) - decimal.Decimal(fh))
    denominator = decimal.Decimal(n**k - decimal.Decimal(1.0))
    return float(numerator/denominator)