Different ways to run the Microsoft Quantum resource estimator

Prerequisites for VS Code

• The latest version of VS Code or open VS Code for the Web.
• The latest version of the QDK extension in VS Code.

Create a new Q# file

1. In VS Code, open the File menu and choose New Text File to create a new file.
2. Save the file as ShorRE.qs. This file is where your write the Q# code for your program.

Create the quantum algorithm

Copy the following code for Shor's algorithm into the ShorRE.qs file:

    import Std.Arrays.*;
    import Std.Canon.*;
    import Std.Convert.*;
    import Std.Diagnostics.*;
    import Std.Math.*;
    import Std.Measurement.*;
    import Std.Arithmetic.*;
    import Std.ResourceEstimation.*;

    operation Main() : Unit {
        let bitsize = 31;

        // When choosing parameters for `EstimateFrequency`, make sure that
        // generator and modules are not co-prime
        let _ = EstimateFrequency(11, 2^bitsize - 1, bitsize);
    }

    // In this sample we concentrate on costing the `EstimateFrequency`
    // operation, which is the core quantum operation in Shor's algorithm, and
    // we omit the classical pre- and post-processing.

    /// # Summary
    /// Estimates the frequency of a generator
    /// in the residue ring Z mod `modulus`.
    ///
    /// # Input
    /// ## generator
    /// The unsigned integer multiplicative order (period)
    /// of which is being estimated. Must be co-prime to `modulus`.
    /// ## modulus
    /// The modulus which defines the residue ring Z mod `modulus`
    /// in which the multiplicative order of `generator` is being estimated.
    /// ## bitsize
    /// Number of bits needed to represent the modulus.
    ///
    /// # Output
    /// The numerator k of dyadic fraction k/2^bitsPrecision
    /// approximating s/r.
    operation EstimateFrequency(
        generator : Int,
        modulus : Int,
        bitsize : Int
    )
    : Int {
        mutable frequencyEstimate = 0;
        let bitsPrecision =  2 * bitsize + 1;

        // Allocate qubits for the superposition of eigenstates of
        // the oracle that is used in period finding.
        use eigenstateRegister = Qubit[bitsize];

        // Initialize eigenstateRegister to 1, which is a superposition of
        // the eigenstates we are estimating the phases of.
        // We first interpret the register as encoding an unsigned integer
        // in little endian encoding.
        ApplyXorInPlace(1, eigenstateRegister);
        let oracle = ApplyOrderFindingOracle(generator, modulus, _, _);

        // Use phase estimation with a semiclassical Fourier transform to
        // estimate the frequency.
        use c = Qubit();
        for idx in bitsPrecision - 1..-1..0 {
            within {
                H(c);
            } apply {
                // `BeginEstimateCaching` and `EndEstimateCaching` are the operations
                // exposed by Microsoft Quantum resource estimator. These will instruct
                // resource counting such that the if-block will be executed
                // only once, its resources will be cached, and appended in
                // every other iteration.
                if BeginEstimateCaching("ControlledOracle", SingleVariant()) {
                    Controlled oracle([c], (1 <<< idx, eigenstateRegister));
                    EndEstimateCaching();
                }
                R1Frac(frequencyEstimate, bitsPrecision - 1 - idx, c);
            }
            if MResetZ(c) == One {
                frequencyEstimate += 1 <<< (bitsPrecision - 1 - idx);
            }
        }

        // Return all the qubits used for oracles eigenstate back to 0 state
        // using Std.Intrinsic.ResetAll
        ResetAll(eigenstateRegister);

        return frequencyEstimate;
    }

    /// # Summary
    /// Interprets `target` as encoding unsigned little-endian integer k
    /// and performs transformation |k⟩ ↦ |gᵖ⋅k mod N ⟩ where
    /// p is `power`, g is `generator` and N is `modulus`.
    ///
    /// # Input
    /// ## generator
    /// The unsigned integer multiplicative order ( period )
    /// of which is being estimated. Must be co-prime to `modulus`.
    /// ## modulus
    /// The modulus which defines the residue ring Z mod `modulus`
    /// in which the multiplicative order of `generator` is being estimated.
    /// ## power
    /// Power of `generator` by which `target` is multiplied.
    /// ## target
    /// Register interpreted as little endian encoded which is multiplied by
    /// given power of the generator. The multiplication is performed modulo
    /// `modulus`.
    internal operation ApplyOrderFindingOracle(
        generator : Int, modulus : Int, power : Int, target : Qubit[]
    )
    : Unit
    is Adj + Ctl {
        // The oracle we use for order finding implements |x⟩ ↦ |x⋅a mod N⟩. We
        // also use `ExpModI` to compute a by which x must be multiplied. Also
        // note that we interpret target as unsigned integer in little-endian
        // encoding.
        ModularMultiplyByConstant(modulus,
                                    ExpModI(generator, power, modulus),
                                    target);
    }

    /// # Summary
    /// Performs modular in-place multiplication by a classical constant.
    ///
    /// # Description
    /// Given the classical constants `c` and `modulus`, and an input
    /// quantum register  |𝑦⟩, this operation
    /// computes `(c*x) % modulus` into |𝑦⟩.
    ///
    /// # Input
    /// ## modulus
    /// Modulus to use for modular multiplication
    /// ## c
    /// Constant by which to multiply |𝑦⟩
    /// ## y
    /// Quantum register of target
    internal operation ModularMultiplyByConstant(modulus : Int, c : Int, y : Qubit[])
    : Unit is Adj + Ctl {
        use qs = Qubit[Length(y)];
        for (idx, yq) in Enumerated(y) {
            let shiftedC = (c <<< idx) % modulus;
            Controlled ModularAddConstant([yq], (modulus, shiftedC, qs));
        }
        ApplyToEachCA(SWAP, Zipped(y, qs));
        let invC = InverseModI(c, modulus);
        for (idx, yq) in Enumerated(y) {
            let shiftedC = (invC <<< idx) % modulus;
            Controlled ModularAddConstant([yq], (modulus, modulus - shiftedC, qs));
        }
    }

    /// # Summary
    /// Performs modular in-place addition of a classical constant into a
    /// quantum register.
    ///
    /// # Description
    /// Given the classical constants `c` and `modulus`, and an input
    /// quantum register  |𝑦⟩, this operation
    /// computes `(x+c) % modulus` into |𝑦⟩.
    ///
    /// # Input
    /// ## modulus
    /// Modulus to use for modular addition
    /// ## c
    /// Constant to add to |𝑦⟩
    /// ## y
    /// Quantum register of target
    internal operation ModularAddConstant(modulus : Int, c : Int, y : Qubit[])
    : Unit is Adj + Ctl {
        body (...) {
            Controlled ModularAddConstant([], (modulus, c, y));
        }
        controlled (ctrls, ...) {
            // We apply a custom strategy to control this operation instead of
            // letting the compiler create the controlled variant for us in which
            // the `Controlled` functor would be distributed over each operation
            // in the body.
            //
            // Here we can use some scratch memory to save ensure that at most one
            // control qubit is used for costly operations such as `AddConstant`
            // and `CompareGreaterThenOrEqualConstant`.
            if Length(ctrls) >= 2 {
                use control = Qubit();
                within {
                    Controlled X(ctrls, control);
                } apply {
                    Controlled ModularAddConstant([control], (modulus, c, y));
                }
            } else {
                use carry = Qubit();
                Controlled AddConstant(ctrls, (c, y + [carry]));
                Controlled Adjoint AddConstant(ctrls, (modulus, y + [carry]));
                Controlled AddConstant([carry], (modulus, y));
                Controlled CompareGreaterThanOrEqualConstant(ctrls, (c, y, carry));
            }
        }
    }

    /// # Summary
    /// Performs in-place addition of a constant into a quantum register.
    ///
    /// # Description
    /// Given a non-empty quantum register |𝑦⟩ of length 𝑛+1 and a positive
    /// constant 𝑐 < 2ⁿ, computes |𝑦 + c⟩ into |𝑦⟩.
    ///
    /// # Input
    /// ## c
    /// Constant number to add to |𝑦⟩.
    /// ## y
    /// Quantum register of second summand and target; must not be empty.
    internal operation AddConstant(c : Int, y : Qubit[]) : Unit is Adj + Ctl {
        // We are using this version instead of the library version that is based
        // on Fourier angles to show an advantage of sparse simulation in this sample.

        let n = Length(y);
        Fact(n > 0, "Bit width must be at least 1");

        Fact(c >= 0, "constant must not be negative");
        Fact(c < 2 ^ n, $"constant must be smaller than {2L ^ n}");

        if c != 0 {
            // If c has j trailing zeroes than the j least significant bits
            // of y won't be affected by the addition and can therefore be
            // ignored by applying the addition only to the other qubits and
            // shifting c accordingly.
            let j = NTrailingZeroes(c);
            use x = Qubit[n - j];
            within {
                ApplyXorInPlace(c >>> j, x);
            } apply {
                IncByLE(x, y[j...]);
            }
        }
    }

    /// # Summary
    /// Performs greater-than-or-equals comparison to a constant.
    ///
    /// # Description
    /// Toggles output qubit `target` if and only if input register `x`
    /// is greater than or equal to `c`.
    ///
    /// # Input
    /// ## c
    /// Constant value for comparison.
    /// ## x
    /// Quantum register to compare against.
    /// ## target
    /// Target qubit for comparison result.
    ///
    /// # Reference
    /// This construction is described in [Lemma 3, arXiv:2201.10200]
    internal operation CompareGreaterThanOrEqualConstant(c : Int, x : Qubit[], target : Qubit)
    : Unit is Adj+Ctl {
        let bitWidth = Length(x);

        if c == 0 {
            X(target);
        } elif c >= 2 ^ bitWidth {
            // do nothing
        } elif c == 2 ^ (bitWidth - 1) {
            ApplyLowTCNOT(Tail(x), target);
        } else {
            // normalize constant
            let l = NTrailingZeroes(c);

            let cNormalized = c >>> l;
            let xNormalized = x[l...];
            let bitWidthNormalized = Length(xNormalized);
            let gates = Rest(IntAsBoolArray(cNormalized, bitWidthNormalized));

            use qs = Qubit[bitWidthNormalized - 1];
            let cs1 = [Head(xNormalized)] + Most(qs);
            let cs2 = Rest(xNormalized);

            within {
                for i in IndexRange(gates) {
                    (gates[i] ? ApplyAnd | ApplyOr)(cs1[i], cs2[i], qs[i]);
                }
            } apply {
                ApplyLowTCNOT(Tail(qs), target);
            }
        }
    }

    /// # Summary
    /// Internal operation used in the implementation of GreaterThanOrEqualConstant.
    internal operation ApplyOr(control1 : Qubit, control2 : Qubit, target : Qubit) : Unit is Adj {
        within {
            ApplyToEachA(X, [control1, control2]);
        } apply {
            ApplyAnd(control1, control2, target);
            X(target);
        }
    }

    internal operation ApplyAnd(control1 : Qubit, control2 : Qubit, target : Qubit)
    : Unit is Adj {
        body (...) {
            CCNOT(control1, control2, target);
        }
        adjoint (...) {
            H(target);
            if (M(target) == One) {
                X(target);
                CZ(control1, control2);
            }
        }
    }


    /// # Summary
    /// Returns the number of trailing zeroes of a number
    ///
    /// ## Example
    /// ```qsharp
    /// let zeroes = NTrailingZeroes(21); // = NTrailingZeroes(0b1101) = 0
    /// let zeroes = NTrailingZeroes(20); // = NTrailingZeroes(0b1100) = 2
    /// ```
    internal function NTrailingZeroes(number : Int) : Int {
        mutable nZeroes = 0;
        mutable copy = number;
        while (copy % 2 == 0) {
            nZeroes += 1;
            copy /= 2;
        }
        return nZeroes;
    }

    /// # Summary
    /// An implementation for `CNOT` that when controlled using a single control uses
    /// a helper qubit and uses `ApplyAnd` to reduce the T-count to 4 instead of 7.
    internal operation ApplyLowTCNOT(a : Qubit, b : Qubit) : Unit is Adj+Ctl {
        body (...) {
            CNOT(a, b);
        }

        adjoint self;

        controlled (ctls, ...) {
            // In this application this operation is used in a way that
            // it is controlled by at most one qubit.
            Fact(Length(ctls) <= 1, "At most one control line allowed");

            if IsEmpty(ctls) {
                CNOT(a, b);
            } else {
                use q = Qubit();
                within {
                    ApplyAnd(Head(ctls), a, q);
                } apply {
                    CNOT(q, b);
                }
            }
        }

        controlled adjoint self;
    }

