From bf6c787306c75f8e5088c9a5aeeb9ca7f52cde0c Mon Sep 17 00:00:00 2001 From: "Stefan J. Wernli" Date: Thu, 11 Jan 2024 14:17:35 -0800 Subject: [PATCH] Add more samples (#1010) This adds some additional samples folks can use for estimation. --- samples/algorithms/Shor.qs | 251 ++++------------------------ samples/estimation/Dynamics.qs | 122 ++++++++++++++ samples/estimation/Precalculated.qs | 25 +++ samples/estimation/ShorRE.qs | 210 +++++++++++++++++++++++ 4 files changed, 385 insertions(+), 223 deletions(-) create mode 100644 samples/estimation/Dynamics.qs create mode 100644 samples/estimation/Precalculated.qs create mode 100644 samples/estimation/ShorRE.qs diff --git a/samples/algorithms/Shor.qs b/samples/algorithms/Shor.qs index a85f83211a..787f168b7b 100644 --- a/samples/algorithms/Shor.qs +++ b/samples/algorithms/Shor.qs @@ -7,6 +7,7 @@ /// /// This Q# program implements Shor's algorithm. namespace Sample { + open Microsoft.Quantum.Convert; open Microsoft.Quantum.Diagnostics; open Microsoft.Quantum.Random; open Microsoft.Quantum.Math; @@ -139,7 +140,7 @@ namespace Sample { Message($"Found factor={factor}"); return (true, (factor, modulus / factor)); } - } + } // Return a flag indicating we hit a trivial case and didn't get // any factors. Message($"Found trivial factors."); @@ -297,64 +298,38 @@ namespace Sample { } /// # 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`. + /// 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`. + /// 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. + /// 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 which is multiplied by given power - /// of the generator. The multiplication is performed modulo `modulus`. - operation ApplyOrderFindingOracle( - generator : Int, - modulus : Int, - power : Int, - target : Qubit[]) - : Unit is Adj + Ctl { - // Check that the parameters satisfy the requirements. - Fact( - GreatestCommonDivisorI(generator, modulus) == 1, - "`generator` and `modulus` must be co-prime"); - - // 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 fromat. + /// Register interpreted as little-endian 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 + // format. ModularMultiplyByConstant( modulus, ExpModI(generator, power, modulus), target); } - // - // Arithmetic helper functions to implement order-finding oracle. - // - - /// # Summary - /// Returns the number of trailing zeroes of a number. - /// - /// ## Example - /// let zeroes = NTrailingZeroes(21); // NTrailingZeroes(0b1101) = 0 - /// let zeroes = NTrailingZeroes(20); // NTrailingZeroes(0b1100) = 2 - function NTrailingZeroes(number : Int) : Int { - Fact(number != 0, "NTrailingZeroes: number cannot be 0."); - mutable nZeroes = 0; - mutable copy = number; - while (copy % 2 == 0) { - set nZeroes += 1; - set copy /= 2; - } - return nZeroes; - } - /// # Summary /// Performs modular in-place multiplication by a classical constant. /// @@ -370,7 +345,7 @@ namespace Sample { /// Constant by which to multiply |𝑦⟩ /// ## y /// Quantum register of target - operation ModularMultiplyByConstant(modulus : Int, c : Int, y : Qubit[]) + internal operation ModularMultiplyByConstant(modulus : Int, c : Int, y : Qubit[]) : Unit is Adj + Ctl { use qs = Qubit[Length(y)]; for idx in IndexRange(y) { @@ -395,7 +370,6 @@ namespace Sample { /// 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 |𝑦⟩ in little-endian format, this operation computes /// `(x+c) % modulus` into |𝑦⟩. @@ -407,7 +381,7 @@ namespace Sample { /// Constant to add to |𝑦⟩ /// ## y /// Quantum register of target - operation ModularAddConstant(modulus : Int, c : Int, y : Qubit[]) + internal operation ModularAddConstant(modulus : Int, c : Int, y : Qubit[]) : Unit is Adj + Ctl { body (...) { Controlled ModularAddConstant([], (modulus, c, y)); @@ -426,184 +400,15 @@ namespace Sample { within { Controlled X(ctrls, control); } apply { - Controlled ModularAddConstant( - [control], - (modulus, c, y)); + 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. - 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 {2^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] - 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)) { - CNOT(Tail(x), target); - } else { - // normalize constant - let l = NTrailingZeroes(c); - let cNormalized = c >>> l; - let xNormalized = x[l...]; - let bitWidthNormalized = Length(xNormalized); - use qs = Qubit[bitWidthNormalized - 1]; - let cs1 = [Head(xNormalized)] + Most(qs); - Fact(Length(cs1) == Length(qs), - "Arrays should be of the same length."); - within { - for i in 0..Length(cs1)-1 { - let op = - cNormalized &&& (1 <<< (i+1)) != 0 ? - ApplyAnd | ApplyOr; - op(cs1[i], xNormalized[i+1], qs[i]); - } - } apply { - CNOT(Tail(qs), target); - } - } - } - - /// # Summary - /// Inverts a given target qubit if and only if both control qubits are in - /// the 1 state, using measurement to perform the adjoint operation. - /// - /// # Description - /// Inverts `target` if and only if both controls are 1, but assumes that - /// `target` is in state 0. The operation has T-count 4, T-depth 2 and - /// requires no helper qubit, and may therefore be preferable to a CCNOT - /// operation, if `target` is known to be 0. - /// The adjoint of this operation is measurement based and requires no T - /// gates. - /// Although the Toffoli gate (CCNOT) will perform faster in in simulations, - /// this version has lower T gate requirements. - /// - /// # Input - /// ## control1 - /// First control qubit - /// ## control2 - /// Second control qubit - /// ## target - /// Target auxiliary qubit; must be in state 0 - /// - /// # References - /// - Cody Jones: "Novel constructions for the fault-tolerant - /// Toffoli gate", - /// Phys. Rev. A 87, 022328, 2013 - /// [arXiv:1212.5069](https://arxiv.org/abs/1212.5069) - /// doi:10.1103/PhysRevA.87.022328 - /// - Craig Gidney: "Halving the cost of quantum addition", - /// Quantum 2, page 74, 2018 - /// [arXiv:1709.06648](https://arxiv.org/abs/1709.06648) - /// doi:10.1103/PhysRevA.85.044302 - /// - Mathias Soeken: "Quantum Oracle Circuits and the Christmas - /// Tree Pattern", - /// [Blog article from December 19, 2019](https://msoeken.github.io/blog_qac.html) - /// (note: explains the multiple controlled construction) - operation ApplyAnd(control1 : Qubit, control2 : Qubit, target : Qubit) - : Unit is Adj { - body (...) { - H(target); - T(target); - CNOT(control1, target); - CNOT(control2, target); - within { - CNOT(target, control1); - CNOT(target, control2); - } - apply { - Adjoint T(control1); - Adjoint T(control2); - T(target); - } - H(target); - S(target); - } - adjoint (...) { - H(target); - if (M(target) == One) { - X(target); - CZ(control1, control2); + Controlled IncByI(ctrls, (c, y + [carry])); + Controlled Adjoint IncByI(ctrls, (modulus, y + [carry])); + Controlled IncByI([carry], (modulus, y)); + Controlled ApplyIfLessOrEqualL(ctrls, (X, IntAsBigInt(c), y, carry)); } } } - - /// # Summary - /// Applies X to the target if any of the controls are 1. - operation ApplyOr(control1 : Qubit, control2 : Qubit, target : Qubit) - : Unit is Adj { - within { - ApplyToEachA(X, [control1, control2]); - } apply { - ApplyAnd(control1, control2, target); - X(target); - } - } } diff --git a/samples/estimation/Dynamics.qs b/samples/estimation/Dynamics.qs new file mode 100644 index 0000000000..fd523bd82d --- /dev/null +++ b/samples/estimation/Dynamics.qs @@ -0,0 +1,122 @@ +/// # Sample +/// Quantum Dynamics +/// +/// # Description +/// This example demonstrates quantum dynamics in a style tailed for +/// resource estimation. The sample is specifically the simulation +/// of an Ising model Hamiltonian on an NxN 2D lattice using a +/// fourth-order Trotter Suzuki product formula, assuming +/// a 2D qubit architecture with nearest-neighbor connectivity. +/// The is an example of a program that is not amenable to simulating +/// classically, but can be run through resource estimation to determine +/// what size of quantum system would be needed to solve the problem. +namespace QuantumDynamics { + open Microsoft.Quantum.Math; + + @EntryPoint() + operation Main() : Unit { + let N = 10; + let J = 1.0; + let g = 1.0; + let totTime = 20.0; + let dt = 0.25; + let eps = 0.010; + + IsingModel2DSim(N, J, g, totTime, dt, eps); + } + + function GetQubitIndex(row : Int, col : Int, n : Int) : Int { + return row % 2 == 0 // if row is even, + ? col + n * row // move from left to right, + | (n - 1 - col) + n * row; // otherwise from right to left. + } + + function SetSequences(len : Int, p : Double, dt : Double, J : Double, g : Double) : (Double[], Double[]) { + // create two arrays of size `len` + mutable seqA = [0.0, size=len]; + mutable seqB = [0.0, size=len]; + + // pre-compute values according to exponents + let values = [ + -J * p * dt, + g * p * dt, + -J * p * dt, + g * p * dt, + -J * (1.0 - 3.0 * p) * dt / 2.0, + g * (1.0 - 4.0 * p) * dt, + -J * (1.0 - 3.0 * p) * dt / 2.0, + g * p * dt, + -J * p * dt, + g * p * dt + ]; + + // assign first and last value of `seqA` + set seqA w/= 0 <- -J * p * dt / 2.0; + set seqA w/= len - 1 <- -J * p * dt / 2.0; + + // assign other values to `seqA` or `seqB` + // in an alternating way + for i in 1..len - 2 { + if i % 2 == 0 { + set seqA w/= i <- values[i % 10]; + } + else { + set seqB w/= i <- values[i % 10]; + } + } + + return (seqA, seqB); + } + + operation ApplyAllX(qs : Qubit[], theta : Double) : Unit { + // This applies `Rx` with an angle of `2.0 * theta` to all qubits in `qs` + // using partial application + ApplyToEach(Rx(2.0 * theta, _), qs); + } + + operation ApplyDoubleZ(n : Int, qs : Qubit[], theta : Double, dir : Bool, grp : Bool) : Unit { + let start = grp ? 0 | 1; // Choose either odd or even indices based on group number + + for i in 0..n - 1 { + for j in start..2..n - 2 { // Iterate through even or odd `j`s based on `grp` + // rows and cols are interchanged depending on direction + let (row, col) = dir ? (i, j) | (j, i); + + // Choose first qubit based on row and col + let ind1 = GetQubitIndex(row, col, n); + // Choose second qubit in column if direction is horizontal and next qubit in row if direction is vertical + let ind2 = dir ? GetQubitIndex(row, col + 1, n) | GetQubitIndex(row + 1, col, n); + + within { + CNOT(qs[ind1], qs[ind2]); + } apply { + Rz(2.0 * theta, qs[ind2]); + } + } + } + } + + operation IsingModel2DSim(N : Int, J : Double, g : Double, totTime : Double, dt : Double, eps : Double) : Unit { + use qs = Qubit[N * N]; + let len = Length(qs); + + let p = 1.0 / (4.0 - (4.0 ^ (1.0 / 3.0))); + let t = Ceiling(totTime / dt); + + let seqLen = 10 * t + 1; + + let (seqA, seqB) = SetSequences(seqLen, p, dt, J, g); + + for i in 0..seqLen - 1 { + // for even indexes + if i % 2 == 0 { + ApplyAllX(qs, seqA[i]); + } else { + // iterate through all possible combinations for `dir` and `grp`. + for (dir, grp) in [(true, true), (true, false), (false, true), (false, false)] { + ApplyDoubleZ(N, qs, seqB[i], dir, grp); + } + } + } + } +} diff --git a/samples/estimation/Precalculated.qs b/samples/estimation/Precalculated.qs new file mode 100644 index 0000000000..5e709fcc81 --- /dev/null +++ b/samples/estimation/Precalculated.qs @@ -0,0 +1,25 @@ +/// # Sample +/// Using pre-calculated logical counts for resource estimation +/// +/// # Description +/// This sample demonstrates how to use the `AccountForEstimates` function to +/// estimate the resources required to run a factoring program from pre-calculated +/// logical counts. The logical counts used in this sample are the ones obtained +/// for a 2048-bit integer factoring application based on the implementation +/// described in [[Quantum 5, 433 (2021)](https://quantum-journal.org/papers/q-2021-04-15-433/)]. +/// Our implementation incorporates all techniques described in the paper, except for +/// carry runways. +namespace PrecalculatedEstimates { + open Microsoft.Quantum.ResourceEstimation; + + @EntryPoint() + operation FactoringFromLogicalCounts() : Unit { + use qubits = Qubit[12581]; + + AccountForEstimates( + [TCount(12), RotationCount(12), RotationDepth(12), + CczCount(3731607428), MeasurementCount(1078154040)], + PSSPCLayout(), qubits); + } + +} diff --git a/samples/estimation/ShorRE.qs b/samples/estimation/ShorRE.qs new file mode 100644 index 0000000000..4f3939ca69 --- /dev/null +++ b/samples/estimation/ShorRE.qs @@ -0,0 +1,210 @@ +/// # Sample +/// Estimating Frequency for Shor's algorithm +/// +/// # Description +/// 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. This makes it ideal for +/// use with the Azure Quantum Resource Estimator. +namespace Shors { + open Microsoft.Quantum.Arrays; + open Microsoft.Quantum.Canon; + open Microsoft.Quantum.Convert; + open Microsoft.Quantum.Diagnostics; + open Microsoft.Quantum.Intrinsic; + open Microsoft.Quantum.Math; + open Microsoft.Quantum.Measurement; + open Microsoft.Quantum.Unstable.Arithmetic; + open Microsoft.Quantum.ResourceEstimation; + + @EntryPoint() + operation RunProgram() : Unit { + let bitsize = 31; + + // When chooseing parameters for `EstimateFrequency`, make sure that + // generator and modules are not co-prime + let _ = EstimateFrequency(11, 2^bitsize - 1, bitsize); + } + + /// # 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 Azure 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 { + set frequencyEstimate += 1 <<< (bitsPrecision - 1 - idx); + } + } + + // Return all the qubits used for oracle's eigenstate back to 0 state + // using Microsoft.Quantum.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 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 + // format. + 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 |𝑦⟩ in little-endian format, 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 in IndexRange(y) { + let shiftedC = (c <<< idx) % modulus; + Controlled ModularAddConstant( + [y[idx]], + (modulus, shiftedC, qs)); + } + for idx in IndexRange(y) { + SWAP(y[idx], qs[idx]); + } + let invC = InverseModI(c, modulus); + for idx in IndexRange(y) { + let shiftedC = (invC <<< idx) % modulus; + Controlled ModularAddConstant( + [y[idx]], + (modulus, modulus - shiftedC, qs)); + } + } + + /// # Summary + /// Performs modular in-place addition of a classical constant into a + /// quantum register. + /// + /// Given the classical constants `c` and `modulus`, and an input quantum + /// register |𝑦⟩ in little-endian format, 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 IncByI(ctrls, (c, y + [carry])); + Controlled Adjoint IncByI(ctrls, (modulus, y + [carry])); + Controlled IncByI([carry], (modulus, y)); + Controlled ApplyIfLessOrEqualL(ctrls, (X, IntAsBigInt(c), y, carry)); + } + } + } +}