def p_correction(n, K=10):
    p0 = n * np.log(n)  # leading term
    logp = np.log(p0)
    sqrtp = np.sqrt(p0)
    corr = 0.0
    for gamma in ZETA_GAMMAS[:K]:
        phi = np.arctan(2 * gamma)
        cos_term = np.cos(gamma * logp - phi)
        denom = np.sqrt(0.25 + gamma*gamma)
        corr += cos_term / denom
    correction = -2.0 * sqrtp / logp * corr
    return p0 + correction