Run the resource estimator

The resource estimator offers six predefined qubit parameters. Four of the parameters have gate-based instruction sets, and the other two have a Majorana instruction set. The resource estimator also offers two quantum error correction codes, surface_code and floquet_code.

In this example, you run the resource estimator with the qubit_gate_us_e3 qubit parameter and the surface_code quantum error correction code.

1. In VS Code, open the View menu and choose Command Palette.

2. Enter QDK: Calculate Resource Estimates and select that command. Or, choose the Estimate code lens from the list of commands that precedes the Main operation. Both options open the Qubit types window for the resource estimator.

3. In the Qubit types window, select only qubit_gate_us_e3 and then choose the OK button.

   

4. In the Error budget window, enter 0.001.

5. In the Friendly name for run window, enter the default run name based on the filename ShorRE.qs. The results from the resource estimator open in a new QDK Estimates tab.

View the resource estimator results

The resource estimator provides multiple estimates for the same algorithm. Compare the estimates to understand the tradeoffs between the number of qubits and the runtime. To view the results and compare estimates, follow these steps:

1. Go to the QDK Estimates tab.

2. The Results dropdown displays a summary of the resource estimation. Choose the menu icon in the first row to select the columns that display. For example, select Estimate type, Qubit type, QEC scheme, Error budget, Code distance, Runtime, and Physical qubits.

   

   The number of optimal combinations of {number of qubits, runtime} for your algorithm is in the Estimate type column of the results table. Each combination appears as a point in the space-time diagram. In this case, there are 14 combinations.

   Note

   If you select more than one qubit parameter and error correction code in the configuration, then the results display in different rows for each selection. Choose a result row from the table to bring up the corresponding space-time diagram and report data.

