def METHOD_NAME(A, V, H, f, p):
    """Extend an existing Arnoldi factorization."""
    k = len(V)
    res = f.norm()[0]
    # the explicit "constant" mode is needed for numpy 1.16
    # mode only gained a default value with numpy 1.17
    H = np.pad(H, ((0, p), (0, p)), mode='constant')
    H[k, k - 1] = res
    v = f * (1 / res)
    V = V.copy()
    V.append(v)
    for i in range(k, k + p):
        v = A.apply(v)
        V.append(v)
        _, R = gram_schmidt(V, return_R=True, atol=0, rtol=0, offset=len(V) - 1, copy=False)
        H[:i + 2, i] = R[:k + p, i + 1]
        v = V[-1]
    return V[:k + p], H, v * R[k + p, k + p]