Extend PurifiedMixedState for additional data preparation

Chemistry/src/Runtime/JordanWigner/JordanWignerOptimizedBlockEncoding.qs

@@ -370,14 +370,14 @@ namespace Microsoft.Quantum.Chemistry.JordanWigner {
     operation _JordanWignerOptimizedBlockEncodingStatePrep(targetError : Double, optimizedBEGeneratorSystem : OptimizedBEGeneratorSystem, qROMIdxRegister : LittleEndian, qROMGarbage : Qubit[], signQubit : Qubit, selectZControlRegisters : Qubit[], OptimizedBEControlRegisters : Qubit[], pauliBases : Qubit[], pauliBasesIdx : LittleEndian, indexRegisters : LittleEndian[]) : Unit {
 
         body (...) {
-            let (nTerms, oneNorm0, intToGenIdx) = optimizedBEGeneratorSystem!;
+            let (nTerms, _, _) = optimizedBEGeneratorSystem!;
             let coefficients = _OptimizedBEGeneratorSystemCoeff_(optimizedBEGeneratorSystem);
-            let (qROMQubitCount, oneNorm, qROMUnitary) = QuantumROM(targetError, coefficients);
+            let purifiedState = PurifiedMixedState(targetError, coefficients);
             let unitaryGenerator = (nTerms, _ToJordanWignerSelectInput_(_, optimizedBEGeneratorSystem));
             let pauliBasesUnitaryGenerator = (5, _ToPauliBases_);
 
             //let multiplexer = MultiplexerFromGenerator;
-            qROMUnitary(qROMIdxRegister, qROMGarbage);
+            purifiedState::Prepare(qROMIdxRegister, new Qubit[0], qROMGarbage);
             MultiplexOperationsFromGenerator(unitaryGenerator, qROMIdxRegister, (signQubit, selectZControlRegisters, OptimizedBEControlRegisters, pauliBasesIdx, indexRegisters));
             MultiplexOperationsFromGenerator(pauliBasesUnitaryGenerator, pauliBasesIdx, pauliBases);
         }
@@ -389,8 +389,8 @@ namespace Microsoft.Quantum.Chemistry.JordanWigner {
 
 
     function _JordanWignerOptimizedBlockEncodingQubitManager_ (targetError : Double, nCoeffs : Int, nZ : Int, nMaj : Int, nIdxRegQubits : Int, ctrlRegister : Qubit[]) : ((LittleEndian, Qubit[], Qubit, Qubit[], Qubit[], Qubit[], LittleEndian, LittleEndian[]), (Qubit, Qubit[], Qubit[], Qubit[], LittleEndian[]), Qubit[]) {
-        let (_, (nIndexRegister, nGarbageQubits)) = QuantumROMQubitCount(targetError, nCoeffs);
-        let parts = Partitioned([nIndexRegister, nGarbageQubits], ctrlRegister);
+        let requirements = PurifiedMixedStateRequirements(targetError, nCoeffs);
+        let parts = Partitioned([requirements::NIndexQubits, requirements::NGarbageQubits], ctrlRegister);
         let ((qROMIdx, qROMGarbage), rest0) = ((LittleEndian(parts[0]), parts[1]), parts[2]);
         let ((signQubit, selectZControlRegisters, optimizedBEControlRegisters, pauliBases, indexRegisters, tmp), rest1) = _JordanWignerSelectQubitManager_(nZ, nMaj, nIdxRegQubits, rest0, new Qubit[0]);
         let registers = Partitioned([3], rest1);
@@ -401,7 +401,7 @@ namespace Microsoft.Quantum.Chemistry.JordanWigner {
     function _JordanWignerOptimizedBlockEncodingQubitCount_ (targetError : Double, nCoeffs : Int, nZ : Int, nMaj : Int, nIdxRegQubits : Int, nTarget : Int) : ((Int, Int), (Int, Int, Int, Int, Int, Int, Int, Int[], Int)) {
 
         let (nSelectTotal, (a0, a1, a2, a3, a4)) = _JordanWignerSelectQubitCount_(nZ, nMaj, nIdxRegQubits);
-        let (nQROMTotal, (b0, b1)) = QuantumROMQubitCount(targetError, nCoeffs);
+        let (nQROMTotal, (b0, b1)) = (PurifiedMixedStateRequirements(targetError, nCoeffs))!;
         let pauliBasesIdx = 3;
         return (((nSelectTotal + nQROMTotal) + pauliBasesIdx, nTarget), (b0, b1, a0, a1, a2, a3, pauliBasesIdx, a4, nTarget));
     }

Chemistry/tests/SamplesTests/DocsInvokingChemistryLibrary.cs

@@ -19,6 +19,7 @@
 // The System namespace provides a number of useful built-in
 // types and methods that we'll use throughout this sample.
 using System;
+using System.IO;
 
 // We use this for convnience functions for manipulation arrays.
 using System.Linq;
@@ -84,7 +85,7 @@ static void LoadFromBroombridgeFile()
             var root = @"Molecules";
 
             // This deserializes Broombridge.
-            var broombridge = Broombridge.Deserializers.DeserializeBroombridge($@"{root}\{filename}");
+            var broombridge = Broombridge.Deserializers.DeserializeBroombridge(Path.Combine(root, filename));
 
             // Note that the deserializer returns a list of `ProblemDescriptions` instances 
             // as the file might describe multiple Hamiltonians. In this example, there is 
@@ -122,7 +123,7 @@ static void LoadFromLiquidFile()
             var root = @"Molecules";
 
             // Deserialize the LiQuiD format.
-            var problem = LiQuiD.Deserialize($@"{root}\{filename}").First();
+            var problem = LiQuiD.Deserialize(Path.Combine(root, filename)).First();
 
             // This extracts the `OrbitalIntegralHamiltonian` from problem
             // description format.
@@ -135,7 +136,7 @@ static void LoadFromLiquidFile()
         static void ResourceEstimate()
         {
             // Filename of Hamiltonian to be loaded.
-            var filename = "Molecules/LiH_0.2.yaml";
+            var filename = Path.Combine("Molecules", "LiH_0.2.yaml");
 
             // This deserializes Broombridge.
             var problem = Broombridge.Deserializers.DeserializeBroombridge(filename).ProblemDescriptions.First();

Standard/src/Arithmetic/Asserts.qs

@@ -7,6 +7,7 @@ namespace Microsoft.Quantum.Arithmetic {
     open Microsoft.Quantum.Convert;
     open Microsoft.Quantum.Arrays;
     open Microsoft.Quantum.Diagnostics;
+    open Microsoft.Quantum.Math;
 
     /// # Summary
     /// Asserts that the probability of a specific state of a quantum register has the
@@ -40,16 +41,23 @@ namespace Microsoft.Quantum.Arithmetic {
     /// - `AssertProbInt(0,0.125,qubits,10e-10);`
     /// - `AssertProbInt(6,0.875,qubits,10e-10);`
     operation AssertProbInt(stateIndex : Int, expected : Double, qubits : LittleEndian, tolerance : Double) : Unit {
-        let nQubits = Length(qubits!);
-        let bits = IntAsBoolArray(stateIndex, nQubits);
-
         using (flag = Qubit()) {
-            (ControlledOnBitString(bits, X))(qubits!, flag);
-            AssertMeasurementProbability([PauliZ], [flag], One, expected, $"AssertProbInt failed on stateIndex {stateIndex}, expected probability {expected}.", tolerance);
+            within {
+                (ControlledOnInt(stateIndex, X))(qubits!, flag);
+            } apply {
+                AssertMeasurementProbability([PauliZ], [flag], One, expected, $"AssertProbInt failed on stateIndex {stateIndex}, expected probability {expected}.", tolerance);
+            }
+        }
+    }
 
-            // Uncompute flag qubit.
-            (ControlledOnBitString(bits, X))(qubits!, flag);
-            Reset(flag);
+    operation AssertSignedProbInt(stateIndex : Int, expected : Double, sign : Qubit, qubits : LittleEndian, tolerance : Double) : Unit {
+        using (flag = Qubit()) {
+            let signOffset = expected < 0.0 ? 1 <<< Length(qubits!) | 0;
+            within {
+                (ControlledOnInt(stateIndex + signOffset, X))(qubits! + [sign], flag);
+            } apply {
+                AssertMeasurementProbability([PauliZ], [flag], One, AbsD(expected), $"AssertSignedProbInt failed on stateIndex {stateIndex}, expected probability {expected}.", tolerance);
+            }
         }
     }
 

Standard/src/Canon/Multiplexer.qs

@@ -5,6 +5,7 @@ namespace Microsoft.Quantum.Canon {
     open Microsoft.Quantum.Intrinsic;
     open Microsoft.Quantum.Arithmetic;
     open Microsoft.Quantum.Arrays;
+    open Microsoft.Quantum.Diagnostics;
     open Microsoft.Quantum.Math;
 
     /// # Summary
@@ -43,29 +44,29 @@ namespace Microsoft.Quantum.Canon {
         }
         if (nUnitaries > 0) {
             let auxiliary = new Qubit[0];
-            Adjoint _MultiplexOperationsFromGenerator(unitaryGeneratorWithOffset, auxiliary, index, target);
+            Adjoint MultiplexOperationsFromGeneratorImpl(unitaryGeneratorWithOffset, auxiliary, index, target);
         }
     }
 
     /// # Summary
     /// Implementation step of `MultiplexOperationsFromGenerator`.
     /// # See Also
     /// - Microsoft.Quantum.Canon.MultiplexOperationsFromGenerator
-    operation _MultiplexOperationsFromGenerator<'T>(unitaryGenerator : (Int, Int, (Int -> ('T => Unit is Adj + Ctl))), auxiliary: Qubit[], index: LittleEndian, target: 'T)
+    internal operation MultiplexOperationsFromGeneratorImpl<'T>(unitaryGenerator : (Int, Int, (Int -> ('T => Unit is Adj + Ctl))), auxiliary: Qubit[], index: LittleEndian, target: 'T)
     : Unit {
         body (...) {
             let nIndex = Length(index!);
             let nStates = 2^nIndex;
 
             let (nUnitaries, unitaryOffset, unitaryFunction) = unitaryGenerator;
 
-            let nUnitariesLeft = MinI(nUnitaries, nStates/2);
+            let nUnitariesLeft = MinI(nUnitaries, nStates / 2);
             let nUnitariesRight = MinI(nUnitaries, nStates);
 
             let leftUnitaries = (nUnitariesLeft, unitaryOffset, unitaryFunction);
-            let rightUnitaries = (nUnitariesRight-nUnitariesLeft, unitaryOffset + nUnitariesLeft, unitaryFunction);
+            let rightUnitaries = (nUnitariesRight - nUnitariesLeft, unitaryOffset + nUnitariesLeft, unitaryFunction);
 
-            let newControls = LittleEndian(index![0..nIndex - 2]);
+            let newControls = LittleEndian(Most(index!));
 
             if (nUnitaries > 0) {
                 if (Length(auxiliary) == 1 and nIndex == 0) {
@@ -74,38 +75,38 @@ namespace Microsoft.Quantum.Canon {
                     (Controlled Adjoint (unitaryFunction(unitaryOffset)))(auxiliary, target);
                 } elif (Length(auxiliary) == 0 and nIndex >= 1) {
                     // Start case
-                    let newauxiliary = [index![Length(index!) - 1]];
-                    if(nUnitariesRight > 0){
-                        _MultiplexOperationsFromGenerator(rightUnitaries, newauxiliary, newControls, target);
+                    let newauxiliary = Tail(index!);
+                    if (nUnitariesRight > 0) {
+                        MultiplexOperationsFromGeneratorImpl(rightUnitaries, [newauxiliary], newControls, target);
                     }
                     within {
-                        X(newauxiliary[0]);
+                        X(newauxiliary);
                     } apply {
-                        _MultiplexOperationsFromGenerator(leftUnitaries, newauxiliary, newControls, target);
+                        MultiplexOperationsFromGeneratorImpl(leftUnitaries, [newauxiliary], newControls, target);
                     }
                 } else {
                     // Recursion that reduces nIndex by 1 & sets Length(auxiliary) to 1.
-                    using (newauxiliary = Qubit[1]) {
-                        let op = LogicalANDMeasAndFix(_, _);
-                        // Naive measurement-free approach uses 4x more T gates with
-                        // let op = (Controlled X);
-                        op([index![Length(index!) - 1]] + auxiliary, newauxiliary[0]);
-                        if (nUnitariesRight > 0) {
-                            _MultiplexOperationsFromGenerator(rightUnitaries, newauxiliary, newControls, target);
-                        }
+                    let controls = [Tail(index!)] + auxiliary;
+                    using ((newauxiliary, andauxiliary) = (Qubit(), Qubit[MaxI(0, Length(controls) - 2)])) {
                         within {
-                            (Controlled X)(auxiliary, newauxiliary[0]);
+                            ApplyAndChain(andauxiliary, controls, newauxiliary);
                         } apply {
-                            _MultiplexOperationsFromGenerator(leftUnitaries, newauxiliary, newControls, target);
+                            if (nUnitariesRight > 0) {
+                                MultiplexOperationsFromGeneratorImpl(rightUnitaries, [newauxiliary], newControls, target);
+                            }
+                            within {
+                                (Controlled X)(auxiliary, newauxiliary);
+                            } apply {
+                                MultiplexOperationsFromGeneratorImpl(leftUnitaries, [newauxiliary], newControls, target);
+                            }
                         }
-                        (Adjoint op)([index![Length(index!) - 1]] + auxiliary, newauxiliary[0]);
                     }
                 }
             }
         }
         adjoint auto;
         controlled (controlRegister, (...)) {
-            _MultiplexOperationsFromGenerator(unitaryGenerator, auxiliary + controlRegister, index, target);
+            MultiplexOperationsFromGeneratorImpl(unitaryGenerator, auxiliary + controlRegister, index, target);
         }
         adjoint controlled auto;
     }
@@ -191,63 +192,24 @@ namespace Microsoft.Quantum.Canon {
     }
 
     /// # Summary
-    /// Computes the logical AND of multiple qubits.
-    /// # Input
-    /// ## ctrlRegister
-    /// Qubits input to the multiple-input AND gate.
-    /// ## target
-    /// Qubit on which output of AND is written to. This is
-    /// assumed to always start in the $\ket{0}$ state.
-    /// # Remarks
-    /// When `ctrlRegister` has exactly $2$ qubits,
-    /// this is equivalent to the `CCNOT` operation but phase with a phase
-    /// $e^{i\Pi/2}$ on $\ket{001}$ and $-e^{i\Pi/2}$ on $\ket{101}$ and $\ket{011}$.
-    /// The Adjoint is also measure-and-fixup approach that is the inverse
-    /// of the original operation only in special cases (see references),
-    /// but uses fewer T-gates.
-    ///
-    /// # References
-    /// - [ *Craig Gidney*, 1709.06648](https://arxiv.org/abs/1709.06648)
-    internal operation LogicalANDMeasAndFix(ctrlRegister : Qubit[], target : Qubit)
-    : Unit {
-        body (...) {
-            if(Length(ctrlRegister) == 2){
-                let c1 = ctrlRegister[0];
-                let c2 = ctrlRegister[1];
-                H(target);
-                T(target);
-                CNOT(c1,target);
-                CNOT(c2,target);
-                CNOT(target,c1);
-                CNOT(target,c2);
-                (Adjoint T)(c1);
-                (Adjoint T)(c2);
-                T(target);
-                CNOT(target,c2);
-                CNOT(target,c1);
-                H(target);
-                S(target);
-            } else {
-                (Controlled X)(ctrlRegister, target);
-            }
-        }
-        adjoint (...)  {
-            if(Length(ctrlRegister) == 2){
-                let c1 = ctrlRegister[0];
-                let c2 = ctrlRegister[1];
-                H(target);
-                let Meas = M(target);
-                if (Meas == One) {
-                    within {
-                        H(c2);
-                    } apply {
-                        CNOT(c1,c2);
-                    }
-                    X(target);
-                }
-            } else {
-                 (Controlled X)(ctrlRegister, target);
-            }
+    /// Computes a chain of AND gates
+    ///
+    /// # Description
+    /// The auxiliary qubits to compute temporary results must be specified explicitly.
+    /// The length of that register is `Length(ctrlRegister) - 2`, if there are at least
+    /// two controls, otherwise the length is 0.
+    internal operation ApplyAndChain(auxRegister : Qubit[], ctrlRegister : Qubit[], target : Qubit)
+    : Unit is Adj {
+        if (Length(ctrlRegister) == 0) {
+            X(target);
+        } elif (Length(ctrlRegister) == 1) {
+            CNOT(Head(ctrlRegister), target);
+        } else {
+            EqualityFactI(Length(auxRegister), Length(ctrlRegister) - 2, "Unexpected number of auxiliary qubits");
+            let controls1 = ctrlRegister[0..0] + auxRegister;
+            let controls2 = Rest(ctrlRegister);
+            let targets = auxRegister + [target];
+            ApplyToEachA(ApplyAnd, Zip3(controls1, controls2, targets));
         }
     }
 }

Standard/src/Preparation/Deprecated.qs

@@ -0,0 +1,95 @@
+// Copyright (c) Microsoft Corporation.
+ // Licensed under the MIT License.
+
+namespace Microsoft.Quantum.Preparation {
+
+    open Microsoft.Quantum.Arithmetic;
+
+    /// # Summary
+    /// Uses the Quantum ROM technique to represent a given density matrix.
+    ///
+    /// Given a list of $N$ coefficients $\alpha_j$, this returns a unitary $U$ that uses the Quantum-ROM
+    /// technique to prepare
+    /// an approximation  $\tilde\rho\sum_{j=0}^{N-1}p_j\ket{j}\bra{j}$ of the purification of the density matrix
+    /// $\rho=\sum_{j=0}^{N-1}\frac{|alpha_j|}{\sum_k |\alpha_k|}\ket{j}\bra{j}$. In this approximation, the
+    /// error $\epsilon$ is such that $|p_j-\frac{|alpha_j|}{\sum_k |\alpha_k|}|\le \epsilon / N$ and
+    /// $\|\tilde\rho - \rho\| \le \epsilon$. In other words,
+    /// $$
+    /// \begin{align}
+    /// U\ket{0}^{\lceil\log_2 N\rceil}\ket{0}^{m}=\sum_{j=0}^{N-1}\sqrt{p_j} \ket{j}\ket{\text{garbage}_j}.
+    /// \end{align}
+    /// $$
+    ///
+    /// # Input
+    /// ## targetError
+    /// The target error $\epsilon$.
+    /// ## coefficients
+    /// Array of $N$ coefficients specifying the probability of basis states.
+    /// Negative numbers $-\alpha_j$ will be treated as positive $|\alpha_j|$.
+    ///
+    /// # Output
+    /// ## First parameter
+    /// A tuple `(x,(y,z))` where `x = y + z` is the total number of qubits allocated,
+    /// `y` is the number of qubits for the `LittleEndian` register, and `z` is the Number
+    /// of garbage qubits.
+    /// ## Second parameter
+    /// The one-norm $\sum_j |\alpha_j|$ of the coefficient array.
+    /// ## Third parameter
+    /// The unitary $U$.
+    ///
+    /// # Remarks
+    /// ## Example
+    /// The following code snippet prepares an purification of the $3$-qubit state
+    /// $\rho=\sum_{j=0}^{4}\frac{|alpha_j|}{\sum_k |\alpha_k|}\ket{j}\bra{j}$, where
+    /// $\vec\alpha=(1.0,2.0,3.0,4.0,5.0)$, and the error is `1e-3`;
+    /// ```qsharp
+    /// let coefficients = [1.0,2.0,3.0,4.0,5.0];
+    /// let targetError = 1e-3;
+    /// let ((nTotalQubits, (nIndexQubits, nGarbageQubits)), oneNorm, op) = QuantumROM(targetError, coefficients);
+    /// using (indexRegister = Qubit[nIndexQubits]) {
+    ///     using (garbageRegister = Qubit[nGarbageQubits]) {
+    ///         op(LittleEndian(indexRegister), garbageRegister);
+    ///     }
+    /// }
+    /// ```
+    ///
+    /// # References
+    /// - Encoding Electronic Spectra in Quantum Circuits with Linear T Complexity
+    ///   Ryan Babbush, Craig Gidney, Dominic W. Berry, Nathan Wiebe, Jarrod McClean, Alexandru Paler, Austin Fowler, Hartmut Neven
+    ///   https://arxiv.org/abs/1805.03662
+    @Deprecated("Microsoft.Quantum.Preparation.PurifiedMixedState")
+    function QuantumROM(targetError: Double, coefficients: Double[])
+    : ((Int, (Int, Int)), Double, ((LittleEndian, Qubit[]) => Unit is Adj + Ctl)) {
+        let preparation = PurifiedMixedState(targetError, coefficients);
+        return (
+            preparation::Requirements!,
+            preparation::Norm,
+            IgnoreDataRegister(preparation::Prepare, _, _)
+        );
+    }
+
+    internal operation IgnoreDataRegister(op : ((LittleEndian, Qubit[], Qubit[]) => Unit is Adj + Ctl), indexRegister : LittleEndian, garbageRegister : Qubit[]) : Unit is Adj + Ctl {
+        op(indexRegister, new Qubit[0], garbageRegister);
+    }
+
+    /// # Summary
+    /// Returns the total number of qubits that must be allocated
+    /// to the operation returned by `QuantumROM`.
+    ///
+    /// # Input
+    /// ## targetError
+    /// The target error $\epsilon$.
+    /// ## nCoeffs
+    /// Number of coefficients specified in `QuantumROM`.
+    ///
+    /// # Output
+    /// ## First parameter
+    /// A tuple `(x,(y,z))` where `x = y + z` is the total number of qubits allocated,
+    /// `y` is the number of qubits for the `LittleEndian` register, and `z` is the Number
+    /// of garbage qubits.
+    @Deprecated("Microsoft.Quantum.Preparation.PurifiedMixedStateRequirements")
+    function QuantumROMQubitCount(targetError: Double, nCoeffs: Int)
+    : (Int, (Int, Int)) {
+        return (PurifiedMixedStateRequirements(targetError, nCoeffs))!;
+    }
+}