3. The Space-time diagram dropdown shows the tradeoffs between the number of physical qubits and the runtime of the algorithm. In this case, the resource estimator finds 14 optimal combinations out of many thousands of possible combinations. To view an estimation summary for a combination, hover over or select the corresponding point in the diagram.

   

   Note

   Select a point in the space-time diagram to see the space diagram and the details of the resource estimation that correspond to that point.

4. The Space diagram dropdown displays the distribution of physical qubits that your algorithm uses and the T factories. Select different points in the space-time diagram to see how the total number of physical qubits, algorithm qubits, and T factory qubits change for each point.

   

5. The Resource Estimates dropdown displays the full list of output data from the resource estimator. To view more information about each resource category, expand the corresponding dropdown. For example, expand the Logical qubit parameters group for the point shown in the previous steps.

   Logical qubit parameter	Value
   QEC scheme	surface_code
   Code distance	23
   Physical qubits	1.06k
   Logical cycle time	14 millisecs
   Logical qubit error rate	3.00e-14
   Crossing prefactor	0.03
   Error correction threshold	0.01
   Logical cycle time formula	(4 * twoQubitGateTime + 2 * oneQubitMeasurementTime) * codeDistance
   Physical qubits formula	2 * codeDistance * codeDistance

   Tip

   Select the Show detailed rows box to display detailed descriptions for each output value in the report data.

   For more information, see Retrieve the output of the Microsoft Quantum resource estimator for the resource estimator.

Change the target parameters

To see resource estimations for other qubit parameter and error correction code values, choose the Estimate code lens again and then select a new value from the Qubit types window.

For example, select only the Majorana-based qubit parameter qubit_maj_ns_e6 + floquet_code. Accept the default error budget and run name values. The resource estimator reruns the estimation with the new target parameters, and the new output appears in the QDK Estimates tab.

For more information, see Target parameters for the resource estimator.

Run multiple configurations of parameters

The Microsoft Quantum resource estimator can run multiple configurations of target parameters at the same time and compare the resource estimation results for each configuration. To run and compare estimates for multiple parameter configurations, follow these steps:

1. Close the QDK Estimates tab.

2. Choose the Estimates code lens from the list of commands that precedes the Main operation.

3. Select qubit_gate_us_e3, qubit_gate_us_e4, qubit_maj_ns_e4 + floquet_code, and qubit_maj_ns_e6 + floquet_code. Then, choose the OK button.

4. Enter the default value of 0.001 in the Error budget window.

5. Enter the default run name value in the Friendly name for run window.

6. In the case of multiple configurations of parameters, the results are displayed in different rows in the Results tab.

7. The Space-time diagram shows the results for all the configurations of parameters. The first column of the results table displays the legend for each configuration of parameters. You can hover over each point to see the details of the resource estimation at that point.

   

8. To explore the resource estimates for different combinations of input parameters and optimizations, select different points on the space-time diagram.

Prerequisites for Jupyter Notebook in VS Code

• A Python environment with Python and Pip installed.

• The latest version of Visual Studio Code or open VS Code on the Web.

• VS Code with the QDK, Python, and Jupyter extensions installed.

• The latest qdk Python library with the optional jupyter extra.

  pip install --upgrade "qdk[jupyter]"
  

Create the quantum algorithm

1. In VS Code, open the View menu and choose Command Palette.

2. Enter and select Create: New Jupyter Notebook.

   VS Code detects and displays the version of Python and the virtual Python environment that you selected for the notebook. If you have multiple Python environments, then use the kernel picker to select that kernel that you want to use. If VS Code doesn't detect a Python environment, then see Jupyter Notebooks in VS Code for setup information.

3. In the first cell of the notebook, import the qsharp package.

   from qdk import qsharp
   

