code: def fidelity(rho0: Density, rho1: Density) -> float:
      
      assert rho0.qubit_nb == rho1.qubit_nb   # FIXME
  
      rho1 = rho1.permute(rho0.qubits)
      op0 = asarray(rho0.asoperator())
      op1 = asarray(rho1.asoperator())
  
      fid = np.real((np.trace(sqrtm(sqrtm(op0) @ op1 @ sqrtm(op0)))) ** 2)
      fid = min(fid, 1.0)
      fid = max(fid, 0.0)     # Correct for rounding errors
  
      return fid
docs: Return the fidelity F(rho0, rho1) between two mixed quantum states.
  
      Note: Fidelity cannot be calculated entirely within the tensor backend.