diff --git a/library/std/src/Std/StatePreparation.qs b/library/std/src/Std/StatePreparation.qs index 0a28b97878..610c4b2f90 100644 --- a/library/std/src/Std/StatePreparation.qs +++ b/library/std/src/Std/StatePreparation.qs @@ -63,8 +63,35 @@ import /// # See Also /// - Std.StatePreparation.ApproximatelyPreparePureStateCP operation PreparePureStateD(coefficients : Double[], qubits : Qubit[]) : Unit is Adj + Ctl { - let coefficientsAsComplexPolar = Mapped(a -> ComplexAsComplexPolar(Complex(a, 0.0)), coefficients); - ApproximatelyPreparePureStateCP(0.0, coefficientsAsComplexPolar, qubits); + let nQubits = Length(qubits); + // pad coefficients at tail length to a power of 2. + let coefficientsPadded = Padded(-2^nQubits, 0.0, coefficients); + let idxTarget = 0; + + // Note we use the reversed qubits array to get the endianness ordering that we expect + // when corresponding qubit state to state vector index. + let qubits = Reversed(qubits); + + // Since we know the coefficients are real, we can optimize the first round of adjoint approximate unpreparation by directly + // computing the disentangling angles and the new coefficients on those doubles without producing intermediate complex numbers. + + // For each 2D block, compute disentangling single-qubit rotation parameters + let (disentanglingY, disentanglingZ, newCoefficients) = StatePreparationSBMComputeCoefficientsD(coefficientsPadded); + + if nQubits == 1 { + let (abs, arg) = newCoefficients[0]!; + if (AbsD(arg) > 0.0) { + Adjoint Exp([PauliI], -1.0 * arg, [qubits[idxTarget]]); + } + } elif (Any(c -> AbsComplexPolar(c) > 0.0, newCoefficients)) { + // Some coefficients are outside tolerance + let newControl = 2..(nQubits - 1); + let newTarget = 1; + Adjoint ApproximatelyUnprepareArbitraryState(0.0, newCoefficients, newControl, newTarget, qubits); + } + + Adjoint ApproximatelyMultiplexPauli(0.0, disentanglingY, PauliY, qubits[1...], qubits[0]); + Adjoint ApproximatelyMultiplexPauli(0.0, disentanglingZ, PauliZ, qubits[1...], qubits[0]); } /// # Summary @@ -148,10 +175,9 @@ operation ApproximatelyUnprepareArbitraryState( ) : Unit is Adj + Ctl { // For each 2D block, compute disentangling single-qubit rotation parameters - let (disentanglingY, disentanglingZ, newCoefficients) = StatePreparationSBMComputeCoefficients(coefficients); + let (disentanglingY, disentanglingZ, newCoefficients) = StatePreparationSBMComputeCoefficientsCP(coefficients); if (AnyOutsideToleranceD(tolerance, disentanglingZ)) { ApproximatelyMultiplexPauli(tolerance, disentanglingZ, PauliZ, register[rngControl], register[idxTarget]); - } if (AnyOutsideToleranceD(tolerance, disentanglingY)) { ApproximatelyMultiplexPauli(tolerance, disentanglingY, PauliY, register[rngControl], register[idxTarget]); @@ -226,7 +252,37 @@ operation ApproximatelyMultiplexPauli( /// # Summary /// Implementation step of arbitrary state preparation procedure. -function StatePreparationSBMComputeCoefficients( +/// This version optimized for purely real coefficients represented by an array of doubles. +function StatePreparationSBMComputeCoefficientsD( + coefficients : Double[] +) : (Double[], Double[], ComplexPolar[]) { + mutable disentanglingZ = []; + mutable disentanglingY = []; + mutable newCoefficients = []; + + for idxCoeff in 0..2..Length(coefficients) - 1 { + let (rt, phi, theta) = { + let abs0 = AbsD(coefficients[idxCoeff]); + let abs1 = AbsD(coefficients[idxCoeff + 1]); + let arg0 = coefficients[idxCoeff] < 0.0 ? PI() | 0.0; + let arg1 = coefficients[idxCoeff + 1] < 0.0 ? PI() | 0.0; + let r = Sqrt(abs0 * abs0 + abs1 * abs1); + let t = 0.5 * (arg0 + arg1); + let phi = arg1 - arg0; + let theta = 2.0 * ArcTan2(abs1, abs0); + (ComplexPolar(r, t), phi, theta) + }; + set disentanglingZ += [0.5 * phi]; + set disentanglingY += [0.5 * theta]; + set newCoefficients += [rt]; + } + + return (disentanglingY, disentanglingZ, newCoefficients); +} + +/// # Summary +/// Implementation step of arbitrary state preparation procedure. +function StatePreparationSBMComputeCoefficientsCP( coefficients : ComplexPolar[] ) : (Double[], Double[], ComplexPolar[]) { @@ -321,25 +377,28 @@ operation ApproximatelyMultiplexZ( target : Qubit ) : Unit is Adj + Ctl { - body (...) { - // pad coefficients length at tail to a power of 2. - let coefficientsPadded = Padded(-2^Length(control), 0.0, coefficients); - - if Length(coefficientsPadded) == 1 { - // Termination case - if AbsD(coefficientsPadded[0]) > tolerance { - Exp([PauliZ], coefficientsPadded[0], [target]); + body ... { + // We separately compute the operation sequence for the multiplex Z steps in a function, which + // provides a performance improvement during partial-evaluation for code generation. + let multiplexZParams = GenerateMultiplexZParams(tolerance, coefficients, control, target); + for (angle, qs) in multiplexZParams { + if Length(qs) == 2 { + CNOT(qs[0], qs[1]); + } elif AbsD(angle) > tolerance { + Exp([PauliZ], angle, qs); } - } else { - // Compute new coefficients. - let (coefficients0, coefficients1) = MultiplexZCoefficients(coefficientsPadded); - ApproximatelyMultiplexZ(tolerance, coefficients0, Most(control), target); - if AnyOutsideToleranceD(tolerance, coefficients1) { - within { - CNOT(Tail(control), target); - } apply { - ApproximatelyMultiplexZ(tolerance, coefficients1, Most(control), target); - } + } + } + + adjoint ... { + // We separately compute the operation sequence for the adjoint multiplex Z steps in a function, which + // provides a performance improvement during partial-evaluation for code generation. + let adjMultiplexZParams = GenerateAdjMultiplexZParams(tolerance, coefficients, control, target); + for (angle, qs) in adjMultiplexZParams { + if Length(qs) == 2 { + CNOT(qs[0], qs[1]); + } elif AbsD(angle) > tolerance { + Exp([PauliZ], -angle, qs); } } } @@ -357,8 +416,70 @@ operation ApproximatelyMultiplexZ( } } } + + controlled adjoint (controlRegister, ...) { + // pad coefficients length to a power of 2. + let coefficientsPadded = Padded(2^(Length(control) + 1), 0.0, Padded(-2^Length(control), 0.0, coefficients)); + let (coefficients0, coefficients1) = MultiplexZCoefficients(coefficientsPadded); + if AnyOutsideToleranceD(tolerance, coefficients1) { + within { + Controlled X(controlRegister, target); + } apply { + Adjoint ApproximatelyMultiplexZ(tolerance, coefficients1, control, target); + } + } + Adjoint ApproximatelyMultiplexZ(tolerance, coefficients0, control, target); + } } +// Provides the sequence of angles or entangling CNOTs to apply for the multiplex Z step of the state preparation procedure, given a set of coefficients and control and target qubits. +function GenerateMultiplexZParams( + tolerance : Double, + coefficients : Double[], + control : Qubit[], + target : Qubit +) : (Double, Qubit[])[] { + // pad coefficients length at tail to a power of 2. + let coefficientsPadded = Padded(-2^Length(control), 0.0, coefficients); + + if Length(coefficientsPadded) == 1 { + // Termination case + [(coefficientsPadded[0], [target])] + } else { + // Compute new coefficients. + let (coefficients0, coefficients1) = MultiplexZCoefficients(coefficientsPadded); + mutable params = GenerateMultiplexZParams(tolerance, coefficients0, Most(control), target); + params += [(0.0, [Tail(control), target])] + GenerateMultiplexZParams(tolerance, coefficients1, Most(control), target); + params += [(0.0, [Tail(control), target])]; + params + } +} + +// Provides the sequence of angles or entangling CNOTs to apply for the adjoint of the multiplex Z step of the state preparation procedure, given a set of coefficients and control and target qubits. +// Note that the adjoint sequence is NOT the reverse of the forward sequence due to the structure of the multiplex Z step, which applies the disentangling rotations in between the two recursive calls. +function GenerateAdjMultiplexZParams( + tolerance : Double, + coefficients : Double[], + control : Qubit[], + target : Qubit +) : (Double, Qubit[])[] { + // pad coefficients length at tail to a power of 2. + let coefficientsPadded = Padded(-2^Length(control), 0.0, coefficients); + + if Length(coefficientsPadded) == 1 { + // Termination case + [(coefficientsPadded[0], [target])] + } else { + // Compute new coefficients. + let (coefficients0, coefficients1) = MultiplexZCoefficients(coefficientsPadded); + mutable params = [(0.0, [Tail(control), target])] + GenerateAdjMultiplexZParams(tolerance, coefficients1, Most(control), target); + params += [(0.0, [Tail(control), target])]; + params += GenerateAdjMultiplexZParams(tolerance, coefficients0, Most(control), target); + params + } +} + + /// # Summary /// Implementation step of multiply-controlled Z rotations. function MultiplexZCoefficients(coefficients : Double[]) : (Double[], Double[]) {