4. Run the following Q# code in a new cell:

   %%qsharp
   
   import Std.Arrays.*;
   import Std.Canon.*;
   import Std.Convert.*;
   import Std.Diagnostics.*;
   import Std.Math.*;
   import Std.Measurement.*;
   import Std.Arithmetic.*;
   import Std.ResourceEstimation.*; 
   
   operation RunProgram() : Unit {
       let bitsize = 31;
   
       // When choosing parameters for `EstimateFrequency`, make sure that
       // generator and modules are not co-prime
       let _ = EstimateFrequency(11, 2^bitsize - 1, bitsize);
   }
   
   
   // In this sample we concentrate on costing the `EstimateFrequency`
   // operation, which is the core quantum operation in Shors algorithm, and
   // we omit the classical pre- and post-processing.
   
   /// # Summary
   /// Estimates the frequency of a generator
   /// in the residue ring Z mod `modulus`.
   ///
   /// # Input
   /// ## generator
   /// The unsigned integer multiplicative order (period)
   /// of which is being estimated. Must be co-prime to `modulus`.
   /// ## modulus
   /// The modulus which defines the residue ring Z mod `modulus`
   /// in which the multiplicative order of `generator` is being estimated.
   /// ## bitsize
   /// Number of bits needed to represent the modulus.
   ///
   /// # Output
   /// The numerator k of dyadic fraction k/2^bitsPrecision
   /// approximating s/r.
   operation EstimateFrequency(
       generator : Int,
       modulus : Int,
       bitsize : Int
   )
   : Int {
       mutable frequencyEstimate = 0;
       let bitsPrecision =  2 * bitsize + 1;
   
       // Allocate qubits for the superposition of eigenstates of
       // the oracle that is used in period finding.
       use eigenstateRegister = Qubit[bitsize];
   
       // Initialize eigenstateRegister to 1, which is a superposition of
       // the eigenstates we are estimating the phases of.
       // We first interpret the register as encoding an unsigned integer
       // in little endian encoding.
       ApplyXorInPlace(1, eigenstateRegister);
       let oracle = ApplyOrderFindingOracle(generator, modulus, _, _);
   
       // Use phase estimation with a semiclassical Fourier transform to
       // estimate the frequency.
       use c = Qubit();
       for idx in bitsPrecision - 1..-1..0 {
           within {
               H(c);
           } apply {
               // `BeginEstimateCaching` and `EndEstimateCaching` are the operations
               // exposed by Microsoft Quantum resource estimator. These will instruct
               // resource counting such that the if-block will be executed
               // only once, its resources will be cached, and appended in
               // every other iteration.
               if BeginEstimateCaching("ControlledOracle", SingleVariant()) {
                   Controlled oracle([c], (1 <<< idx, eigenstateRegister));
                   EndEstimateCaching();
               }
               R1Frac(frequencyEstimate, bitsPrecision - 1 - idx, c);
           }
           if MResetZ(c) == One {
               frequencyEstimate += 1 <<< (bitsPrecision - 1 - idx);
           }
       }
   
       // Return all the qubits used for oracle eigenstate back to 0 state
       // using Std.Intrinsic.ResetAll
       ResetAll(eigenstateRegister);
   
       return frequencyEstimate;
   }
   
   /// # Summary
   /// Interprets `target` as encoding unsigned little-endian integer k
   /// and performs transformation |k⟩ ↦ |gᵖ⋅k mod N ⟩ where
   /// p is `power`, g is `generator` and N is `modulus`.
   ///
   /// # Input
   /// ## generator
   /// The unsigned integer multiplicative order ( period )
   /// of which is being estimated. Must be co-prime to `modulus`.
   /// ## modulus
   /// The modulus which defines the residue ring Z mod `modulus`
   /// in which the multiplicative order of `generator` is being estimated.
   /// ## power
   /// Power of `generator` by which `target` is multiplied.
   /// ## target
   /// Register interpreted as little endian encoded which is multiplied by
   /// given power of the generator. The multiplication is performed modulo
   /// `modulus`.
   internal operation ApplyOrderFindingOracle(
       generator : Int, modulus : Int, power : Int, target : Qubit[]
   )
   : Unit
   is Adj + Ctl {
       // The oracle we use for order finding implements |x⟩ ↦ |x⋅a mod N⟩. We
       // also use `ExpModI` to compute a by which x must be multiplied. Also
       // note that we interpret target as unsigned integer in little-endian
       // encoding.
       ModularMultiplyByConstant(modulus,
                                   ExpModI(generator, power, modulus),
                                   target);
   }
   
   /// # Summary
   /// Performs modular in-place multiplication by a classical constant.
   ///
   /// # Description
   /// Given the classical constants `c` and `modulus`, and an input
   /// quantum register |𝑦⟩, this operation
   /// computes `(c*x) % modulus` into |𝑦⟩.
   ///
   /// # Input
   /// ## modulus
   /// Modulus to use for modular multiplication
   /// ## c
   /// Constant by which to multiply |𝑦⟩
   /// ## y
   /// Quantum register of target
   internal operation ModularMultiplyByConstant(modulus : Int, c : Int, y : Qubit[])
   : Unit is Adj + Ctl {
       use qs = Qubit[Length(y)];
       for (idx, yq) in Enumerated(y) {
           let shiftedC = (c <<< idx) % modulus;
           Controlled ModularAddConstant([yq], (modulus, shiftedC, qs));
       }
       ApplyToEachCA(SWAP, Zipped(y, qs));
       let invC = InverseModI(c, modulus);
       for (idx, yq) in Enumerated(y) {
           let shiftedC = (invC <<< idx) % modulus;
           Controlled ModularAddConstant([yq], (modulus, modulus - shiftedC, qs));
       }
   }
   
   /// # Summary
   /// Performs modular in-place addition of a classical constant into a
   /// quantum register.
   ///
   /// # Description
   /// Given the classical constants `c` and `modulus`, and an input
   /// quantum register  |𝑦⟩, this operation
   /// computes `(x+c) % modulus` into |𝑦⟩.
   ///
   /// # Input
   /// ## modulus
   /// Modulus to use for modular addition
   /// ## c
   /// Constant to add to |𝑦⟩
   /// ## y
   /// Quantum register of target
   internal operation ModularAddConstant(modulus : Int, c : Int, y : Qubit[])
   : Unit is Adj + Ctl {
       body (...) {
           Controlled ModularAddConstant([], (modulus, c, y));
       }
       controlled (ctrls, ...) {
           // We apply a custom strategy to control this operation instead of
           // letting the compiler create the controlled variant for us in which
           // the `Controlled` functor would be distributed over each operation
           // in the body.
           //
           // Here we can use some scratch memory to save ensure that at most one
           // control qubit is used for costly operations such as `AddConstant`
           // and `CompareGreaterThenOrEqualConstant`.
           if Length(ctrls) >= 2 {
               use control = Qubit();
               within {
                   Controlled X(ctrls, control);
               } apply {
                   Controlled ModularAddConstant([control], (modulus, c, y));
               }
           } else {
               use carry = Qubit();
               Controlled AddConstant(ctrls, (c, y + [carry]));
               Controlled Adjoint AddConstant(ctrls, (modulus, y + [carry]));
               Controlled AddConstant([carry], (modulus, y));
               Controlled CompareGreaterThanOrEqualConstant(ctrls, (c, y, carry));
           }
       }
   }
   
   /// # Summary
   /// Performs in-place addition of a constant into a quantum register.
   ///
   /// # Description
   /// Given a non-empty quantum register |𝑦⟩ of length 𝑛+1 and a positive
   /// constant 𝑐 < 2ⁿ, computes |𝑦 + c⟩ into |𝑦⟩.
   ///
   /// # Input
   /// ## c
   /// Constant number to add to |𝑦⟩.
   /// ## y
   /// Quantum register of second summand and target; must not be empty.
   internal operation AddConstant(c : Int, y : Qubit[]) : Unit is Adj + Ctl {
       // We are using this version instead of the library version that is based
       // on Fourier angles to show an advantage of sparse simulation in this sample.
   
       let n = Length(y);
       Fact(n > 0, "Bit width must be at least 1");
   
       Fact(c >= 0, "constant must not be negative");
       Fact(c < 2 ^ n, $"constant must be smaller than {2L ^ n}");
   
       if c != 0 {
           // If c has j trailing zeroes than the j least significant bits
           // of y will not be affected by the addition and can therefore be
           // ignored by applying the addition only to the other qubits and
           // shifting c accordingly.
           let j = NTrailingZeroes(c);
           use x = Qubit[n - j];
           within {
               ApplyXorInPlace(c >>> j, x);
           } apply {
               IncByLE(x, y[j...]);
           }
       }
   }
   
   /// # Summary
   /// Performs greater-than-or-equals comparison to a constant.
   ///
   /// # Description
   /// Toggles output qubit `target` if and only if input register `x`
   /// is greater than or equal to `c`.
   ///
   /// # Input
   /// ## c
   /// Constant value for comparison.
   /// ## x
   /// Quantum register to compare against.
   /// ## target
   /// Target qubit for comparison result.
   ///
   /// # Reference
   /// This construction is described in [Lemma 3, arXiv:2201.10200]
   internal operation CompareGreaterThanOrEqualConstant(c : Int, x : Qubit[], target : Qubit)
   : Unit is Adj+Ctl {
       let bitWidth = Length(x);
   
       if c == 0 {
           X(target);
       } elif c >= 2 ^ bitWidth {
           // do nothing
       } elif c == 2 ^ (bitWidth - 1) {
           ApplyLowTCNOT(Tail(x), target);
       } else {
           // normalize constant
           let l = NTrailingZeroes(c);
   
           let cNormalized = c >>> l;
           let xNormalized = x[l...];
           let bitWidthNormalized = Length(xNormalized);
           let gates = Rest(IntAsBoolArray(cNormalized, bitWidthNormalized));
   
           use qs = Qubit[bitWidthNormalized - 1];
           let cs1 = [Head(xNormalized)] + Most(qs);
           let cs2 = Rest(xNormalized);
   
           within {
               for i in IndexRange(gates) {
                   (gates[i] ? ApplyAnd | ApplyOr)(cs1[i], cs2[i], qs[i]);
               }
           } apply {
               ApplyLowTCNOT(Tail(qs), target);
           }
       }
   }
   
   /// # Summary
   /// Internal operation used in the implementation of GreaterThanOrEqualConstant.
   internal operation ApplyOr(control1 : Qubit, control2 : Qubit, target : Qubit) : Unit is Adj {
       within {
           ApplyToEachA(X, [control1, control2]);
       } apply {
           ApplyAnd(control1, control2, target);
           X(target);
       }
   }
   
   internal operation ApplyAnd(control1 : Qubit, control2 : Qubit, target : Qubit)
   : Unit is Adj {
       body (...) {
           CCNOT(control1, control2, target);
       }
       adjoint (...) {
           H(target);
           if (M(target) == One) {
               X(target);
               CZ(control1, control2);
           }
       }
   }
   
   
   /// # Summary
   /// Returns the number of trailing zeroes of a number
   ///
   /// ## Example
   /// ```qsharp
   /// let zeroes = NTrailingZeroes(21); // = NTrailingZeroes(0b1101) = 0
   /// let zeroes = NTrailingZeroes(20); // = NTrailingZeroes(0b1100) = 2
   /// ```
   internal function NTrailingZeroes(number : Int) : Int {
       mutable nZeroes = 0;
       mutable copy = number;
       while (copy % 2 == 0) {
           nZeroes += 1;
           copy /= 2;
       }
       return nZeroes;
   }
   
   /// # Summary
   /// An implementation for `CNOT` that when controlled using a single control uses
   /// a helper qubit and uses `ApplyAnd` to reduce the T-count to 4 instead of 7.
   internal operation ApplyLowTCNOT(a : Qubit, b : Qubit) : Unit is Adj+Ctl {
       body (...) {
           CNOT(a, b);
       }
   
       adjoint self;
   
       controlled (ctls, ...) {
           // In this application this operation is used in a way that
           // it is controlled by at most one qubit.
           Fact(Length(ctls) <= 1, "At most one control line allowed");
   
           if IsEmpty(ctls) {
               CNOT(a, b);
           } else {
               use q = Qubit();
               within {
                   ApplyAnd(Head(ctls), a, q);
               } apply {
                   CNOT(q, b);
               }
           }
       }
   
       controlled adjoint self;
   }
   

