From fadccd5c7544314c379570e3092d2ff2cd1a8105 Mon Sep 17 00:00:00 2001 From: Daniel Strano Date: Wed, 26 Dec 2018 10:12:57 -0500 Subject: [PATCH] WIP: Adding logarithm of gate method (#117) * Adding logarithm of gate method * Adding 2x2 exp/log convenience methods * Fixing imaginary Exp convention * Anti-controlled Exp options * Renaming convenience exp2x2/log2x2/mul2x2 * HamiltonianOp control bit toggles * Basic HamiltonianOp toggle unit test * Commentary on HamiltonianOp toggles --- include/common/qrack_types.hpp | 8 ++- include/hamiltonian.hpp | 32 +++++++++- include/qinterface.hpp | 15 ++++- src/qinterface/operators.cpp | 18 +++++- src/qinterface/protected.cpp | 103 ++++++++++++++++++++++++++++++++- src/qinterface/rotational.cpp | 98 +++++++------------------------ test/tests.cpp | 48 ++++++++++++++- 7 files changed, 232 insertions(+), 90 deletions(-) diff --git a/include/common/qrack_types.hpp b/include/common/qrack_types.hpp index 2f424d2af..c801dee5b 100644 --- a/include/common/qrack_types.hpp +++ b/include/common/qrack_types.hpp @@ -12,6 +12,8 @@ #pragma once +#include + #define bitLenInt uint8_t #define bitCapInt uint64_t #define bitsInByte 8 @@ -39,4 +41,8 @@ namespace Qrack { typedef std::shared_ptr BitOp; -} + +void mul2x2(complex* left, complex* right, complex* out); +void exp2x2(complex* matrix2x2, complex* outMatrix2x2); +void log2x2(complex* matrix2x2, complex* outMatrix2x2); +} // namespace Qrack diff --git a/include/hamiltonian.hpp b/include/hamiltonian.hpp index 5060170ae..44d85101e 100644 --- a/include/hamiltonian.hpp +++ b/include/hamiltonian.hpp @@ -27,22 +27,34 @@ struct HamiltonianOp { BitOp matrix; bitLenInt* controls; bitLenInt controlLen; + bool anti; + bool* toggles; HamiltonianOp(bitLenInt target, BitOp mtrx) : targetBit(target) , matrix(mtrx) , controls(NULL) , controlLen(0) + , anti(false) + , toggles(NULL) { } - HamiltonianOp(bitLenInt* ctrls, bitLenInt ctrlLen, bitLenInt target, BitOp mtrx) + HamiltonianOp(bitLenInt* ctrls, bitLenInt ctrlLen, bitLenInt target, BitOp mtrx, bool antiCtrled = false, + bool* ctrlToggles = NULL) : targetBit(target) , matrix(mtrx) , controls(new bitLenInt[ctrlLen]) , controlLen(ctrlLen) + , anti(antiCtrled) + , toggles(NULL) { std::copy(ctrls, ctrls + ctrlLen, controls); + + if (ctrlToggles) { + toggles = new bool[ctrlLen]; + std::copy(ctrlToggles, ctrlToggles + ctrlLen, toggles); + } } ~HamiltonianOp() @@ -50,6 +62,10 @@ struct HamiltonianOp { if (controls) { delete[] controls; } + + if (toggles) { + delete[] toggles; + } } }; @@ -57,7 +73,19 @@ struct HamiltonianOp { * To define a Hamiltonian, give a vector of controlled single bit gates ("HamiltonianOp" instances) that are * applied by left-multiplication in low-to-high vector index order on the state vector. * - * \warning Hamiltonian components might not commute, and observe the component factor of 2 * pi. + * To specify Hamiltonians with interaction terms, arbitrary sets of control bits may be specified for each term, in + * which case the term is acted only if the (superposable) control bits are true. The "antiCtrled" bool flips the + * overall convention, such that the term is acted only if all control bits are false. Additionally, for a combination + * of control bits and "anti-control" bits, an array of booleans, "ctrlToggles," of length "ctrlLen" may be specified + * that flips the activation state for each control bit, (relative the global anti- on/off convention,) without altering + * the state of the control bits. + * + * The point of this "toggle" behavior is to allow enumeration of arbitrary local Hamiltonian terms with permutations of + * a set of control bits. For example, a Hamiltonian might represent an array of local electromagnetic potential wells. + * If there are 4 wells, each with independent potentials, control "toggles" could be used on two control bits, to + * enumerate all four permutations of two control bits with four different local Hamiltonian terms. + * + * \warning Hamiltonian components might not commute. * * As a general point of linear algebra, where A and B are linear operators, e^{i * (A + B) * t} = e^{i * A * t} * * e^{i * B * t} might NOT hold, if the operators A and B do not commute. As a rule of thumb, A will commute with B diff --git a/include/qinterface.hpp b/include/qinterface.hpp index 304b7204e..ff48ce664 100644 --- a/include/qinterface.hpp +++ b/include/qinterface.hpp @@ -634,11 +634,20 @@ class QInterface { virtual void Exp(real1 radians, bitLenInt qubitIndex); /** - * Exponentiation of arbitrary 2x2 gate + * Imaginary exponentiation of arbitrary 2x2 gate * - * Applies \f$ e^{-i*Op*I} \f$, where "Op" is a 2x2 matrix, (with controls on the application of the gate). + * Applies \f$ e^{-i*Op} \f$, where "Op" is a 2x2 matrix, (with controls on the application of the gate). */ - virtual void Exp(bitLenInt* controls, bitLenInt controlLen, bitLenInt qubit, complex* matrix2x2); + virtual void Exp( + bitLenInt* controls, bitLenInt controlLen, bitLenInt qubit, complex* matrix2x2, bool antiCtrled = false); + + /** + * Logarithm of arbitrary 2x2 gate + * + * Applies \f$ log(Op) \f$, where "Op" is a 2x2 matrix, (with controls on the application of the gate). + */ + virtual void Log( + bitLenInt* controls, bitLenInt controlLen, bitLenInt qubit, complex* matrix2x2, bool antiCtrled = false); /** * Dyadic fraction (identity) exponentiation gate diff --git a/src/qinterface/operators.cpp b/src/qinterface/operators.cpp index 3754b586a..38c83308a 100644 --- a/src/qinterface/operators.cpp +++ b/src/qinterface/operators.cpp @@ -171,7 +171,23 @@ void QInterface::TimeEvolve(Hamiltonian h, real1 timeDiff) for (int j = 0; j < 4; j++) { mtrx[j] = opMtrx[j] * (-timeDiff); } - Exp(op->controls, op->controlLen, op->targetBit, mtrx); + if (op->toggles) { + for (bitLenInt j = 0; j < op->controlLen; j++) { + if (op->toggles[j]) { + X(op->controls[j]); + } + } + } + + Exp(op->controls, op->controlLen, op->targetBit, mtrx, op->anti); + + if (op->toggles) { + for (bitLenInt j = 0; j < op->controlLen; j++) { + if (op->toggles[j]) { + X(op->controls[j]); + } + } + } } } } // namespace Qrack diff --git a/src/qinterface/protected.cpp b/src/qinterface/protected.cpp index ba07cd5a3..cb5ccb632 100644 --- a/src/qinterface/protected.cpp +++ b/src/qinterface/protected.cpp @@ -10,9 +10,9 @@ // See LICENSE.md in the project root or https://www.gnu.org/licenses/lgpl-3.0.en.html // for details. +#include "common/qrack_types.hpp" #include - -#include "qinterface.hpp" +#include namespace Qrack { @@ -36,4 +36,103 @@ void rotate(BidirectionalIterator first, BidirectionalIterator middle, Bidirecti template void rotate(complex* first, complex* middle, complex* last, bitCapInt stride); +void mul2x2(complex* left, complex* right, complex* out) +{ + out[0] = (left[0] * right[0]) + (left[1] * right[2]); + out[1] = (left[0] * right[1]) + (left[1] * right[3]); + out[2] = (left[2] * right[0]) + (left[3] * right[2]); + out[3] = (left[2] * right[1]) + (left[3] * right[3]); +} + +void _expLog2x2(complex* matrix2x2, complex* outMatrix2x2, bool isExp) +{ + // Solve for the eigenvalues and eigenvectors of a 2x2 matrix, diagonalize, exponentiate, return to the original + // basis, and apply. + + // Diagonal matrices are a special case. + bool isDiag = true; + if (norm(matrix2x2[1]) > min_norm) { + isDiag = false; + } else if (norm(matrix2x2[2]) > min_norm) { + isDiag = false; + } + + complex expOfGate[4]; + complex jacobian[4]; + complex inverseJacobian[4]; + complex tempMatrix2x2[4]; + + // Diagonalize the matrix, if it is not already diagonal. Otherwise, copy it into the temporary matrix. + if (!isDiag) { + complex trace = matrix2x2[0] + matrix2x2[3]; + complex determinant = (matrix2x2[0] * matrix2x2[3]) - (matrix2x2[1] * matrix2x2[2]); + complex quadraticRoot = + sqrt((matrix2x2[0] - matrix2x2[3]) * (matrix2x2[0] - matrix2x2[3]) - (real1)(4.0) * determinant); + complex eigenvalue1 = (trace + quadraticRoot) / (real1)2.0; + complex eigenvalue2 = (trace - quadraticRoot) / (real1)2.0; + + if (norm(matrix2x2[1]) > min_norm) { + jacobian[0] = matrix2x2[1]; + jacobian[2] = eigenvalue1 - matrix2x2[0]; + + jacobian[1] = matrix2x2[1]; + jacobian[3] = eigenvalue2 - matrix2x2[0]; + } else { + jacobian[0] = eigenvalue1 - matrix2x2[3]; + jacobian[2] = matrix2x2[2]; + + jacobian[1] = eigenvalue2 - matrix2x2[3]; + jacobian[3] = matrix2x2[2]; + } + + real1 nrm = std::sqrt(norm(jacobian[0]) + norm(jacobian[2])); + jacobian[0] /= nrm; + jacobian[2] /= nrm; + + nrm = std::sqrt(norm(jacobian[1]) + norm(jacobian[3])); + jacobian[1] /= nrm; + jacobian[3] /= nrm; + + determinant = (jacobian[0] * jacobian[3]) - (jacobian[1] * jacobian[2]); + inverseJacobian[0] = jacobian[3] / determinant; + inverseJacobian[1] = -jacobian[1] / determinant; + inverseJacobian[2] = -jacobian[2] / determinant; + inverseJacobian[3] = jacobian[0] / determinant; + + mul2x2(matrix2x2, jacobian, tempMatrix2x2); + mul2x2(inverseJacobian, tempMatrix2x2, expOfGate); + } else { + std::copy(matrix2x2, matrix2x2 + 4, expOfGate); + } + + if (isExp) { + // In this branch, we calculate e^(matrix2x2). + + // Note: For a (2x2) hermitian input gate, this theoretically produces a unitary output transformation. + expOfGate[0] = ((real1)std::exp(real(expOfGate[0]))) * + complex((real1)cos(imag(expOfGate[0])), (real1)sin(imag(expOfGate[0]))); + expOfGate[1] = complex(ZERO_R1, ZERO_R1); + expOfGate[2] = complex(ZERO_R1, ZERO_R1); + expOfGate[3] = ((real1)std::exp(real(expOfGate[3]))) * + complex((real1)cos(imag(expOfGate[3])), (real1)sin(imag(expOfGate[3]))); + } else { + // In this branch, we calculate log(matrix2x2). + expOfGate[0] = complex(std::log(abs(expOfGate[0])), arg(expOfGate[0])); + expOfGate[1] = complex(ZERO_R1, ZERO_R1); + expOfGate[2] = complex(ZERO_R1, ZERO_R1); + expOfGate[3] = complex(std::log(abs(expOfGate[3])), arg(expOfGate[3])); + } + + if (!isDiag) { + mul2x2(expOfGate, inverseJacobian, tempMatrix2x2); + mul2x2(jacobian, tempMatrix2x2, expOfGate); + } + + std::copy(expOfGate, expOfGate + 4, outMatrix2x2); +} + +void exp2x2(complex* matrix2x2, complex* outMatrix2x2) { _expLog2x2(matrix2x2, outMatrix2x2, true); } + +void log2x2(complex* matrix2x2, complex* outMatrix2x2) { _expLog2x2(matrix2x2, outMatrix2x2, false); } + } // namespace Qrack diff --git a/src/qinterface/rotational.cpp b/src/qinterface/rotational.cpp index 45a2c8abb..940f22a78 100644 --- a/src/qinterface/rotational.cpp +++ b/src/qinterface/rotational.cpp @@ -64,90 +64,32 @@ void QInterface::Exp(real1 radians, bitLenInt qubit) ApplySingleBit(expIdentity, true, qubit); } -void matrix2x2Mul(complex* left, complex* right, complex* out) +/// Imaginary exponentiate of arbitrary single bit gate +void QInterface::Exp(bitLenInt* controls, bitLenInt controlLen, bitLenInt qubit, complex* matrix2x2, bool antiCtrled) { - out[0] = (left[0] * right[0]) + (left[1] * right[2]); - out[1] = (left[0] * right[1]) + (left[1] * right[3]); - out[2] = (left[2] * right[0]) + (left[3] * right[2]); - out[3] = (left[2] * right[1]) + (left[3] * right[3]); -} - -/// Exponentiate of arbitrary single bit gate -void QInterface::Exp(bitLenInt* controls, bitLenInt controlLen, bitLenInt qubit, complex* matrix2x2) -{ - // Solve for the eigenvalues and eigenvectors of a 2x2 matrix, diagonalize, exponentiate, return to the original - // basis, and apply. - - // Diagonal matrices are a special case. - bool isDiag = true; - if (norm(matrix2x2[1]) > min_norm) { - isDiag = false; - } else if (norm(matrix2x2[2]) > min_norm) { - isDiag = false; + complex timesI[4]; + complex toApply[4]; + for (bitLenInt i = 0; i < 4; i++) { + timesI[i] = complex(ZERO_R1, ONE_R1) * matrix2x2[i]; } - - complex expOfGate[4]; - complex jacobian[4]; - complex inverseJacobian[4]; - complex tempMatrix2x2[4]; - - // Diagonalize the matrix, if it is not already diagonal. Otherwise, copy it into the temporary matrix. - if (!isDiag) { - complex trace = matrix2x2[0] + matrix2x2[3]; - complex determinant = (matrix2x2[0] * matrix2x2[3]) - (matrix2x2[1] * matrix2x2[2]); - complex quadraticRoot = - sqrt((matrix2x2[0] - matrix2x2[3]) * (matrix2x2[0] - matrix2x2[3]) - (real1)(4.0) * determinant); - complex eigenvalue1 = (trace + quadraticRoot) / (real1)2.0; - complex eigenvalue2 = (trace - quadraticRoot) / (real1)2.0; - - if (norm(matrix2x2[1]) > min_norm) { - jacobian[0] = matrix2x2[1]; - jacobian[2] = eigenvalue1 - matrix2x2[0]; - - jacobian[1] = matrix2x2[1]; - jacobian[3] = eigenvalue2 - matrix2x2[0]; - } else { - jacobian[0] = eigenvalue1 - matrix2x2[3]; - jacobian[2] = matrix2x2[2]; - - jacobian[1] = eigenvalue2 - matrix2x2[3]; - jacobian[3] = matrix2x2[2]; - } - - real1 nrm = std::sqrt(norm(jacobian[0]) + norm(jacobian[2])); - jacobian[0] /= nrm; - jacobian[2] /= nrm; - - nrm = std::sqrt(norm(jacobian[1]) + norm(jacobian[3])); - jacobian[1] /= nrm; - jacobian[3] /= nrm; - - determinant = (jacobian[0] * jacobian[3]) - (jacobian[1] * jacobian[2]); - inverseJacobian[0] = jacobian[3] / determinant; - inverseJacobian[1] = -jacobian[1] / determinant; - inverseJacobian[2] = -jacobian[2] / determinant; - inverseJacobian[3] = jacobian[0] / determinant; - - matrix2x2Mul(matrix2x2, jacobian, tempMatrix2x2); - matrix2x2Mul(inverseJacobian, tempMatrix2x2, expOfGate); + Qrack::exp2x2(timesI, toApply); + if (antiCtrled) { + ApplyAntiControlledSingleBit(controls, controlLen, qubit, toApply); } else { - std::copy(matrix2x2, matrix2x2 + 4, expOfGate); + ApplyControlledSingleBit(controls, controlLen, qubit, toApply); } +} - // Note: For a (2x2) hermitian input gate, this theoretically produces a unitary output transformation. - expOfGate[0] = - ((real1)exp(-imag(expOfGate[0]))) * complex((real1)cos(real(expOfGate[0])), (real1)sin(real(expOfGate[0]))); - expOfGate[1] = complex(ZERO_R1, ZERO_R1); - expOfGate[2] = complex(ZERO_R1, ZERO_R1); - expOfGate[3] = - ((real1)exp(-imag(expOfGate[3]))) * complex((real1)cos(real(expOfGate[3])), (real1)sin(real(expOfGate[3]))); - - if (!isDiag) { - matrix2x2Mul(expOfGate, inverseJacobian, tempMatrix2x2); - matrix2x2Mul(jacobian, tempMatrix2x2, expOfGate); +/// Logarithm of arbitrary single bit gate +void QInterface::Log(bitLenInt* controls, bitLenInt controlLen, bitLenInt qubit, complex* matrix2x2, bool antiCtrled) +{ + complex toApply[4]; + Qrack::log2x2(matrix2x2, toApply); + if (antiCtrled) { + ApplyAntiControlledSingleBit(controls, controlLen, qubit, toApply); + } else { + ApplyControlledSingleBit(controls, controlLen, qubit, toApply); } - - ApplyControlledSingleBit(controls, controlLen, qubit, expOfGate); } /// Exponentiate Pauli X operator diff --git a/test/tests.cpp b/test/tests.cpp index aab6018a2..e8d1082d9 100644 --- a/test/tests.cpp +++ b/test/tests.cpp @@ -2564,9 +2564,51 @@ TEST_CASE_METHOD(QInterfaceTestFixture, "test_timeevolve") qftReg->SetPermutation(0); qftReg->TimeEvolve(h, tDiff); - // std::cout << qftReg->Prob(0) << std::endl; - // std::cout << sin(aParam * tDiff) * sin(aParam * tDiff) << std::endl; - // std::cout << abs(qftReg->Prob(0) - sin(aParam * tDiff) * sin(aParam * tDiff)) << std::endl; + REQUIRE_FLOAT(abs((ONE_R1 - qftReg->Prob(0)) - sin(aParam * tDiff) * sin(aParam * tDiff)), 0); + REQUIRE_FLOAT(abs(qftReg->Prob(0) - cos(aParam * tDiff) * cos(aParam * tDiff)), 0); + + bitLenInt controls[1] = { 1 }; + bool controlToggles[1] = { false }; + + HamiltonianOpPtr h1 = std::make_shared(controls, 1, 0, o2neg1, false, controlToggles); + h[0] = h1; + + // The point of this "toggle" behavior is to allow enumeration of arbitrary local Hamiltonian terms with + // permutations of a set of control bits. For example, a Hamiltonian might represent an array of local + // electromagnetic potential wells. If there are 4 wells, each with independent potentials, control "toggles" could + // be used on two control bits, to enumerate all four permutations of two control bits with four different local + // Hamiltonian terms. + + qftReg->SetPermutation(2); + qftReg->TimeEvolve(h, tDiff); + + REQUIRE_FLOAT(abs((ONE_R1 - qftReg->Prob(0)) - sin(aParam * tDiff) * sin(aParam * tDiff)), 0); + REQUIRE_FLOAT(abs(qftReg->Prob(0) - cos(aParam * tDiff) * cos(aParam * tDiff)), 0); + + controlToggles[0] = true; + HamiltonianOpPtr h2 = std::make_shared(controls, 1, 0, o2neg1, false, controlToggles); + h[0] = h2; + + qftReg->SetPermutation(2); + qftReg->TimeEvolve(h, tDiff); + + REQUIRE_FLOAT(qftReg->Prob(0), ZERO_R1); + + controlToggles[0] = false; + HamiltonianOpPtr h3 = std::make_shared(controls, 1, 0, o2neg1, true, controlToggles); + h[0] = h3; + + qftReg->SetPermutation(2); + qftReg->TimeEvolve(h, tDiff); + + REQUIRE_FLOAT(qftReg->Prob(0), ZERO_R1); + + controlToggles[0] = true; + HamiltonianOpPtr h4 = std::make_shared(controls, 1, 0, o2neg1, true, controlToggles); + h[0] = h4; + + qftReg->SetPermutation(2); + qftReg->TimeEvolve(h, tDiff); REQUIRE_FLOAT(abs((ONE_R1 - qftReg->Prob(0)) - sin(aParam * tDiff) * sin(aParam * tDiff)), 0); REQUIRE_FLOAT(abs(qftReg->Prob(0) - cos(aParam * tDiff) * cos(aParam * tDiff)), 0);