練習 - 估計實際問題的資源

10 分鐘

Shor 的分解演算法是最著名的量子演算法之一。它提供超過任何已知傳統分解演算法的指數型加速。

傳統密碼編譯會使用實體或數學方式 (例如計算困難) 來完成工作。熱門的密碼編譯通訊協定是 Rivest-Shamir–Adleman (RSA) 配置，其以假設使用傳統電腦難以分解質數為基礎。

Shor 演算法意味著只要有足夠龐大的量子電腦，便可以破解公開金鑰加密。評估 Shor 演算法所需的資源，對於評估這些密碼編譯配置類型的弱點很重要。

在下列練習中，您將計算分解 2,048 位元整數的資源估計值。針對此應用程式，您將直接從預先計算的邏輯資源估計值計算實體資源估計值。針對容許的錯誤預算，您將會使用 $\epsilon = 1/3$。

撰寫 Shor 演算法

在 Visual Studio Code 中，選取 [檢視]>[命令選擇區]，然後選取 [建立：新增 Jupyter Notebook]。
將筆記本儲存為 ShorRE.ipynb。
在第一個儲存格中，匯入 qsharp 套件：
```
import qsharp
```

使用 Microsoft.Quantum.ResourceEstimation 命名空間來定義 Shor 整數分解演算法的快取最佳化版本。新增儲存格並複製貼上下列程式碼：

%%qsharp
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;

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 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 oracles 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 LittleEndian 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 by using the `LittleEndian` type.
    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 (as LittleEndian) |𝑦⟩, 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 (as LittleEndian) |𝑦⟩, 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) {
        set nZeroes += 1;
        set 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;
}

估計 Shor 演算法

現在，使用預設假設來估計 RunProgram 作業的實體資源。新增儲存格並複製貼上下列程式碼：

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

qsharp.estimate 函式會建立結果物件，您可以用來顯示具有整體實體資源計數的資料表。您可以藉由折迭群組來檢查成本詳細資料，這些群組具有詳細資訊。

例如，摺疊 [邏輯量子位元參數] 群組，以查看程式碼距離為 21，而實體量子位元的數目為 882。

邏輯量子位元參數	值
QEC 配置	surface_code
程式碼距離	21
實際量子位元	882
邏輯週期時間	8 毫秒
邏輯量子位元錯誤率	3.00E-13
交叉前置要素	0.03
錯誤修正臨界值	0.01
邏輯週期時間公式	(4 * `twoQubitGateTime` + 2 * `oneQubitMeasurementTime`) * `codeDistance`
實體量子位元公式	2 * `codeDistance` * `codeDistance`

提示

如需更精簡的輸出資料表版本，您可以使用 result.summary。

空間圖

演算法和 T Factory 所使用的實體量子位元分佈是可能會影響演算法設計的因素。您可以使用 qsharp-widgets 套件將此分佈視覺化，以進一步了解演算法的估計空間需求。

from qsharp_widgets import SpaceChart
SpaceChart(result)

在此範例中，執行演算法所需的實體量子位元數目 829766，其中 196686 為演算法量子位元，而其中 633080 為 T Factory 量子位元。

顯示資源估算器空間圖的螢幕擷取畫面。

比較不同量子位元技術的資源估計值

Azure Quantum 資源估算器可讓您執行多個目標參數的組態，並比較結果。當您想要比較不同量子位元模型、QEC 配置或錯誤預算的成本時，這會很有用。

您也可以使用 EstimatorParams 物件來建構估計參數的清單。

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(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)

量子位元模型	邏輯量子位元	邏輯深度	T 狀態	程式碼距離	T Factory	T Factory 分數	實際量子位元	rQOPS	實際執行時間
閘道型 µs, 10⁻³	223 3.64M	4.70M	17	13	40.54 %	216.77k	21.86k	10 個小時
閘道型 µs, 10⁻⁴	223	3.64M	4.70M	9	14	43.17 %	63.57k	41.30k	5 個小時
閘道型 ns, 10⁻³	223 3.64M	4.70M	17	16	69.08 %	416.89k	32.79M	25 secs
閘道型 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

從邏輯資源計數擷取資源估計值

如果您已經知道作業的部分估計值，資源估算器可讓您將已知的估計值納入整體程式成本，以減少執行時間。您可以使用 LogicalCounts 類別，從預先計算的資源估計值中擷取邏輯資源估計值。

選取 [程式碼] 以新增儲存格，然後輸入並執行下列程式碼：

logical_counts = LogicalCounts({
    'numQubits': 12581,
    'tCount': 12,
    'rotationCount': 12,
    'rotationDepth': 12,
    'cczCount': 3731607428,
    'measurementCount': 1078154040})

logical_counts.estimate(params).summary_data_frame(labels=labels)

量子位元模型	邏輯量子位元	邏輯深度	T 狀態	程式碼距離	T Factory	T Factory 分數	實際量子位元	實際執行時間
閘道型 µs, 10⁻³	25481	1.2e+10	1.5e+10	27	13	0.6%	37.38M	6 年
閘道型 µs, 10⁻⁴	25481	1.2e+10	1.5e+10	13	14	0.8%	8.68M	3 年
閘道型 ns, 10⁻³	25481	1.2e+10	1.5e+10	27	15	1.3%	37.65M	2 天
閘道型 ns, 10⁻⁴	25481	1.2e+10	1.5e+10	13	18	1.2%	8.72M	18 小時
Majorana ns, 10⁻⁴	25481	1.2e+10	1.5e+10	15	15	1.3%	26.11M	15 小時
Majorana ns, 10⁻⁶	25481	1.2e+10	1.5e+10	7	13	0.5%	6.25M	7 小時

推論

在最不理想的案例中，使用閘道型 µs 量子位元 (在奈秒制中具有作業時間的量子位元，例如超導量子位元) 與表面 QEC 程式碼的量子電腦，需要 6 年和 37.38M 個量子位元，才能使用 Shor 的演算法來分解 2,048 位元整數。

如果您使用不同的量子位元技術，例如閘道型 ns 離子量子位元及相同的表面碼，則量子位元的數目不會變更太多，但執行時間在最不理想的情況下會變成 2 天，而理想的案例中為 18 小時。例如，如果您使用 Majorana 型量子位元變更量子位元技術和 QEC 程式碼，則使用 Shor 演算法來分解 2,048 位元整數時，可以在最佳案例中以 6.25M 個量子位元的陣列來完成。

從實驗中，您可以判斷使用 Majorana 量子位元和 Floquet QEC 程式碼是執行 Shor 演算法及分解 2,048 位元整數的最佳選擇。

練習 - 估計實際問題的資源

撰寫 Shor 演算法

估計 Shor 演算法

空間圖

比較不同量子位元技術的資源估計值

從邏輯資源計數擷取資源估計值

推論

意見反應