Estimate the required resources to run your quantum algorithm

Estimate the physical resources required for the RunProgram operation with the default resource estimator input parameter values. Run the following code in a new cell:

result = qsharp.estimate("RunProgram()")
result

The qsharp.estimate function creates a result object that displays a table with the estimated physical resource requirements. To inspect the results, expand the dropdown groups in the output. For more information, see Retrieve the output of the Microsoft Quantum resource estimator.

For example, expand the Logical qubit parameters group to see the following estimation results:

Space diagram

Your algorithm design might depend on how the physical qubits are distributed between algorithm qubits and T factory qubits. To visualize this distribution and better understand the estimated space requirements for your algorithm, use tools from the qdk.widgets module. Run the following code in a new cell:

from qdk.widgets import SpaceChart, EstimateDetails

SpaceChart(result)

The space diagram shows the proportion of algorithm qubits and T factory qubits. The number of T factory copies, in this case 19, contributes to the number of physical qubits that are used for T factories according to the following equation:

$$\text{T factories} \cdot \text{physical qubit per T factory}= 19 \cdot 33,320 = 633,080$$$

For more information, see T factory physical estimation.

Change the input parameters and compare resource estimates

When you submit a resource estimate request for your program, you can specify some optional parameters. To see all the parameters that you can pass to the job run, along with the default values that are used for each parameter, use the jobParams field.

