Repository: microsoft/qsharp
Pull request #1215 (2024-02-28 16:13:02)

QSC: display the path when emitting HIR/QIR fails
I wasted a little time on this error as I could not figure out why the emit was failing for me, so I hope the change makes sense to everyone.

When emitting HIR/QIR fails (e.g. directory permissions or output path is invalid) it is not obvious where it failed because the code swallows the path:

```
Error:   × could not emit QIR
  ╰─▶ Permission denied (os error 13)
```

or 

```
Error:   × could not emit QIR
  ╰─▶ No such file or directory (os error 2)
```

With this change the path is also shown:

```
Error:   × could not emit QIR file `./../qir/qir.ll`
  ╰─▶ No such file or directory (os error 2)
```

This is symmetrical to how the path is already shown for source file errors

```
Error:   × could not read source file `../source/foo.qs`
  ╰─▶ No such file or directory (os error 2)
```


Patch:
diff --git a/Cargo.lock b/Cargo.lock
index 9c95bae8..f4853e39 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -971,6 +971,7 @@ dependencies = [
  "miette",
  "num-bigint",
  "num-complex",
+ "num-traits",
  "qsc_data_structures",
  "qsc_fir",
  "qsc_frontend",
diff --git a/compiler/qsc/src/bin/qsc.rs b/compiler/qsc/src/bin/qsc.rs
index b1f3fdd4..18ff3a78 100644
--- a/compiler/qsc/src/bin/qsc.rs
+++ b/compiler/qsc/src/bin/qsc.rs
@@ -165,12 +165,12 @@ fn read_source(path: impl AsRef<Path>) -> miette::Result<(SourceName, SourceCont
 fn emit_hir(package: &Package, dir: impl AsRef<Path>) -> miette::Result<()> {
     let path = dir.as_ref().join("hir.txt");
     info!(
-        "Writing hir output file to: {}",
+        "Writing HIR output file to: {}",
         path.to_str().unwrap_or_default()
     );
-    fs::write(path, package.to_string())
+    fs::write(&path, package.to_string())
         .into_diagnostic()
-        .context("could not emit HIR")
+        .with_context(|| format!("could not emit HIR file `{}`", path.display()))
 }
 
 fn emit_qir(out_dir: &Path, store: &PackageStore, package_id: PackageId) -> Result<(), Report> {
@@ -179,13 +179,12 @@ fn emit_qir(out_dir: &Path, store: &PackageStore, package_id: PackageId) -> Resu
     match result {
         Ok(qir) => {
             info!(
-                "Writing qir output file to: {}",
+                "Writing QIR output file to: {}",
                 path.to_str().unwrap_or_default()
             );
-            fs::write(path, qir)
+            fs::write(&path, qir)
                 .into_diagnostic()
-                .context("could not emit QIR")?;
-            Ok(())
+                .with_context(|| format!("could not emit QIR file `{}`", path.display()))
         }
         Err((error, _)) => {
             let unit = store.get(package_id).expect("package should be in store");
diff --git a/compiler/qsc_eval/Cargo.toml b/compiler/qsc_eval/Cargo.toml
index 6f4aaf29..8bd73470 100644
--- a/compiler/qsc_eval/Cargo.toml
+++ b/compiler/qsc_eval/Cargo.toml
@@ -12,6 +12,7 @@ license.workspace = true
 miette = { workspace = true }
 num-bigint = { workspace = true }
 num-complex = { workspace = true }
+num-traits = { workspace = true }
 quantum-sparse-sim = { workspace = true }
 qsc_data_structures = { path = "../qsc_data_structures" }
 qsc_fir = { path = "../qsc_fir" }
diff --git a/compiler/qsc_eval/src/intrinsic.rs b/compiler/qsc_eval/src/intrinsic.rs
index d817dd84..6d63de31 100644
--- a/compiler/qsc_eval/src/intrinsic.rs
+++ b/compiler/qsc_eval/src/intrinsic.rs
@@ -1,6 +1,8 @@
 // Copyright (c) Microsoft Corporation.
 // Licensed under the MIT License.
 
+mod utils;
+
 #[cfg(test)]
 mod tests;
 
@@ -13,6 +15,7 @@ use crate::{
 };
 use num_bigint::BigInt;
 use rand::{rngs::StdRng, Rng};
+use rustc_hash::FxHashSet;
 use std::array;
 
 #[allow(clippy::too_many_lines)]
@@ -40,6 +43,23 @@ pub(crate) fn call(
                 Err(_) => Err(Error::OutputFail(name_span)),
             }
         }
+        "DumpRegister" => {
+            let qubits = arg.unwrap_array();
+            let qubits = qubits
+                .iter()
+                .map(|q| q.clone().unwrap_qubit().0)
+                .collect::<Vec<_>>();
+            if qubits.len() != qubits.iter().collect::<FxHashSet<_>>().len() {
+                return Err(Error::QubitUniqueness(arg_span));
+            }
+            let (state, qubit_count) = sim.capture_quantum_state();
+            let state = utils::split_state(&qubits, state, qubit_count)
+                .map_err(|()| Error::QubitsNotSeparable(arg_span))?;
+            match out.state(state, qubits.len()) {
+                Ok(()) => Ok(Value::unit()),
+                Err(_) => Err(Error::OutputFail(name_span)),
+            }
+        }
         "Message" => match out.message(&arg.unwrap_string()) {
             Ok(()) => Ok(Value::unit()),
             Err(_) => Err(Error::OutputFail(name_span)),
diff --git a/compiler/qsc_eval/src/intrinsic/utils.rs b/compiler/qsc_eval/src/intrinsic/utils.rs
new file mode 100644
index 00000000..6730593b
--- /dev/null
+++ b/compiler/qsc_eval/src/intrinsic/utils.rs
@@ -0,0 +1,177 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT License.
+
+use std::collections::hash_map::Entry;
+
+use num_bigint::BigUint;
+use num_complex::{Complex, Complex64};
+use num_traits::{One, Zero};
+use rustc_hash::FxHashMap;
+
+/// Given a state and a set of qubits, split the state into two parts: the qubits to dump and the remaining qubits.
+/// This function will return an error if the state is not separable using the provided qubit identifiers.
+pub fn split_state(
+    qubits: &[usize],
+    state: Vec<(BigUint, Complex64)>,
+    qubit_count: usize,
+) -> Result<Vec<(BigUint, Complex64)>, ()> {
+    let state = state.into_iter().collect::<FxHashMap<_, _>>();
+    let mut dump_state = FxHashMap::default();
+    let mut other_state = FxHashMap::default();
+
+    // Compute the mask for the qubits to dump and the mask for the other qubits.
+    let (dump_mask, other_mask) = compute_mask(qubit_count, qubits);
+
+    // Try to split out the state for the given qubits from the whole state, detecting any entanglement
+    // and returning an error if the qubits are not separable.
+    let dump_norm = collect_split_state(
+        &state,
+        &dump_mask,
+        &other_mask,
+        &mut dump_state,
+        &mut other_state,
+    )?;
+
+    // If the product of the collected states is not equal to the total number of input states, then that
+    // implies some states are zero amplitude that would have to be non-zero for the state to be separable.
+    if state.len() != dump_state.len() * other_state.len() {
+        return Err(());
+    }
+
+    let dump_norm = 1.0 / dump_norm.sqrt();
+    Ok(dump_state
+        .into_iter()
+        .filter_map(|(label, val)| {
+            normalize_and_reorder(val, dump_norm, qubits, &label, qubit_count)
+        })
+        .collect())
+}
+
+/// From the qubit identifiers provided, compute the bit masks for the qubits to dump and the remaining qubits.
+/// These masks can be applied to the state labels to separate the label into the two parts needed.
+fn compute_mask(qubit_count: usize, qubits: &[usize]) -> (BigUint, BigUint) {
+    let mut dump_mask = BigUint::zero();
+    let mut other_mask = BigUint::zero();
+    for q in 0..qubit_count {
+        // Note that the qubit order is reversed to match the order of the qubits in the state.
+        if qubits.contains(&q) {
+            dump_mask.set_bit((qubit_count - q - 1) as u64, true);
+        } else {
+            other_mask.set_bit((qubit_count - q - 1) as u64, true);
+        }
+    }
+    (dump_mask, other_mask)
+}
+
+/// Iterates through the given state and for each entry uses the mask to calculate what the separated labels would be
+/// and finds the amplitude for each separated state. If the state is not separable, returns an error.
+/// On success, the `dump_state` and `other_state` maps will be populated with the separated states, and the
+/// function returns the accumulated norm of the dump state.
+fn collect_split_state(
+    state: &FxHashMap<BigUint, Complex64>,
+    dump_mask: &BigUint,
+    other_mask: &BigUint,
+    dump_state: &mut FxHashMap<BigUint, Complex64>,
+    other_state: &mut FxHashMap<BigUint, Complex64>,
+) -> Result<f64, ()> {
+    let mut state_iter = state.iter();
+    let (base_label, base_val) = state_iter.next().expect("state should never be empty");
+    let dump_base_label = base_label & dump_mask;
+    let other_base_label = base_label & other_mask;
+    let mut dump_norm = 1.0_f64;
+
+    // Start with an amplitude of 1 in the first expected split states. This becomes the basis
+    // of the split later, but will get normalized as part of collecting the remaining states.
+    dump_state.insert(dump_base_label.clone(), Complex64::one());
+    other_state.insert(other_base_label.clone(), Complex64::one());
+
+    for (curr_label, curr_val) in state_iter {
+        let dump_label = curr_label & dump_mask;
+        let other_label = curr_label & other_mask;
+
+        // If either the state identified by the dump mask or the state identified by the other mask
+        // is None, that means it has zero amplitude and we can conclude the state is not separable.
+        let Some(dump_val) = state.get(&(&dump_label | &other_base_label)) else {
+            return Err(());
+        };
+        let Some(other_val) = state.get(&(&dump_base_label | &other_label)) else {
+            return Err(());
+        };
+
+        if !(dump_val * other_val - base_val * curr_val)
+            .norm()
+            .is_nearly_zero()
+        {
+            // Coefficients are not equal, so the state is not separable.
+            return Err(());
+        }
+
+        if let Entry::Vacant(entry) = dump_state.entry(dump_label) {
+            // When capturing the amplitude for the dump state, we must divide out the amplitude for the other
+            // state, and vice-versa below.
+            let amplitude = curr_val / other_val;
+            let norm = amplitude.norm();
+            if !norm.is_nearly_zero() {
+                entry.insert(amplitude);
+                dump_norm += norm;
+            }
+        }
+        if let Entry::Vacant(entry) = other_state.entry(other_label) {
+            let amplitude = curr_val / dump_val;
+            let norm = amplitude.norm();
+            if !norm.is_nearly_zero() {
+                entry.insert(amplitude);
+            }
+        }
+    }
+    Ok(dump_norm)
+}
+
+/// Given a dump state amplitude, the normalization factor, the qubits to dump, the label, and the qubit count,
+/// normalize the amplitude and reorder the label to match the provided qubit order.
+/// Specifically, qubits in the requested array may not be in the same allocation order that is used in the state
+/// labels, so the bits in the label must be reordered to match the qubit order.
+fn normalize_and_reorder(
+    val: Complex64,
+    dump_norm: f64,
+    qubits: &[usize],
+    label: &BigUint,
+    qubit_count: usize,
+) -> Option<(BigUint, Complex64)> {
+    // Normalize the dump state by the collected factor.
+    let new_val = val * dump_norm;
+    // Drop any zero amplitude states.
+    if new_val.is_nearly_zero() {
+        None
+    } else {
+        // Reorder the bits in the label to match the provided qubit order.
+        let mut new_label = BigUint::zero();
+        for (i, q) in qubits.iter().enumerate() {
+            // Note that the qubit order is reversed to match the order of the qubits in the state.
+            if label.bit((qubit_count - *q - 1) as u64) {
+                new_label.set_bit((qubits.len() - i - 1) as u64, true);
+            }
+        }
+
+        Some((new_label, new_val))
+    }
+}
+
+trait NearlyZero {
+    fn is_nearly_zero(&self) -> bool;
+}
+
+impl NearlyZero for f64 {
+    fn is_nearly_zero(&self) -> bool {
+        self.abs() <= 1e-10
+    }
+}
+
+impl<T> NearlyZero for Complex<T>
+where
+    T: NearlyZero,
+{
+    fn is_nearly_zero(&self) -> bool {
+        self.re.is_nearly_zero() && self.im.is_nearly_zero()
+    }
+}
diff --git a/compiler/qsc_eval/src/lib.rs b/compiler/qsc_eval/src/lib.rs
index 05091b92..4811ba8c 100644
--- a/compiler/qsc_eval/src/lib.rs
+++ b/compiler/qsc_eval/src/lib.rs
@@ -85,10 +85,15 @@ pub enum Error {
     #[diagnostic(code("Qsc.Eval.OutputFail"))]
     OutputFail(#[label("failed to generate output")] PackageSpan),
 
-    #[error("qubits in gate invocation are not unique")]
+    #[error("qubits in invocation are not unique")]
     #[diagnostic(code("Qsc.Eval.QubitUniqueness"))]
     QubitUniqueness(#[label] PackageSpan),
 
+    #[error("qubits are not separable")]
+    #[diagnostic(help("subset of qubits provided as arguments must not be entangled with any qubits outside of the subset"))]
+    #[diagnostic(code("Qsc.Eval.QubitsNotSeparable"))]
+    QubitsNotSeparable(#[label] PackageSpan),
+
     #[error("range with step size of zero")]
     #[diagnostic(code("Qsc.Eval.RangeStepZero"))]
     RangeStepZero(#[label("invalid range")] PackageSpan),
@@ -134,6 +139,7 @@ impl Error {
             | Error::InvalidNegativeInt(_, span)
             | Error::OutputFail(span)
             | Error::QubitUniqueness(span)
+            | Error::QubitsNotSeparable(span)
             | Error::RangeStepZero(span)
             | Error::ReleasedQubitNotZero(_, span)
             | Error::UnboundName(span)
diff --git a/library/std/diagnostics.qs b/library/std/diagnostics.qs
index d5b66a26..1a976b0a 100644
--- a/library/std/diagnostics.qs
+++ b/library/std/diagnostics.qs
@@ -4,10 +4,71 @@
 namespace Microsoft.Quantum.Diagnostics {
     open QIR.Intrinsic;
 
+    /// # Summary
+    /// Dumps the current target machine's status.
+    ///
+    /// # Description
+    /// This method allows you to dump information about the current quantum state.
+    /// The actual information generated and the semantics are specific to each target machine.
+    ///
+    /// For the local sparse-state simulator distributed as part of the
+    /// Quantum Development Kit, this method will write the wave function as a
+    /// one-dimensional array of pairs of state indices and complex numbers, in which each element represents
+    /// the amplitudes of the probability of measuring the corresponding state.
+    ///
+    /// # Example
+    /// When run on the sparse-state simulator, the following snippet dumps
+    /// the Bell state (|00〉 + |11〉 ) / √2 to the console:
+    /// ```qsharp
+    /// use left = Qubit();
+    /// use right = Qubit();
+    /// within {
+    ///     H(left);
+    ///     CNOT(left, right);
+    /// } apply {
+    ///     DumpMachine();
+    /// }
+    /// ```
     function DumpMachine() : Unit {
         body intrinsic;
     }
 
+    /// # Summary
+    /// Dumps the current target machine's status associated with the given qubits.
+    ///
+    /// # Input
+    /// ## qubits
+    /// The list of qubits to report.
+    ///
+    /// # Remarks
+    /// This method allows you to dump the information associated with the state of the
+    /// given qubits.
+    ///
+    /// For the local sparse-state simulator distributed as part of the
+    /// Quantum Development Kit, this method will write the
+    /// state of the given qubits (i.e. the wave function of the corresponding subsystem) as a
+    /// one-dimensional array of pairs of state indices and complex numbers, in which each element represents
+    /// the amplitudes of the probability of measuring the corresponding state.
+    /// If the given qubits are entangled with some other qubit and their
+    /// state can't be separated, it fails with a runtime error indicating that the qubits are entangled.
+    ///
+    /// # Example
+    /// When run on the sparse-state simulator, the following snippet dumps
+    /// the Bell state (|00〉 + |11〉 ) / √2 to the console:
+    /// ```qsharp
+    /// use left = Qubit();
+    /// use right = Qubit();
+    /// within {
+    ///     H(left);
+    ///     CNOT(left, right);
+    /// } apply {
+    ///     DumpRegister([left, right]);
+    /// }
+    /// ```
+    function DumpRegister(register : Qubit[]) : Unit {
+        body intrinsic;
+    }
+
     @Config(Unrestricted)
     operation CheckZero(qubit : Qubit) : Bool {
         body intrinsic;
diff --git a/resource_estimator/src/estimates/physical_estimation.rs b/resource_estimator/src/estimates/physical_estimation.rs
index 34b2265a..b27c4f10 100644
--- a/resource_estimator/src/estimates/physical_estimation.rs
+++ b/resource_estimator/src/estimates/physical_estimation.rs
@@ -57,6 +57,13 @@ pub trait ErrorCorrection {
         &self,
         lower_bound: Option<&Self::Parameter>,
     ) -> impl Iterator<Item = Self::Parameter>;
+
+    fn code_parameter_cmp(
+        &self,
+        qubit: &Self::Qubit,
+        p1: &Self::Parameter,
+        p2: &Self::Parameter,
+    ) -> Ordering;
 }
 
 pub trait FactoryBuilder<E: ErrorCorrection> {
@@ -87,7 +94,7 @@ pub trait Factory {
     fn max_code_parameter(&self) -> Self::Parameter;
 }
 
-pub struct PhysicalResourceEstimationResult<E: ErrorCorrection, F: Factory, L: Overhead + Clone> {
+pub struct PhysicalResourceEstimationResult<E: ErrorCorrection, F, L> {
     logical_qubit: LogicalQubit<E>,
     num_cycles: u64,
     factory: Option<F>,
@@ -104,10 +111,11 @@ pub struct PhysicalResourceEstimationResult<E: ErrorCorrection, F: Factory, L: O
     error_budget: ErrorBudget,
 }
 
-impl<E: ErrorCorrection, F: Factory<Parameter = E::Parameter> + Clone, L: Overhead + Clone>
-    PhysicalResourceEstimationResult<E, F, L>
-where
-    E::Parameter: Clone + Ord + Default,
+impl<
+        E: ErrorCorrection<Parameter = impl Clone>,
+        F: Factory<Parameter = E::Parameter> + Clone,
+        L: Overhead + Clone,
+    > PhysicalResourceEstimationResult<E, F, L>
 {
     pub fn new(
         estimation: &PhysicalResourceEstimation<E, impl FactoryBuilder<E, Factory = F>, L>,
@@ -118,7 +126,7 @@ where
         required_logical_qubit_error_rate: f64,
         required_logical_magic_state_error_rate: Option<f64>,
     ) -> Self {
-        // Compute statistics for single T-factory
+        // Compute statistics for single factory
         let magic_states_per_run = factory
             .as_ref()
             .map_or(0, |factory| num_factories * factory.num_output_states());
@@ -139,7 +147,7 @@ where
         };
         let physical_qubits_for_single_factory = factory.as_ref().map_or(0, F::physical_qubits);
 
-        // Compute statistics for all T-factories and total overhead
+        // Compute statistics for all factories and total overhead
         let physical_qubits_for_factories = num_factories * physical_qubits_for_single_factory;
         let physical_qubits_for_algorithm =
             estimation.layout_overhead.logical_qubits() * logical_qubit.physical_qubits();
@@ -170,6 +178,23 @@ where
         }
     }
 
+    pub fn without_factories(
+        estimation: &PhysicalResourceEstimation<E, impl FactoryBuilder<E, Factory = F>, L>,
+        logical_qubit: LogicalQubit<E>,
+        num_cycles: u64,
+        required_logical_qubit_error_rate: f64,
+    ) -> Self {
+        Self::new(
+            estimation,
+            logical_qubit,
+            num_cycles,
+            None,
+            0,
+            required_logical_qubit_error_rate,
+            None,
+        )
+    }
+
     pub fn logical_qubit(&self) -> &LogicalQubit<E> {
         &self.logical_qubit
     }
@@ -247,10 +272,7 @@ where
     }
 }
 
-pub struct PhysicalResourceEstimation<E: ErrorCorrection, Builder: FactoryBuilder<E>, L: Overhead>
-where
-    Builder::Factory: Factory + Clone,
-{
+pub struct PhysicalResourceEstimation<E: ErrorCorrection, Builder, L> {
     // required parameters
     ftp: E,
     qubit: Rc<E::Qubit>,
@@ -264,11 +286,11 @@ where
     max_physical_qubits: Option<u64>,
 }
 
-impl<E: ErrorCorrection, Builder: FactoryBuilder<E>, L: Overhead + Clone>
-    PhysicalResourceEstimation<E, Builder, L>
-where
-    Builder::Factory: Factory<Parameter = E::Parameter> + Clone,
-    E::Parameter: Clone + Ord + Default,
+impl<
+        E: ErrorCorrection<Parameter = impl Clone>,
+        Builder: FactoryBuilder<E, Factory = impl Factory<Parameter = E::Parameter> + Clone>,
+        L: Overhead + Clone,
+    > PhysicalResourceEstimation<E, Builder, L>
 {
     pub fn new(
         ftp: E,
@@ -334,9 +356,9 @@ where
     ) -> Result<Vec<PhysicalResourceEstimationResult<E, Builder::Factory, L>>, Error> {
         let num_cycles_required_by_layout_overhead = self.compute_num_cycles()?;
 
-        // The required T-state error rate is computed by dividing the total
-        // error budget for T states by the number of T-states required for the
-        // algorithm.
+        // The required magic state error rate is computed by dividing the total
+        // error budget for magic states by the number of magic states required
+        // for the algorithm.
         let num_magic_states_per_rotation = self
             .layout_overhead
             .num_magic_states_per_rotation(self.error_budget.rotations());
@@ -412,7 +434,11 @@ where
             // state factory.
             if last_code_parameter
                 .as_ref()
-                .map(|d| *d > code_parameter)
+                .map(|d| {
+                    self.ftp
+                        .code_parameter_cmp(self.qubit.as_ref(), d, &code_parameter)
+                        .is_gt()
+                })
                 .unwrap_or(true)
             {
                 last_factories = self.factory_builder.find_factories(
@@ -422,7 +448,7 @@ where
                     &code_parameter,
                 );
 
-                last_code_parameter = Some(Self::find_highest_code_parameter(&last_factories));
+                last_code_parameter = self.find_highest_code_parameter(&last_factories);
             }
 
             if let Some((factory, _)) = Self::try_pick_factory_with_num_cycles(
@@ -430,20 +456,17 @@ where
                 &logical_qubit,
                 max_num_cycles_allowed,
             ) {
-                // Here we compute the number of T-factories required limited by the
+                // Here we compute the number of factories required limited by the
                 // maximum number of cycles allowed by the duration constraint (and the error rate).
-                let min_num_tfactories =
+                let min_num_factories =
                     self.num_factories(&logical_qubit, &factory, max_num_cycles_allowed);
 
-                let mut num_tfactories = min_num_tfactories;
+                let mut num_factories = min_num_factories;
 
                 loop {
-                    // Based on the num_tfactories we compute the number of cycles required
-                    // which must be smaller than the maximum number of cycles allowed by the
-                    // duration constraint (and the error rate).
                     let num_cycles_required_for_magic_states = self
                         .compute_num_cycles_required_for_magic_states(
-                            num_tfactories,
+                            num_factories,
                             &factory,
                             &logical_qubit,
                         );
@@ -459,25 +482,25 @@ where
                         LogicalQubit::new(&self.ftp, code_parameter.clone(), self.qubit.clone())?,
                         num_cycles,
                         Some(factory.clone()),
-                        num_tfactories,
+                        num_factories,
                         required_logical_qubit_error_rate,
                         Some(required_logical_magic_state_error_rate),
                     );
 
                     let value1 = result.runtime() as f64;
                     let value2 = result.physical_qubits();
-                    let num_t_factory_runs = result.num_factory_runs();
+                    let num_factory_runs = result.num_factory_runs();
                     let point = Point2D::new(result, value1, value2);
                     best_estimation_results.push(point);
 
                     if num_cycles_required_for_magic_states
                         <= num_cycles_required_by_layout_overhead
-                        || num_t_factory_runs <= 1
+                        || num_factory_runs <= 1
                     {
                         break;
                     }
 
-                    num_tfactories += 1;
+                    num_factories += 1;
                 }
             }
         }
@@ -531,9 +554,10 @@ where
                     None,
                 );
             }
-            // The required T-state error rate is computed by dividing the total
-            // error budget for T states by the number of T-states required for the
-            // algorithm.
+
+            // The required magic state error rate is computed by dividing the total
+            // error budget for magic states by the number of magic states required
+            // for the algorithm.
             let required_logical_magic_state_error_rate = self.error_budget.magic_states()
                 / (self
                     .layout_overhead
@@ -603,8 +627,8 @@ where
                 num_cycles,
             )
         {
-            let num_tfactories = self.num_factories(logical_qubit, &factory, num_cycles);
-            return Some((factory, num_cycles, num_tfactories));
+            let num_factories = self.num_factories(logical_qubit, &factory, num_cycles);
+            return Some((factory, num_cycles, num_factories));
         }
         if let Some((factory, num_cycles_required)) = self
             .try_find_factory_for_code_parameter_duration_and_max_factories(
@@ -614,25 +638,23 @@ where
             )
         {
             if num_cycles_required <= max_allowed_num_cycles_for_code_parameter {
-                let num_tfactories =
+                let num_factories =
                     self.num_factories(logical_qubit, &factory, num_cycles_required);
-                return Some((factory, num_cycles_required, num_tfactories));
+                return Some((factory, num_cycles_required, num_factories));
             }
         }
 
         None
     }
 
-    #[allow(clippy::too_many_lines)]
-    pub fn estimate_with_max_duration(
+    fn compute_initial_optimization_values(
         &self,
-        max_duration_in_nanoseconds: u64,
-    ) -> Result<PhysicalResourceEstimationResult<E, Builder::Factory, L>, Error> {
+    ) -> Result<InitialOptimizationValues<E::Parameter>, Error> {
         let num_cycles_required_by_layout_overhead = self.compute_num_cycles()?;
 
-        // The required T-state error rate is computed by dividing the total
-        // error budget for T states by the number of T-states required for the
-        // algorithm.
+        // The required magic state error rate is computed by dividing the total
+        // error budget for magic states by the number of magic states required
+        // for the algorithm.
         let num_magic_states_per_rotation = self
             .layout_overhead
             .num_magic_states_per_rotation(self.error_budget.rotations());
@@ -652,6 +674,28 @@ where
             .compute_code_parameter(&self.qubit, required_logical_qubit_error_rate)
             .map_err(Error::CodeParameterComputationFailed)?;
 
+        Ok(InitialOptimizationValues {
+            num_magic_states_per_rotation,
+            min_code_parameter,
+            num_cycles_required_by_layout_overhead,
+            required_logical_qubit_error_rate,
+            required_logical_magic_state_error_rate,
+        })
+    }
+
+    #[allow(clippy::too_many_lines)]
+    pub fn estimate_with_max_duration(
+        &self,
+        max_duration_in_nanoseconds: u64,
+    ) -> Result<PhysicalResourceEstimationResult<E, Builder::Factory, L>, Error> {
+        let InitialOptimizationValues {
+            num_magic_states_per_rotation,
+            min_code_parameter,
+            num_cycles_required_by_layout_overhead,
+            required_logical_qubit_error_rate,
+            required_logical_magic_state_error_rate,
+        } = self.compute_initial_optimization_values()?;
+
         if self
             .layout_overhead
             .num_magic_states(num_magic_states_per_rotation.unwrap_or_default())
@@ -660,17 +704,14 @@ where
             let logical_qubit =
                 LogicalQubit::new(&self.ftp, min_code_parameter, self.qubit.clone())?;
 
-            if num_cycles_required_by_layout_overhead * (logical_qubit.logical_cycle_time())
-                <= (max_duration_in_nanoseconds)
+            if num_cycles_required_by_layout_overhead * logical_qubit.logical_cycle_time()
+                <= max_duration_in_nanoseconds
             {
-                return Ok(PhysicalResourceEstimationResult::new(
+                return Ok(PhysicalResourceEstimationResult::without_factories(
                     self,
                     logical_qubit,
                     num_cycles_required_by_layout_overhead,
-                    None,
-                    0,
                     required_logical_qubit_error_rate,
-                    None,
                 ));
             }
             return Err(Error::MaxDurationTooSmall);
@@ -722,7 +763,11 @@ where
             // state factory.
             if last_code_parameter
                 .as_ref()
-                .map(|d| *d > code_parameter)
+                .map(|d| {
+                    self.ftp
+                        .code_parameter_cmp(self.qubit.as_ref(), d, &code_parameter)
+                        .is_gt()
+                })
                 .unwrap_or(true)
             {
                 last_factories = self.factory_builder.find_factories(
@@ -732,7 +777,7 @@ where
                     &code_parameter,
                 );
 
-                last_code_parameter = Some(Self::find_highest_code_parameter(&last_factories));
+                last_code_parameter = self.find_highest_code_parameter(&last_factories);
             }
 
             if let Some((factory, _)) = Self::try_pick_factory_with_num_cycles(
@@ -740,14 +785,9 @@ where
                 &logical_qubit,
                 max_num_cycles_allowed,
             ) {
-                // Here we compute the number of T-factories required limited by the
-                // maximum number of cycles allowed by the duration constraint (and the error rate).
                 let num_factories =
                     self.num_factories(&logical_qubit, &factory, max_num_cycles_allowed);
 
-                // Based on the num_tfactories we compute the number of cycles required
-                // which must be smaller than the maximum number of cycles allowed by the
-                // duration constraint (and the error rate).
                 let num_cycles_required_for_magic_states = self
                     .compute_num_cycles_required_for_magic_states(
                         num_factories,
@@ -761,8 +801,8 @@ where
                 let num_cycles = num_cycles_required_for_magic_states
                     .max(num_cycles_required_by_layout_overhead);
 
-                if let Some(max_tfactories) = self.max_factories {
-                    if num_factories > max_tfactories {
+                if let Some(max_factories) = self.max_factories {
+                    if num_factories > max_factories {
                         continue;
                     }
                 }
@@ -794,29 +834,13 @@ where
         &self,
         max_num_qubits: u64,
     ) -> Result<PhysicalResourceEstimationResult<E, Builder::Factory, L>, Error> {
-        let min_num_cycles_required_by_layout_overhead = self.compute_num_cycles()?;
-
-        // The required T-state error rate is computed by dividing the total
-        // error budget for T states by the number of T-states required for the
-        // algorithm.
-        let num_magic_states_per_rotation = self
-            .layout_overhead
-            .num_magic_states_per_rotation(self.error_budget.rotations());
-        let required_logical_magic_state_error_rate = self.error_budget.magic_states()
-            / (self
-                .layout_overhead
-                .num_magic_states(num_magic_states_per_rotation.unwrap_or_default())
-                as f64);
-
-        // Required logical error rate (\eps_{\log} / (Q * C) in the paper)
-        let required_logical_qubit_error_rate = self.error_budget.logical()
-            / ((self.layout_overhead.logical_qubits()) * min_num_cycles_required_by_layout_overhead)
-                as f64;
-
-        let min_code_parameter = self
-            .ftp
-            .compute_code_parameter(&self.qubit, required_logical_qubit_error_rate)
-            .map_err(Error::CodeParameterComputationFailed)?;
+        let InitialOptimizationValues {
+            num_magic_states_per_rotation,
+            min_code_parameter,
+            num_cycles_required_by_layout_overhead,
+            required_logical_qubit_error_rate,
+            required_logical_magic_state_error_rate,
+        } = self.compute_initial_optimization_values()?;
 
         if self
             .layout_overhead
@@ -828,14 +852,11 @@ where
             if self.layout_overhead.logical_qubits() * logical_qubit.physical_qubits()
                 <= max_num_qubits
             {
-                return Ok(PhysicalResourceEstimationResult::new(
+                return Ok(PhysicalResourceEstimationResult::without_factories(
                     self,
                     logical_qubit,
-                    min_num_cycles_required_by_layout_overhead,
-                    None,
-                    0,
+                    num_cycles_required_by_layout_overhead,
                     required_logical_qubit_error_rate,
-                    None,
                 ));
             }
             return Err(Error::MaxPhysicalQubitsTooSmall);
@@ -875,7 +896,7 @@ where
                     * min_allowed_logical_qubit_error_rate))
                 .floor() as u64;
 
-            if max_num_cycles_allowed_by_error_rate < min_num_cycles_required_by_layout_overhead {
+            if max_num_cycles_allowed_by_error_rate < num_cycles_required_by_layout_overhead {
                 continue;
             }
 
@@ -885,7 +906,11 @@ where
             // state factory.
             if last_code_parameter
                 .as_ref()
-                .map(|d| *d > code_parameter)
+                .map(|d| {
+                    self.ftp
+                        .code_parameter_cmp(self.qubit.as_ref(), d, &code_parameter)
+                        .is_gt()
+                })
                 .unwrap_or(true)
             {
                 last_factories = self.factory_builder.find_factories(
@@ -895,7 +920,7 @@ where
                     &code_parameter,
                 );
 
-                last_code_parameter = Some(Self::find_highest_code_parameter(&last_factories));
+                last_code_parameter = self.find_highest_code_parameter(&last_factories);
             }
 
             if let Some(factory) = Self::try_pick_factory_below_or_equal_num_qubits(
@@ -918,7 +943,7 @@ where
                     );
 
                 let num_cycles = num_cycles_required_for_magic_states
-                    .max(min_num_cycles_required_by_layout_overhead);
+                    .max(num_cycles_required_by_layout_overhead);
 
                 if num_cycles > max_num_cycles_allowed_by_error_rate {
                     continue;
@@ -952,6 +977,9 @@ where
         best_estimation_result.ok_or(Error::MaxPhysicalQubitsTooSmall)
     }
 
+    /// Based on num_factories, we compute the number of cycles required which
+    /// must be smaller than the maximum number of cycles allowed by the
+    /// duration constraint (and the error rate).
     fn compute_num_cycles_required_for_magic_states(
         &self,
         num_factories: u64,
@@ -982,7 +1010,7 @@ where
             .min_by(|&p, &q| {
                 p.normalized_volume()
                     .partial_cmp(&q.normalized_volume())
-                    .expect("Could not compare T-factories normalized volume")
+                    .expect("Could not compare factories normalized volume")
             })
             .cloned()
     }
@@ -1023,7 +1051,7 @@ where
             .min_by(|&p, &q| {
                 p.normalized_volume()
                     .partial_cmp(&q.normalized_volume())
-                    .expect("Could not compare T-factories normalized volume")
+                    .expect("Could not compare factories normalized volume")
             })
             .cloned()
     }
@@ -1055,12 +1083,12 @@ where
         factories: &[Builder::Factory],
         logical_qubit: &LogicalQubit<E>,
         max_allowed_num_cycles_for_code_parameter: u64,
-        max_tfactories: u64,
+        max_factories: u64,
     ) -> Option<(Builder::Factory, u64)> {
         factories
             .iter()
             .map(|factory| {
-                let magic_states_per_run = max_tfactories * factory.num_output_states();
+                let magic_states_per_run = max_factories * factory.num_output_states();
                 let num_magic_states_per_rotation = self
                     .layout_overhead
                     .num_magic_states_per_rotation(self.error_budget.rotations());
@@ -1080,11 +1108,11 @@ where
                 let comp1 = p
                     .normalized_volume()
                     .partial_cmp(&q.normalized_volume())
-                    .expect("Could not compare T-factories normalized volume");
+                    .expect("Could not compare factories normalized volume");
                 if comp1 == Ordering::Equal {
                     num_p
                         .partial_cmp(num_q)
-                        .expect("Could not compare T-factories num cycles")
+                        .expect("Could not compare factories num cycles")
                 } else {
                     comp1
                 }
@@ -1107,16 +1135,15 @@ where
             .min_by(|(p, _), (q, _)| {
                 p.normalized_volume()
                     .partial_cmp(&q.normalized_volume())
-                    .expect("Could not compare T-factories normalized volume")
+                    .expect("Could not compare factories normalized volume")
             })
     }
 
-    fn find_highest_code_parameter(factories: &[Builder::Factory]) -> E::Parameter {
+    fn find_highest_code_parameter(&self, factories: &[Builder::Factory]) -> Option<E::Parameter> {
         factories
             .iter()
             .map(|p| p.max_code_parameter())
-            .max()
-            .unwrap_or_default()
+            .max_by(|a, b| self.ftp.code_parameter_cmp(self.qubit.as_ref(), a, b))
     }
 
     // Possibly adjusts number of cycles C from initial starting point C_min
@@ -1135,7 +1162,7 @@ where
             num_cycles = ((num_cycles as f64) * logical_depth_scaling).ceil() as u64;
         }
 
-        // We cannot perform resource estimation when there are neither T-states nor cycles
+        // We cannot perform resource estimation when there are neither magic states nor cycles
         if self
             .layout_overhead
             .num_magic_states(num_magic_states_per_rotation.unwrap_or_default())
@@ -1148,9 +1175,9 @@ where
         Ok(num_cycles)
     }
 
-    // Choose number of T factories to use; we can safely use unwrap on
-    // the t_count here because the algorithm only finds T-factories
-    // that provide this number
+    // Choose number of factories to use; we can safely use unwrap on the number
+    // of magic states here, because the algorithm only finds factories that
+    // provide this number
     fn num_factories(
         &self,
         logical_qubit: &LogicalQubit<E>,
@@ -1185,3 +1212,11 @@ where
         }
     }
 }
+
+struct InitialOptimizationValues<Parameter> {
+    num_magic_states_per_rotation: Option<u64>,
+    min_code_parameter: Parameter,
+    num_cycles_required_by_layout_overhead: u64,
+    required_logical_qubit_error_rate: f64,
+    required_logical_magic_state_error_rate: f64,
+}
diff --git a/resource_estimator/src/system/modeling/fault_tolerance.rs b/resource_estimator/src/system/modeling/fault_tolerance.rs
index 754f3bc7..5ca63222 100644
--- a/resource_estimator/src/system/modeling/fault_tolerance.rs
+++ b/resource_estimator/src/system/modeling/fault_tolerance.rs
@@ -556,4 +556,8 @@ impl ErrorCorrection for Protocol {
     fn code_parameter_range(&self, lower_bound: Option<&u64>) -> impl Iterator<Item = u64> {
         (lower_bound.copied().unwrap_or(1)..=self.max_code_distance).step_by(2)
     }
+
+    fn code_parameter_cmp(&self, _qubit: &PhysicalQubit, p1: &u64, p2: &u64) -> std::cmp::Ordering {
+        p1.cmp(p2)
+    }
 }


Test patch:
diff --git a/compiler/qsc/src/interpret/tests.rs b/compiler/qsc/src/interpret/tests.rs
index ff721388..34c89fb4 100644
--- a/compiler/qsc/src/interpret/tests.rs
+++ b/compiler/qsc/src/interpret/tests.rs
@@ -471,7 +471,7 @@ mod given_interpreter {
                 &result,
                 &output,
                 &expect![[r#"
-                    runtime error: qubits in gate invocation are not unique
+                    runtime error: qubits in invocation are not unique
                        [intrinsic.qs] [(control, target)]
                 "#]],
             );
diff --git a/compiler/qsc_eval/src/intrinsic/tests.rs b/compiler/qsc_eval/src/intrinsic/tests.rs
index b34b6481..295ad60c 100644
--- a/compiler/qsc_eval/src/intrinsic/tests.rs
+++ b/compiler/qsc_eval/src/intrinsic/tests.rs
@@ -293,6 +293,149 @@ fn dump_machine_endianness() {
     );
 }
 
+#[test]
+fn dump_register_all_qubits() {
+    check_intrinsic_output(
+        "",
+        indoc! {"{
+            use qs = Qubit[4];
+            X(qs[1]);
+            Microsoft.Quantum.Diagnostics.DumpRegister(qs);
+            X(qs[1]);
+        }"},
+        &expect![[r#"
+            STATE:
+            |0100⟩: 1.0000+0.0000𝑖
+        "#]],
+    );
+}
+
+#[test]
+fn dump_register_subset_qubits() {
+    check_intrinsic_output(
+        "",
+        indoc! {"{
+            use qs = Qubit[4];
+            X(qs[1]);
+            Microsoft.Quantum.Diagnostics.DumpRegister([qs[1], qs[2]]);
+            X(qs[1]);
+        }"},
+        &expect![[r#"
+            STATE:
+            |10⟩: 1.0000+0.0000𝑖
+        "#]],
+    );
+}
+
+#[test]
+fn dump_register_subset_entangled_within_subset_is_separable() {
+    check_intrinsic_output(
+        "",
+        indoc! {"{
+            use (q1, q2, q3) = (Qubit(), Qubit(), Qubit());
+            H(q1);
+            CNOT(q1, q3);
+            Microsoft.Quantum.Diagnostics.DumpRegister([q1, q3]);
+            Reset(q1);
+            Reset(q2);
+            Reset(q3);
+        }"},
+        &expect![[r#"
+            STATE:
+            |00⟩: 0.7071+0.0000𝑖
+            |11⟩: 0.7071+0.0000𝑖
+        "#]],
+    );
+}
+
+#[test]
+fn dump_register_subset_entangled_with_other_qubits_not_separable() {
+    check_intrinsic_result(
+        "",
+        indoc! {"{
+            use (q1, q2, q3) = (Qubit(), Qubit(), Qubit());
+            H(q1);
+            CNOT(q1, q3);
+            Microsoft.Quantum.Diagnostics.DumpRegister([q1, q2]);
+        }"},
+        &expect!["qubits are not separable"],
+    );
+}
+
+#[test]
+fn dump_register_other_qubits_superposition_is_separable() {
+    check_intrinsic_output(
+        "",
+        indoc! {"{
+            use qs = Qubit[3];
+            H(qs[0]);
+            H(qs[2]);
+            Microsoft.Quantum.Diagnostics.DumpRegister(qs[...1]);
+            ResetAll(qs);
+        }"},
+        &expect![[r#"
+            STATE:
+            |00⟩: 0.7071+0.0000𝑖
+            |10⟩: 0.7071+0.0000𝑖
+        "#]],
+    );
+}
+
+#[test]
+fn dump_register_other_qubits_one_state_is_separable() {
+    check_intrinsic_output(
+        "",
+        indoc! {"{
+            use qs = Qubit[3];
+            H(qs[0]);
+            X(qs[2]);
+            Microsoft.Quantum.Diagnostics.DumpRegister(qs[...1]);
+            ResetAll(qs);
+        }"},
+        &expect![[r#"
+            STATE:
+            |00⟩: 0.7071+0.0000𝑖
+            |10⟩: 0.7071+0.0000𝑖
+        "#]],
+    );
+}
+
+#[test]
+fn dump_register_qubits_reorder_output() {
+    check_intrinsic_output(
+        "",
+        indoc! {"{
+            use qs = Qubit[5];
+            H(qs[0]);
+            X(qs[2]);
+            Microsoft.Quantum.Diagnostics.DumpMachine();
+            Microsoft.Quantum.Diagnostics.DumpRegister(qs[2..-1...]);
+            ResetAll(qs);
+        }"},
+        &expect![[r#"
+            STATE:
+            |00100⟩: 0.7071+0.0000𝑖
+            |10100⟩: 0.7071+0.0000𝑖
+            STATE:
+            |100⟩: 0.7071+0.0000𝑖
+            |101⟩: 0.7071+0.0000𝑖
+        "#]],
+    );
+}
+
+#[test]
+fn dump_register_qubits_not_unique_fails() {
+    check_intrinsic_result(
+        "",
+        indoc! {"{
+            use qs = Qubit[3];
+            H(qs[0]);
+            Microsoft.Quantum.Diagnostics.DumpRegister([qs[0], qs[0]]);
+        }"},
+        &expect!["qubits in invocation are not unique"],
+    );
+}
+
 #[test]
 fn message() {
     check_intrinsic_output(
@@ -1077,7 +1220,7 @@ fn qubit_not_unique_two_qubit_error() {
             use q = Qubit();
             CNOT(q , q);
         }"},
-        &expect!["qubits in gate invocation are not unique"],
+        &expect!["qubits in invocation are not unique"],
     );
 }
 
@@ -1089,7 +1232,7 @@ fn qubit_not_unique_two_qubit_rotation_error() {
             use q = Qubit();
             Rxx(0.1, q, q);
         }"},
-        &expect!["qubits in gate invocation are not unique"],
+        &expect!["qubits in invocation are not unique"],
     );
 }
 
@@ -1102,7 +1245,7 @@ fn qubit_not_unique_three_qubit_error_first_second() {
             use a = Qubit();
             CCNOT(q , q, a);
         }"},
-        &expect!["qubits in gate invocation are not unique"],
+        &expect!["qubits in invocation are not unique"],
     );
 }
 
@@ -1115,7 +1258,7 @@ fn qubit_not_unique_three_qubit_error_first_third() {
             use a = Qubit();
             CCNOT(q , a, q);
         }"},
-        &expect!["qubits in gate invocation are not unique"],
+        &expect!["qubits in invocation are not unique"],
     );
 }
 
@@ -1128,7 +1271,7 @@ fn qubit_not_unique_three_qubit_error_second_third() {
             use a = Qubit();
             CCNOT(a , q, q);
         }"},
-        &expect!["qubits in gate invocation are not unique"],
+        &expect!["qubits in invocation are not unique"],
     );
 }
 
diff --git a/wasm/src/tests.rs b/wasm/src/tests.rs
index c2706f4a..9ea48d73 100644
--- a/wasm/src/tests.rs
+++ b/wasm/src/tests.rs
@@ -242,9 +242,9 @@ fn test_run_error_program_multiple_shots() {
     )
     .expect("code should compile and run");
     expect![[r#"
-        {"result":{"code":"Qsc.Eval.QubitUniqueness","message":"runtime error: qubits in gate invocation are not unique","range":{"end":{"character":1,"line":0},"start":{"character":0,"line":0}},"severity":"error"},"success":false,"type":"Result"}
-        {"result":{"code":"Qsc.Eval.QubitUniqueness","message":"runtime error: qubits in gate invocation are not unique","range":{"end":{"character":1,"line":0},"start":{"character":0,"line":0}},"severity":"error"},"success":false,"type":"Result"}
-        {"result":{"code":"Qsc.Eval.QubitUniqueness","message":"runtime error: qubits in gate invocation are not unique","range":{"end":{"character":1,"line":0},"start":{"character":0,"line":0}},"severity":"error"},"success":false,"type":"Result"}"#]]
+        {"result":{"code":"Qsc.Eval.QubitUniqueness","message":"runtime error: qubits in invocation are not unique","range":{"end":{"character":1,"line":0},"start":{"character":0,"line":0}},"severity":"error"},"success":false,"type":"Result"}
+        {"result":{"code":"Qsc.Eval.QubitUniqueness","message":"runtime error: qubits in invocation are not unique","range":{"end":{"character":1,"line":0},"start":{"character":0,"line":0}},"severity":"error"},"success":false,"type":"Result"}
+        {"result":{"code":"Qsc.Eval.QubitUniqueness","message":"runtime error: qubits in invocation are not unique","range":{"end":{"character":1,"line":0},"start":{"character":0,"line":0}},"severity":"error"},"success":false,"type":"Result"}"#]]
     .assert_eq(&output.join("\n"));
 }