microsoft/qdk

Public

mirrored fromhttps://github.com/microsoft/qdkAvailable

CodeCommitsIssuesPull requestsActionsInsightsSecurity
iadavis/python-cli

Branches

Tags

  • No tags available.
0Branches0Tags
Go to file
Add file
Code

Clone

HTTPS

Download ZIP

samples/estimation/df-chemistry/src/prepare.qs

222lines · modecode

1// Copyright (c) Microsoft Corporation. All rights reserved.
2// Licensed under the MIT License.
3namespace Microsoft.Quantum.Applications.Chemistry {
4 open Microsoft.Quantum.Arrays;
5 open Microsoft.Quantum.Convert;
6 open Microsoft.Quantum.Diagnostics;
7 open Microsoft.Quantum.Intrinsic;
8 open Microsoft.Quantum.Math;
9 open Microsoft.Quantum.Unstable.Arithmetic;
10 open Microsoft.Quantum.Unstable.TableLookup;
11
12 // ------------------------------------- //
13 // State preparation (public operations) //
14 // ------------------------------------- //
15
16 operation PrepareSingleQubit(p0 : Double, p1 : Double, target : Qubit) : Unit is Adj + Ctl {
17 let oneNorm = p0 + p1;
18 let alpha = ArcCos(Sqrt(p0 / oneNorm));
19
20 Ry(2.0 * alpha, target);
21 }
22
23 operation PrepareUniformSuperposition(numStates : Int, qs : Qubit[]) : Unit is Adj + Ctl {
24 Fact(numStates >= 1, "numStates must be positive");
25 Fact(numStates <= 2^Length(qs), $"numStates must be smaller or equal to {2^Length(qs)}");
26
27 let qsAdjusted = qs[...Ceiling(Lg(IntAsDouble(numStates))) - 1];
28
29 let (factor, pow) = DecomposePowerOf2(numStates);
30
31 if factor == 1 {
32 ApplyToEachCA(H, qsAdjusted[0..pow - 1]);
33 } else {
34 use tgt = Qubit();
35
36 let sqrt = Sqrt(IntAsDouble(1 <<< Length(qsAdjusted)) / IntAsDouble(numStates));
37 let angle = 2.0 * ArcSin(0.5 * sqrt);
38
39 ApplyToEachCA(H, qsAdjusted);
40
41 ApplyIfGreaterL(Ry(2.0 * angle, _), IntAsBigInt(numStates), qsAdjusted, tgt);
42
43 within {
44 ApplyToEachA(H, qsAdjusted[pow...]);
45 } apply {
46 ReflectAboutInteger(0, qsAdjusted[pow...] + [tgt]);
47 Ry(-angle, tgt);
48 }
49
50 X(tgt);
51 }
52 }
53
54 newtype PrepareArbitrarySuperposition = (
55 NIndexQubits : Int,
56 NGarbageQubits : Int,
57 Prepare : (Qubit[], Qubit[], Qubit[]) => Unit is Adj + Ctl,
58 PrepareWithSelect : ((Bool[][], Qubit[], Qubit[]) => Unit is Adj + Ctl, Qubit[], Qubit[], Qubit[]) => Unit is Adj + Ctl
59 );
60
61 function MakePrepareArbitrarySuperposition(targetError : Double, coefficients : Double[]) : PrepareArbitrarySuperposition {
62 let nBitsPrecision = -Ceiling(Lg(0.5 * targetError)) + 1;
63 let positiveCoefficients = Mapped(AbsD, coefficients);
64 let (keepCoeff, altIndex) = DiscretizedProbabilityDistribution(nBitsPrecision, positiveCoefficients);
65 let nCoeffs = Length(positiveCoefficients);
66 let nBitsIndices = Ceiling(Lg(IntAsDouble(nCoeffs)));
67
68 let op = PrepareQuantumROMState(nBitsPrecision, nCoeffs, nBitsIndices, keepCoeff, altIndex, [], Select, _, _, _);
69 let opWithSelect = PrepareQuantumROMState(nBitsPrecision, nCoeffs, nBitsIndices, keepCoeff, altIndex, [], _, _, _, _);
70 let (nIndexQubits, nGarbageQubits) = ArbitrarySuperpositionRegisterLengths(targetError, nCoeffs);
71 return PrepareArbitrarySuperposition(nIndexQubits, nGarbageQubits, op, opWithSelect);
72 }
73
74 function MakePrepareArbitrarySuperpositionWithData(targetError : Double, coefficients : Double[], data : Bool[][]) : PrepareArbitrarySuperposition {
75 let nBitsPrecision = -Ceiling(Lg(0.5 * targetError)) + 1;
76 let positiveCoefficients = Mapped(AbsD, coefficients);
77 let (keepCoeff, altIndex) = DiscretizedProbabilityDistribution(nBitsPrecision, positiveCoefficients);
78 let nCoeffs = Length(positiveCoefficients);
79 let nBitsIndices = Ceiling(Lg(IntAsDouble(nCoeffs)));
80
81 let op = PrepareQuantumROMState(nBitsPrecision, nCoeffs, nBitsIndices, keepCoeff, altIndex, data, Select, _, _, _);
82 let opWithSelect = PrepareQuantumROMState(nBitsPrecision, nCoeffs, nBitsIndices, keepCoeff, altIndex, data, _, _, _, _);
83 let (nIndexQubits, nGarbageQubits) = ArbitrarySuperpositionRegisterLengths(targetError, nCoeffs);
84 return PrepareArbitrarySuperposition(nIndexQubits, nGarbageQubits + Length(data[0]), op, opWithSelect);
85 }
86
87 // -------------------------------------- //
88 // State preparation (private operations) //
89 // -------------------------------------- //
90
91 internal function DecomposePowerOf2(number : Int) : (Int, Int) {
92 mutable pow = 0;
93 mutable factor = number;
94
95 while factor % 2 == 0 {
96 set factor /= 2;
97 set pow += 1;
98 }
99
100 (factor, pow)
101 }
102
103 internal function ArbitrarySuperpositionRegisterLengths(targetError : Double, nCoefficients : Int) : (Int, Int) {
104 Fact(targetError > 0.0, "targetError must be positive");
105 Fact(nCoefficients > 0, "nCoefficients must be positive");
106
107 let nBitsPrecision = -Ceiling(Lg(0.5 * targetError)) + 1;
108 let nIndexQubits = Ceiling(Lg(IntAsDouble(nCoefficients)));
109 let nGarbageQubits = nIndexQubits + 2 * nBitsPrecision + 1;
110 (nIndexQubits, nGarbageQubits)
111 }
112
113 // Computes discretized probability distribution as described in Section 3
114 // and Fig. 13 in [arXiv:1805.03662](https://arxiv.org/pdf/1805.03662.pdf)
115 internal function DiscretizedProbabilityDistribution(bitsPrecision : Int, coefficients : Double[]) : (Int[], Int[]) {
116 let oneNorm = PNorm(1.0, coefficients);
117 let nCoefficients = Length(coefficients);
118 Fact(bitsPrecision <= 31, $"Bits of precision {bitsPrecision} unsupported. Max is 31.");
119 Fact(nCoefficients > 1, "Cannot prepare state with less than 2 coefficients.");
120 Fact(oneNorm != 0.0, "State must have at least one coefficient > 0");
121
122 let barHeight = 2^bitsPrecision - 1;
123
124 mutable altIndex = SequenceI(0, nCoefficients - 1);
125 mutable keepCoeff = Mapped(
126 coefficient -> Round((AbsD(coefficient) / oneNorm) * IntAsDouble(nCoefficients) * IntAsDouble(barHeight)),
127 coefficients
128 );
129
130 // Calculate difference between number of discretized bars vs. maximum
131 let bars = Fold((state, value) -> state + value - barHeight, 0, keepCoeff);
132
133 // Uniformly distribute excess bars across coefficients.
134 for idx in 0..AbsI(bars) - 1 {
135 set keepCoeff w/= idx <- keepCoeff[idx] + (bars > 0 ? -1 | + 1);
136 }
137
138 mutable barSink = [];
139 mutable barSource = [];
140
141 for idxCoeff in IndexRange(keepCoeff) {
142 if keepCoeff[idxCoeff] > barHeight {
143 set barSource += [idxCoeff];
144 } elif keepCoeff[idxCoeff] < barHeight {
145 set barSink += [idxCoeff];
146 }
147 }
148
149 for rep in 0..nCoefficients * 10 {
150 if Length(barSink) > 0 and Length(barSource) > 0 {
151 let idxSink = Tail(barSink);
152 let idxSource = Tail(barSource);
153 set barSink = Most(barSink);
154 set barSource = Most(barSource);
155
156 set keepCoeff w/= idxSource <- keepCoeff[idxSource] - barHeight + keepCoeff[idxSink];
157 set altIndex w/= idxSink <- idxSource;
158
159 if keepCoeff[idxSource] < barHeight {
160 set barSink += [idxSource];
161 } elif keepCoeff[idxSource] > barHeight {
162 set barSource += [idxSource];
163 }
164 } elif Length(barSource) > 0 {
165 let idxSource = Tail(barSource);
166 set barSource = Most(barSource);
167 set keepCoeff w/= idxSource <- barHeight;
168 } else {
169 return (keepCoeff, altIndex);
170 }
171 }
172
173 return (keepCoeff, altIndex);
174 }
175
176 // Used in QuantumROM implementation.
177 internal operation PrepareQuantumROMState(
178 nBitsPrecision : Int,
179 nCoeffs : Int,
180 nBitsIndices : Int,
181 keepCoeff : Int[],
182 altIndex : Int[],
183 data : Bool[][],
184 selectOperation : (Bool[][], Qubit[], Qubit[]) => Unit is Adj + Ctl,
185 indexRegister : Qubit[],
186 dataQubits : Qubit[],
187 garbageRegister : Qubit[]
188 ) : Unit is Adj + Ctl {
189 let garbageIdx0 = nBitsIndices;
190 let garbageIdx1 = garbageIdx0 + nBitsPrecision;
191 let garbageIdx2 = garbageIdx1 + nBitsPrecision;
192 let garbageIdx3 = garbageIdx2 + 1;
193
194 let altIndexRegister = garbageRegister[0..garbageIdx0 - 1];
195 let keepCoeffRegister = garbageRegister[garbageIdx0..garbageIdx1 - 1];
196 let uniformKeepCoeffRegister = garbageRegister[garbageIdx1..garbageIdx2 - 1];
197 let flagQubit = garbageRegister[garbageIdx3 - 1];
198 let dataRegister = dataQubits;
199 let altDataRegister = garbageRegister[garbageIdx3...];
200
201 // Create uniform superposition over index and alt coeff register.
202 PrepareUniformSuperposition(nCoeffs, indexRegister);
203 ApplyToEachCA(H, uniformKeepCoeffRegister);
204
205 // Write bitstrings to altIndex and keepCoeff register.
206 let target = keepCoeffRegister + altIndexRegister + dataRegister + altDataRegister;
207 let selectData = MappedOverRange(idx -> IntAsBoolArray(keepCoeff[idx], Length(keepCoeffRegister)) + IntAsBoolArray(altIndex[idx], Length(altIndexRegister)) + (IsEmpty(data) ? [] | data[idx] + data[altIndex[idx]]), 0..nCoeffs - 1);
208 selectOperation(selectData, indexRegister, target);
209
210 // Perform comparison
211 ApplyIfGreaterLE(X, uniformKeepCoeffRegister, keepCoeffRegister, flagQubit);
212
213 let indexRegisterSize = Length(indexRegister);
214
215 // Swap in register based on comparison
216 let lhs = indexRegister + dataRegister;
217 let rhs = altIndexRegister + altDataRegister;
218 for i in IndexRange(lhs) {
219 Controlled SWAP([flagQubit], (lhs[i], rhs[i]));
220 }
221 }
222}
223