Run the following code in a new cell:

result['jobParams']

{'qecScheme': {'name': 'surface_code',
  'errorCorrectionThreshold': 0.01,
  'crossingPrefactor': 0.03,
  'distanceCoefficientPower': 0,
  'logicalCycleTime': '(4 * twoQubitGateTime + 2 * oneQubitMeasurementTime) * codeDistance',
  'physicalQubitsPerLogicalQubit': '2 * codeDistance * codeDistance',
  'maxCodeDistance': 50},
 'errorBudget': 0.001,
 'qubitParams': {'instructionSet': 'GateBased',
  'name': 'qubit_gate_ns_e3',
  'oneQubitMeasurementTime': '100 ns',
  'oneQubitGateTime': '50 ns',
  'twoQubitGateTime': '50 ns',
  'tGateTime': '50 ns',
  'oneQubitMeasurementErrorRate': 0.001,
  'oneQubitGateErrorRate': 0.001,
  'twoQubitGateErrorRate': 0.001,
  'tGateErrorRate': 0.001,
  'idleErrorRate': 0.001},
 'constraints': {'maxDistillationRounds': 3},
 'estimateType': 'singlePoint'}

The output is a dictionary of parameter names and values. For example, the resource estimator takes the qubit_gate_ns_e3 qubit model, the surface_code error correction scheme, and 0.001 error budget as default values for the estimation.

You can pass custom values for the following target parameters:

For more information on resource estimator parameters, see Target parameters and the qdk.qsharp.estimator API reference.

Change the qubit model

Estimate the cost for the same algorithm, but use the Majorana-based qubit parameter qubit_maj_ns_e6 instead. Run the following code in a new cell:

inputParams = {"qubitParams": {"name": "qubit_maj_ns_e6"}}

result_maj = qsharp.estimate("RunProgram()", params=inputParams)

EstimateDetails(result_maj)

Change quantum error correction scheme

Rerun the resource estimation job for the same example on the Majorana-based qubit parameter, but with a Floquet QEC scheme instead. Run the following code in a new cell:

inputParams = {"qubitParams": {"name": "qubit_maj_ns_e6"},
               "qecScheme": {"name": "floquet_code"}
               }

result_maj = qsharp.estimate("RunProgram()", params=inputParams)

EstimateDetails(result_maj)

Change the error budget

Rerun the same quantum circuit with an errorBudget of 10%.

inputParams = {"qubitParams": {"name": "qubit_maj_ns_e6"},
               "qecScheme": {"name": "floquet_code"},
                "errorBudget": 0.1
               }

result_maj = qsharp.estimate("RunProgram()", params=inputParams)

EstimateDetails(result_maj)

Run batch resource estimates

The Microsoft Quantum resource estimator allows you to run multiple configurations of target parameters and compare the results for each configuration. To run a batch resource estimate, follow these steps:

1. Pass a list of target parameters to the params parameter of the qsharp.estimate function. For example, run your algorithm with the default parameters and with the Majorana-based qubit parameter with a Floquet QEC scheme. Run the following code in a new cell:

   inputParams = {"qubitParams": {"name": "qubit_maj_ns_e6"},
              "qecScheme": {"name": "floquet_code"}
              }
   
   # An empty dictionary sets all parameters to their default values
   result_batch = qsharp.estimate("RunProgram()", params=[{}, inputParams])
   
   result_batch.summary_data_frame(labels=["Gate-based ns, 10⁻³", "Majorana ns, 10⁻⁶"])
   

   Model	Logical qubits	Logical depth	T states	Code distance	T factories	T factory fraction	Physical qubits	rQOPS	Physical runtime
   Gate-based ns, 10⁻³	223	3.64M	4.70M	21	19	76.30 %	829.77k	26.55M	31 secs
   Majorana ns, 10⁻⁶	223	3.64M	4.70M	5	19	63.02 %	79.60k	148.67M	5 secs

2. You can also construct a list of estimation parameters by using the EstimatorParams class. Run the following code in a new cell:

   from qsharp.estimator import EstimatorParams, QubitParams, QECScheme, LogicalCounts
   
   labels = ["Gate-based µs, 10⁻³", "Gate-based µs, 10⁻⁴", "Gate-based ns, 10⁻³", "Gate-based ns, 10⁻⁴", "Majorana ns, 10⁻⁴", "Majorana ns, 10⁻⁶"]
   
   params = EstimatorParams(num_items=6)
   params.error_budget = 0.333
   params.items[0].qubit_params.name = QubitParams.GATE_US_E3
   params.items[1].qubit_params.name = QubitParams.GATE_US_E4
   params.items[2].qubit_params.name = QubitParams.GATE_NS_E3
   params.items[3].qubit_params.name = QubitParams.GATE_NS_E4
   params.items[4].qubit_params.name = QubitParams.MAJ_NS_E4
   params.items[4].qec_scheme.name = QECScheme.FLOQUET_CODE
   params.items[5].qubit_params.name = QubitParams.MAJ_NS_E6
   params.items[5].qec_scheme.name = QECScheme.FLOQUET_CODE
   
   qsharp.estimate("RunProgram()", params=params).summary_data_frame(labels=labels)
   

   Model	Logical qubits	Logical depth	T states	Code distance	T factories	T factory fraction	Physical qubits	rQOPS	Physical runtime
   Gate-based µs, 10⁻³	223	3.64M	4.70M	17	13	40.54 %	216.77k	21.86k	10 hours
   Gate-based µs, 10⁻⁴	223	3.64M	4.70M	9	14	43.17 %	63.57k	41.30k	5 hours
   Gate-based ns, 10⁻³	223	3.64M	4.70M	17	16	69.08 %	416.89k	32.79M	25 secs
   Gate-based ns, 10⁻⁴	223	3.64M	4.70M	9	14	43.17 %	63.57k	61.94M	13 secs
   Majorana ns, 10⁻⁴	223	3.64M	4.70M	9	19	82.75 %	501.48k	82.59M	10 secs
   Majorana ns, 10⁻⁶	223	3.64M	4.70M	5	13	31.47 %	42.96k	148.67M	5 secs

