GitVita
New repositoryImport repositoryNew issue
Your repositories Your dashboard Sign in

microsoft/qdk

Public

mirrored from https://github.com/microsoft/qdkAvailable

Watch0Fork0Star0
CodeCommitsIssuesPull requestsActionsInsightsSecurity
alex/testHarnessIntegTests

Branches

Tags

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

Clone

HTTPS

Download ZIP

qdk/samples/estimation/df-chemistry/src

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

220lines · modecode

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