microsoft/qdk
Publicmirrored fromhttps://github.com/microsoft/qdkAvailable
samples/estimation/EkeraHastadFactoring.qs
195lines · modecode
| 1 | /// # Sample |
| 2 | /// Resource Estimation for Integer Factoring |
| 3 | /// |
| 4 | /// # Description |
| 5 | /// In this sample we concentrate on costing quantum part in the algorithm for |
| 6 | /// factoring RSA integers based on Ekerå and Håstad |
| 7 | /// [ia.cr/2017/077](https://eprint.iacr.org/2017/077) based on the |
| 8 | /// implementation described in |
| 9 | /// [arXiv:1905.09749](https://arxiv.org/abs/1905.09749). This makes it ideal |
| 10 | /// for use with the Azure Quantum Resource Estimator. |
| 11 | import Std.Convert.*; |
| 12 | import Std.Math.*; |
| 13 | import Std.ResourceEstimation.*; |
| 14 | import Std.Arrays.*; |
| 15 | import Std.Arithmetic.*; |
| 16 | import Std.TableLookup.*; |
| 17 | |
| 18 | // !!! IMPORTANT !!! |
| 19 | // When computing resource estimtes from the VS Code plugin directly on this |
| 20 | // file, make sure that you set the error budget to 0.333. |
| 21 | |
| 22 | operation Main() : Unit { |
| 23 | // Try different instances of the algorithm by commenting in and out |
| 24 | // the following lines. You can find more RSA numbers at |
| 25 | // https://en.wikipedia.org/wiki/RSA_numbers |
| 26 | |
| 27 | // RSA-100 (330 bits) |
| 28 | EkeraHastad(330, 1522605027922533360535618378132637429718068114961380688657908494580122963258952897654000350692006139L, 7L); |
| 29 | |
| 30 | // RSA-1024 (1024 bits) |
| 31 | // EkeraHastad(1024, 135066410865995223349603216278805969938881475605667027524485143851526510604859533833940287150571909441798207282164471551373680419703964191743046496589274256239341020864383202110372958725762358509643110564073501508187510676594629205563685529475213500852879416377328533906109750544334999811150056977236890927563L, 7L); |
| 32 | |
| 33 | // RSA-2048 (2048 bits) |
| 34 | // EkeraHastad(2048, 25195908475657893494027183240048398571429282126204032027777137836043662020707595556264018525880784406918290641249515082189298559149176184502808489120072844992687392807287776735971418347270261896375014971824691165077613379859095700097330459748808428401797429100642458691817195118746121515172654632282216869987549182422433637259085141865462043576798423387184774447920739934236584823824281198163815010674810451660377306056201619676256133844143603833904414952634432190114657544454178424020924616515723350778707749817125772467962926386356373289912154831438167899885040445364023527381951378636564391212010397122822120720357L, 7L); |
| 35 | } |
| 36 | |
| 37 | /// # Summary |
| 38 | /// Main algorithm based on quantum phase estimation |
| 39 | /// |
| 40 | /// # Reference |
| 41 | /// [ia.cr/2017/077, Section 4.3](https://eprint.iacr.org/2017/077) |
| 42 | operation EkeraHastad(numBits : Int, N : BigInt, g : BigInt) : Unit { |
| 43 | let x = ExpModL(g, ((N - 1L) / 2L), N); |
| 44 | let xinv = InverseModL(x, N); |
| 45 | |
| 46 | let m = numBits / 2; |
| 47 | use c1 = Qubit[2 * m]; |
| 48 | use c2 = Qubit[m]; |
| 49 | use target = Qubit[numBits]; |
| 50 | |
| 51 | let ne = 3 * m; |
| 52 | let cpad = Ceiling(2.0 * Lg(IntAsDouble(numBits)) + Lg(IntAsDouble(ne)) + 10.0); |
| 53 | |
| 54 | // The algorithm uses the coset representation to replace modular |
| 55 | // addition inside MultiplyExpMod with regular addition. |
| 56 | use padding = Qubit[cpad]; |
| 57 | |
| 58 | ApplyToEach(H, c1 + c2); |
| 59 | InitCoset(N, target, padding); |
| 60 | MultiplyExpMod(g, N, c1, target + padding); |
| 61 | MultiplyExpMod(xinv, N, c2, target + padding); |
| 62 | |
| 63 | Adjoint ApplyQFT(c1); |
| 64 | Adjoint ApplyQFT(c2); |
| 65 | } |
| 66 | |
| 67 | // ------------------------------ // |
| 68 | // Modular arithmetic (constants) // |
| 69 | // ------------------------------ // |
| 70 | |
| 71 | /// Window size for exponentiation (c_exp) |
| 72 | internal function ExponentWindowLength_() : Int { 5 } |
| 73 | |
| 74 | /// Window size for multiplication (c_mul) |
| 75 | internal function MultiplicationWindowLength_() : Int { 5 } |
| 76 | |
| 77 | // ------------------------------- // |
| 78 | // Modular arithmetic (operations) // |
| 79 | // ------------------------------- // |
| 80 | |
| 81 | /// # Summary |
| 82 | /// Encodes register in coset representation |
| 83 | /// |
| 84 | /// # Reference |
| 85 | /// [arXiv:quant-ph/0601097, Section 4.1](https://arxiv.org/abs/quant-ph/0601097) |
| 86 | internal operation InitCoset(mod : BigInt, xs : Qubit[], padding : Qubit[]) : Unit { |
| 87 | use helper = Qubit(); |
| 88 | let cpad = Length(padding); |
| 89 | let n = Length(xs); |
| 90 | |
| 91 | let combined = xs + padding; |
| 92 | |
| 93 | for j in 0..cpad - 1 { |
| 94 | Controlled IncByLUsingIncByLE([helper], (RippleCarryCGIncByLE, mod, combined[j..j + n - 1])); |
| 95 | |
| 96 | ApplyIfLessOrEqualL(X, mod, combined[j..j + n - 1], helper); |
| 97 | } |
| 98 | } |
| 99 | |
| 100 | /// # Summary |
| 101 | /// Computes zs *= (base ^ xs) % mod (for a large register xs) |
| 102 | /// |
| 103 | /// # Reference |
| 104 | /// [arXiv:1905.07682, Fig. 7](https://arxiv.org/abs/1905.07682) |
| 105 | internal operation MultiplyExpMod( |
| 106 | base : BigInt, |
| 107 | mod : BigInt, |
| 108 | xs : Qubit[], |
| 109 | zs : Qubit[] |
| 110 | ) : Unit { |
| 111 | let expWindows = Chunks(ExponentWindowLength_(), xs); |
| 112 | |
| 113 | within { |
| 114 | RepeatEstimates(Length(expWindows)); |
| 115 | } apply { |
| 116 | let i = 0; // in simulation this i must be iterated over IndexRange(expWindows) |
| 117 | |
| 118 | let adjustedBase = ExpModL(base, 1L <<< (i * ExponentWindowLength_()), mod); |
| 119 | MultiplyExpModWindowed(adjustedBase, mod, expWindows[i], zs); |
| 120 | } |
| 121 | } |
| 122 | |
| 123 | /// # Summary |
| 124 | /// Computes zs *= (base ^ xs) % mod (for a small register xs) |
| 125 | /// |
| 126 | /// # Reference |
| 127 | /// [arXiv:1905.07682, Fig. 6](https://arxiv.org/abs/1905.07682) |
| 128 | internal operation MultiplyExpModWindowed( |
| 129 | base : BigInt, |
| 130 | mod : BigInt, |
| 131 | xs : Qubit[], |
| 132 | zs : Qubit[] |
| 133 | ) : Unit { |
| 134 | let n = Length(zs); |
| 135 | |
| 136 | use qs = Qubit[n]; |
| 137 | AddExpModWindowed(base, mod, 1, xs, zs, qs); |
| 138 | AddExpModWindowed(InverseModL(base, mod), mod, -1, xs, qs, zs); |
| 139 | for i in IndexRange(zs) { |
| 140 | SWAP(zs[i], qs[i]); |
| 141 | } |
| 142 | } |
| 143 | |
| 144 | /// # Summary |
| 145 | /// Computes zs += ys * (base ^ xs) % mod (for small registers xs and ys) |
| 146 | /// |
| 147 | /// # Reference |
| 148 | /// [arXiv:1905.07682, Fig. 5](https://arxiv.org/abs/1905.07682) |
| 149 | /// |
| 150 | /// # Remark |
| 151 | /// Unlike in the reference, this implementation uses regular addition |
| 152 | /// instead of modular addition because the target register is encoded |
| 153 | /// using the coset representation. |
| 154 | internal operation AddExpModWindowed( |
| 155 | base : BigInt, |
| 156 | mod : BigInt, |
| 157 | sign : Int, |
| 158 | xs : Qubit[], |
| 159 | ys : Qubit[], |
| 160 | zs : Qubit[] |
| 161 | ) : Unit { |
| 162 | // split factor into parts |
| 163 | let factorWindows = Chunks(MultiplicationWindowLength_(), ys); |
| 164 | |
| 165 | within { |
| 166 | RepeatEstimates(Length(factorWindows)); |
| 167 | } apply { |
| 168 | let i = 0; // in simulation this i must be iterated over IndexRange(factorWindows) |
| 169 | |
| 170 | // compute data for table lookup |
| 171 | let factorValue = ExpModL(2L, IntAsBigInt(i * MultiplicationWindowLength_()), mod); |
| 172 | let data = LookupData(factorValue, Length(xs), Length(factorWindows[i]), base, mod, sign, Length(zs)); |
| 173 | |
| 174 | use output = Qubit[Length(data[0])]; |
| 175 | |
| 176 | within { |
| 177 | Select(data, xs + factorWindows[i], output); |
| 178 | } apply { |
| 179 | RippleCarryCGIncByLE(output, zs); |
| 180 | } |
| 181 | } |
| 182 | } |
| 183 | |
| 184 | internal function LookupData(factor : BigInt, expLength : Int, mulLength : Int, base : BigInt, mod : BigInt, sign : Int, numBits : Int) : Bool[][] { |
| 185 | mutable data = [[false, size = numBits], size = 2^(expLength + mulLength)]; |
| 186 | for b in 0..2^mulLength - 1 { |
| 187 | for a in 0..2^expLength - 1 { |
| 188 | let idx = b * 2^expLength + a; |
| 189 | let value = ModulusL(factor * IntAsBigInt(b) * IntAsBigInt(sign) * (base^a), mod); |
| 190 | set data w/= idx <- BigIntAsBoolArray(value, numBits); |
| 191 | } |
| 192 | } |
| 193 | |
| 194 | data |
| 195 | } |
| 196 | |