Run a Pareto frontier estimation

When you estimate the resources required to run an algorithm, it's important to consider the tradeoff between the number of physical qubits and the runtime of the algorithm. If you use more qubits, then you can probably reduce the runtime of your algorithm. However, the number of physical qubits that you can use is limited by the the hardware design.

The Pareto frontier estimation provides multiple estimates for the same algorithm, each with a tradeoff between the number of qubits and the runtime. To perform a Pareto frontier estimation, follow these steps:

1. Pass the value "frontier" to the "estimateType" target parameter when you call the qsharp.estimate function. For example, perform a Pareto frontier estimate for an algorithm with the Majorana-based qubit parameter with a surface code.

   inputParams = {"qubitParams": {"name": "qubit_maj_ns_e4"},
              "qecScheme": {"name": "surface_code"},
               "estimateType": "frontier"
              }
   
   result = qsharp.estimate("RunProgram()", params=inputParams)
   

2. To display a table with the estimation results, use the EstimatesOverview function from the qdk.widgets module. Run the following code in a new cell:

   from qdk.widgets import EstimatesOverview
   
   EstimatesOverview(result)
   

3. Choose the menu icon to select the columns that you want to display.

   

In the Estimate type column of the results table, you can see the number of different combinations of {number of qubits, runtime} for your algorithm. In this case, the resource estimator finds 22 different optimal combinations out of many thousands of possible combinations.

Space-time diagram

The EstimatesOverview function also displays the space-time diagram of the resource estimator.

The space-time diagram shows the number of physical qubits and the runtime of the algorithm for each optimal pair of {number of qubits, runtime}. Hover over a point to see a summary of the resource estimation results for that point.

Run a batch Pareto frontier estimation

1. To estimate and compare multiple configurations of target parameters with frontier estimation, include "frontier" for the "estimateType" parameter in each parameter set. For example, run the following code in a new cell to perform a batch frontier estimation for two different parameter sets:

   inputParams1 = {"qubitParams": {"name": "qubit_maj_ns_e4"},
                   "qecScheme": {"name": "surface_code"},
                   "estimateType": "frontier" # Pareto frontier estimation
               }   
   inputParams2 = {"qubitParams": {"name": "qubit_maj_ns_e6"},
                   "qecScheme": {"name": "floquet_code"},
                   "estimateType": "frontier" # Pareto frontier estimation
               }
   
   result = qsharp.estimate("RunProgram()", params=[inputParams1, inputParams2])
   
   EstimatesOverview(result, colors=["#1f77b4", "#ff7f0e"], runNames=["e4 Surface Code", "e6 Floquet Code"])
   

   

   Note

   To define colors and run names for the qubit-time diagram, use the EstimatesOverview function.

2. When you run multiple configurations of target parameters for Pareto frontier estimations, you can programmatically access the resource estimates for a specific point in the space-time diagram. For example, run the following code in a new cell to show the estimation results for the second (estimate index = 0) run and the fourth (point index = 3) shortest runtime.

   EstimateDetails(result[1], 4)
   

3. You can also see the space diagram for a specific point of the space-time diagram. For example, run the following code in a new cell to show the space diagram for the first run of combinations (estimate index = 0) and the third shortest runtime (point index = 2).

   SpaceChart(result[0], 2)
   

Prerequisites for Qiskit in VS Code

• A Python environment with Python and Pip installed.

• The latest version of VS Code or open VS Code for the Web.

• Install the QDK, Python, and Jupyter extensions in VS Code.

• The latest qdk Python library with the jupyter and qiskit extras.

  pip install --upgrade "qdk[jupyter,qiskit]"
  

Create a new Jupyter Notebook

1. In VS Code, open the View menu and choose Command Palette.
2. Enter and select Create: New Jupyter Notebook.

VS Code detects and displays the version of Python and the virtual Python environment that you selected for the notebook. If you have multiple Python environments, then use the kernel picker to select that kernel that you want to use. If VS Code doesn't detect a Python environment, then see Jupyter Notebooks in VS Code for setup information.

Create the quantum algorithm

In this example, you create a quantum circuit for a multiplier that's based on the construction presented in Ruiz-Perez and Garcia-Escartin (arXiv:1411.5949), which uses the Quantum Fourier Transform to implement arithmetic.

To adjust the size of the multiplier, change the bitwidth variable. The circuit generation is wrapped in a function that you can call with the bitwidth value of the multiplier. The operation has two input registers, each the size of bitwidth, and one output register that's twice the size of bitwidth. The function also prints some logical resource counts for the multiplier extracted directly from the quantum circuit.

Run the following code in the first cell of your notebook:

from qiskit.circuit.library import RGQFTMultiplier 

def create_algorithm(bitwidth):
    print(f"[INFO] Create a QFT-based multiplier with bitwidth {bitwidth}")

    circ = RGQFTMultiplier(num_state_qubits=bitwidth)

    return circ

Estimate the quantum algorithm

To create an instance of your algorithm, use the create_algorithm function. To adjust the size of the multiplier, change the value of the bitwidth variable.

For example, run the following code in a new cell to create an instance of your algorithm with a bitwidth of 4:

