code: def phase_estimation(U: np.ndarray, accuracy: int, reg_offset: int = 0) -> Program:
      
      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
docs: 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.