def phase_estimation(U: np.ndarray, accuracy: int, reg_offset: int = 0) -> Program:
    """
    Generate a circuit for quantum phase estimation.

    :param U: A unitary matrix.
    :param accuracy: Number of bits of accuracy desired.
    :param reg_offset: Where to start writing measurements (default 0).
    :return: A Quil program to perform phase estimation.
    """
    assert isinstance(accuracy, int)
    rows, cols = U.shape
    m = int(log2(rows))
    output_qubits = range(0, accuracy)
    U_qubits = range(accuracy, accuracy + m)
    p = Program()
    ro = p.declare('ro', 'BIT', len(output_qubits))

    # Hadamard initialization
    for i in output_qubits:
        p.inst(H(i))
    # Controlled unitaries
    for i in output_qubits:
        if i > 0:
            U = np.dot(U, U)
        cU = controlled(U)
        name = "CONTROLLED-U{0}".format(2 ** i)
        # define the gate
        p.defgate(name, cU)
        # apply it
        p.inst((name, i) + tuple(U_qubits))
    # Compute the QFT
    p = p + inverse_qft(output_qubits)
    # Perform the measurements
    for i in output_qubits:
        p.measure(i, ro[reg_offset + i])

    return p
Generate a circuit for quantum phase estimation.

    :param U: A unitary matrix.
    :param accuracy: Number of bits of accuracy desired.
    :param reg_offset: Where to start writing measurements (default 0).
    :return: A Quil program to perform phase estimation.