bitwidth = 4

circ = create_algorithm(bitwidth)

Use the default resource estimator input parameters to estimate the physical resources needed to run this operation. To get the output from the resource estimator, call the estimate method, which is overloaded to accept a QuantumCircuit object from Qiskit.

Run the following code in a new cell:

from qdk import qsharp
from qsharp.estimator import EstimatorParams
from qsharp.interop.qiskit import estimate

params = EstimatorParams()
result = estimate(circ, params)

The result object contains the output of the resource estimation job. Use the EstimateDetails function to display the results in a readable format.

Run the following code in a new cell:

from qdk.widgets import EstimateDetails

EstimateDetails(result)

The EstimateDetails function displays a list of dropdowns with the estimates for different physical resource categories. To get more information on each resource, expand the dropdown groups in the output. For more information, see Retrieve the output of the Microsoft Quantum resource estimator.

For example, expand the Logical qubit parameter group to see a table of values for that category, such as 15 for the error correction code distance.

Space diagram

Your algorithm design might depend on how the physical qubits are distributed between algorithm qubits and T factory qubits. To visualize this distribution and better understand the estimated space requirements for your algorithm, run the following code in a new cell:

from qdk.widgets import SpaceChart

SpaceChart(result)

The space diagram shows the proportion of algorithm qubits and T factory qubits. The number of T factory copies, in this case 19, contributes to the number of physical qubits that are used for T factories according to the following equation:

$$\text{T factories} \cdot \text{physical qubit per T factory}= 19 \cdot 18,000 = 342,000$$$

For more information, see T factory physical estimation.

Change the input parameters and compare resource estimates

When you submit a resource estimate request for your program, you can specify some optional parameters. To see all the parameters that you can pass to the job run, along with the default values that are used for each parameter, use the jobParams field.

Run the following code in a new cell:

result.data()["jobParams"]

{'qecScheme': {'name': 'surface_code',
  'errorCorrectionThreshold': 0.01,
  'crossingPrefactor': 0.03,
  'distanceCoefficientPower': 0,
  'logicalCycleTime': '(4 * twoQubitGateTime + 2 * oneQubitMeasurementTime) * codeDistance',
  'physicalQubitsPerLogicalQubit': '2 * codeDistance * codeDistance',
  'maxCodeDistance': 50},
 'errorBudget': 0.001,
 'qubitParams': {'instructionSet': 'GateBased',
  'name': 'qubit_gate_ns_e3',
  'oneQubitMeasurementTime': '100 ns',
  'oneQubitGateTime': '50 ns',
  'twoQubitGateTime': '50 ns',
  'tGateTime': '50 ns',
  'oneQubitMeasurementErrorRate': 0.001,
  'oneQubitGateErrorRate': 0.001,
  'twoQubitGateErrorRate': 0.001,
  'tGateErrorRate': 0.001,
  'idleErrorRate': 0.001},
 'constraints': {'maxDistillationRounds': 3},
 'estimateType': 'singlePoint'}

You can pass custom values for the following target parameters:

For more information on resource estimator parameters, see Target parameters and the qdk.qsharp.estimator API reference.

Change the qubit model

Estimate the cost for the same algorithm, but use the Majorana-based qubit parameter qubit_maj_ns_e6 instead. Run the following code in a new cell:

qubitParams = {"name": "qubit_maj_ns_e6"}

result = estimate(circ, qubitParams)

Call the data method to programmatically inspect resource estimates. For example, run the following code in a new cell to explore details about the T factory that's created to run the program.

result.data()["tfactory"]

{'physicalQubits': 18000,
 'runtime': 78000,
 'numTstates': 1,
 'numInputTstates': 30,
 'numRounds': 1,
 'numUnitsPerRound': [2],
 'unitNamePerRound': ['15-to-1 space efficient'],
 'codeDistancePerRound': [15],
 'physicalQubitsPerRound': [18000],
 'runtimePerRound': [78000],
 'logicalErrorRate': 3.713e-08}

You can use this data to explain how the T factories produce the required T states. For example, run the following code in a new cell:

data = result.data()
tfactory = data["tfactory"]
breakdown = data["physicalCounts"]["breakdown"]
producedTstates = breakdown["numTfactories"] * breakdown["numTfactoryRuns"] * tfactory["numTstates"]

print(f"""A single T factory produces {tfactory["logicalErrorRate"]:.2e} T states with an error rate of (required T state error rate is {breakdown["requiredLogicalTstateErrorRate"]:.2e}).""")
print(f"""{breakdown["numTfactories"]} copie(s) of a T factory are executed {breakdown["numTfactoryRuns"]} time(s) to produce {producedTstates} T states ({breakdown["numTstates"]} are required by the algorithm).""")
print(f"""A single T factory is composed of {tfactory["numRounds"]} rounds of distillation:""")
for round in range(tfactory["numRounds"]):
    print(f"""- {tfactory["numUnitsPerRound"][round]} {tfactory["unitNamePerRound"][round]} unit(s)""")

A single T factory produces 3.71e-08 T states with an error rate of (required T state error rate is 3.99e-08).
19 copie(s) of a T factory are executed 441 time(s) to produce 8379 T states (8362 are required by the algorithm).
A single T factory is composed of 1 rounds of distillation:
- 2 15-to-1 space efficient unit(s)

Change the quantum error correction scheme

Rerun the resource estimation job for the same example on the Majorana-based qubit parameter, but with a floquet QEC scheme instead. Run the following code in a new cell:

params = {
    "qubitParams": {"name": "qubit_maj_ns_e6"},
    "qecScheme": {"name": "floquet_code"}
}

result_maj_floquet = estimate(circ, params)
EstimateDetails(result_maj_floquet)

Change the error budget

Rerun the same quantum circuit with an errorBudget of 10%.

params = {
    "errorBudget": 0.01,
    "qubitParams": {"name": "qubit_maj_ns_e6"},
    "qecScheme": {"name": "floquet_code"},
}
result_maj_floquet_e1 = backend.run(circ, params).result()
EstimateDetails(result_maj_floquet_e1)

Compare the output to see how the different combinations of qubit model, error correction scheme, and error budget parameters affect the estimation of the resources needed to run your algorithm on a real quantum computer.