test_operators.cpp
Go to the documentation of this file.
1 
2 #include "catch.hpp"
3 #include "QuEST.h"
4 #include "utilities.hpp"
5 
10 #define PREPARE_TEST(quregVec, quregMatr, refVec, refMatr) \
11  Qureg quregVec = createQureg(NUM_QUBITS, QUEST_ENV); \
12  Qureg quregMatr = createDensityQureg(NUM_QUBITS, QUEST_ENV); \
13  initDebugState(quregVec); \
14  initDebugState(quregMatr); \
15  QVector refVec = toQVector(quregVec); \
16  QMatrix refMatr = toQMatrix(quregMatr);
17 
19 #define CLEANUP_TEST(quregVec, quregMatr) \
20  destroyQureg(quregVec, QUEST_ENV); \
21  destroyQureg(quregMatr, QUEST_ENV);
22 
23 /* allows concise use of Contains in catch's REQUIRE_THROWS_WITH */
24 using Catch::Matchers::Contains;
25 
26 
27 
32 TEST_CASE( "applyDiagonalOp", "[operators]" ) {
33 
34  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
35 
36  SECTION( "correctness" ) {
37 
38  // try 10 random operators
39  GENERATE( range(0,10) );
40 
41  // make a totally random (non-Hermitian) diagonal oeprator
42  DiagonalOp op = createDiagonalOp(NUM_QUBITS, QUEST_ENV);
43  for (long long int i=0; i<op.numElemsPerChunk; i++) {
44  op.real[i] = getRandomReal(-5, 5);
45  op.imag[i] = getRandomReal(-5, 5);
46  }
47  syncDiagonalOp(op);
48 
49  SECTION( "state-vector" ) {
50 
51  QVector ref = toQMatrix(op) * refVec;
52  applyDiagonalOp(quregVec, op);
53  REQUIRE( areEqual(quregVec, ref) );
54  }
55  SECTION( "density-matrix" ) {
56 
57  QMatrix ref = toQMatrix(op) * refMatr;
58  applyDiagonalOp(quregMatr, op);
59  REQUIRE( areEqual(quregMatr, ref, 100*REAL_EPS) );
60  }
61 
62  destroyDiagonalOp(op, QUEST_ENV);
63  }
64  SECTION( "input validation" ) {
65 
66  SECTION( "mismatching size" ) {
67 
68  DiagonalOp op = createDiagonalOp(NUM_QUBITS + 1, QUEST_ENV);
69 
70  REQUIRE_THROWS_WITH( applyDiagonalOp(quregVec, op), Contains("equal number of qubits"));
71  REQUIRE_THROWS_WITH( applyDiagonalOp(quregMatr, op), Contains("equal number of qubits"));
72 
73  destroyDiagonalOp(op, QUEST_ENV);
74  }
75  }
76  CLEANUP_TEST( quregVec, quregMatr );
77 }
78 
83 TEST_CASE( "applyFullQFT", "[operators]" ) {
84 
85  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
86 
87  SECTION( "correctness" ) {
88 
89  // try 10 times for each kind of input
90  GENERATE( range(0,10) );
91 
92  SECTION( "state-vector" ) {
93 
94  SECTION( "normalised" ) {
95 
96  refVec = getRandomStateVector(NUM_QUBITS);
97  toQureg(quregVec, refVec);
98  refVec = getDFT(refVec);
99 
100  applyFullQFT(quregVec);
101  REQUIRE( areEqual(quregVec, refVec) );
102  }
103  SECTION( "unnormalised" ) {
104 
105  refVec = getRandomQVector(1 << NUM_QUBITS);
106  toQureg(quregVec, refVec);
107  refVec = getDFT(refVec);
108 
109  applyFullQFT(quregVec);
110  REQUIRE( areEqual(quregVec, refVec) );
111  }
112  }
113  SECTION( "density-matrix" ) {
114 
115  SECTION( "pure" ) {
116 
117  /* a pure density matrix should be mapped to a pure state
118  * corresponding to the state-vector DFT
119  */
120 
121  refVec = getRandomStateVector(NUM_QUBITS);
122  refMatr = getPureDensityMatrix(refVec);
123 
124  toQureg(quregMatr, refMatr);
125  applyFullQFT(quregMatr);
126 
127  refVec = getDFT(refVec);
128  refMatr = getPureDensityMatrix(refVec);
129 
130  REQUIRE( areEqual(quregMatr, refMatr) );
131  }
132  SECTION( "mixed" ) {
133 
134  /* a mixed density matrix, conceptualised as a mixture of orthogonal
135  * state-vectors, should be mapped to an equally weighted mixture
136  * of DFTs of each state-vector (because QFT is unitary and hence
137  * maintains state orthogonality)
138  */
139 
140  int numStates = (1 << NUM_QUBITS)/4; // quarter as many states as possible
141  std::vector<QVector> states = getRandomOrthonormalVectors(NUM_QUBITS, numStates);
142  std::vector<qreal> probs = getRandomProbabilities(numStates);
143 
144  // set qureg to random mixture
145  refMatr = getMixedDensityMatrix(probs, states);
146  toQureg(quregMatr, refMatr);
147 
148  // apply QFT to mixture
149  applyFullQFT(quregMatr);
150 
151  // compute dft of mixture, via dft of each state
152  refMatr = getZeroMatrix(1 << NUM_QUBITS);
153  for (int i=0; i<numStates; i++)
154  refMatr += probs[i] * getPureDensityMatrix(getDFT(states[i]));
155 
156  REQUIRE( areEqual(quregMatr, refMatr) );
157  }
158  SECTION( "unnormalised" ) {
159 
160  /* repeat method above, except that we use unnormalised vectors,
161  * and mix them with arbitrary complex numbers instead of probabilities,
162  * yielding an unnormalised density matrix
163  */
164 
165  int numVecs = (1 << NUM_QUBITS)/4; // quarter as many states as possible
166  std::vector<QVector> vecs;
167  std::vector<qcomp> coeffs;
168  for (int i=0; i<numVecs; i++) {
169  vecs.push_back(getRandomQVector(1 << NUM_QUBITS));
170  coeffs.push_back(getRandomComplex());
171  }
172 
173  // produce unnormalised matrix
174  refMatr = getZeroMatrix(1 << NUM_QUBITS);
175  for (int i=0; i<numVecs; i++)
176  refMatr += coeffs[i] * getPureDensityMatrix(vecs[i]);
177 
178  toQureg(quregMatr, refMatr);
179  applyFullQFT(quregMatr);
180 
181  // compute target matrix via dft of each unnormalised vector
182  refMatr = getZeroMatrix(1 << NUM_QUBITS);
183  for (int i=0; i<numVecs; i++)
184  refMatr += coeffs[i] * getPureDensityMatrix(getDFT(vecs[i]));
185 
186  REQUIRE( areEqual(quregMatr, refMatr) );
187  }
188  }
189  }
190  SECTION( "input validation" ) {
191 
192  SUCCEED( );
193  }
194  CLEANUP_TEST( quregVec, quregMatr );
195 }
196 
197 
198 
203 TEST_CASE( "applyMatrix2", "[operators]" ) {
204 
205  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
206 
207  // every test will use a unique random matrix
208  QMatrix op = getRandomQMatrix(2); // 2-by-2
209  ComplexMatrix2 matr = toComplexMatrix2(op);
210 
211  SECTION( "correctness" ) {
212 
213  int target = GENERATE( range(0,NUM_QUBITS) );
214 
215  // reference boilerplate
216  int* ctrls = NULL;
217  int numCtrls = 0;
218  int targs[] = {target};
219  int numTargs = 1;
220 
221  SECTION( "state-vector" ) {
222 
223  applyMatrix2(quregVec, target, matr);
224  applyReferenceMatrix(refVec, ctrls, numCtrls, targs, numTargs, op);
225 
226  REQUIRE( areEqual(quregVec, refVec) );
227  }
228  SECTION( "density-matrix" ) {
229 
230  applyMatrix2(quregMatr, target, matr);
231  applyReferenceMatrix(refMatr, ctrls, numCtrls, targs, numTargs, op);
232 
233  REQUIRE( areEqual(quregMatr, refMatr, 10*REAL_EPS) );
234  }
235  }
236  SECTION( "input validation" ) {
237 
238  SECTION( "qubit indices" ) {
239 
240  int target = GENERATE( -1, NUM_QUBITS );
241  REQUIRE_THROWS_WITH( applyMatrix2(quregVec, target, matr), Contains("Invalid target") );
242  }
243  }
244  CLEANUP_TEST( quregVec, quregMatr );
245 }
246 
247 
248 
253 TEST_CASE( "applyMatrix4", "[operators]" ) {
254 
255  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
256 
257  // in distributed mode, each node must be able to fit all amps modified by matrix
258  REQUIRE( quregVec.numAmpsPerChunk >= 4 );
259 
260  // every test will use a unique random matrix
261  QMatrix op = getRandomQMatrix(4); // 4-by-4
262  ComplexMatrix4 matr = toComplexMatrix4(op);
263 
264  SECTION( "correctness" ) {
265 
266  int targ1 = GENERATE( range(0,NUM_QUBITS) );
267  int targ2 = GENERATE_COPY( filter([=](int t){ return t!=targ1; }, range(0,NUM_QUBITS)) );
268 
269  // reference boilerplate
270  int* ctrls = NULL;
271  int numCtrls = 0;
272  int targs[] = {targ1, targ2};
273  int numTargs = 2;
274 
275  SECTION( "state-vector" ) {
276 
277  applyMatrix4(quregVec, targ1, targ2, matr);
278  applyReferenceMatrix(refVec, ctrls, numCtrls, targs, numTargs, op);
279  REQUIRE( areEqual(quregVec, refVec) );
280  }
281  SECTION( "density-matrix" ) {
282 
283  applyMatrix4(quregMatr, targ1, targ2, matr);
284  applyReferenceMatrix(refMatr, ctrls, numCtrls, targs, numTargs, op);
285  REQUIRE( areEqual(quregMatr, refMatr, 10*REAL_EPS) );
286  }
287  }
288  SECTION( "input validation" ) {
289 
290  SECTION( "qubit indices" ) {
291 
292  int targ1 = GENERATE( -1, NUM_QUBITS );
293  int targ2 = 0;
294  REQUIRE_THROWS_WITH( applyMatrix4(quregVec, targ1, targ2, matr), Contains("Invalid target") );
295  REQUIRE_THROWS_WITH( applyMatrix4(quregVec, targ2, targ1, matr), Contains("Invalid target") );
296  }
297  SECTION( "repetition of targets" ) {
298 
299  int qb = 0;
300  REQUIRE_THROWS_WITH( applyMatrix4(quregVec, qb, qb, matr), Contains("target") && Contains("unique") );
301  }
302  SECTION( "matrix fits in node" ) {
303 
304  // pretend we have a very limited distributed memory
305  quregVec.numAmpsPerChunk = 1;
306  REQUIRE_THROWS_WITH( applyMatrix4(quregVec, 0, 1, matr), Contains("targets too many qubits"));
307  }
308  }
309  CLEANUP_TEST( quregVec, quregMatr );
310 }
311 
312 
313 
318 TEST_CASE( "applyMatrixN", "[operators]" ) {
319 
320  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
321 
322  // figure out max-num (inclusive) targs allowed by hardware backend
323  int maxNumTargs = calcLog2(quregVec.numAmpsPerChunk);
324 
325  SECTION( "correctness" ) {
326 
327  // generate all possible qubit arrangements
328  int numTargs = GENERATE_COPY( range(1,maxNumTargs+1) ); // inclusive upper bound
329  int* targs = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numTargs) );
330 
331  // for each qubit arrangement, use a new random matrix
332  QMatrix op = getRandomQMatrix(1 << numTargs);
333  ComplexMatrixN matr = createComplexMatrixN(numTargs);
334  toComplexMatrixN(op, matr);
335 
336  // reference boilerplate
337  int* ctrls = NULL;
338  int numCtrls = 0;
339 
340  SECTION( "state-vector" ) {
341 
342  applyMatrixN(quregVec, targs, numTargs, matr);
343  applyReferenceMatrix(refVec, ctrls, numCtrls, targs, numTargs, op);
344  REQUIRE( areEqual(quregVec, refVec) );
345  }
346  SECTION( "density-matrix" ) {
347 
348  applyMatrixN(quregMatr, targs, numTargs, matr);
349  applyReferenceMatrix(refMatr, ctrls, numCtrls, targs, numTargs, op);
350  REQUIRE( areEqual(quregMatr, refMatr, 100*REAL_EPS) );
351  }
352  destroyComplexMatrixN(matr);
353  }
354  SECTION( "input validation" ) {
355 
356  SECTION( "number of targets" ) {
357 
358  // there cannot be more targets than qubits in register
359  int numTargs = GENERATE( -1, 0, NUM_QUBITS+1 );
360  int targs[NUM_QUBITS+1]; // prevents seg-fault if validation doesn't trigger
361  ComplexMatrixN matr = createComplexMatrixN(NUM_QUBITS+1); // prevent seg-fault
362 
363  REQUIRE_THROWS_WITH( applyMatrixN(quregVec, targs, numTargs, matr), Contains("Invalid number of target"));
364  destroyComplexMatrixN(matr);
365  }
366  SECTION( "repetition in targets" ) {
367 
368  int numTargs = 3;
369  int targs[] = {1,2,2};
370  ComplexMatrixN matr = createComplexMatrixN(numTargs); // prevents seg-fault if validation doesn't trigger
371 
372  REQUIRE_THROWS_WITH( applyMatrixN(quregVec, targs, numTargs, matr), Contains("target") && Contains("unique"));
373  destroyComplexMatrixN(matr);
374  }
375  SECTION( "qubit indices" ) {
376 
377  int numTargs = 3;
378  int targs[] = {1,2,3};
379  ComplexMatrixN matr = createComplexMatrixN(numTargs); // prevents seg-fault if validation doesn't trigger
380 
381  int inv = GENERATE( -1, NUM_QUBITS );
382  targs[GENERATE_COPY( range(0,numTargs) )] = inv; // make invalid target
383  REQUIRE_THROWS_WITH( applyMatrixN(quregVec, targs, numTargs, matr), Contains("Invalid target") );
384 
385  destroyComplexMatrixN(matr);
386  }
387  SECTION( "matrix creation" ) {
388 
389  int numTargs = 3;
390  int targs[] = {1,2,3};
391 
392  /* compilers don't auto-initialise to NULL; the below circumstance
393  * only really occurs when 'malloc' returns NULL in createComplexMatrixN,
394  * which actually triggers its own validation. Hence this test is useless
395  * currently.
396  */
397  ComplexMatrixN matr;
398  matr.real = NULL;
399  matr.imag = NULL;
400  REQUIRE_THROWS_WITH( applyMatrixN(quregVec, targs, numTargs, matr), Contains("created") );
401  }
402  SECTION( "matrix dimensions" ) {
403 
404  int targs[2] = {1,2};
405  ComplexMatrixN matr = createComplexMatrixN(3); // intentionally wrong size
406 
407  REQUIRE_THROWS_WITH( applyMatrixN(quregVec, targs, 2, matr), Contains("matrix size"));
408  destroyComplexMatrixN(matr);
409  }
410  SECTION( "matrix fits in node" ) {
411 
412  // pretend we have a very limited distributed memory (judged by matr size)
413  quregVec.numAmpsPerChunk = 1;
414  int qb[] = {1,2};
415  ComplexMatrixN matr = createComplexMatrixN(2); // prevents seg-fault if validation doesn't trigger
416  REQUIRE_THROWS_WITH( applyMatrixN(quregVec, qb, 2, matr), Contains("targets too many qubits"));
417  destroyComplexMatrixN(matr);
418  }
419  }
420  CLEANUP_TEST( quregVec, quregMatr );
421 }
422 
423 
424 
429 TEST_CASE( "applyMultiControlledMatrixN", "[operators]" ) {
430 
431  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
432 
433  // figure out max-num targs (inclusive) allowed by hardware backend
434  int maxNumTargs = calcLog2(quregVec.numAmpsPerChunk);
435  if (maxNumTargs >= NUM_QUBITS)
436  maxNumTargs = NUM_QUBITS - 1; // leave room for min-number of control qubits
437 
438  SECTION( "correctness" ) {
439 
440  // try all possible numbers of targets and controls
441  int numTargs = GENERATE_COPY( range(1,maxNumTargs+1) );
442  int maxNumCtrls = NUM_QUBITS - numTargs;
443  int numCtrls = GENERATE_COPY( range(1,maxNumCtrls+1) );
444 
445  // generate all possible valid qubit arrangements
446  int* targs = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numTargs) );
447  int* ctrls = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numCtrls, targs, numTargs) );
448 
449  // for each qubit arrangement, use a new random unitary
450  QMatrix op = getRandomQMatrix(1 << numTargs);
451  ComplexMatrixN matr = createComplexMatrixN(numTargs);
452  toComplexMatrixN(op, matr);
453 
454  SECTION( "state-vector" ) {
455 
456  applyMultiControlledMatrixN(quregVec, ctrls, numCtrls, targs, numTargs, matr);
457  applyReferenceMatrix(refVec, ctrls, numCtrls, targs, numTargs, op);
458  REQUIRE( areEqual(quregVec, refVec) );
459  }
460  SECTION( "density-matrix" ) {
461 
462  applyMultiControlledMatrixN(quregMatr, ctrls, numCtrls, targs, numTargs, matr);
463  applyReferenceMatrix(refMatr, ctrls, numCtrls, targs, numTargs, op);
464  REQUIRE( areEqual(quregMatr, refMatr, 100*REAL_EPS) );
465  }
466  destroyComplexMatrixN(matr);
467  }
468  SECTION( "input validation" ) {
469 
470  SECTION( "number of targets" ) {
471 
472  // there cannot be more targets than qubits in register
473  // (numTargs=NUM_QUBITS is caught elsewhere, because that implies ctrls are invalid)
474  int numTargs = GENERATE( -1, 0, NUM_QUBITS+1 );
475  int targs[NUM_QUBITS+1]; // prevents seg-fault if validation doesn't trigger
476  int ctrls[] = {0};
477  ComplexMatrixN matr = createComplexMatrixN(NUM_QUBITS+1); // prevent seg-fault
478  toComplexMatrixN(getRandomQMatrix( 1 << (NUM_QUBITS+1)), matr);
479 
480  REQUIRE_THROWS_WITH( applyMultiControlledMatrixN(quregVec, ctrls, 1, targs, numTargs, matr), Contains("Invalid number of target"));
481  destroyComplexMatrixN(matr);
482  }
483  SECTION( "repetition in targets" ) {
484 
485  int ctrls[] = {0};
486  int numTargs = 3;
487  int targs[] = {1,2,2};
488  ComplexMatrixN matr = createComplexMatrixN(numTargs); // prevents seg-fault if validation doesn't trigger
489  toComplexMatrixN(getRandomQMatrix(1 << numTargs), matr);
490 
491  REQUIRE_THROWS_WITH( applyMultiControlledMatrixN(quregVec, ctrls, 1, targs, numTargs, matr), Contains("target") && Contains("unique"));
492  destroyComplexMatrixN(matr);
493  }
494  SECTION( "number of controls" ) {
495 
496  int numCtrls = GENERATE( -1, 0, NUM_QUBITS, NUM_QUBITS+1 );
497  int ctrls[NUM_QUBITS+1]; // avoids seg-fault if validation not triggered
498  int targs[1] = {0};
499  ComplexMatrixN matr = createComplexMatrixN(1);
500  toComplexMatrixN(getRandomQMatrix(1 << 1), matr);
501 
502  REQUIRE_THROWS_WITH( applyMultiControlledMatrixN(quregVec, ctrls, numCtrls, targs, 1, matr), Contains("Invalid number of control"));
503  destroyComplexMatrixN(matr);
504  }
505  SECTION( "repetition in controls" ) {
506 
507  int ctrls[] = {0,1,1};
508  int targs[] = {3};
509  ComplexMatrixN matr = createComplexMatrixN(1);
510  toComplexMatrixN(getRandomQMatrix(1 << 1), matr);
511 
512  REQUIRE_THROWS_WITH( applyMultiControlledMatrixN(quregVec, ctrls, 3, targs, 1, matr), Contains("control") && Contains("unique"));
513  destroyComplexMatrixN(matr);
514  }
515  SECTION( "control and target collision" ) {
516 
517  int ctrls[] = {0,1,2};
518  int targs[] = {3,1,4};
519  ComplexMatrixN matr = createComplexMatrixN(3);
520  toComplexMatrixN(getRandomQMatrix(1 << 3), matr);
521 
522  REQUIRE_THROWS_WITH( applyMultiControlledMatrixN(quregVec, ctrls, 3, targs, 3, matr), Contains("Control") && Contains("target") && Contains("disjoint"));
523  destroyComplexMatrixN(matr);
524  }
525  SECTION( "qubit indices" ) {
526 
527  // valid inds
528  int numQb = 2;
529  int qb1[2] = {0,1};
530  int qb2[2] = {2,3};
531  ComplexMatrixN matr = createComplexMatrixN(numQb);
532  toComplexMatrixN(getRandomQMatrix(1 << numQb), matr);
533 
534  // make qb1 invalid
535  int inv = GENERATE( -1, NUM_QUBITS );
536  qb1[GENERATE_COPY(range(0,numQb))] = inv;
537 
538  REQUIRE_THROWS_WITH( applyMultiControlledMatrixN(quregVec, qb1, numQb, qb2, numQb, matr), Contains("Invalid control") );
539  REQUIRE_THROWS_WITH( applyMultiControlledMatrixN(quregVec, qb2, numQb, qb1, numQb, matr), Contains("Invalid target") );
540  destroyComplexMatrixN(matr);
541  }
542  SECTION( "matrix creation" ) {
543 
544  int ctrls[1] = {0};
545  int targs[3] = {1,2,3};
546 
547  /* compilers don't auto-initialise to NULL; the below circumstance
548  * only really occurs when 'malloc' returns NULL in createComplexMatrixN,
549  * which actually triggers its own validation. Hence this test is useless
550  * currently.
551  */
552  ComplexMatrixN matr;
553  matr.real = NULL;
554  matr.imag = NULL;
555  REQUIRE_THROWS_WITH( applyMultiControlledMatrixN(quregVec, ctrls, 1, targs, 3, matr), Contains("created") );
556  }
557  SECTION( "matrix dimensions" ) {
558 
559  int ctrls[1] = {0};
560  int targs[2] = {1,2};
561  ComplexMatrixN matr = createComplexMatrixN(3); // intentionally wrong size
562  toComplexMatrixN(getRandomQMatrix(1 << 3), matr);
563 
564  REQUIRE_THROWS_WITH( applyMultiControlledMatrixN(quregVec, ctrls, 1, targs, 2, matr), Contains("matrix size"));
565  destroyComplexMatrixN(matr);
566  }
567  SECTION( "matrix fits in node" ) {
568 
569  // pretend we have a very limited distributed memory (judged by matr size)
570  quregVec.numAmpsPerChunk = 1;
571  int ctrls[1] = {0};
572  int targs[2] = {1,2};
573  ComplexMatrixN matr = createComplexMatrixN(2);
574  toComplexMatrixN(getRandomQMatrix(1 << 2), matr);
575 
576  REQUIRE_THROWS_WITH( applyMultiControlledMatrixN(quregVec, ctrls, 1, targs, 2, matr), Contains("targets too many qubits"));
577  destroyComplexMatrixN(matr);
578  }
579  }
580  CLEANUP_TEST( quregVec, quregMatr );
581 }
582 
583 
584 
589 TEST_CASE( "applyMultiVarPhaseFunc", "[operators]" ) {
590 
591  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
592 
593  SECTION( "correctness" ) {
594 
595  // try every kind of binary encodings
596  enum bitEncoding encoding = GENERATE( UNSIGNED,TWOS_COMPLEMENT );
597 
598  // try every possible number of registers
599  // (between #qubits containing 1, and 1 containing #qubits)
600  int numRegs;
601  int maxNumRegs = 0;
602  if (encoding == UNSIGNED)
603  maxNumRegs = NUM_QUBITS;
604  if (encoding == TWOS_COMPLEMENT)
605  maxNumRegs = NUM_QUBITS/2; // floors
606  numRegs = GENERATE_COPY( range(1, maxNumRegs+1) );
607 
608  // try every possible total number of involed qubits
609  int totalNumQubits;
610  int minTotalQubits = 0;
611  if (encoding == UNSIGNED)
612  // each register must contain at least 1 qubit
613  minTotalQubits = numRegs;
614  if (encoding == TWOS_COMPLEMENT)
615  // each register must contain at least 2 qubits
616  minTotalQubits = 2*numRegs;
617  totalNumQubits = GENERATE_COPY( range(minTotalQubits,NUM_QUBITS+1) );
618 
619  // try every qubits subset and ordering
620  int* regs = GENERATE_COPY( sublists(range(0,NUM_QUBITS), totalNumQubits) );
621 
622  // assign each sub-reg its minimum length
623  int unallocQubits = totalNumQubits;
624  int numQubitsPerReg[numRegs];
625  for (int i=0; i<numRegs; i++)
626  if (encoding == UNSIGNED) {
627  numQubitsPerReg[i] = 1;
628  unallocQubits -= 1;
629  }
630  else if (encoding == TWOS_COMPLEMENT) {
631  numQubitsPerReg[i] = 2;
632  unallocQubits -= 2;
633  }
634  // and randomly allocate the remaining qubits between the registers
635  while (unallocQubits > 0) {
636  numQubitsPerReg[getRandomInt(0,numRegs)] += 1;
637  unallocQubits--;
638  }
639 
640 
641  // each register gets a random number of terms in the full phase function
642  int numTermsPerReg[numRegs];
643  int numTermsTotal = 0;
644  for (int r=0; r<numRegs; r++) {
645  int numTerms = getRandomInt(1,5);
646  numTermsPerReg[r] = numTerms;
647  numTermsTotal += numTerms;
648  }
649 
650  // populate the multi-var phase function with random but POSITIVE-power terms,
651  // which must further be integers in two's complement
652  int flatInd = 0;
653  qreal coeffs[numTermsTotal];
654  qreal expons[numTermsTotal];
655  for (int r=0; r<numRegs; r++) {
656  for (int t=0; t<numTermsPerReg[r]; t++) {
657  coeffs[flatInd] = getRandomReal(-10,10);
658  if (encoding == TWOS_COMPLEMENT)
659  expons[flatInd] = getRandomInt(0, 3+1);
660  else if (encoding == UNSIGNED)
661  expons[flatInd] = getRandomReal(0, 3);
662 
663  flatInd++;
664  }
665  }
666 
667  /* To perform this calculation more distinctly from the QuEST
668  * core method, we can exploit that
669  * exp(i (f1[x1] + f2[x2]))|x2>|x1> = exp(i f2[x2])|x2> (x) exp(i f1[x1])|x1>
670  * and so compute a separate diagonal matrix for each sub-register with
671  * phases determined only by its own variable, and Kronecker multiply
672  * them together. The target qubits of the Kronecker product is then
673  * just the full register, stored in *regs.
674  */
675  QMatrix allRegMatr{{1}};
676  int startInd = 0;
677 
678  for (int r=0; r<numRegs; r++) {
679 
680  QMatrix singleRegMatr = getZeroMatrix( 1 << numQubitsPerReg[r] );
681 
682  for (size_t i=0; i<singleRegMatr.size(); i++) {
683 
684  long long int ind = 0;
685  if (encoding == UNSIGNED)
686  ind = i;
687  if (encoding == TWOS_COMPLEMENT)
688  ind = getTwosComplement(i, numQubitsPerReg[r]);
689 
690  qreal phase = 0;
691  for (int t=0; t<numTermsPerReg[r]; t++)
692  phase += coeffs[t+startInd] * pow(ind, expons[t+startInd]);
693 
694  singleRegMatr[i][i] = expI(phase);
695  }
696  allRegMatr = getKroneckerProduct(singleRegMatr, allRegMatr);
697  startInd += numTermsPerReg[r];
698  }
699 
700  SECTION( "state-vector" ) {
701 
702  applyMultiVarPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, coeffs, expons, numTermsPerReg);
703  applyReferenceOp(refVec, regs, totalNumQubits, allRegMatr);
704  REQUIRE( areEqual(quregVec, refVec, 1E4*REAL_EPS) );
705  }
706  SECTION( "density-matrix" ) {
707 
708  applyMultiVarPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, coeffs, expons, numTermsPerReg);
709  applyReferenceOp(refMatr, regs, totalNumQubits, allRegMatr);
710  REQUIRE( areEqual(quregMatr, refMatr, 1E6*REAL_EPS) );
711  }
712  }
713  SECTION( "input validation" ) {
714 
715  int numRegs = 2;
716  int numQubitsPerReg[] = {2,3};
717  int qubits[] = {0,1,2,3,4};
718 
719  SECTION( "number of registers" ) {
720 
721  numRegs = GENERATE_COPY( -1, 0, 1+MAX_NUM_REGS_APPLY_ARBITRARY_PHASE );
722  REQUIRE_THROWS_WITH( applyMultiVarPhaseFunc(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, NULL, NULL, NULL), Contains("Invalid number of qubit subregisters") );
723  }
724  SECTION( "number of qubits" ) {
725 
726  numQubitsPerReg[GENERATE_COPY(range(0,numRegs))] = GENERATE( -1, 0, 1+NUM_QUBITS );
727  REQUIRE_THROWS_WITH( applyMultiVarPhaseFunc(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, NULL, NULL, NULL), Contains("Invalid number of qubits") );
728  }
729  SECTION( "repetition of qubits" ) {
730 
731  qubits[GENERATE(2,3,4)] = qubits[1];
732  REQUIRE_THROWS_WITH( applyMultiVarPhaseFunc(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, NULL, NULL, NULL), Contains("The qubits must be unique") );
733  }
734  SECTION( "qubit indices" ) {
735 
736  qubits[GENERATE(range(0,NUM_QUBITS))] = GENERATE( -1, NUM_QUBITS );
737  REQUIRE_THROWS_WITH( applyMultiVarPhaseFunc(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, NULL, NULL, NULL), Contains("Invalid qubit index") );
738  }
739  SECTION( "number of terms" ) {
740 
741  int numTermsPerReg[] = {3, 3};
742 
743  numTermsPerReg[GENERATE_COPY(range(0,numRegs))] = GENERATE( -1, 0 );
744  REQUIRE_THROWS_WITH( applyMultiVarPhaseFunc(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, NULL, NULL, numTermsPerReg), Contains("Invalid number of terms in the phase function") );
745  }
746  SECTION( "bit encoding name" ) {
747 
748  enum bitEncoding enc = (enum bitEncoding) GENERATE( -1, 2 );
749  REQUIRE_THROWS_WITH( applyMultiVarPhaseFunc(quregVec, qubits, numQubitsPerReg, numRegs, enc, NULL, NULL, NULL), Contains("Invalid bit encoding") );
750  }
751  SECTION( "two's complement register" ) {
752 
753  numQubitsPerReg[GENERATE_COPY(range(0,numRegs))] = 1;
754  REQUIRE_THROWS_WITH( applyMultiVarPhaseFunc(quregVec, qubits, numQubitsPerReg, numRegs, TWOS_COMPLEMENT, NULL, NULL, NULL), Contains("A sub-register contained too few qubits to employ TWOS_COMPLEMENT encoding") );
755  }
756  SECTION( "fractional exponent" ) {
757 
758  int numTermsPerReg[] = {3, 3};
759  qreal coeffs[] = {0,0,0, 0,0,0};
760  qreal expos[] = {1,2,3, 1,2,3};
761 
762  expos[GENERATE(range(0,6))] = GENERATE( 0.5, 1.999, 5.0001 );
763  REQUIRE_THROWS_WITH( applyMultiVarPhaseFunc(quregVec, qubits, numQubitsPerReg, numRegs, TWOS_COMPLEMENT, coeffs, expos, numTermsPerReg), Contains("The phase function contained a fractional exponent, which is illegal in TWOS_COMPLEMENT") );
764 
765  // ensure fractional exponents are valid in unsigned mode however
766  REQUIRE_NOTHROW( applyMultiVarPhaseFunc(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, coeffs, expos, numTermsPerReg) );
767 
768  }
769  SECTION( "negative exponent" ) {
770 
771  int numTermsPerReg[] = {3, 3};
772  qreal coeffs[] = {0,0,0, 0,0,0};
773  qreal expos[] = {1,2,3, 1,2,3};
774 
775  expos[GENERATE(range(0,6))] = GENERATE( -1, -2, -2.5 );
776  enum bitEncoding enc = GENERATE( UNSIGNED, TWOS_COMPLEMENT );
777 
778  REQUIRE_THROWS_WITH( applyMultiVarPhaseFunc(quregVec, qubits, numQubitsPerReg, numRegs, enc, coeffs, expos, numTermsPerReg), Contains("The phase function contained an illegal negative exponent") );
779  }
780  }
781  CLEANUP_TEST( quregVec, quregMatr );
782 }
783 
784 
785 
790 TEST_CASE( "applyMultiVarPhaseFuncOverrides", "[operators]" ) {
791 
792  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
793 
794  SECTION( "correctness" ) {
795 
796  // try every kind of binary encodings
797  enum bitEncoding encoding = GENERATE( UNSIGNED,TWOS_COMPLEMENT );
798 
799  // try every possible number of registers
800  // (between #qubits containing 1, and 1 containing #qubits)
801  int numRegs;
802  int maxNumRegs = 0;
803  if (encoding == UNSIGNED)
804  maxNumRegs = NUM_QUBITS;
805  if (encoding == TWOS_COMPLEMENT)
806  maxNumRegs = NUM_QUBITS/2; // floors
807  numRegs = GENERATE_COPY( range(1, maxNumRegs+1) );
808 
809  // try every possible total number of involed qubits
810  int totalNumQubits;
811  int minTotalQubits = 0;
812  if (encoding == UNSIGNED)
813  // each register must contain at least 1 qubit
814  minTotalQubits = numRegs;
815  if (encoding == TWOS_COMPLEMENT)
816  // each register must contain at least 2 qubits
817  minTotalQubits = 2*numRegs;
818  totalNumQubits = GENERATE_COPY( range(minTotalQubits,NUM_QUBITS+1) );
819 
820  // try every qubits subset and ordering
821  int* regs = GENERATE_COPY( sublists(range(0,NUM_QUBITS), totalNumQubits) );
822 
823  // assign each sub-reg its minimum length
824  int unallocQubits = totalNumQubits;
825  int numQubitsPerReg[numRegs];
826  for (int i=0; i<numRegs; i++)
827  if (encoding == UNSIGNED) {
828  numQubitsPerReg[i] = 1;
829  unallocQubits -= 1;
830  }
831  else if (encoding == TWOS_COMPLEMENT) {
832  numQubitsPerReg[i] = 2;
833  unallocQubits -= 2;
834  }
835  // and randomly allocate the remaining qubits between the registers
836  while (unallocQubits > 0) {
837  numQubitsPerReg[getRandomInt(0,numRegs)] += 1;
838  unallocQubits--;
839  }
840 
841 
842  // each register gets a random number of terms in the full phase function
843  int numTermsPerReg[numRegs];
844  int numTermsTotal = 0;
845  for (int r=0; r<numRegs; r++) {
846  int numTerms = getRandomInt(1,5);
847  numTermsPerReg[r] = numTerms;
848  numTermsTotal += numTerms;
849  }
850 
851  // populate the multi-var phase function with random but POSITIVE-power terms,
852  // which must further be integers in two's complement
853  int flatInd = 0;
854  qreal coeffs[numTermsTotal];
855  qreal expons[numTermsTotal];
856  for (int r=0; r<numRegs; r++) {
857  for (int t=0; t<numTermsPerReg[r]; t++) {
858  coeffs[flatInd] = getRandomReal(-10,10);
859  if (encoding == TWOS_COMPLEMENT)
860  expons[flatInd] = getRandomInt(0, 3+1);
861  else if (encoding == UNSIGNED)
862  expons[flatInd] = getRandomReal(0, 3);
863 
864  flatInd++;
865  }
866  }
867 
868 
869  // choose a random number of overrides (even overriding every amplitude)
870  int numOverrides = getRandomInt(0, (1<<totalNumQubits) + 1);
871 
872  // randomise each override index (uniqueness isn't checked)
873  long long int overrideInds[numOverrides*numRegs];
874  flatInd = 0;
875  for (int v=0; v<numOverrides; v++) {
876  for (int r=0; r<numRegs; r++) {
877  if (encoding == UNSIGNED)
878  overrideInds[flatInd] = getRandomInt(0, 1<<numQubitsPerReg[r]);
879  else if (encoding == TWOS_COMPLEMENT)
880  overrideInds[flatInd] = getRandomInt(-(1<<(numQubitsPerReg[r]-1)), (1<<(numQubitsPerReg[r]-1))-1);
881  flatInd++;
882  }
883  }
884 
885  // override to a random phase
886  qreal overridePhases[numOverrides];
887  for (int v=0; v<numOverrides; v++)
888  overridePhases[v] = getRandomReal(-4, 4); // periodic in [-pi, pi]
889 
890  /* To perform this calculation more distinctly from the QuEST
891  * core method, we can exploit that
892  * exp(i (f1[x1] + f2[x2]))|x2>|x1> = exp(i f2[x2])|x2> (x) exp(i f1[x1])|x1>
893  * and so compute a separate diagonal matrix for each sub-register with
894  * phases determined only by its own variable, and Kronecker multiply
895  * them together. The target qubits of the Kronecker product is then
896  * just the full register, stored in *regs. We do this, then iterate
897  * the list of overrides directly to overwrite the ultimate diagonal
898  * entries
899  */
900  QMatrix allRegMatr{{1}};
901  int startInd = 0;
902  for (int r=0; r<numRegs; r++) {
903 
904  QMatrix singleRegMatr = getZeroMatrix( 1 << numQubitsPerReg[r] );
905 
906  for (size_t i=0; i<singleRegMatr.size(); i++) {
907 
908  long long int ind = 0;
909  if (encoding == UNSIGNED)
910  ind = i;
911  if (encoding == TWOS_COMPLEMENT)
912  ind = getTwosComplement(i, numQubitsPerReg[r]);
913 
914  qreal phase = 0;
915  for (int t=0; t<numTermsPerReg[r]; t++)
916  phase += coeffs[t+startInd] * pow(ind, expons[t+startInd]);
917 
918  singleRegMatr[i][i] = expI(phase);
919  }
920  allRegMatr = getKroneckerProduct(singleRegMatr, allRegMatr);
921  startInd += numTermsPerReg[r];
922  }
923  setDiagMatrixOverrides(allRegMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
924 
925  SECTION( "state-vector" ) {
926 
927  applyMultiVarPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, coeffs, expons, numTermsPerReg, overrideInds, overridePhases, numOverrides);
928  applyReferenceOp(refVec, regs, totalNumQubits, allRegMatr);
929  REQUIRE( areEqual(quregVec, refVec, 1E4*REAL_EPS) );
930  }
931  SECTION( "density-matrix" ) {
932 
933  applyMultiVarPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, coeffs, expons, numTermsPerReg, overrideInds, overridePhases, numOverrides);
934  applyReferenceOp(refMatr, regs, totalNumQubits, allRegMatr);
935  REQUIRE( areEqual(quregMatr, refMatr, 1E6*REAL_EPS) );
936  }
937  }
938  SECTION( "input validation" ) {
939 
940  int numRegs = 2;
941  int numQubitsPerReg[] = {2,3};
942  int qubits[] = {0,1,2,3,4};
943 
944  SECTION( "number of registers" ) {
945 
946  numRegs = GENERATE_COPY( -1, 0, 1+MAX_NUM_REGS_APPLY_ARBITRARY_PHASE );
947  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, NULL, NULL, NULL, NULL, NULL, 0), Contains("Invalid number of qubit subregisters") );
948  }
949  SECTION( "number of qubits" ) {
950 
951  numQubitsPerReg[GENERATE_COPY(range(0,numRegs))] = GENERATE( -1, 0, 1+NUM_QUBITS );
952  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, NULL, NULL, NULL, NULL, NULL, 0), Contains("Invalid number of qubits") );
953  }
954  SECTION( "repetition of qubits" ) {
955 
956  qubits[GENERATE(2,3,4)] = qubits[1];
957  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, NULL, NULL, NULL, NULL, NULL, 0), Contains("The qubits must be unique") );
958  }
959  SECTION( "qubit indices" ) {
960 
961  qubits[GENERATE(range(0,NUM_QUBITS))] = GENERATE( -1, NUM_QUBITS );
962  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, NULL, NULL, NULL, NULL, NULL, 0), Contains("Invalid qubit index") );
963  }
964  SECTION( "number of terms" ) {
965 
966  int numTermsPerReg[] = {3, 3};
967 
968  numTermsPerReg[GENERATE_COPY(range(0,numRegs))] = GENERATE( -1, 0 );
969  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, NULL, NULL, numTermsPerReg, NULL, NULL, 0), Contains("Invalid number of terms in the phase function") );
970  }
971  SECTION( "bit encoding name" ) {
972 
973  enum bitEncoding enc = (enum bitEncoding) GENERATE( -1, 2 );
974  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, enc, NULL, NULL, NULL, NULL, NULL, 0), Contains("Invalid bit encoding") );
975  }
976  SECTION( "two's complement register" ) {
977 
978  numQubitsPerReg[GENERATE_COPY(range(0,numRegs))] = 1;
979  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, TWOS_COMPLEMENT, NULL, NULL, NULL, NULL, NULL, 0), Contains("A sub-register contained too few qubits to employ TWOS_COMPLEMENT encoding") );
980  }
981  SECTION( "fractional exponent" ) {
982 
983  int numTermsPerReg[] = {3, 3};
984  qreal coeffs[] = {0,0,0, 0,0,0};
985  qreal expos[] = {1,2,3, 1,2,3};
986 
987  expos[GENERATE(range(0,6))] = GENERATE( 0.5, 1.999, 5.0001 );
988  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, TWOS_COMPLEMENT, coeffs, expos, numTermsPerReg, NULL, NULL, 0), Contains("The phase function contained a fractional exponent, which is illegal in TWOS_COMPLEMENT") );
989 
990  // ensure fractional exponents are valid in unsigned mode however
991  REQUIRE_NOTHROW( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, coeffs, expos, numTermsPerReg, NULL, NULL, 0) );
992  }
993  SECTION( "negative exponent" ) {
994 
995  int numTermsPerReg[] = {3, 3};
996  qreal coeffs[] = {0,0,0, 0,0,0};
997  qreal expos[] = {1,2,3, 1,2,3};
998 
999  expos[GENERATE(range(0,6))] = GENERATE( -1, -2, -2.5 );
1000  enum bitEncoding enc = GENERATE( UNSIGNED, TWOS_COMPLEMENT );
1001 
1002  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, enc, coeffs, expos, numTermsPerReg, NULL, NULL, 0), Contains("The phase function contained an illegal negative exponent") );
1003  }
1004  SECTION( "number of overrides" ) {
1005 
1006  int numTermsPerReg[] = {3, 3};
1007  qreal coeffs[] = {0,0,0, 0,0,0};
1008  qreal expos[] = {1,2,3, 1,2,3};
1009 
1010  int numOverrides = -1;
1011  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, UNSIGNED, coeffs, expos, numTermsPerReg, NULL, NULL, numOverrides), Contains("Invalid number of phase function overrides specified") );
1012  }
1013  SECTION( "override indices" ) {
1014 
1015  int numTermsPerReg[] = {3, 3};
1016  qreal coeffs[] = {0,0,0, 0,0,0};
1017  qreal expos[] = {1,2,3, 1,2,3};
1018 
1019  // numQubitsPerReg = {2, 3}
1020  int numOverrides = 3;
1021  long long int overrideInds[] = {0,0, 0,0, 0,0}; // repetition not checked
1022  qreal overridePhases[] = {.1, .1, .1};
1023 
1024  // first element of overrideInds coordinate is a 2 qubit register
1025  enum bitEncoding enc = UNSIGNED;
1026  int badInd = GENERATE(0, 2, 4);
1027  overrideInds[badInd] = GENERATE( -1, (1<<2) );
1028  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, enc, coeffs, expos, numTermsPerReg, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the UNSIGNED encoding") );
1029  overrideInds[badInd] = 0;
1030 
1031  // second element of overrideInds coordinate is a 3 qubit register
1032  badInd += 1;
1033  overrideInds[badInd] = GENERATE( -1, (1<<3) );
1034  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, enc, coeffs, expos, numTermsPerReg, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the UNSIGNED encoding") );
1035  overrideInds[badInd] = 0;
1036  badInd -= 1;
1037 
1038  enc = TWOS_COMPLEMENT;
1039  int minInd = -(1<<(numQubitsPerReg[0]-1));
1040  int maxInd = (1<<(numQubitsPerReg[0]-1)) - 1;
1041  overrideInds[badInd] = GENERATE_COPY( minInd-1, maxInd+1 );
1042  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, enc, coeffs, expos, numTermsPerReg, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the TWOS_COMPLEMENT encoding") );
1043  overrideInds[badInd] = 0;
1044 
1045  badInd++;
1046  minInd = -(1<<(numQubitsPerReg[1]-1));
1047  maxInd = (1<<(numQubitsPerReg[1]-1)) -1;
1048  overrideInds[badInd] = GENERATE_COPY( minInd-1, maxInd+1 );
1049  REQUIRE_THROWS_WITH( applyMultiVarPhaseFuncOverrides(quregVec, qubits, numQubitsPerReg, numRegs, enc, coeffs, expos, numTermsPerReg, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the TWOS_COMPLEMENT encoding") );
1050  }
1051  }
1052  CLEANUP_TEST( quregVec, quregMatr );
1053 }
1054 
1055 
1056 
1061 TEST_CASE( "applyNamedPhaseFunc", "[operators]" ) {
1062 
1063  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
1064 
1065  SECTION( "correctness" ) {
1066 
1067  // try every kind of binary encoding
1068  enum bitEncoding encoding = GENERATE( UNSIGNED,TWOS_COMPLEMENT );
1069 
1070  // try every possible number of registers
1071  // (between #qubits containing 1, and 1 containing #qubits)
1072  int numRegs;
1073  int maxNumRegs = 0;
1074  if (encoding == UNSIGNED)
1075  maxNumRegs = NUM_QUBITS;
1076  if (encoding == TWOS_COMPLEMENT)
1077  maxNumRegs = NUM_QUBITS/2; // floors
1078  numRegs = GENERATE_COPY( range(1, maxNumRegs+1) );
1079 
1080  // try every possible total number of involved qubits
1081  int totalNumQubits;
1082  int minTotalQubits = 0;
1083  if (encoding == UNSIGNED)
1084  // each register must contain at least 1 qubit
1085  minTotalQubits = numRegs;
1086  if (encoding == TWOS_COMPLEMENT)
1087  // each register must contain at least 2 qubits
1088  minTotalQubits = 2*numRegs;
1089  totalNumQubits = GENERATE_COPY( range(minTotalQubits,NUM_QUBITS+1) );
1090 
1091  // try every qubits subset and ordering
1092  int* regs = GENERATE_COPY( sublists(range(0,NUM_QUBITS), totalNumQubits) );
1093 
1094  // assign each sub-reg its minimum length
1095  int unallocQubits = totalNumQubits;
1096  int numQubitsPerReg[numRegs];
1097  for (int i=0; i<numRegs; i++)
1098  if (encoding == UNSIGNED) {
1099  numQubitsPerReg[i] = 1;
1100  unallocQubits -= 1;
1101  }
1102  else if (encoding == TWOS_COMPLEMENT) {
1103  numQubitsPerReg[i] = 2;
1104  unallocQubits -= 2;
1105  }
1106  // and randomly allocate the remaining qubits between the registers
1107  while (unallocQubits > 0) {
1108  numQubitsPerReg[getRandomInt(0,numRegs)] += 1;
1109  unallocQubits--;
1110  }
1111 
1112  // for reference, determine the values corresponding to each register for all basis states
1113  qreal regVals[1<<totalNumQubits][numRegs];
1114  for (long long int i=0; i<(1<<totalNumQubits); i++) {
1115 
1116  long long int bits = i;
1117  for (int r=0; r<numRegs; r++) {
1118  regVals[i][r] = bits % (1 << numQubitsPerReg[r]);
1119  bits = bits >> numQubitsPerReg[r];
1120 
1121  if (encoding == TWOS_COMPLEMENT)
1122  regVals[i][r] = getTwosComplement(regVals[i][r], numQubitsPerReg[r]);
1123  }
1124  }
1125 
1126  /* the reference diagonal matrix which assumes the qubits are
1127  * contiguous and strictly increasing between the registers, and hence
1128  * only depends on the number of qubits in each register.
1129  */
1130  QMatrix diagMatr = getZeroMatrix(1 << totalNumQubits);
1131 
1132  SECTION( "NORM" ) {
1133 
1134  for (size_t i=0; i<diagMatr.size(); i++) {
1135  qreal phase = 0;
1136  for (int r=0; r<numRegs; r++)
1137  phase += pow(regVals[i][r], 2);
1138  phase = sqrt(phase);
1139  diagMatr[i][i] = expI(phase);
1140  }
1141 
1142  SECTION( "state-vector" ) {
1143 
1144  applyNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, NORM);
1145  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1146  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1147  }
1148  SECTION( "density-matrix" ) {
1149 
1150  applyNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, NORM);
1151  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1152  REQUIRE( areEqual(quregMatr, refMatr, 142*REAL_EPS) );
1153  }
1154  }
1155  SECTION( "PRODUCT" ) {
1156 
1157  for (size_t i=0; i<diagMatr.size(); i++) {
1158  qreal phase = 1;
1159  for (int r=0; r<numRegs; r++)
1160  phase *= regVals[i][r];
1161  diagMatr[i][i] = expI(phase);
1162  }
1163 
1164  SECTION( "state-vector" ) {
1165 
1166  applyNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, PRODUCT);
1167  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1168  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1169  }
1170  SECTION( "density-matrix" ) {
1171 
1172  applyNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, PRODUCT);
1173  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1174  REQUIRE( areEqual(quregMatr, refMatr, 142*REAL_EPS) );
1175  }
1176  }
1177  SECTION( "DISTANCE" ) {
1178 
1179  // test only if there are an even number of registers
1180  if (numRegs%2 == 0) {
1181 
1182  for (size_t i=0; i<diagMatr.size(); i++) {
1183  qreal phase = 0;
1184  for (int r=0; r<numRegs; r+=2)
1185  phase += pow(regVals[i][r+1]-regVals[i][r], 2);
1186  phase = sqrt(phase);
1187  diagMatr[i][i] = expI(phase);
1188  }
1189  }
1190 
1191  SECTION( "state-vector" ) {
1192 
1193  if (numRegs%2 == 0) {
1194  applyNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, DISTANCE);
1195  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1196  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1197  }
1198  }
1199  SECTION( "density-matrix" ) {
1200 
1201  if (numRegs%2 == 0) {
1202  applyNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, DISTANCE);
1203  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1204  REQUIRE( areEqual(quregMatr, refMatr, 142*REAL_EPS) );
1205  }
1206  }
1207  }
1208  }
1209  SECTION( "input validation" ) {
1210 
1211  int numRegs = 2;
1212  int numQubitsPerReg[] = {2,3};
1213  int regs[] = {0,1,2,3,4};
1214 
1215  SECTION( "number of registers" ) {
1216 
1217  numRegs = GENERATE_COPY( -1, 0, 1+MAX_NUM_REGS_APPLY_ARBITRARY_PHASE );
1218  REQUIRE_THROWS_WITH( applyNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM), Contains("Invalid number of qubit subregisters") );
1219  }
1220  SECTION( "number of qubits" ) {
1221 
1222  numQubitsPerReg[GENERATE_COPY(range(0,numRegs))] = GENERATE( -1, 0, 1+NUM_QUBITS );
1223  REQUIRE_THROWS_WITH( applyNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM), Contains("Invalid number of qubits") );
1224  }
1225  SECTION( "repetition of qubits" ) {
1226 
1227  regs[GENERATE(2,3,4)] = regs[1];
1228  REQUIRE_THROWS_WITH( applyNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM), Contains("The qubits must be unique") );
1229  }
1230  SECTION( "qubit indices" ) {
1231 
1232  regs[GENERATE(range(0,NUM_QUBITS))] = GENERATE( -1, NUM_QUBITS );
1233  REQUIRE_THROWS_WITH( applyNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM), Contains("Invalid qubit index") );
1234  }
1235  SECTION( "bit encoding name" ) {
1236 
1237  enum bitEncoding enc = (enum bitEncoding) GENERATE( -1, 2 );
1238  REQUIRE_THROWS_WITH( applyNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, enc, NORM), Contains("Invalid bit encoding") );
1239  }
1240  SECTION( "two's complement register" ) {
1241 
1242  numQubitsPerReg[GENERATE_COPY(range(0,numRegs))] = 1;
1243  REQUIRE_THROWS_WITH( applyNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, TWOS_COMPLEMENT, NORM), Contains("A sub-register contained too few qubits to employ TWOS_COMPLEMENT encoding") );
1244  }
1245  SECTION( "phase function name" ) {
1246 
1247  enum phaseFunc func = (enum phaseFunc) GENERATE( -1, 14 );
1248  REQUIRE_THROWS_WITH( applyNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func), Contains("Invalid named phase function") );
1249  }
1250  SECTION( "phase function parameters" ) {
1251 
1252  enum phaseFunc func = GENERATE( SCALED_NORM, INVERSE_NORM, SCALED_INVERSE_NORM, SCALED_INVERSE_SHIFTED_NORM, SCALED_PRODUCT, INVERSE_PRODUCT, SCALED_INVERSE_PRODUCT, SCALED_DISTANCE, INVERSE_DISTANCE, SCALED_INVERSE_DISTANCE, SCALED_INVERSE_SHIFTED_DISTANCE );
1253  REQUIRE_THROWS_WITH( applyNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func), Contains("Invalid number of parameters") );
1254  }
1255  SECTION( "distance pair registers" ) {
1256 
1257  int numQb[] = {1,1,1,1,1};
1258  int qb[] = {0,1,2,3,4};
1259 
1260  numRegs = GENERATE( 1, 3, 5 );
1261  REQUIRE_THROWS_WITH( applyNamedPhaseFunc(quregVec, qb, numQb, numRegs, UNSIGNED, DISTANCE), Contains("Phase functions DISTANCE") && Contains("even number of sub-registers") );
1262  }
1263  }
1264  CLEANUP_TEST( quregVec, quregMatr );
1265 }
1266 
1267 
1268 
1273 TEST_CASE( "applyNamedPhaseFuncOverrides", "[operators]" ) {
1274 
1275  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
1276 
1277  SECTION( "correctness" ) {
1278 
1279  // try every kind of binary encoding
1280  enum bitEncoding encoding = GENERATE( UNSIGNED,TWOS_COMPLEMENT );
1281 
1282  // try every possible number of registers
1283  // (between #qubits containing 1, and 1 containing #qubits)
1284  int numRegs;
1285  int maxNumRegs = 0;
1286  if (encoding == UNSIGNED)
1287  maxNumRegs = NUM_QUBITS;
1288  if (encoding == TWOS_COMPLEMENT)
1289  maxNumRegs = NUM_QUBITS/2; // floors
1290  numRegs = GENERATE_COPY( range(1, maxNumRegs+1) );
1291 
1292  // try every possible total number of involved qubits
1293  int totalNumQubits;
1294  int minTotalQubits = 0;
1295  if (encoding == UNSIGNED)
1296  // each register must contain at least 1 qubit
1297  minTotalQubits = numRegs;
1298  if (encoding == TWOS_COMPLEMENT)
1299  // each register must contain at least 2 qubits
1300  minTotalQubits = 2*numRegs;
1301  totalNumQubits = GENERATE_COPY( range(minTotalQubits,NUM_QUBITS+1) );
1302 
1303  // try every qubits subset and ordering
1304  int* regs = GENERATE_COPY( sublists(range(0,NUM_QUBITS), totalNumQubits) );
1305 
1306  // assign each sub-reg its minimum length
1307  int unallocQubits = totalNumQubits;
1308  int numQubitsPerReg[numRegs];
1309  for (int i=0; i<numRegs; i++)
1310  if (encoding == UNSIGNED) {
1311  numQubitsPerReg[i] = 1;
1312  unallocQubits -= 1;
1313  }
1314  else if (encoding == TWOS_COMPLEMENT) {
1315  numQubitsPerReg[i] = 2;
1316  unallocQubits -= 2;
1317  }
1318  // and randomly allocate the remaining qubits between the registers
1319  while (unallocQubits > 0) {
1320  numQubitsPerReg[getRandomInt(0,numRegs)] += 1;
1321  unallocQubits--;
1322  }
1323 
1324 
1325  // choose a random number of overrides (even overriding every amplitude)
1326  int numOverrides = getRandomInt(0, (1<<totalNumQubits) + 1);
1327 
1328  // randomise each override index (uniqueness isn't checked)
1329  long long int overrideInds[numOverrides*numRegs];
1330  int flatInd = 0;
1331  for (int v=0; v<numOverrides; v++) {
1332  for (int r=0; r<numRegs; r++) {
1333  if (encoding == UNSIGNED)
1334  overrideInds[flatInd] = getRandomInt(0, 1<<numQubitsPerReg[r]);
1335  else if (encoding == TWOS_COMPLEMENT)
1336  overrideInds[flatInd] = getRandomInt(-(1<<(numQubitsPerReg[r]-1)), (1<<(numQubitsPerReg[r]-1))-1);
1337  flatInd++;
1338  }
1339  }
1340 
1341  // override to a random phase
1342  qreal overridePhases[numOverrides];
1343  for (int v=0; v<numOverrides; v++)
1344  overridePhases[v] = getRandomReal(-4, 4); // periodic in [-pi, pi]
1345 
1346 
1347  // determine the values corresponding to each register for all basis states
1348  qreal regVals[1<<totalNumQubits][numRegs];
1349  for (long long int i=0; i<(1<<totalNumQubits); i++) {
1350 
1351  long long int bits = i;
1352  for (int r=0; r<numRegs; r++) {
1353  regVals[i][r] = bits % (1 << numQubitsPerReg[r]);
1354  bits = bits >> numQubitsPerReg[r];
1355 
1356  if (encoding == TWOS_COMPLEMENT)
1357  regVals[i][r] = getTwosComplement(regVals[i][r], numQubitsPerReg[r]);
1358  }
1359  }
1360 
1361  /* a reference diagonal matrix which assumes the qubits are
1362  * contiguous and strictly increasing between the registers, and hence
1363  * only depends on the number of qubits in each register.
1364  */
1365  QMatrix diagMatr = getZeroMatrix(1 << totalNumQubits);
1366 
1367  SECTION( "NORM" ) {
1368 
1369  for (size_t i=0; i<diagMatr.size(); i++) {
1370  qreal phase = 0;
1371  for (int r=0; r<numRegs; r++)
1372  phase += pow(regVals[i][r], 2);
1373  phase = sqrt(phase);
1374  diagMatr[i][i] = expI(phase);
1375  }
1376  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
1377 
1378  SECTION( "state-vector" ) {
1379 
1380  applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, NORM, overrideInds, overridePhases, numOverrides);
1381  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1382  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1383  }
1384  SECTION( "density-matrix" ) {
1385 
1386  applyNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, NORM, overrideInds, overridePhases, numOverrides);
1387  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1388  REQUIRE( areEqual(quregMatr, refMatr, 1E4*REAL_EPS) );
1389  }
1390  }
1391  SECTION( "PRODUCT" ) {
1392 
1393  for (size_t i=0; i<diagMatr.size(); i++) {
1394  qreal phase = 1;
1395  for (int r=0; r<numRegs; r++)
1396  phase *= regVals[i][r];
1397  diagMatr[i][i] = expI(phase);
1398  }
1399 
1400  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
1401 
1402  SECTION( "state-vector" ) {
1403 
1404  applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, PRODUCT, overrideInds, overridePhases, numOverrides);
1405  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1406  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1407  }
1408  SECTION( "density-matrix" ) {
1409 
1410  applyNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, PRODUCT, overrideInds, overridePhases, numOverrides);
1411  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1412  REQUIRE( areEqual(quregMatr, refMatr, 1E4*REAL_EPS) );
1413  }
1414  }
1415  SECTION( "DISTANCE" ) {
1416 
1417  // test only if there are an even number of registers
1418  if (numRegs%2 == 0) {
1419 
1420  for (size_t i=0; i<diagMatr.size(); i++) {
1421  qreal phase = 0;
1422  for (int r=0; r<numRegs; r+=2)
1423  phase += pow(regVals[i][r+1]-regVals[i][r], 2);
1424  phase = sqrt(phase);
1425  diagMatr[i][i] = expI(phase);
1426  }
1427 
1428  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
1429  }
1430 
1431  SECTION( "state-vector" ) {
1432 
1433  if (numRegs%2 == 0) {
1434  applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, DISTANCE, overrideInds, overridePhases, numOverrides);
1435  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1436  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1437  }
1438  }
1439  SECTION( "density-matrix" ) {
1440 
1441  if (numRegs%2 == 0) {
1442  applyNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, DISTANCE, overrideInds, overridePhases, numOverrides);
1443  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1444  REQUIRE( areEqual(quregMatr, refMatr, 1E4*REAL_EPS) );
1445  }
1446  }
1447  }
1448  }
1449  SECTION( "input validation" ) {
1450 
1451  int numRegs = 2;
1452  int numQubitsPerReg[] = {2,3};
1453  int regs[] = {0,1,2,3,4};
1454 
1455  SECTION( "number of registers" ) {
1456 
1457  numRegs = GENERATE_COPY( -1, 0, 1+MAX_NUM_REGS_APPLY_ARBITRARY_PHASE );
1458  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, NULL, 0), Contains("Invalid number of qubit subregisters") );
1459  }
1460  SECTION( "number of qubits" ) {
1461 
1462  numQubitsPerReg[GENERATE_COPY(range(0,numRegs))] = GENERATE( -1, 0, 1+NUM_QUBITS );
1463  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, NULL, 0), Contains("Invalid number of qubits") );
1464  }
1465  SECTION( "repetition of qubits" ) {
1466 
1467  regs[GENERATE(2,3,4)] = regs[1];
1468  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, NULL, 0), Contains("The qubits must be unique") );
1469  }
1470  SECTION( "qubit indices" ) {
1471 
1472  regs[GENERATE(range(0,NUM_QUBITS))] = GENERATE( -1, NUM_QUBITS );
1473  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, NULL, 0), Contains("Invalid qubit index") );
1474  }
1475  SECTION( "bit encoding name" ) {
1476 
1477  enum bitEncoding enc = (enum bitEncoding) GENERATE( -1, 2 );
1478  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, enc, NORM, NULL, NULL, 0), Contains("Invalid bit encoding") );
1479  }
1480  SECTION( "two's complement register" ) {
1481 
1482  numQubitsPerReg[GENERATE_COPY(range(0,numRegs))] = 1;
1483  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, TWOS_COMPLEMENT, NORM, NULL, NULL, 0), Contains("A sub-register contained too few qubits to employ TWOS_COMPLEMENT encoding") );
1484  }
1485  SECTION( "phase function name" ) {
1486 
1487  enum phaseFunc func = (enum phaseFunc) GENERATE( -1, 14 );
1488  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, NULL, NULL, 0), Contains("Invalid named phase function") );
1489  }
1490  SECTION( "phase function parameters" ) {
1491 
1492  enum phaseFunc func = GENERATE( SCALED_NORM, INVERSE_NORM, SCALED_INVERSE_NORM, SCALED_INVERSE_SHIFTED_NORM, SCALED_PRODUCT, INVERSE_PRODUCT, SCALED_INVERSE_PRODUCT, SCALED_DISTANCE, INVERSE_DISTANCE, SCALED_INVERSE_DISTANCE, SCALED_INVERSE_SHIFTED_DISTANCE );
1493  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, NULL, NULL, 0), Contains("Invalid number of parameters") );
1494  }
1495  SECTION( "distance pair registers" ) {
1496 
1497  int numQb[] = {1,1,1,1,1};
1498  int qb[] = {0,1,2,3,4};
1499 
1500  numRegs = GENERATE( 1, 3, 5 );
1501  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, qb, numQb, numRegs, UNSIGNED, DISTANCE, NULL, NULL, 0), Contains("Phase functions DISTANCE") && Contains("even number of sub-registers") );
1502  }
1503  SECTION( "number of overrides" ) {
1504 
1505  int numOverrides = -1;
1506  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, NULL, numOverrides), Contains("Invalid number of phase function overrides specified") );
1507  }
1508  SECTION( "override indices" ) {
1509 
1510  // numQubitsPerReg = {2, 3}
1511  int numOverrides = 3;
1512  long long int overrideInds[] = {0,0, 0,0, 0,0}; // repetition not checked
1513  qreal overridePhases[] = {.1, .1, .1};
1514 
1515  // first element of overrideInds coordinate is a 2 qubit register
1516  enum bitEncoding enc = UNSIGNED;
1517  int badInd = GENERATE(0, 2, 4);
1518  overrideInds[badInd] = GENERATE( -1, (1<<2) );
1519  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, enc, NORM, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the UNSIGNED encoding") );
1520  overrideInds[badInd] = 0;
1521 
1522  // second element of overrideInds coordinate is a 3 qubit register
1523  badInd += 1;
1524  overrideInds[badInd] = GENERATE( -1, (1<<3) );
1525  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, enc, NORM, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the UNSIGNED encoding") );
1526  overrideInds[badInd] = 0;
1527  badInd -= 1;
1528 
1529  enc = TWOS_COMPLEMENT;
1530  int minInd = -(1<<(numQubitsPerReg[0]-1));
1531  int maxInd = (1<<(numQubitsPerReg[0]-1)) - 1;
1532  overrideInds[badInd] = GENERATE_COPY( minInd-1, maxInd+1 );
1533  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, enc, NORM, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the TWOS_COMPLEMENT encoding") );
1534  overrideInds[badInd] = 0;
1535 
1536  badInd++;
1537  minInd = -(1<<(numQubitsPerReg[1]-1));
1538  maxInd = (1<<(numQubitsPerReg[1]-1)) -1;
1539  overrideInds[badInd] = GENERATE_COPY( minInd-1, maxInd+1 );
1540  REQUIRE_THROWS_WITH( applyNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, enc, NORM, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the TWOS_COMPLEMENT encoding") );
1541  }
1542  }
1543  CLEANUP_TEST( quregVec, quregMatr );
1544 }
1545 
1546 
1547 
1552 TEST_CASE( "applyParamNamedPhaseFunc", "[operators]" ) {
1553 
1554  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
1555 
1556  SECTION( "correctness" ) {
1557 
1558  // try every kind of binary encoding
1559  enum bitEncoding encoding = GENERATE( UNSIGNED,TWOS_COMPLEMENT );
1560 
1561  // try every possible number of registers
1562  // (between #qubits containing 1, and 1 containing #qubits)
1563  int numRegs;
1564  int maxNumRegs = 0;
1565  if (encoding == UNSIGNED)
1566  maxNumRegs = NUM_QUBITS;
1567  if (encoding == TWOS_COMPLEMENT)
1568  maxNumRegs = NUM_QUBITS/2; // floors
1569  numRegs = GENERATE_COPY( range(1, maxNumRegs+1) );
1570 
1571  // try every possible total number of involved qubits
1572  int totalNumQubits;
1573  int minTotalQubits = 0;
1574  if (encoding == UNSIGNED)
1575  // each register must contain at least 1 qubit
1576  minTotalQubits = numRegs;
1577  if (encoding == TWOS_COMPLEMENT)
1578  // each register must contain at least 2 qubits
1579  minTotalQubits = 2*numRegs;
1580  totalNumQubits = GENERATE_COPY( range(minTotalQubits,NUM_QUBITS+1) );
1581 
1582  // try every qubits subset and ordering
1583  int* regs = GENERATE_COPY( sublists(range(0,NUM_QUBITS), totalNumQubits) );
1584 
1585  // assign each sub-reg its minimum length
1586  int unallocQubits = totalNumQubits;
1587  int numQubitsPerReg[numRegs];
1588  for (int i=0; i<numRegs; i++)
1589  if (encoding == UNSIGNED) {
1590  numQubitsPerReg[i] = 1;
1591  unallocQubits -= 1;
1592  }
1593  else if (encoding == TWOS_COMPLEMENT) {
1594  numQubitsPerReg[i] = 2;
1595  unallocQubits -= 2;
1596  }
1597  // and randomly allocate the remaining qubits between the registers
1598  while (unallocQubits > 0) {
1599  numQubitsPerReg[getRandomInt(0,numRegs)] += 1;
1600  unallocQubits--;
1601  }
1602 
1603  /* produce a reference diagonal matrix which assumes the qubits are
1604  * contiguous and strictly increasing between the registers, and hence
1605  * only depends on the number of qubits in each register.
1606  */
1607  QMatrix diagMatr = getZeroMatrix(1 << totalNumQubits);
1608 
1609  // determine the values corresponding to each register for all basis states
1610  qreal regVals[1<<totalNumQubits][numRegs];
1611  for (long long int i=0; i<(1<<totalNumQubits); i++) {
1612 
1613  long long int bits = i;
1614  for (int r=0; r<numRegs; r++) {
1615  regVals[i][r] = bits % (1 << numQubitsPerReg[r]);
1616  bits = bits >> numQubitsPerReg[r];
1617 
1618  if (encoding == TWOS_COMPLEMENT)
1619  regVals[i][r] = getTwosComplement(regVals[i][r], numQubitsPerReg[r]);
1620  }
1621  }
1622 
1623  SECTION( "INVERSE_NORM" ) {
1624 
1625  enum phaseFunc func = INVERSE_NORM;
1626  qreal divPhase = getRandomReal(-4, 4);
1627 
1628  for (size_t i=0; i<diagMatr.size(); i++) {
1629  qreal phase = 0;
1630  for (int r=0; r<numRegs; r++)
1631  phase += pow(regVals[i][r], 2);
1632  phase = (phase == 0.)? divPhase : 1/sqrt(phase);
1633  diagMatr[i][i] = expI(phase);
1634  }
1635 
1636  qreal params[] = {divPhase};
1637  int numParams = 1;
1638 
1639  SECTION( "state-vector" ) {
1640 
1641  applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1642  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1643  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1644  }
1645  SECTION( "density-matrix" ) {
1646 
1647  applyParamNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1648  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1649  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
1650  }
1651 
1652  }
1653  SECTION( "SCALED_NORM" ) {
1654 
1655  enum phaseFunc func = SCALED_NORM;
1656  qreal coeff = getRandomReal(-10, 10);
1657 
1658  for (size_t i=0; i<diagMatr.size(); i++) {
1659  qreal phase = 0;
1660  for (int r=0; r<numRegs; r++)
1661  phase += pow(regVals[i][r], 2);
1662  phase = coeff * sqrt(phase);
1663  diagMatr[i][i] = expI(phase);
1664  }
1665 
1666  qreal params[] = {coeff};
1667  int numParams = 1;
1668 
1669  SECTION( "state-vector" ) {
1670 
1671  applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1672  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1673  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1674  }
1675  SECTION( "density-matrix" ) {
1676 
1677  applyParamNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1678  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1679  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
1680  }
1681  }
1682  SECTION( "SCALED_INVERSE_NORM" ) {
1683 
1684  enum phaseFunc func = SCALED_INVERSE_NORM;
1685  qreal coeff = getRandomReal(-10, 10);
1686  qreal divPhase = getRandomReal(-4, 4);
1687 
1688  for (size_t i=0; i<diagMatr.size(); i++) {
1689  qreal phase = 0;
1690  for (int r=0; r<numRegs; r++)
1691  phase += pow(regVals[i][r], 2);
1692  phase = (phase == 0.)? divPhase : coeff/sqrt(phase);
1693  diagMatr[i][i] = expI(phase);
1694  }
1695 
1696  qreal params[] = {coeff, divPhase};
1697  int numParams = 2;
1698 
1699  SECTION( "state-vector" ) {
1700 
1701  applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1702  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1703  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1704  }
1705  SECTION( "density-matrix" ) {
1706 
1707  applyParamNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1708  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1709  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
1710  }
1711  }
1712  SECTION( "SCALED_INVERSE_SHIFTED_NORM" ) {
1713 
1714  enum phaseFunc func = SCALED_INVERSE_SHIFTED_NORM;
1715  int numParams = 2 + numRegs;
1716  qreal params[numParams];
1717  params[0] = getRandomReal(-10, 10); // scaling
1718  params[1] = getRandomReal(-4, 4); // divergence override
1719  for (int r=0; r<numRegs; r++)
1720  params[2+r] = getRandomReal(-8, 8); // shifts
1721 
1722  for (size_t i=0; i<diagMatr.size(); i++) {
1723  qreal phase = 0;
1724  for (int r=0; r<numRegs; r++)
1725  phase += pow(regVals[i][r] - params[2+r], 2);
1726  phase = sqrt(phase);
1727  phase = (phase <= REAL_EPS)? params[1] : params[0]/phase;
1728  diagMatr[i][i] = expI(phase);
1729  }
1730 
1731  SECTION( "state-vector" ) {
1732 
1733  applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1734  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1735  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1736  }
1737  SECTION( "density-matrix" ) {
1738 
1739  applyParamNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1740  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1741  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
1742  }
1743  }
1744  SECTION( "INVERSE_PRODUCT" ) {
1745 
1746  enum phaseFunc func = INVERSE_PRODUCT;
1747  qreal divPhase = getRandomReal(-4, 4);
1748 
1749  for (size_t i=0; i<diagMatr.size(); i++) {
1750  qreal phase = 1;
1751  for (int r=0; r<numRegs; r++)
1752  phase *= regVals[i][r];
1753  phase = (phase == 0.)? divPhase : 1. / phase;
1754  diagMatr[i][i] = expI(phase);
1755  }
1756 
1757  qreal params[] = {divPhase};
1758  int numParams = 1;
1759 
1760  SECTION( "state-vector" ) {
1761 
1762  applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1763  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1764  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1765  }
1766  SECTION( "density-matrix" ) {
1767 
1768  applyParamNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1769  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1770  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
1771  }
1772  }
1773  SECTION( "SCALED_PRODUCT" ) {
1774 
1775  enum phaseFunc func = SCALED_PRODUCT;
1776  qreal coeff = getRandomReal(-10, 10);
1777 
1778  for (size_t i=0; i<diagMatr.size(); i++) {
1779  qreal phase = 1;
1780  for (int r=0; r<numRegs; r++)
1781  phase *= regVals[i][r];
1782  diagMatr[i][i] = expI(coeff * phase);
1783  }
1784 
1785  qreal params[] = {coeff};
1786  int numParams = 1;
1787 
1788  SECTION( "state-vector" ) {
1789 
1790  applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1791  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1792  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1793  }
1794  SECTION( "density-matrix" ) {
1795 
1796  applyParamNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1797  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1798  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
1799  }
1800  }
1801  SECTION( "SCALED_INVERSE_PRODUCT" ) {
1802 
1803  enum phaseFunc func = SCALED_INVERSE_PRODUCT;
1804  qreal coeff = getRandomReal(-10, 10);
1805  qreal divPhase = getRandomReal(-4, 4);
1806  qreal params[] = {coeff, divPhase};
1807  int numParams = 2;
1808 
1809  for (size_t i=0; i<diagMatr.size(); i++) {
1810  qreal phase = 1;
1811  for (int r=0; r<numRegs; r++)
1812  phase *= regVals[i][r];
1813  phase = (phase == 0)? divPhase : coeff / phase;
1814  diagMatr[i][i] = expI(phase);
1815  }
1816 
1817  SECTION( "state-vector" ) {
1818 
1819  applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1820  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1821  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1822  }
1823  SECTION( "density-matrix" ) {
1824 
1825  applyParamNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1826  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1827  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
1828  }
1829  }
1830  SECTION( "INVERSE_DISTANCE" ) {
1831 
1832  enum phaseFunc func = INVERSE_DISTANCE;
1833  qreal divPhase = getRandomReal( -4, 4 );
1834  qreal params[] = {divPhase};
1835  int numParams = 1;
1836 
1837  // test only if there are an even number of registers
1838  if (numRegs%2 == 0) {
1839 
1840  for (size_t i=0; i<diagMatr.size(); i++) {
1841  qreal phase = 0;
1842  for (int r=0; r<numRegs; r+=2)
1843  phase += pow(regVals[i][r+1]-regVals[i][r], 2);
1844  phase = (phase == 0.)? divPhase : 1./sqrt(phase);
1845  diagMatr[i][i] = expI(phase);
1846  }
1847  }
1848 
1849  SECTION( "state-vector" ) {
1850 
1851  if (numRegs%2 == 0) {
1852  applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1853  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1854  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1855  }
1856  }
1857  SECTION( "density-matrix" ) {
1858 
1859  if (numRegs%2 == 0) {
1860  applyParamNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1861  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1862  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
1863  }
1864  }
1865  }
1866  SECTION( "SCALED_DISTANCE" ) {
1867 
1868  enum phaseFunc func = SCALED_DISTANCE;
1869  qreal coeff = getRandomReal( -10, 10 );
1870  qreal params[] = {coeff};
1871  int numParams = 1;
1872 
1873  // test only if there are an even number of registers
1874  if (numRegs%2 == 0) {
1875 
1876  for (size_t i=0; i<diagMatr.size(); i++) {
1877  qreal phase = 0;
1878  for (int r=0; r<numRegs; r+=2)
1879  phase += pow(regVals[i][r+1]-regVals[i][r], 2);
1880  phase = coeff * sqrt(phase);
1881  diagMatr[i][i] = expI(phase);
1882  }
1883  }
1884 
1885  SECTION( "state-vector" ) {
1886 
1887  if (numRegs%2 == 0) {
1888  applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1889  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1890  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1891  }
1892  }
1893  SECTION( "density-matrix" ) {
1894 
1895  if (numRegs%2 == 0) {
1896  applyParamNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1897  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1898  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
1899  }
1900  }
1901  }
1902  SECTION( "SCALED_INVERSE_DISTANCE" ) {
1903 
1904  enum phaseFunc func = SCALED_INVERSE_DISTANCE;
1905  qreal coeff = getRandomReal( -10, 10 );
1906  qreal divPhase = getRandomReal( -4, 4 );
1907  qreal params[] = {coeff, divPhase};
1908  int numParams = 2;
1909 
1910  // test only if there are an even number of registers
1911  if (numRegs%2 == 0) {
1912 
1913  for (size_t i=0; i<diagMatr.size(); i++) {
1914  qreal phase = 0;
1915  for (int r=0; r<numRegs; r+=2)
1916  phase += pow(regVals[i][r+1]-regVals[i][r], 2);
1917  phase = (phase == 0.)? divPhase : coeff/sqrt(phase);
1918  diagMatr[i][i] = expI(phase);
1919  }
1920  }
1921 
1922  SECTION( "state-vector" ) {
1923 
1924  if (numRegs%2 == 0) {
1925  applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1926  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1927  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1928  }
1929  }
1930  SECTION( "density-matrix" ) {
1931 
1932  if (numRegs%2 == 0) {
1933  applyParamNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1934  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1935  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
1936  }
1937  }
1938  }
1939  SECTION( "SCALED_INVERSE_SHIFTED_DISTANCE" ) {
1940 
1941  enum phaseFunc func = SCALED_INVERSE_SHIFTED_DISTANCE;
1942  int numParams = 2 + numRegs/2;
1943  qreal params[numParams];
1944 
1945  // test only if there are an even number of registers
1946  if (numRegs%2 == 0) {
1947 
1948  params[0] = getRandomReal( -10, 10 ); // scaling
1949  params[1] = getRandomReal( -4, 4 ); // divergence override
1950  for (int r=0; r<numRegs/2; r++)
1951  params[2+r] = getRandomReal( -8, 8 ); // shifts
1952 
1953  for (size_t i=0; i<diagMatr.size(); i++) {
1954  qreal phase = 0;
1955  for (int r=0; r<numRegs; r+=2)
1956  phase += pow(regVals[i][r]-regVals[i][r+1]-params[2+r/2], 2);
1957  phase = sqrt(phase);
1958  phase = (phase <= REAL_EPS)? params[1] : params[0]/phase;
1959  diagMatr[i][i] = expI(phase);
1960  }
1961  }
1962 
1963  SECTION( "state-vector" ) {
1964 
1965  if (numRegs%2 == 0) {
1966  applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1967  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
1968  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
1969  }
1970  }
1971  SECTION( "density-matrix" ) {
1972 
1973  if (numRegs%2 == 0) {
1974  applyParamNamedPhaseFunc(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams);
1975  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
1976  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
1977  }
1978  }
1979  }
1980  }
1981  SECTION( "input validation" ) {
1982 
1983  int numRegs = 2;
1984  int numQubitsPerReg[] = {2,3};
1985  int regs[] = {0,1,2,3,4};
1986 
1987  SECTION( "number of registers" ) {
1988 
1989  numRegs = GENERATE_COPY( -1, 0, 1+MAX_NUM_REGS_APPLY_ARBITRARY_PHASE );
1990  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, 0), Contains("Invalid number of qubit subregisters") );
1991  }
1992  SECTION( "number of qubits" ) {
1993 
1994  numQubitsPerReg[GENERATE_COPY(range(0,numRegs))] = GENERATE( -1, 0, 1+NUM_QUBITS );
1995  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, 0), Contains("Invalid number of qubits") );
1996  }
1997  SECTION( "repetition of qubits" ) {
1998 
1999  regs[GENERATE(2,3,4)] = regs[1];
2000  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, 0), Contains("The qubits must be unique") );
2001  }
2002  SECTION( "qubit indices" ) {
2003 
2004  regs[GENERATE(range(0,NUM_QUBITS))] = GENERATE( -1, NUM_QUBITS );
2005  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, 0), Contains("Invalid qubit index") );
2006  }
2007  SECTION( "bit encoding name" ) {
2008 
2009  enum bitEncoding enc = (enum bitEncoding) GENERATE( -1, 2 );
2010  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, enc, NORM, NULL, 0), Contains("Invalid bit encoding") );
2011  }
2012  SECTION( "two's complement register" ) {
2013 
2014  numQubitsPerReg[GENERATE_COPY(range(0,numRegs))] = 1;
2015  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, TWOS_COMPLEMENT, NORM, NULL, 0), Contains("A sub-register contained too few qubits to employ TWOS_COMPLEMENT encoding") );
2016  }
2017  SECTION( "phase function name" ) {
2018 
2019  enum phaseFunc func = (enum phaseFunc) GENERATE( -1, 14 );
2020  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, NULL, 0), Contains("Invalid named phase function") );
2021  }
2022  SECTION( "phase function parameters" ) {
2023 
2024  qreal params[numRegs + 3];
2025  for (int r=0; r<numRegs + 3; r++)
2026  params[r] = 0;
2027 
2028  SECTION( "no parameter functions" ) {
2029 
2030  enum phaseFunc func = GENERATE( NORM, PRODUCT, DISTANCE );
2031  int numParams = GENERATE( -1, 1, 2 );
2032  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, params, numParams), Contains("Invalid number of parameters") );
2033  }
2034  SECTION( "single parameter functions" ) {
2035 
2036  enum phaseFunc func = GENERATE( SCALED_NORM, INVERSE_NORM, SCALED_PRODUCT, INVERSE_PRODUCT, SCALED_DISTANCE, INVERSE_DISTANCE );
2037  int numParams = GENERATE( -1, 0, 2 );
2038  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, params, numParams), Contains("Invalid number of parameters") );
2039  }
2040  SECTION( "two parameter functions" ) {
2041 
2042  enum phaseFunc func = GENERATE( SCALED_INVERSE_NORM, SCALED_INVERSE_PRODUCT, SCALED_INVERSE_DISTANCE );
2043  int numParams = GENERATE( 0, 1, 3 );
2044  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, params, numParams), Contains("Invalid number of parameters") );
2045  }
2046  SECTION( "shifted distance" ) {
2047 
2048  if (numRegs%2 == 0) {
2049  enum phaseFunc func = SCALED_INVERSE_SHIFTED_DISTANCE;
2050  int numParams = GENERATE_COPY( 0, 1, numRegs/2 - 1, numRegs/2, numRegs/2 + 1, numRegs/2 + 3 );
2051  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, params, numParams), Contains("Invalid number of parameters") );
2052  }
2053  }
2054  SECTION( "shifted norm" ) {
2055 
2056  enum phaseFunc func = SCALED_INVERSE_SHIFTED_NORM;
2057  int numParams = GENERATE_COPY( 0, 1, numRegs-1, numRegs, numRegs+1, numRegs+3 );
2058  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, params, numParams), Contains("Invalid number of parameters") );
2059  }
2060  }
2061  SECTION( "distance pair registers" ) {
2062 
2063  int numQb[] = {1,1,1,1,1};
2064  int qb[] = {0,1,2,3,4};
2065 
2066  numRegs = GENERATE( 1, 3, 5 );
2067  REQUIRE_THROWS_WITH( applyParamNamedPhaseFunc(quregVec, qb, numQb, numRegs, UNSIGNED, DISTANCE, NULL, 0), Contains("Phase functions DISTANCE") && Contains("even number of sub-registers") );
2068  }
2069  }
2070  CLEANUP_TEST( quregVec, quregMatr );
2071 }
2072 
2073 
2074 
2079 TEST_CASE( "applyParamNamedPhaseFuncOverrides", "[operators]" ) {
2080 
2081  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
2082 
2083  SECTION( "correctness" ) {
2084 
2085  // try every kind of binary encoding
2086  enum bitEncoding encoding = GENERATE( UNSIGNED,TWOS_COMPLEMENT );
2087 
2088  // try every possible number of registers
2089  // (between #qubits containing 1, and 1 containing #qubits)
2090  int numRegs;
2091  int maxNumRegs = 0;
2092  if (encoding == UNSIGNED)
2093  maxNumRegs = NUM_QUBITS;
2094  if (encoding == TWOS_COMPLEMENT)
2095  maxNumRegs = NUM_QUBITS/2; // floors
2096  numRegs = GENERATE_COPY( range(1, maxNumRegs+1) );
2097 
2098  // try every possible total number of involved qubits
2099  int totalNumQubits;
2100  int minTotalQubits = 0;
2101  if (encoding == UNSIGNED)
2102  // each register must contain at least 1 qubit
2103  minTotalQubits = numRegs;
2104  if (encoding == TWOS_COMPLEMENT)
2105  // each register must contain at least 2 qubits
2106  minTotalQubits = 2*numRegs;
2107  totalNumQubits = GENERATE_COPY( range(minTotalQubits,NUM_QUBITS+1) );
2108 
2109  // try every qubits subset and ordering
2110  int* regs = GENERATE_COPY( sublists(range(0,NUM_QUBITS), totalNumQubits) );
2111 
2112  // assign each sub-reg its minimum length
2113  int unallocQubits = totalNumQubits;
2114  int numQubitsPerReg[numRegs];
2115  for (int i=0; i<numRegs; i++)
2116  if (encoding == UNSIGNED) {
2117  numQubitsPerReg[i] = 1;
2118  unallocQubits -= 1;
2119  }
2120  else if (encoding == TWOS_COMPLEMENT) {
2121  numQubitsPerReg[i] = 2;
2122  unallocQubits -= 2;
2123  }
2124  // and randomly allocate the remaining qubits between the registers
2125  while (unallocQubits > 0) {
2126  numQubitsPerReg[getRandomInt(0,numRegs)] += 1;
2127  unallocQubits--;
2128  }
2129 
2130 
2131  // choose a random number of overrides (even overriding every amplitude)
2132  int numOverrides = getRandomInt(0, (1<<totalNumQubits) + 1);
2133 
2134  // randomise each override index (uniqueness isn't checked)
2135  long long int overrideInds[numOverrides*numRegs];
2136  int flatInd = 0;
2137  for (int v=0; v<numOverrides; v++) {
2138  for (int r=0; r<numRegs; r++) {
2139  if (encoding == UNSIGNED)
2140  overrideInds[flatInd] = getRandomInt(0, 1<<numQubitsPerReg[r]);
2141  else if (encoding == TWOS_COMPLEMENT)
2142  overrideInds[flatInd] = getRandomInt(-(1<<(numQubitsPerReg[r]-1)), (1<<(numQubitsPerReg[r]-1))-1);
2143  flatInd++;
2144  }
2145  }
2146 
2147  // override to a random phase
2148  qreal overridePhases[numOverrides];
2149  for (int v=0; v<numOverrides; v++)
2150  overridePhases[v] = getRandomReal(-4, 4); // periodic in [-pi, pi]
2151 
2152 
2153  // determine the values corresponding to each register for all basis states
2154  qreal regVals[1<<totalNumQubits][numRegs];
2155  for (long long int i=0; i<(1<<totalNumQubits); i++) {
2156 
2157  long long int bits = i;
2158  for (int r=0; r<numRegs; r++) {
2159  regVals[i][r] = bits % (1 << numQubitsPerReg[r]);
2160  bits = bits >> numQubitsPerReg[r];
2161 
2162  if (encoding == TWOS_COMPLEMENT)
2163  regVals[i][r] = getTwosComplement(regVals[i][r], numQubitsPerReg[r]);
2164  }
2165  }
2166 
2167  /* produce a reference diagonal matrix which assumes the qubits are
2168  * contiguous and strictly increasing between the registers, and hence
2169  * only depends on the number of qubits in each register.
2170  */
2171  QMatrix diagMatr = getZeroMatrix(1 << totalNumQubits);
2172 
2173 
2174  SECTION( "INVERSE_NORM" ) {
2175 
2176  enum phaseFunc func = INVERSE_NORM;
2177  qreal divPhase = getRandomReal(-4, 4);
2178  qreal params[] = {divPhase};
2179  int numParams = 1;
2180 
2181  for (size_t i=0; i<diagMatr.size(); i++) {
2182  qreal phase = 0;
2183  for (int r=0; r<numRegs; r++)
2184  phase += pow(regVals[i][r], 2);
2185  phase = (phase == 0.)? divPhase : 1/sqrt(phase);
2186  diagMatr[i][i] = expI(phase);
2187  }
2188  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
2189 
2190  SECTION( "state-vector" ) {
2191 
2192  applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2193  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
2194  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
2195  }
2196  SECTION( "density-matrix" ) {
2197 
2198  applyParamNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2199  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
2200  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
2201  }
2202  }
2203  SECTION( "SCALED_NORM" ) {
2204 
2205  enum phaseFunc func = SCALED_NORM;
2206  qreal coeff = getRandomReal(-10, 10);
2207  qreal params[] = {coeff};
2208  int numParams = 1;
2209 
2210  for (size_t i=0; i<diagMatr.size(); i++) {
2211  qreal phase = 0;
2212  for (int r=0; r<numRegs; r++)
2213  phase += pow(regVals[i][r], 2);
2214  phase = coeff * sqrt(phase);
2215  diagMatr[i][i] = expI(phase);
2216  }
2217  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
2218 
2219  SECTION( "state-vector" ) {
2220 
2221  applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2222  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
2223  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
2224  }
2225  SECTION( "density-matrix" ) {
2226 
2227  applyParamNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2228  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
2229  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
2230  }
2231  }
2232  SECTION( "SCALED_INVERSE_NORM" ) {
2233 
2234  enum phaseFunc func = SCALED_INVERSE_NORM;
2235  qreal coeff = getRandomReal(-10, 10);
2236  qreal divPhase = getRandomReal(-4, 4);
2237  qreal params[] = {coeff, divPhase};
2238  int numParams = 2;
2239 
2240  for (size_t i=0; i<diagMatr.size(); i++) {
2241  qreal phase = 0;
2242  for (int r=0; r<numRegs; r++)
2243  phase += pow(regVals[i][r], 2);
2244  phase = (phase == 0.)? divPhase : coeff/sqrt(phase);
2245  diagMatr[i][i] = expI(phase);
2246  }
2247  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
2248 
2249  SECTION( "state-vector" ) {
2250 
2251  applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2252  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
2253  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
2254  }
2255  SECTION( "density-matrix" ) {
2256 
2257  applyParamNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2258  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
2259  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
2260  }
2261  }
2262  SECTION( "SCALED_INVERSE_SHIFTED_NORM" ) {
2263 
2264  enum phaseFunc func = SCALED_INVERSE_SHIFTED_NORM;
2265  int numParams = 2 + numRegs;
2266  qreal params[numParams];
2267  params[0] = getRandomReal(-10, 10); // scaling
2268  params[1] = getRandomReal(-4, 4); // divergence override
2269  for (int r=0; r<numRegs; r++)
2270  params[2+r] = getRandomReal(-8, 8); // shifts
2271 
2272  for (size_t i=0; i<diagMatr.size(); i++) {
2273  qreal phase = 0;
2274  for (int r=0; r<numRegs; r++)
2275  phase += pow(regVals[i][r] - params[2+r], 2);
2276  phase = sqrt(phase);
2277  phase = (phase <= REAL_EPS)? params[1] : params[0]/phase;
2278  diagMatr[i][i] = expI(phase);
2279  }
2280  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
2281 
2282  SECTION( "state-vector" ) {
2283 
2284  applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2285  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
2286  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
2287  }
2288  SECTION( "density-matrix" ) {
2289 
2290  applyParamNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2291  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
2292  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
2293  }
2294  }
2295  SECTION( "INVERSE_PRODUCT" ) {
2296 
2297  enum phaseFunc func = INVERSE_PRODUCT;
2298  qreal divPhase = getRandomReal(-4, 4);
2299  qreal params[] = {divPhase};
2300  int numParams = 1;
2301 
2302  for (size_t i=0; i<diagMatr.size(); i++) {
2303  qreal phase = 1;
2304  for (int r=0; r<numRegs; r++)
2305  phase *= regVals[i][r];
2306  phase = (phase == 0.)? divPhase : 1. / phase;
2307  diagMatr[i][i] = expI(phase);
2308  }
2309  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
2310 
2311  SECTION( "state-vector" ) {
2312 
2313  applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2314  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
2315  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
2316  }
2317  SECTION( "density-matrix" ) {
2318 
2319  applyParamNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2320  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
2321  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
2322  }
2323  }
2324  SECTION( "SCALED_PRODUCT" ) {
2325 
2326  enum phaseFunc func = SCALED_PRODUCT;
2327  qreal coeff = getRandomReal(-10, 10);
2328  qreal params[] = {coeff};
2329  int numParams = 1;
2330 
2331  for (size_t i=0; i<diagMatr.size(); i++) {
2332  qreal phase = 1;
2333  for (int r=0; r<numRegs; r++)
2334  phase *= regVals[i][r];
2335  diagMatr[i][i] = expI(coeff * phase);
2336  }
2337  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
2338 
2339  SECTION( "state-vector" ) {
2340 
2341  applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2342  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
2343  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
2344  }
2345  SECTION( "density-matrix" ) {
2346 
2347  applyParamNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2348  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
2349  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
2350  }
2351  }
2352  SECTION( "SCALED_INVERSE_PRODUCT" ) {
2353 
2354  enum phaseFunc func = SCALED_INVERSE_PRODUCT;
2355  qreal coeff = getRandomReal(-10, 10);
2356  qreal divPhase = getRandomReal(-4, 4);
2357  qreal params[] = {coeff, divPhase};
2358  int numParams = 2;
2359 
2360  for (size_t i=0; i<diagMatr.size(); i++) {
2361  qreal phase = 1;
2362  for (int r=0; r<numRegs; r++)
2363  phase *= regVals[i][r];
2364  phase = (phase == 0)? divPhase : coeff / phase;
2365  diagMatr[i][i] = expI(phase);
2366  }
2367  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
2368 
2369  SECTION( "state-vector" ) {
2370 
2371  applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2372  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
2373  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
2374  }
2375  SECTION( "density-matrix" ) {
2376 
2377  applyParamNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2378  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
2379  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
2380  }
2381  }
2382  SECTION( "INVERSE_DISTANCE" ) {
2383 
2384  enum phaseFunc func = INVERSE_DISTANCE;
2385  qreal divPhase = getRandomReal( -4, 4 );
2386  qreal params[] = {divPhase};
2387  int numParams = 1;
2388 
2389  // test only if there are an even number of registers
2390  if (numRegs%2 == 0) {
2391 
2392  for (size_t i=0; i<diagMatr.size(); i++) {
2393  qreal phase = 0;
2394  for (int r=0; r<numRegs; r+=2)
2395  phase += pow(regVals[i][r+1]-regVals[i][r], 2);
2396  phase = (phase == 0.)? divPhase : 1./sqrt(phase);
2397  diagMatr[i][i] = expI(phase);
2398  }
2399  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
2400  }
2401 
2402  SECTION( "state-vector" ) {
2403 
2404  if (numRegs%2 == 0) {
2405  applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2406  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
2407  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
2408  }
2409  }
2410  SECTION( "density-matrix" ) {
2411 
2412  if (numRegs%2 == 0) {
2413  applyParamNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2414  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
2415  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
2416  }
2417  }
2418  }
2419  SECTION( "SCALED_DISTANCE" ) {
2420 
2421  enum phaseFunc func = SCALED_DISTANCE;
2422  qreal coeff = getRandomReal( -10, 10 );
2423  qreal params[] = {coeff};
2424  int numParams = 1;
2425 
2426  // test only if there are an even number of registers
2427  if (numRegs%2 == 0) {
2428 
2429  for (size_t i=0; i<diagMatr.size(); i++) {
2430  qreal phase = 0;
2431  for (int r=0; r<numRegs; r+=2)
2432  phase += pow(regVals[i][r+1]-regVals[i][r], 2);
2433  phase = coeff * sqrt(phase);
2434  diagMatr[i][i] = expI(phase);
2435  }
2436  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
2437  }
2438 
2439  SECTION( "state-vector" ) {
2440 
2441  if (numRegs%2 == 0) {
2442  applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2443  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
2444  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
2445  }
2446  }
2447  SECTION( "density-matrix" ) {
2448 
2449  if (numRegs%2 == 0) {
2450  applyParamNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2451  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
2452  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
2453  }
2454  }
2455  }
2456  SECTION( "SCALED_INVERSE_DISTANCE" ) {
2457 
2458  enum phaseFunc func = SCALED_INVERSE_DISTANCE;
2459  qreal coeff = getRandomReal( -10, 10 );
2460  qreal divPhase = getRandomReal( -4, 4 );
2461  qreal params[] = {coeff, divPhase};
2462  int numParams = 2;
2463 
2464  // test only if there are an even number of registers
2465  if (numRegs%2 == 0) {
2466 
2467  for (size_t i=0; i<diagMatr.size(); i++) {
2468  qreal phase = 0;
2469  for (int r=0; r<numRegs; r+=2)
2470  phase += pow(regVals[i][r+1]-regVals[i][r], 2);
2471  phase = (phase == 0.)? divPhase : coeff/sqrt(phase);
2472  diagMatr[i][i] = expI(phase);
2473  }
2474  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
2475  }
2476 
2477  SECTION( "state-vector" ) {
2478 
2479  if (numRegs%2 == 0) {
2480  applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2481  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
2482  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
2483  }
2484  }
2485  SECTION( "density-matrix" ) {
2486 
2487  if (numRegs%2 == 0) {
2488  applyParamNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2489  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
2490  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
2491  }
2492  }
2493  }
2494  SECTION( "SCALED_INVERSE_SHIFTED_DISTANCE" ) {
2495 
2496  enum phaseFunc func = SCALED_INVERSE_SHIFTED_DISTANCE;
2497  int numParams = 2 + numRegs/2;
2498  qreal params[numParams];
2499 
2500  // test only if there are an even number of registers
2501  if (numRegs%2 == 0) {
2502 
2503  params[0] = getRandomReal( -10, 10 ); // scaling
2504  params[1] = getRandomReal( -4, 4 ); // divergence override
2505  for (int r=0; r<numRegs/2; r++)
2506  params[2+r] = getRandomReal( -8, 8 ); // shifts
2507 
2508  for (size_t i=0; i<diagMatr.size(); i++) {
2509  qreal phase = 0;
2510  for (int r=0; r<numRegs; r+=2)
2511  phase += pow(regVals[i][r]-regVals[i][r+1]-params[2+r/2], 2);
2512  phase = sqrt(phase);
2513  phase = (phase <= REAL_EPS)? params[1] : params[0]/phase;
2514  diagMatr[i][i] = expI(phase);
2515  }
2516 
2517  setDiagMatrixOverrides(diagMatr, numQubitsPerReg, numRegs, encoding, overrideInds, overridePhases, numOverrides);
2518  }
2519 
2520  SECTION( "state-vector" ) {
2521 
2522  if (numRegs%2 == 0) {
2523  applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2524  applyReferenceOp(refVec, regs, totalNumQubits, diagMatr);
2525  REQUIRE( areEqual(quregVec, refVec, 1E2*REAL_EPS) );
2526  }
2527  }
2528  SECTION( "density-matrix" ) {
2529 
2530  if (numRegs%2 == 0) {
2531  applyParamNamedPhaseFuncOverrides(quregMatr, regs, numQubitsPerReg, numRegs, encoding, func, params, numParams, overrideInds, overridePhases, numOverrides);
2532  applyReferenceOp(refMatr, regs, totalNumQubits, diagMatr);
2533  REQUIRE( areEqual(quregMatr, refMatr, 1E2*REAL_EPS) );
2534  }
2535  }
2536  }
2537  }
2538  SECTION( "input validation" ) {
2539 
2540  int numRegs = 2;
2541  int numQubitsPerReg[] = {2,3};
2542  int regs[] = {0,1,2,3,4};
2543 
2544  SECTION( "number of registers" ) {
2545 
2546  numRegs = GENERATE_COPY( -1, 0, 1+MAX_NUM_REGS_APPLY_ARBITRARY_PHASE );
2547  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, 0, NULL, NULL, 0), Contains("Invalid number of qubit subregisters") );
2548  }
2549  SECTION( "number of qubits" ) {
2550 
2551  numQubitsPerReg[GENERATE_COPY(range(0,numRegs))] = GENERATE( -1, 0, 1+NUM_QUBITS );
2552  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, 0, NULL, NULL, 0), Contains("Invalid number of qubits") );
2553  }
2554  SECTION( "repetition of qubits" ) {
2555 
2556  regs[GENERATE(2,3,4)] = regs[1];
2557  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, 0, NULL, NULL, 0), Contains("The qubits must be unique") );
2558  }
2559  SECTION( "qubit indices" ) {
2560 
2561  regs[GENERATE(range(0,NUM_QUBITS))] = GENERATE( -1, NUM_QUBITS );
2562  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, 0, NULL, NULL, 0), Contains("Invalid qubit index") );
2563  }
2564  SECTION( "bit encoding name" ) {
2565 
2566  enum bitEncoding enc = (enum bitEncoding) GENERATE( -1, 2 );
2567  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, enc, NORM, NULL, 0, NULL, NULL, 0), Contains("Invalid bit encoding") );
2568  }
2569  SECTION( "two's complement register" ) {
2570 
2571  numQubitsPerReg[GENERATE_COPY(range(0,numRegs))] = 1;
2572  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, TWOS_COMPLEMENT, NORM, NULL, 0, NULL, NULL, 0), Contains("A sub-register contained too few qubits to employ TWOS_COMPLEMENT encoding") );
2573  }
2574  SECTION( "phase function name" ) {
2575 
2576  enum phaseFunc func = (enum phaseFunc) GENERATE( -1, 14 );
2577  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, NULL, 0, NULL, NULL, 0), Contains("Invalid named phase function") );
2578  }
2579  SECTION( "phase function parameters" ) {
2580 
2581  qreal params[numRegs + 3];
2582  for (int r=0; r<numRegs + 3; r++)
2583  params[r] = 0;
2584 
2585  SECTION( "no parameter functions" ) {
2586 
2587  enum phaseFunc func = GENERATE( NORM, PRODUCT, DISTANCE );
2588  int numParams = GENERATE( -1, 1, 2 );
2589  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, params, numParams, NULL, NULL, 0), Contains("Invalid number of parameters") );
2590  }
2591  SECTION( "single parameter functions" ) {
2592 
2593  enum phaseFunc func = GENERATE( SCALED_NORM, INVERSE_NORM, SCALED_PRODUCT, INVERSE_PRODUCT, SCALED_DISTANCE, INVERSE_DISTANCE );
2594  int numParams = GENERATE( -1, 0, 2 );
2595  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, params, numParams, NULL, NULL, 0), Contains("Invalid number of parameters") );
2596  }
2597  SECTION( "two parameter functions" ) {
2598 
2599  enum phaseFunc func = GENERATE( SCALED_INVERSE_NORM, SCALED_INVERSE_PRODUCT, SCALED_INVERSE_DISTANCE );
2600  int numParams = GENERATE( 0, 1, 3 );
2601  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, params, numParams, NULL, NULL, 0), Contains("Invalid number of parameters") );
2602  }
2603  SECTION( "shifted distance" ) {
2604 
2605  if (numRegs%2 == 0) {
2606  enum phaseFunc func = SCALED_INVERSE_SHIFTED_DISTANCE;
2607  int numParams = GENERATE_COPY( 0, 1, numRegs/2 - 1, numRegs/2, numRegs/2 + 1, numRegs/2 + 3 );
2608  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, params, numParams, NULL, NULL, 0), Contains("Invalid number of parameters") );
2609  }
2610  }
2611  SECTION( "shifted norm" ) {
2612 
2613  enum phaseFunc func = SCALED_INVERSE_SHIFTED_NORM;
2614  int numParams = GENERATE_COPY( 0, 1, numRegs-1, numRegs, numRegs+1, numRegs+3 );
2615  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, func, params, numParams, NULL, NULL, 0), Contains("Invalid number of parameters") );
2616  }
2617  }
2618  SECTION( "distance pair registers" ) {
2619 
2620  int numQb[] = {1,1,1,1,1};
2621  int qb[] = {0,1,2,3,4};
2622 
2623  numRegs = GENERATE( 1, 3, 5 );
2624  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, qb, numQb, numRegs, UNSIGNED, DISTANCE, NULL, 0, NULL, NULL, 0), Contains("Phase functions DISTANCE") && Contains("even number of sub-registers") );
2625  }
2626  SECTION( "number of overrides" ) {
2627 
2628  int numOverrides = -1;
2629  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, UNSIGNED, NORM, NULL, 0, NULL, NULL, numOverrides), Contains("Invalid number of phase function overrides specified") );
2630  }
2631  SECTION( "override indices" ) {
2632 
2633  // numQubitsPerReg = {2, 3}
2634  int numOverrides = 3;
2635  long long int overrideInds[] = {0,0, 0,0, 0,0}; // repetition not checked
2636  qreal overridePhases[] = {.1, .1, .1};
2637 
2638  // first element of overrideInds coordinate is a 2 qubit register
2639  enum bitEncoding enc = UNSIGNED;
2640  int badInd = GENERATE(0, 2, 4);
2641  overrideInds[badInd] = GENERATE( -1, (1<<2) );
2642  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, enc, NORM, NULL, 0, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the UNSIGNED encoding") );
2643  overrideInds[badInd] = 0;
2644 
2645  // second element of overrideInds coordinate is a 3 qubit register
2646  badInd += 1;
2647  overrideInds[badInd] = GENERATE( -1, (1<<3) );
2648  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, enc, NORM, NULL, 0, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the UNSIGNED encoding") );
2649  overrideInds[badInd] = 0;
2650  badInd -= 1;
2651 
2652  enc = TWOS_COMPLEMENT;
2653  int minInd = -(1<<(numQubitsPerReg[0]-1));
2654  int maxInd = (1<<(numQubitsPerReg[0]-1)) - 1;
2655  overrideInds[badInd] = GENERATE_COPY( minInd-1, maxInd+1 );
2656  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, enc, NORM, NULL, 0, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the TWOS_COMPLEMENT encoding") );
2657  overrideInds[badInd] = 0;
2658 
2659  badInd++;
2660  minInd = -(1<<(numQubitsPerReg[1]-1));
2661  maxInd = (1<<(numQubitsPerReg[1]-1)) -1;
2662  overrideInds[badInd] = GENERATE_COPY( minInd-1, maxInd+1 );
2663  REQUIRE_THROWS_WITH( applyParamNamedPhaseFuncOverrides(quregVec, regs, numQubitsPerReg, numRegs, enc, NORM, NULL, 0, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the TWOS_COMPLEMENT encoding") );
2664  }
2665  }
2666  CLEANUP_TEST( quregVec, quregMatr );
2667 }
2668 
2669 
2670 
2675 TEST_CASE( "applyPauliHamil", "[operators]" ) {
2676 
2677  Qureg vecIn = createQureg(NUM_QUBITS, QUEST_ENV);
2678  Qureg vecOut = createQureg(NUM_QUBITS, QUEST_ENV);
2679  Qureg matIn = createDensityQureg(NUM_QUBITS, QUEST_ENV);
2680  Qureg matOut = createDensityQureg(NUM_QUBITS, QUEST_ENV);
2681 
2682  initDebugState(vecIn);
2683  initDebugState(matIn);
2684 
2685  SECTION( "correctness" ) {
2686 
2687  /* it's too expensive to try every possible Pauli configuration, so
2688  * we'll try 10 random codes, and for each, random coefficients
2689  */
2690  GENERATE( range(0,10) );
2691 
2692  int numTerms = GENERATE( 1, 2, 10, 15 );
2693  PauliHamil hamil = createPauliHamil(NUM_QUBITS, numTerms);
2694  setRandomPauliSum(hamil);
2695  QMatrix refHamil = toQMatrix(hamil);
2696 
2697  SECTION( "state-vector" ) {
2698 
2699  QVector vecRef = toQVector(vecIn);
2700  applyPauliHamil(vecIn, hamil, vecOut);
2701 
2702  // ensure vecIn barely changes under precision
2703  REQUIRE( areEqual(vecIn, vecRef) );
2704 
2705  // ensure vecOut changed correctly
2706  REQUIRE( areEqual(vecOut, refHamil * vecRef) );
2707  }
2708  SECTION( "density-matrix" ) {
2709 
2710  QMatrix matRef = toQMatrix(matIn);
2711  applyPauliHamil(matIn, hamil, matOut);
2712 
2713  // ensure matIn barely changes under precision
2714  REQUIRE( areEqual(matIn, matRef) );
2715 
2716  // ensure matOut changed correctly
2717  REQUIRE( areEqual(matOut, refHamil * matRef, 1E2*REAL_EPS) );
2718  }
2719 
2720  destroyPauliHamil(hamil);
2721  }
2722  SECTION( "input validation" ) {
2723 
2724  SECTION( "pauli codes" ) {
2725 
2726  int numTerms = 3;
2727  PauliHamil hamil = createPauliHamil(NUM_QUBITS, numTerms);
2728 
2729  // make one pauli code wrong
2730  hamil.pauliCodes[GENERATE_COPY( range(0,numTerms*NUM_QUBITS) )] = (pauliOpType) GENERATE( -1, 4 );
2731  REQUIRE_THROWS_WITH( applyPauliHamil(vecIn, hamil, vecOut), Contains("Invalid Pauli code") );
2732 
2733  destroyPauliHamil(hamil);
2734  }
2735  SECTION( "qureg dimensions" ) {
2736 
2737  Qureg badVec = createQureg(NUM_QUBITS+1, QUEST_ENV);
2738  PauliHamil hamil = createPauliHamil(NUM_QUBITS, 1);
2739 
2740  REQUIRE_THROWS_WITH( applyPauliHamil(vecIn, hamil, badVec), Contains("Dimensions of the qubit registers don't match") );
2741 
2742  destroyQureg(badVec, QUEST_ENV);
2743  destroyPauliHamil(hamil);
2744  }
2745  SECTION( "qureg types" ) {
2746 
2747  PauliHamil hamil = createPauliHamil(NUM_QUBITS, 1);
2748 
2749  REQUIRE_THROWS_WITH( applyPauliHamil(vecIn, hamil, matOut), Contains("Registers must both be state-vectors or both be density matrices") );
2750 
2751  destroyPauliHamil(hamil);
2752  }
2753  SECTION( "matching hamiltonian qubits" ) {
2754 
2755  PauliHamil hamil = createPauliHamil(NUM_QUBITS + 1, 1);
2756 
2757  REQUIRE_THROWS_WITH( applyPauliHamil(vecIn, hamil, vecOut), Contains("same number of qubits") );
2758  REQUIRE_THROWS_WITH( applyPauliHamil(matIn, hamil, matOut), Contains("same number of qubits") );
2759 
2760  destroyPauliHamil(hamil);
2761  }
2762  }
2763  destroyQureg(vecIn, QUEST_ENV);
2764  destroyQureg(vecOut, QUEST_ENV);
2765  destroyQureg(matIn, QUEST_ENV);
2766  destroyQureg(matOut, QUEST_ENV);
2767 }
2768 
2769 
2770 
2775 TEST_CASE( "applyPauliSum", "[operators]" ) {
2776 
2777  Qureg vecIn = createQureg(NUM_QUBITS, QUEST_ENV);
2778  Qureg vecOut = createQureg(NUM_QUBITS, QUEST_ENV);
2779  Qureg matIn = createDensityQureg(NUM_QUBITS, QUEST_ENV);
2780  Qureg matOut = createDensityQureg(NUM_QUBITS, QUEST_ENV);
2781 
2782  initDebugState(vecIn);
2783  initDebugState(matIn);
2784 
2785  SECTION( "correctness" ) {
2786 
2787  /* it's too expensive to try ALL Pauli sequences, via
2788  * pauliOpType* paulis = GENERATE_COPY( pauliseqs(numPaulis) );.
2789  * Furthermore, take(10, pauliseqs(numTargs)) will try the same pauli codes.
2790  * Hence, we instead opt to repeatedly randomly generate pauliseqs
2791  */
2792  GENERATE( range(0,10) ); // gen 10 random pauli-codes
2793 
2794  int numTerms = GENERATE( 1, 2, 10, 15);
2795  int numPaulis = numTerms * NUM_QUBITS;
2796 
2797  // each test will use random coefficients
2798  qreal coeffs[numTerms];
2799  pauliOpType paulis[numPaulis];
2800  setRandomPauliSum(coeffs, paulis, NUM_QUBITS, numTerms);
2801  QMatrix pauliSum = toQMatrix(coeffs, paulis, NUM_QUBITS, numTerms);
2802 
2803  SECTION( "state-vector" ) {
2804 
2805  QVector vecRef = toQVector(vecIn);
2806  applyPauliSum(vecIn, paulis, coeffs, numTerms, vecOut);
2807 
2808  // ensure vecIn barely changes under precision
2809  REQUIRE( areEqual(vecIn, vecRef) );
2810 
2811  // ensure vecOut changed correctly
2812  REQUIRE( areEqual(vecOut, pauliSum * vecRef) );
2813  }
2814  SECTION( "density-matrix" ) {
2815 
2816  QMatrix matRef = toQMatrix(matIn);
2817  applyPauliSum(matIn, paulis, coeffs, numTerms, matOut);
2818 
2819  // ensure matIn barely changes under precision
2820  REQUIRE( areEqual(matIn, matRef) );
2821 
2822  // ensure matOut changed correctly
2823  REQUIRE( areEqual(matOut, pauliSum * matRef, 1E2*REAL_EPS) );
2824  }
2825  }
2826  SECTION( "input validation" ) {
2827 
2828  SECTION( "number of terms" ) {
2829 
2830  int numTerms = GENERATE( -1, 0 );
2831  REQUIRE_THROWS_WITH( applyPauliSum(vecIn, NULL, NULL, numTerms, vecOut), Contains("Invalid number of terms") );
2832  }
2833  SECTION( "pauli codes" ) {
2834 
2835  // valid codes
2836  int numTerms = 3;
2837  int numPaulis = numTerms*NUM_QUBITS;
2838  pauliOpType paulis[numPaulis];
2839  for (int i=0; i<numPaulis; i++)
2840  paulis[i] = PAULI_I;
2841 
2842  // make one invalid
2843  paulis[GENERATE_COPY( range(0,numPaulis) )] = (pauliOpType) GENERATE( -1, 4 );
2844  REQUIRE_THROWS_WITH( applyPauliSum(vecIn, paulis, NULL, numTerms, vecOut), Contains("Invalid Pauli code") );
2845  }
2846  SECTION( "qureg dimensions" ) {
2847 
2848  Qureg badVec = createQureg(NUM_QUBITS+1, QUEST_ENV);
2849  pauliOpType paulis[NUM_QUBITS];
2850  REQUIRE_THROWS_WITH( applyPauliSum(vecIn, paulis, NULL, 1, badVec), Contains("Dimensions of the qubit registers don't match") );
2851  destroyQureg(badVec, QUEST_ENV);
2852  }
2853  SECTION( "qureg types" ) {
2854 
2855  pauliOpType paulis[NUM_QUBITS];
2856  REQUIRE_THROWS_WITH( applyPauliSum(vecIn, paulis, NULL, 1, matOut), Contains("Registers must both be state-vectors or both be density matrices") );
2857  }
2858  }
2859  destroyQureg(vecIn, QUEST_ENV);
2860  destroyQureg(vecOut, QUEST_ENV);
2861  destroyQureg(matIn, QUEST_ENV);
2862  destroyQureg(matOut, QUEST_ENV);
2863 }
2864 
2865 
2866 
2871 TEST_CASE( "applyPhaseFunc", "[operators]" ) {
2872 
2873  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
2874 
2875  SECTION( "correctness" ) {
2876 
2877  // try every kind of binary encodings
2878  enum bitEncoding encoding = GENERATE( UNSIGNED,TWOS_COMPLEMENT );
2879 
2880  // try every sub-register size
2881  int numQubits = GENERATE_COPY( range(1,NUM_QUBITS+1) );
2882 
2883  // force at least 2 qubits in two's compement though
2884  if (encoding == TWOS_COMPLEMENT && numQubits == 1)
2885  numQubits++;
2886 
2887  // try every possible sub-register
2888  int* qubits = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numQubits) );
2889 
2890  // choose a random number of terms in the phase function
2891  int numTerms = getRandomInt(1, 5);
2892 
2893  // populate the phase function with random but POSITIVE-power terms,
2894  // and in two's complement mode, strictly integer powers
2895  qreal coeffs[numTerms];
2896  qreal expons[numTerms];
2897  for (int t=0; t<numTerms; t++) {
2898  coeffs[t] = getRandomReal(-10,10);
2899  if (encoding == TWOS_COMPLEMENT)
2900  expons[t] = getRandomInt(0, 3+1); // repetition of power is ok
2901  else
2902  expons[t] = getRandomReal(0, 3);
2903  }
2904 
2905  // build a reference diagonal matrix, on the reduced Hilbert space
2906  QMatrix matr = getZeroMatrix( 1 << numQubits );
2907  for (size_t i=0; i<matr.size(); i++) {
2908 
2909  long long int ind = 0;
2910  if (encoding == UNSIGNED)
2911  ind = i;
2912  if (encoding == TWOS_COMPLEMENT)
2913  ind = getTwosComplement(i, numQubits);
2914 
2915  qreal phase = 0;
2916  for (int t=0; t<numTerms; t++)
2917  phase += coeffs[t] * pow(ind, expons[t]);
2918 
2919  matr[i][i] = expI(phase);
2920  }
2921 
2922  SECTION( "state-vector" ) {
2923 
2924  applyPhaseFunc(quregVec, qubits, numQubits, encoding, coeffs, expons, numTerms);
2925  applyReferenceOp(refVec, qubits, numQubits, matr);
2926  REQUIRE( areEqual(quregVec, refVec, 1E4*REAL_EPS) );
2927  }
2928  SECTION( "density-matrix" ) {
2929 
2930  applyPhaseFunc(quregMatr, qubits, numQubits, encoding, coeffs, expons, numTerms);
2931  applyReferenceOp(refMatr, qubits, numQubits, matr);
2932  REQUIRE( areEqual(quregMatr, refMatr, 1E6*REAL_EPS) );
2933  }
2934  }
2935  SECTION( "input validation" ) {
2936 
2937  int numQubits = 3;
2938  int qubits[] = {0,1,2};
2939 
2940  SECTION( "number of qubits" ) {
2941 
2942  numQubits = GENERATE_COPY( -1, 0, NUM_QUBITS+1 );
2943  REQUIRE_THROWS_WITH( applyPhaseFunc(quregVec, NULL, numQubits, UNSIGNED, NULL, NULL, 1), Contains("Invalid number of qubits") );
2944  }
2945  SECTION( "repetition of qubits" ) {
2946 
2947  qubits[GENERATE(1,2)] = qubits[0];
2948  REQUIRE_THROWS_WITH( applyPhaseFunc(quregVec, qubits, numQubits, UNSIGNED, NULL, NULL, 1), Contains("The qubits must be unique") );
2949  }
2950  SECTION( "qubit indices" ) {
2951 
2952  int inv = GENERATE( -1, NUM_QUBITS );
2953  qubits[ GENERATE_COPY( range(0,numQubits) )] = inv;
2954  REQUIRE_THROWS_WITH( applyPhaseFunc(quregVec, qubits, numQubits, UNSIGNED, NULL, NULL, 1), Contains("Invalid qubit index") );
2955  }
2956  SECTION( "number of terms" ) {
2957 
2958  int numTerms = GENERATE( -1, 0 );
2959  REQUIRE_THROWS_WITH( applyPhaseFunc(quregVec, qubits, numQubits, UNSIGNED, NULL, NULL, numTerms), Contains("Invalid number of terms in the phase function") );
2960  }
2961  SECTION( "bit encoding name" ) {
2962 
2963  enum bitEncoding enc = (enum bitEncoding) GENERATE( -1,2 );
2964  REQUIRE_THROWS_WITH( applyPhaseFunc(quregVec, qubits, numQubits, enc, NULL, NULL, 1), Contains("Invalid bit encoding") );
2965  }
2966  SECTION( "two's complement register" ) {
2967 
2968  numQubits = 1;
2969  REQUIRE_THROWS_WITH( applyPhaseFunc(quregVec, qubits, numQubits, TWOS_COMPLEMENT, NULL, NULL, 1), Contains("too few qubits to employ TWOS_COMPLEMENT") );
2970  }
2971  SECTION( "fractional exponent" ) {
2972 
2973  int numTerms = 3;
2974  qreal coeffs[] = {0,0,0};
2975  qreal expos[] = {1,2,3};
2976  expos[GENERATE_COPY( range(0,numTerms) )] = GENERATE( 0.5, 1.999, 5.0001 );
2977 
2978  REQUIRE_THROWS_WITH( applyPhaseFunc(quregVec, qubits, numQubits, TWOS_COMPLEMENT, coeffs, expos, numTerms), Contains("fractional exponent") && Contains("TWOS_COMPLEMENT") && Contains("negative indices were not overriden") );
2979  }
2980  SECTION( "negative exponent" ) {
2981 
2982  int numTerms = 3;
2983  qreal coeffs[] = {0,0,0};
2984  qreal expos[] = {1,2,3};
2985  expos[GENERATE_COPY( range(0,numTerms) )] = GENERATE( -1, -2, -1.5 );
2986 
2987  enum bitEncoding encoding = GENERATE( UNSIGNED, TWOS_COMPLEMENT );
2988 
2989  REQUIRE_THROWS_WITH( applyPhaseFunc(quregVec, qubits, numQubits, encoding, coeffs, expos, numTerms), Contains("The phase function contained a negative exponent which would diverge at zero, but the zero index was not overriden") );
2990  }
2991  }
2992  CLEANUP_TEST( quregVec, quregMatr );
2993 }
2994 
2995 
2996 
3001 TEST_CASE( "applyPhaseFuncOverrides", "[operators]" ) {
3002 
3003  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
3004 
3005  SECTION( "correctness" ) {
3006 
3007  // try every kind of binary encodings
3008  enum bitEncoding encoding = GENERATE( UNSIGNED,TWOS_COMPLEMENT );
3009 
3010  // try every sub-register size
3011  int numQubits = GENERATE_COPY( range(1,NUM_QUBITS+1) );
3012 
3013  // force at least 2 qubits in two's compement though
3014  if (encoding == TWOS_COMPLEMENT && numQubits == 1)
3015  numQubits++;
3016 
3017  // try every possible sub-register
3018  int* qubits = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numQubits) );
3019 
3020  // choose a random number of terms in the phase function
3021  int numTerms = getRandomInt(1, 21);
3022 
3023  // populate the phase function with random powers.
3024  // this includes negative powers, so we must always override 0.
3025  // in two's complement, we must use only integer powers
3026  qreal coeffs[numTerms];
3027  qreal expons[numTerms];
3028  for (int t=0; t<numTerms; t++) {
3029  coeffs[t] = getRandomReal(-10,10);
3030  if (encoding == TWOS_COMPLEMENT)
3031  expons[t] = getRandomInt(-3, 3+1); // note we COULD do getRandomReal(), and override all negatives
3032  else
3033  expons[t] = getRandomReal(-3, 3);
3034  }
3035 
3036  // choose a random number of overrides (even overriding every amplitude)
3037  int numOverrides = getRandomInt(1, (1<<numQubits) + 1);
3038 
3039  // randomise each override index (uniqueness isn't checked)
3040  long long int overrideInds[numOverrides];
3041  overrideInds[0] = 0LL;
3042  for (int i=1; i<numOverrides; i++)
3043  if (encoding == UNSIGNED)
3044  overrideInds[i] = getRandomInt(0, 1<<numQubits);
3045  else if (encoding == TWOS_COMPLEMENT)
3046  overrideInds[i] = getRandomInt(-(1<<(numQubits-1)), (1<<(numQubits-1))-1);
3047 
3048  // override to a random phase
3049  qreal overridePhases[numOverrides];
3050  for (int i=0; i<numOverrides; i++)
3051  overridePhases[i] = getRandomReal(-4, 4); // periodic in [-pi, pi]
3052 
3053  // build a reference diagonal matrix, on the reduced Hilbert space
3054  QMatrix matr = getZeroMatrix( 1 << numQubits );
3055  for (size_t i=0; i<matr.size(); i++) {
3056 
3057  long long int ind = 0;
3058  if (encoding == UNSIGNED)
3059  ind = i;
3060  if (encoding == TWOS_COMPLEMENT)
3061  ind = getTwosComplement(i, numQubits);
3062 
3063  // reference diagonal matrix incorporates overriden phases
3064  qreal phase;
3065  bool overriden = false;
3066  for (int v=0; v<numOverrides; v++) {
3067  if (ind == overrideInds[v]) {
3068  phase = overridePhases[v];
3069  overriden = true;
3070  break;
3071  }
3072  }
3073 
3074  if (!overriden) {
3075  phase = 0;
3076  for (int t=0; t<numTerms; t++)
3077  phase += coeffs[t] * pow(ind, expons[t]);
3078  }
3079 
3080  matr[i][i] = expI(phase);
3081  }
3082 
3083  SECTION( "state-vector" ) {
3084 
3085  applyPhaseFuncOverrides(quregVec, qubits, numQubits, encoding, coeffs, expons, numTerms, overrideInds, overridePhases, numOverrides);
3086  applyReferenceOp(refVec, qubits, numQubits, matr);
3087  REQUIRE( areEqual(quregVec, refVec, 1E4*REAL_EPS) );
3088  }
3089  SECTION( "density-matrix" ) {
3090 
3091  applyPhaseFuncOverrides(quregMatr, qubits, numQubits, encoding, coeffs, expons, numTerms, overrideInds, overridePhases, numOverrides);
3092  applyReferenceOp(refMatr, qubits, numQubits, matr);
3093  REQUIRE( areEqual(quregMatr, refMatr, 1E6*REAL_EPS) );
3094  }
3095  }
3096  SECTION( "input validation" ) {
3097 
3098  int numQubits = 3;
3099  int qubits[] = {0,1,2};
3100 
3101  SECTION( "number of qubits" ) {
3102 
3103  int numQubits = GENERATE_COPY( -1, 0, NUM_QUBITS+1 );
3104  REQUIRE_THROWS_WITH( applyPhaseFuncOverrides(quregVec, NULL, numQubits, UNSIGNED, NULL, NULL, 1, NULL, NULL, 0), Contains("Invalid number of qubits") );
3105  }
3106  SECTION( "repetition qubits" ) {
3107 
3108  qubits[GENERATE(1,2)] = qubits[0];
3109  REQUIRE_THROWS_WITH( applyPhaseFuncOverrides(quregVec, qubits, numQubits, UNSIGNED, NULL, NULL, 1, NULL, NULL, 0), Contains("The qubits must be unique") );
3110  }
3111  SECTION( "qubit indices" ) {
3112 
3113  int inv = GENERATE( -1, NUM_QUBITS );
3114  qubits[ GENERATE_COPY( range(0,numQubits) )] = inv;
3115  REQUIRE_THROWS_WITH( applyPhaseFuncOverrides(quregVec, qubits, numQubits, UNSIGNED, NULL, NULL, 1, NULL, NULL, 0), Contains("Invalid qubit index") );
3116  }
3117  SECTION( "number of terms" ) {
3118 
3119  int numTerms = GENERATE( -1, 0 );
3120  REQUIRE_THROWS_WITH( applyPhaseFuncOverrides(quregVec, qubits, numQubits, UNSIGNED, NULL, NULL, numTerms, NULL, NULL, 0), Contains("Invalid number of terms in the phase function") );
3121  }
3122  SECTION( "bit encoding name" ) {
3123 
3124  enum bitEncoding enc = (enum bitEncoding) GENERATE( -1,2 );
3125  REQUIRE_THROWS_WITH( applyPhaseFuncOverrides(quregVec, qubits, numQubits, enc, NULL, NULL, 1, NULL, NULL, 0), Contains("Invalid bit encoding") );
3126  }
3127  SECTION( "two's complement register" ) {
3128 
3129  numQubits = 1;
3130  REQUIRE_THROWS_WITH( applyPhaseFuncOverrides(quregVec, qubits, numQubits, TWOS_COMPLEMENT, NULL, NULL, 1, NULL, NULL, 0), Contains("too few qubits to employ TWOS_COMPLEMENT") );
3131  }
3132  SECTION( "number of overrides" ) {
3133 
3134  qreal dummyTerms[] = {0};
3135 
3136  int numOverrides = GENERATE_COPY( -1, 1 + (1<<numQubits) );
3137  REQUIRE_THROWS_WITH( applyPhaseFuncOverrides(quregVec, qubits, numQubits, UNSIGNED, dummyTerms, dummyTerms, 1, NULL, NULL, numOverrides), Contains("Invalid number of phase function overrides") );
3138  }
3139  SECTION( "override indices" ) {
3140 
3141  int numOverrides = 3;
3142  long long int overrideInds[] = {0,1,2};
3143  qreal overridePhases[] = {.1,.1,.1};
3144  qreal dummyTerms[] = {0};
3145 
3146  enum bitEncoding encoding = UNSIGNED;
3147  overrideInds[GENERATE(0,1,2)] = GENERATE_COPY( -1, (1<<numQubits) );
3148  REQUIRE_THROWS_WITH( applyPhaseFuncOverrides(quregVec, qubits, numQubits, encoding, dummyTerms, dummyTerms, 1, overrideInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the UNSIGNED encoding") );
3149 
3150  encoding = TWOS_COMPLEMENT;
3151  long long int newInds[] = {0,1,2};
3152  int minInd = -(1<<(numQubits-1));
3153  int maxInd = (1<<(numQubits-1)) -1;
3154  newInds[GENERATE(0,1,2)] = GENERATE_COPY( minInd-1, maxInd+1 );
3155  REQUIRE_THROWS_WITH( applyPhaseFuncOverrides(quregVec, qubits, numQubits, encoding, dummyTerms, dummyTerms, 1, newInds, overridePhases, numOverrides), Contains("Invalid phase function override index, in the TWOS_COMPLEMENT encoding") );
3156  }
3157  SECTION( "fractional exponent" ) {
3158 
3159  int numTerms = 3;
3160  qreal coeffs[] = {0,0,0};
3161  qreal expos[] = {1,2,3};
3162 
3163  // make one exponent fractional, thereby requiring negative overrides
3164  expos[GENERATE_COPY( range(0,numTerms) )] = GENERATE( 0.5, 1.999, 5.0001 );
3165 
3166  // catch when no negative indices are overridden
3167  REQUIRE_THROWS_WITH( applyPhaseFuncOverrides(quregVec, qubits, numQubits, TWOS_COMPLEMENT, coeffs, expos, numTerms, NULL, NULL, 0), Contains("fractional exponent") && Contains("TWOS_COMPLEMENT") && Contains("negative indices were not overriden") );
3168 
3169  int numNegs = 1 << (numQubits-1);
3170  long long int overrideInds[numNegs];
3171  qreal overridePhases[numNegs];
3172  for (int i=0; i<numNegs; i++) {
3173  overrideInds[i] = -(i+1);
3174  overridePhases[i] = 0;
3175  }
3176 
3177  // ensure no throw when all are overriden
3178  REQUIRE_NOTHROW( applyPhaseFuncOverrides(quregVec, qubits, numQubits, TWOS_COMPLEMENT, coeffs, expos, numTerms, overrideInds, overridePhases, numNegs) );
3179 
3180  // catch when at least one isn't overriden
3181  overrideInds[GENERATE_COPY( range(0,numNegs) )] = 0; // override a non-negative
3182  REQUIRE_THROWS_WITH( applyPhaseFuncOverrides(quregVec, qubits, numQubits, TWOS_COMPLEMENT, coeffs, expos, numTerms, overrideInds, overridePhases, numNegs), Contains("fractional exponent") && Contains("TWOS_COMPLEMENT") && Contains("negative indices were not overriden") );
3183  }
3184  SECTION( "negative exponent" ) {
3185 
3186  int numTerms = 3;
3187  qreal coeffs[] = {0,0,0};
3188  qreal expos[] = {1,2,3};
3189  expos[GENERATE_COPY( range(0,numTerms) )] = GENERATE( -1, -2 );
3190 
3191  enum bitEncoding encoding = GENERATE( UNSIGNED, TWOS_COMPLEMENT );
3192 
3193  // test both when giving no overrides, and giving all non-zero overrides
3194  int numOverrides = GENERATE( 0, 3 );
3195  long long int overrideInds[] = {1,2,3};
3196  qreal overridePhases[] = {0,0,0};
3197  REQUIRE_THROWS_WITH( applyPhaseFuncOverrides(quregVec, qubits, numQubits, encoding, coeffs, expos, numTerms, overrideInds, overridePhases, numOverrides), Contains("The phase function contained a negative exponent which would diverge at zero, but the zero index was not overriden") );
3198 
3199  // but ensure that when the zero IS overriden (anywhere), there's no error
3200  numOverrides = 3;
3201  overrideInds[GENERATE_COPY(range(0,numOverrides))] = 0;
3202  REQUIRE_NOTHROW( applyPhaseFuncOverrides(quregVec, qubits, numQubits, encoding, coeffs, expos, numTerms, overrideInds, overridePhases, numOverrides) );
3203  }
3204  }
3205  CLEANUP_TEST( quregVec, quregMatr );
3206 }
3207 
3208 
3209 
3214 TEST_CASE( "applyProjector", "[operators]" ) {
3215 
3216  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
3217  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
3218 
3219  SECTION( "correctness" ) {
3220 
3221  int qubit = GENERATE( range(0,NUM_QUBITS) );
3222  int outcome = GENERATE( 0, 1 );
3223 
3224  // repeat these random tests 10 times on every qubit, and for both outcomes
3225  GENERATE( range(0,10) );
3226 
3227  SECTION( "state-vector" ) {
3228 
3229  SECTION( "normalised" ) {
3230 
3231  // use a random L2 state for every qubit & outcome
3232  QVector vecRef = getRandomStateVector(NUM_QUBITS);
3233  toQureg(vec, vecRef);
3234 
3235  // zero non-outcome reference amps
3236  for (size_t ind=0; ind<vecRef.size(); ind++) {
3237  int bit = (ind >> qubit) & 1; // target-th bit
3238  if (bit != outcome)
3239  vecRef[ind] = 0;
3240  }
3241 
3242  applyProjector(vec, qubit, outcome);
3243  REQUIRE( areEqual(vec, vecRef) );
3244  }
3245  SECTION( "unnormalised" ) {
3246 
3247  // use a random non-physical state for every qubit & outcome
3248  QVector vecRef = getRandomQVector(1 << NUM_QUBITS);
3249  toQureg(vec, vecRef);
3250 
3251  // zero non-outcome reference amps
3252  for (size_t ind=0; ind<vecRef.size(); ind++) {
3253  int bit = (ind >> qubit) & 1; // target-th bit
3254  if (bit != outcome)
3255  vecRef[ind] = 0;
3256  }
3257 
3258  applyProjector(vec, qubit, outcome);
3259  REQUIRE( areEqual(vec, vecRef) );
3260  }
3261  }
3262  SECTION( "density-matrix" ) {
3263 
3264  SECTION( "pure" ) {
3265 
3266  QVector vecRef = getRandomStateVector(NUM_QUBITS);
3267  QMatrix matRef = getPureDensityMatrix(vecRef);
3268 
3269  toQureg(mat, matRef);
3270  applyProjector(mat, qubit, outcome);
3271 
3272  // zero any amplitudes that aren't |outcome><outcome|
3273  for (size_t r=0; r<matRef.size(); r++) {
3274  for (size_t c=0; c<matRef.size(); c++) {
3275  int ketBit = (c >> qubit) & 1;
3276  int braBit = (r >> qubit) & 1;
3277  if (!(ketBit == outcome && braBit == outcome))
3278  matRef[r][c] = 0;
3279  }
3280  }
3281 
3282  REQUIRE( areEqual(mat, matRef) );
3283  }
3284  SECTION( "mixed" ) {
3285 
3286  QMatrix matRef = getRandomDensityMatrix(NUM_QUBITS);
3287 
3288  toQureg(mat, matRef);
3289  applyProjector(mat, qubit, outcome);
3290 
3291  // zero any amplitudes that aren't |outcome><outcome|
3292  for (size_t r=0; r<matRef.size(); r++) {
3293  for (size_t c=0; c<matRef.size(); c++) {
3294  int ketBit = (c >> qubit) & 1;
3295  int braBit = (r >> qubit) & 1;
3296  if (!(ketBit == outcome && braBit == outcome))
3297  matRef[r][c] = 0;
3298  }
3299  }
3300 
3301  REQUIRE( areEqual(mat, matRef) );
3302  }
3303  SECTION( "unnormalised" ) {
3304 
3305  QMatrix matRef = getRandomQMatrix(1 << NUM_QUBITS);
3306 
3307  toQureg(mat, matRef);
3308  applyProjector(mat, qubit, outcome);
3309 
3310  // zero any amplitudes that aren't |outcome><outcome|
3311  for (size_t r=0; r<matRef.size(); r++) {
3312  for (size_t c=0; c<matRef.size(); c++) {
3313  int ketBit = (c >> qubit) & 1;
3314  int braBit = (r >> qubit) & 1;
3315  if (!(ketBit == outcome && braBit == outcome))
3316  matRef[r][c] = 0;
3317  }
3318  }
3319 
3320  REQUIRE( areEqual(mat, matRef) );
3321  }
3322  }
3323  }
3324  SECTION( "input validation" ) {
3325 
3326  SECTION( "qubit index" ) {
3327 
3328  int qubit = GENERATE( -1, NUM_QUBITS );
3329  int outcome = 0;
3330  REQUIRE_THROWS_WITH( applyProjector(mat, qubit, outcome), Contains("Invalid target qubit") );
3331  }
3332  SECTION( "outcome value" ) {
3333 
3334  int qubit = 0;
3335  int outcome = GENERATE( -1, 2 );
3336  REQUIRE_THROWS_WITH( applyProjector(mat, qubit, outcome), Contains("Invalid measurement outcome") );
3337  }
3338  }
3339  destroyQureg(vec, QUEST_ENV);
3340  destroyQureg(mat, QUEST_ENV);
3341 }
3342 
3343 
3344 
3349 TEST_CASE( "applyQFT", "[operators]" ) {
3350 
3351  PREPARE_TEST( quregVec, quregMatr, refVec, refMatr );
3352 
3353  SECTION( "correctness" ) {
3354 
3355  // try every sub-register size
3356  int numQubits = GENERATE_COPY( range(1,NUM_QUBITS+1) );
3357 
3358  // try every possible sub-register
3359  int* qubits = GENERATE_COPY( sublists(range(0,NUM_QUBITS), numQubits) );
3360 
3361  SECTION( "state-vector" ) {
3362 
3363  SECTION( "normalised" ) {
3364 
3365  QVector refVec = getRandomStateVector(NUM_QUBITS);
3366  toQureg(quregVec, refVec);
3367 
3368  applyQFT(quregVec, qubits, numQubits);
3369  refVec = getDFT(refVec, qubits, numQubits);
3370 
3371  REQUIRE( areEqual(quregVec, refVec) );
3372  }
3373  SECTION( "unnormalised" ) {
3374 
3375  QVector refVec = getRandomQVector(1 << NUM_QUBITS);
3376  toQureg(quregVec, refVec);
3377 
3378  applyQFT(quregVec, qubits, numQubits);
3379  refVec = getDFT(refVec, qubits, numQubits);
3380 
3381  REQUIRE( areEqual(quregVec, refVec) );
3382  }
3383  }
3384  SECTION( "density-matrix" ) {
3385 
3386  SECTION( "pure" ) {
3387 
3388  /* a pure density matrix should be mapped to a pure state
3389  * corresponding to the state-vector DFT
3390  */
3391 
3392  refVec = getRandomStateVector(NUM_QUBITS);
3393  refMatr = getPureDensityMatrix(refVec);
3394  toQureg(quregMatr, refMatr);
3395 
3396  applyQFT(quregMatr, qubits, numQubits);
3397  refVec = getDFT(refVec, qubits, numQubits);
3398  refMatr = getPureDensityMatrix(refVec);
3399 
3400  REQUIRE( areEqual(quregMatr, refMatr) );
3401  }
3402  SECTION( "mixed" ) {
3403 
3404  /* a mixed density matrix, conceptualised as a mixture of orthogonal
3405  * state-vectors, should be mapped to an equally weighted mixture
3406  * of DFTs of each state-vector (because QFT is unitary and hence
3407  * maintains state orthogonality)
3408  */
3409 
3410  int numStates = (1 << NUM_QUBITS)/4; // a quarter of as many states as are possible
3411  std::vector<QVector> states = getRandomOrthonormalVectors(NUM_QUBITS, numStates);
3412  std::vector<qreal> probs = getRandomProbabilities(numStates);
3413 
3414  // set qureg to random mixture of state-vectors
3415  refMatr = getMixedDensityMatrix(probs, states);
3416  toQureg(quregMatr, refMatr);
3417 
3418  // apply QFT to mixture
3419  applyQFT(quregMatr, qubits, numQubits);
3420 
3421  // compute dft of mixture, via dft of each state
3422  refMatr = getZeroMatrix(1 << NUM_QUBITS);
3423  for (int i=0; i<numStates; i++) {
3424  QVector dft = getDFT(states[i], qubits, numQubits);
3425  refMatr += probs[i] * getPureDensityMatrix(dft);
3426  }
3427 
3428  REQUIRE( areEqual(quregMatr, refMatr) );
3429  }
3430  SECTION( "unnormalised" ) {
3431 
3432  /* repeat method above, except that we use unnormalised vectors,
3433  * and mix them with arbitrary complex numbers instead of probabilities,
3434  * yielding an unnormalised density matrix
3435  */
3436 
3437  int numVecs = (1 << NUM_QUBITS)/4; // a quarter of as many states as are possible
3438  std::vector<QVector> vecs;
3439  std::vector<qcomp> coeffs;
3440  for (int i=0; i<numVecs; i++) {
3441  vecs.push_back(getRandomQVector(1 << NUM_QUBITS));
3442  coeffs.push_back(getRandomComplex());
3443  }
3444 
3445  // produce unnormalised matrix via random complex sum of random unnormalised vectors
3446  refMatr = getZeroMatrix(1 << NUM_QUBITS);
3447  for (int i=0; i<numVecs; i++)
3448  refMatr += coeffs[i] * getPureDensityMatrix(vecs[i]);
3449 
3450  toQureg(quregMatr, refMatr);
3451  applyQFT(quregMatr, qubits, numQubits);
3452 
3453  // compute target matrix via dft of each unnormalised vector
3454  refMatr = getZeroMatrix(1 << NUM_QUBITS);
3455  for (int i=0; i<numVecs; i++) {
3456  QVector dft = getDFT(vecs[i], qubits, numQubits);
3457  refMatr += coeffs[i] * getPureDensityMatrix(dft);
3458  }
3459 
3460  REQUIRE( areEqual(quregMatr, refMatr) );
3461  }
3462  }
3463  }
3464  SECTION( "input validation" ) {
3465 
3466  SECTION( "number of targets" ) {
3467 
3468  // there cannot be more targets than qubits in register
3469  int numQubits = GENERATE( -1, 0, NUM_QUBITS+1 );
3470  int qubits[NUM_QUBITS+1];
3471 
3472  REQUIRE_THROWS_WITH( applyQFT(quregVec, qubits, numQubits), Contains("Invalid number of target"));
3473  }
3474  SECTION( "repetition in targets" ) {
3475 
3476  int numQubits = 3;
3477  int qubits[] = {1,2,2};
3478 
3479  REQUIRE_THROWS_WITH( applyQFT(quregVec, qubits, numQubits), Contains("target") && Contains("unique"));
3480  }
3481  SECTION( "qubit indices" ) {
3482 
3483  int numQubits = 3;
3484  int qubits[] = {1,2,3};
3485 
3486  int inv = GENERATE( -1, NUM_QUBITS );
3487  qubits[GENERATE_COPY( range(0,numQubits) )] = inv; // make invalid target
3488  REQUIRE_THROWS_WITH( applyQFT(quregVec, qubits, numQubits), Contains("Invalid target") );
3489  }
3490  }
3491  CLEANUP_TEST( quregVec, quregMatr );
3492 }
3493 
3494 
3495 
3500 TEST_CASE( "applyTrotterCircuit", "[operators]" ) {
3501 
3502  Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
3503  Qureg mat = createDensityQureg(NUM_QUBITS, QUEST_ENV);
3504  initDebugState(vec);
3505  initDebugState(mat);
3506 
3507  Qureg vecRef = createCloneQureg(vec, QUEST_ENV);
3508  Qureg matRef = createCloneQureg(mat, QUEST_ENV);
3509 
3510  SECTION( "correctness" ) {
3511 
3512  SECTION( "one term" ) {
3513 
3514  // a Hamiltonian with one term has an exact (trivial) Trotterisation
3515  PauliHamil hamil = createPauliHamil(NUM_QUBITS, 1);
3516 
3517  // H = coeff X Y Z (on qubits 0,1,2)
3518  qreal coeff = getRandomReal(-5, 5);
3519  int numTargs = 3;
3520  int targs[] = {0, 1, 2};
3521  pauliOpType codes[] = {PAULI_X, PAULI_Y, PAULI_Z};
3522  hamil.termCoeffs[0] = coeff;
3523  for (int i=0; i<numTargs; i++)
3524  hamil.pauliCodes[targs[i]] = codes[i];
3525 
3526  // time can be negative
3527  qreal time = getRandomReal(-2,2);
3528 
3529  // by commutation, all reps & orders yield the same total unitary
3530  int reps = GENERATE( range(1,5) );
3531 
3532  // applyTrotter(t, 1, 1) = exp(-i t coeff paulis)
3533  // multiRotatePauli(a) = exp(- i a / 2 paulis)
3534 
3535  SECTION( "state-vector" ) {
3536 
3537  int order = GENERATE( 1, 2, 4 );
3538 
3539  applyTrotterCircuit(vec, hamil, time, order, reps);
3540  multiRotatePauli(vecRef, targs, codes, numTargs, 2*time*coeff);
3541  REQUIRE( areEqual(vec, vecRef, 10*REAL_EPS) );
3542  }
3543  SECTION( "density-matrix" ) {
3544 
3545  int order = GENERATE( 1, 2 ); // precision bites density-matrices quickly
3546 
3547  applyTrotterCircuit(mat, hamil, time, order, reps);
3548  multiRotatePauli(matRef, targs, codes, numTargs, 2*time*coeff);
3549  REQUIRE( areEqual(mat, matRef, 1E2*REAL_EPS) );
3550  }
3551 
3552  destroyPauliHamil(hamil);
3553  }
3554  SECTION( "commuting terms" ) {
3555 
3556  // a Hamiltonian of commuting terms, Trotterises exactly
3557  PauliHamil hamil = createPauliHamil(NUM_QUBITS, 2);
3558 
3559  // H = c0 X Y I + c1 X I Z (on qubits 0,1,2)
3560  int targs[] = {0, 1, 2};
3561  hamil.pauliCodes[0] = PAULI_X;
3562  hamil.pauliCodes[0 + NUM_QUBITS] = PAULI_X;
3563  hamil.pauliCodes[1] = PAULI_Y;
3564  hamil.pauliCodes[1 + NUM_QUBITS] = PAULI_I;
3565  hamil.pauliCodes[2] = PAULI_I;
3566  hamil.pauliCodes[2 + NUM_QUBITS] = PAULI_Z;
3567  for (int i=0; i<hamil.numSumTerms; i++)
3568  hamil.termCoeffs[i] = getRandomReal(-5,5);
3569 
3570  // time can be negative
3571  qreal time = getRandomReal(-2,2);
3572 
3573  // applyTrotter(t, 1, 1) = exp(-i t c0 paulis0) exp(-i t c1 paulis1)
3574  // multiRotatePauli(a) = exp(- i a / 2 paulis)
3575 
3576  SECTION( "state-vector" ) {
3577 
3578  int reps = GENERATE( range(1,5) );
3579  int order = GENERATE( 1, 2, 4 );
3580 
3581  applyTrotterCircuit(vec, hamil, time, order, reps);
3582  multiRotatePauli(vecRef, targs, hamil.pauliCodes, 3, 2*time*hamil.termCoeffs[0]);
3583  multiRotatePauli(vecRef, targs, &(hamil.pauliCodes[NUM_QUBITS]), 3, 2*time*hamil.termCoeffs[1]);
3584  REQUIRE( areEqual(vec, vecRef, 10*REAL_EPS) );
3585  }
3586  SECTION( "density-matrix" ) {
3587 
3588  int reps = GENERATE( range(1,5) );
3589  int order = GENERATE( 1, 2 ); // precision hurts density matrices quickly
3590 
3591  applyTrotterCircuit(mat, hamil, time, order, reps);
3592  multiRotatePauli(matRef, targs, hamil.pauliCodes, 3, 2*time*hamil.termCoeffs[0]);
3593  multiRotatePauli(matRef, targs, &(hamil.pauliCodes[NUM_QUBITS]), 3, 2*time*hamil.termCoeffs[1]);
3594  REQUIRE( areEqual(mat, matRef, 1E2*REAL_EPS) );
3595  }
3596 
3597  destroyPauliHamil(hamil);
3598  }
3599  SECTION( "general" ) {
3600 
3601  /* We'll consider an analytic time-evolved state, so that we can avoid
3602  * comparing applyTrotterCircuit to other numerical approximations.
3603  * We can construct such a state, by using a Hamiltonian with known
3604  * analytic eigenvalues, and hence a known period. Time evolution of the
3605  * period will just yield the input state.
3606  *
3607  * E.g. H = pi sqrt(2) X Y Z X Y + pi Y Z X Y Z + pi Z X Y Z X
3608  * has (degenerate) eigenvalues +- 2 pi, so the period
3609  * of the Hamiltonian is t=1.
3610  */
3611 
3612  // hardcoded 5 qubits here in the Pauli codes
3613  REQUIRE( NUM_QUBITS == 5 );
3614 
3615  PauliHamil hamil = createPauliHamil(NUM_QUBITS, 3);
3616  qreal coeffs[] = {(qreal) (M_PI * sqrt(2.0)), M_PI, M_PI};
3617  enum pauliOpType codes[] = {
3618  PAULI_X, PAULI_Y, PAULI_Z, PAULI_X, PAULI_Y,
3619  PAULI_Y, PAULI_Z, PAULI_X, PAULI_Y, PAULI_Z,
3620  PAULI_Z, PAULI_X, PAULI_Y, PAULI_Z, PAULI_X};
3621  initPauliHamil(hamil, coeffs, codes);
3622 
3623  // evolving to t=1 should leave the input state unchanged
3624  qreal time = 1;
3625 
3626  // since unnormalised (initDebugState), max fid is 728359.8336
3627  qreal fidNorm = 728359.8336;
3628 
3629  SECTION( "absolute" ) {
3630 
3631  // such a high order and reps should yield precise solution
3632  int order = 4;
3633  int reps = 20;
3634  applyTrotterCircuit(vec, hamil, time, 4, 20);
3635  qreal fid = calcFidelity(vec, vecRef) / fidNorm;
3636 
3637  REQUIRE( fid == Approx(1).epsilon(1E-8) );
3638  }
3639  SECTION( "repetitions scaling" ) {
3640 
3641  // exclude order 1; too few reps for monotonic increase of accuracy
3642  int order = GENERATE( 2, 4, 6 );
3643 
3644  // accuracy should increase with increasing repetitions
3645  int reps[] = {1, 5, 10};
3646 
3647  qreal prevFid = 0;
3648  for (int i=0; i<3; i++) {
3649  initDebugState(vec);
3650  applyTrotterCircuit(vec, hamil, time, order, reps[i]);
3651  qreal fid = calcFidelity(vec, vecRef) / fidNorm;
3652 
3653  REQUIRE( fid >= prevFid );
3654  prevFid = fid;
3655  }
3656  }
3657  SECTION( "order scaling" ) {
3658 
3659  // exclude order 1; too few reps for monotonic increase of accuracy
3660  int reps = GENERATE( 5, 10 );
3661 
3662  // accuracy should increase with increasing repetitions
3663  int orders[] = {1, 2, 4, 6};
3664 
3665  qreal prevFid = 0;
3666  for (int i=0; i<4; i++) {
3667  initDebugState(vec);
3668  applyTrotterCircuit(vec, hamil, time, orders[i], reps);
3669  qreal fid = calcFidelity(vec, vecRef) / fidNorm;
3670 
3671  REQUIRE( fid >= prevFid );
3672  prevFid = fid;
3673  }
3674  }
3675 
3676  destroyPauliHamil(hamil);
3677  }
3678  }
3679  SECTION( "input validation" ) {
3680 
3681  SECTION( "repetitions" ) {
3682 
3683  PauliHamil hamil = createPauliHamil(NUM_QUBITS, 1);
3684  int reps = GENERATE( -1, 0 );
3685 
3686  REQUIRE_THROWS_WITH( applyTrotterCircuit(vec, hamil, 1, 1, reps), Contains("repetitions must be >=1") );
3687 
3688  destroyPauliHamil(hamil);
3689  }
3690  SECTION( "order" ) {
3691 
3692  PauliHamil hamil = createPauliHamil(NUM_QUBITS, 1);
3693  int order = GENERATE( -1, 0, 3, 5, 7 );
3694 
3695  REQUIRE_THROWS_WITH( applyTrotterCircuit(vec, hamil, 1, order, 1), Contains("order must be 1, or an even number") );
3696 
3697  destroyPauliHamil(hamil);
3698  }
3699  SECTION( "pauli codes" ) {
3700 
3701  int numTerms = 3;
3702  PauliHamil hamil = createPauliHamil(NUM_QUBITS, numTerms);
3703 
3704  // make one pauli code wrong
3705  hamil.pauliCodes[GENERATE_COPY( range(0,numTerms*NUM_QUBITS) )] = (pauliOpType) GENERATE( -1, 4 );
3706  REQUIRE_THROWS_WITH( applyTrotterCircuit(vec, hamil, 1, 1, 1), Contains("Invalid Pauli code") );
3707 
3708  destroyPauliHamil(hamil);
3709  }
3710  SECTION( "matching hamiltonian qubits" ) {
3711 
3712  PauliHamil hamil = createPauliHamil(NUM_QUBITS + 1, 1);
3713 
3714  REQUIRE_THROWS_WITH( applyTrotterCircuit(vec, hamil, 1, 1, 1), Contains("same number of qubits") );
3715  REQUIRE_THROWS_WITH( applyTrotterCircuit(mat, hamil, 1, 1, 1), Contains("same number of qubits") );
3716 
3717  destroyPauliHamil(hamil);
3718  }
3719  }
3720 
3721  destroyQureg(vec, QUEST_ENV);
3722  destroyQureg(mat, QUEST_ENV);
3723  destroyQureg(vecRef, QUEST_ENV);
3724  destroyQureg(matRef, QUEST_ENV);
3725 }
3726 
INVERSE_PRODUCT
@ INVERSE_PRODUCT
Definition: QuEST.h:233
pauliOpType
pauliOpType
Codes for specifying Pauli operators.
Definition: QuEST.h:96
applyPhaseFunc
void applyPhaseFunc(Qureg qureg, int *qubits, int numQubits, enum bitEncoding encoding, qreal *coeffs, qreal *exponents, int numTerms)
Induces a phase change upon each amplitude of qureg, determined by the passed exponential polynomial ...
Definition: QuEST.c:726
applyMultiControlledMatrixN
void applyMultiControlledMatrixN(Qureg qureg, int *ctrls, int numCtrls, int *targs, int numTargs, ComplexMatrixN u)
Apply a general N-by-N matrix, which may be non-unitary, with additional controlled qubits.
Definition: QuEST.c:1114
applyNamedPhaseFuncOverrides
void applyNamedPhaseFuncOverrides(Qureg qureg, int *qubits, int *numQubitsPerReg, int numRegs, enum bitEncoding encoding, enum phaseFunc functionNameCode, long long int *overrideInds, qreal *overridePhases, int numOverrides)
Induces a phase change upon each amplitude of qureg, determined by a named (and potentially multi-var...
Definition: QuEST.c:813
initPauliHamil
void initPauliHamil(PauliHamil hamil, qreal *coeffs, enum pauliOpType *codes)
Initialise PauliHamil instance hamil with the given term coefficients and Pauli codes (one for every ...
Definition: QuEST.c:1504
QUEST_ENV
QuESTEnv QUEST_ENV
The global QuESTEnv instance, to be created and destroyed once in this main(), so that the MPI enviro...
Definition: main.cpp:20
applyProjector
void applyProjector(Qureg qureg, int qubit, int outcome)
Force the target qubit of qureg into the given classical outcome, via a non-renormalising projection.
Definition: QuEST.c:888
PAULI_Z
@ PAULI_Z
Definition: QuEST.h:96
getKroneckerProduct
QVector getKroneckerProduct(QVector b, QVector a)
Returns b (otimes) a.
Definition: utilities.cpp:143
DISTANCE
@ DISTANCE
Definition: QuEST.h:234
destroyComplexMatrixN
void destroyComplexMatrixN(ComplexMatrixN m)
Destroy a ComplexMatrixN instance created with createComplexMatrixN()
Definition: QuEST.c:1369
applyNamedPhaseFunc
void applyNamedPhaseFunc(Qureg qureg, int *qubits, int *numQubitsPerReg, int numRegs, enum bitEncoding encoding, enum phaseFunc functionNameCode)
Induces a phase change upon each amplitude of qureg, determined by a named (and potentially multi-var...
Definition: QuEST.c:796
PAULI_I
@ PAULI_I
Definition: QuEST.h:96
getRandomInt
int getRandomInt(int min, int max)
Returns a random integer between min (inclusive) and max (exclusive), from the uniform distribution.
Definition: utilities.cpp:526
createComplexMatrixN
ComplexMatrixN createComplexMatrixN(int numQubits)
Allocate dynamic memory for a square complex matrix of any size, which can be passed to functions lik...
Definition: QuEST.c:1348
getMixedDensityMatrix
QMatrix getMixedDensityMatrix(std::vector< qreal > probs, std::vector< QVector > states)
Returns a mixed density matrix formed from mixing the given pure states, which are assumed normalised...
Definition: utilities.cpp:640
syncDiagonalOp
void syncDiagonalOp(DiagonalOp op)
Update the GPU memory with the current values in op.real and op.imag.
Definition: QuEST.c:1531
destroyDiagonalOp
void destroyDiagonalOp(DiagonalOp op, QuESTEnv env)
Destroys a DiagonalOp created with createDiagonalOp(), freeing its memory.
Definition: QuEST.c:1524
NUM_QUBITS
#define NUM_QUBITS
The default number of qubits in the registers created for unit testing (both statevectors and density...
Definition: utilities.hpp:36
CLEANUP_TEST
#define CLEANUP_TEST(quregVec, quregMatr)
Destroys the data structures made by PREPARE_TEST.
Definition: test_operators.cpp:19
TWOS_COMPLEMENT
@ TWOS_COMPLEMENT
Definition: QuEST.h:269
calcFidelity
qreal calcFidelity(Qureg qureg, Qureg pureState)
Calculates the fidelity of qureg (a state-vector or density matrix) against a reference pure state (n...
Definition: QuEST.c:1191
getRandomReal
qreal getRandomReal(qreal min, qreal max)
Returns a random real between min (inclusive) and max (exclusive), from the uniform distribution.
Definition: utilities.cpp:421
NORM
@ NORM
Definition: QuEST.h:232
setDiagMatrixOverrides
void setDiagMatrixOverrides(QMatrix &matr, int *numQubitsPerReg, int numRegs, enum bitEncoding encoding, long long int *overrideInds, qreal *overridePhases, int numOverrides)
Modifies the given diagonal matrix such that the diagonal elements which correspond to the coordinate...
Definition: utilities.cpp:1316
getRandomOrthonormalVectors
std::vector< QVector > getRandomOrthonormalVectors(int numQb, int numStates)
Returns a list of random orthonormal complex vectors, from an undisclosed distribution.
Definition: utilities.cpp:613
ComplexMatrix4
Represents a 4x4 matrix of complex numbers.
Definition: QuEST.h:175
SCALED_INVERSE_DISTANCE
@ SCALED_INVERSE_DISTANCE
Definition: QuEST.h:234
applyPauliHamil
void applyPauliHamil(Qureg inQureg, PauliHamil hamil, Qureg outQureg)
Modifies outQureg to be the result of applying PauliHamil (a Hermitian but not necessarily unitary op...
Definition: QuEST.c:1059
applyQFT
void applyQFT(Qureg qureg, int *qubits, int numQubits)
Applies the quantum Fourier transform (QFT) to a specific subset of qubits of the register qureg.
Definition: QuEST.c:866
UNSIGNED
@ UNSIGNED
Definition: QuEST.h:269
applyPhaseFuncOverrides
void applyPhaseFuncOverrides(Qureg qureg, int *qubits, int numQubits, enum bitEncoding encoding, qreal *coeffs, qreal *exponents, int numTerms, long long int *overrideInds, qreal *overridePhases, int numOverrides)
Induces a phase change upon each amplitude of qureg, determined by the passed exponential polynomial ...
Definition: QuEST.c:743
ComplexMatrixN
Represents a general 2^N by 2^N matrix of complex numbers.
Definition: QuEST.h:186
INVERSE_DISTANCE
@ INVERSE_DISTANCE
Definition: QuEST.h:234
qreal
#define qreal
multiRotatePauli
void multiRotatePauli(Qureg qureg, int *targetQubits, enum pauliOpType *targetPaulis, int numTargets, qreal angle)
Apply a multi-qubit multi-Pauli rotation, also known as a Pauli gadget, on a selected number of qubit...
Definition: QuEST.c:685
toQMatrix
QMatrix toQMatrix(ComplexMatrix2 src)
Returns a copy of the given 2-by-2 matrix.
Definition: utilities.cpp:1044
toQureg
void toQureg(Qureg qureg, QVector vec)
Initialises the state-vector qureg to have the same amplitudes as vec.
Definition: utilities.cpp:1201
calcLog2
unsigned int calcLog2(long unsigned int num)
returns log2 of numbers which must be gauranteed to be 2^n
Definition: QuEST_validation.c:328
PAULI_X
@ PAULI_X
Definition: QuEST.h:96
TEST_CASE
TEST_CASE("applyDiagonalOp", "[operators]")
Definition: test_operators.cpp:32
setRandomPauliSum
void setRandomPauliSum(qreal *coeffs, pauliOpType *codes, int numQubits, int numTerms)
Populates the coeffs array with random qreals in (-5, 5), and populates codes with random Pauli codes...
Definition: utilities.cpp:1229
toComplexMatrixN
void toComplexMatrixN(QMatrix qm, ComplexMatrixN cm)
Initialises cm with the values of qm.
Definition: utilities.cpp:1033
expI
qcomp expI(qreal phase)
Returns the unit-norm complex number exp(i*phase).
Definition: utilities.cpp:417
QVector
std::vector< qcomp > QVector
A complex vector, which can be zero-initialised with QVector(numAmps).
Definition: utilities.hpp:60
DiagonalOp::imag
qreal * imag
The imaginary values of the 2^numQubits complex elements.
Definition: QuEST.h:310
phaseFunc
phaseFunc
Flags for specifying named phase functions.
Definition: QuEST.h:231
toQVector
QVector toQVector(Qureg qureg)
Returns an equal-size copy of the given state-vector qureg.
Definition: utilities.cpp:1113
applyDiagonalOp
void applyDiagonalOp(Qureg qureg, DiagonalOp op)
Apply a diagonal operator, which is possibly non-unitary and non-Hermitian, to the entire qureg.
Definition: QuEST.c:1127
PauliHamil::termCoeffs
qreal * termCoeffs
The real coefficient of each Pauli product. This is an array of length PauliHamil....
Definition: QuEST.h:283
applyMatrix4
void applyMatrix4(Qureg qureg, int targetQubit1, int targetQubit2, ComplexMatrix4 u)
Apply a general 4-by-4 matrix, which may be non-unitary.
Definition: QuEST.c:1093
PauliHamil::pauliCodes
enum pauliOpType * pauliCodes
The Pauli operators acting on each qubit, flattened over every operator.
Definition: QuEST.h:281
toComplexMatrix4
ComplexMatrix4 toComplexMatrix4(QMatrix qm)
Returns a ComplexMatrix4 copy of QMatix qm.
Definition: utilities.cpp:1027
SCALED_PRODUCT
@ SCALED_PRODUCT
Definition: QuEST.h:233
MAX_NUM_REGS_APPLY_ARBITRARY_PHASE
#define MAX_NUM_REGS_APPLY_ARBITRARY_PHASE
Definition: QuEST_precision.h:72
PauliHamil::numSumTerms
int numSumTerms
The number of terms in the weighted sum, or the number of Pauli products.
Definition: QuEST.h:285
DiagonalOp
Represents a diagonal complex operator on the full Hilbert state of a Qureg.
Definition: QuEST.h:297
PAULI_Y
@ PAULI_Y
Definition: QuEST.h:96
PauliHamil
A Pauli Hamiltonian, expressed as a real-weighted sum of pauli products, and which can hence represen...
Definition: QuEST.h:277
destroyQureg
void destroyQureg(Qureg qureg, QuESTEnv env)
Deallocate a Qureg, freeing its memory.
Definition: QuEST.c:77
applyParamNamedPhaseFuncOverrides
void applyParamNamedPhaseFuncOverrides(Qureg qureg, int *qubits, int *numQubitsPerReg, int numRegs, enum bitEncoding encoding, enum phaseFunc functionNameCode, qreal *params, int numParams, long long int *overrideInds, qreal *overridePhases, int numOverrides)
Induces a phase change upon each amplitude of qureg, determined by a named, parameterised (and potent...
Definition: QuEST.c:848
SCALED_INVERSE_SHIFTED_NORM
@ SCALED_INVERSE_SHIFTED_NORM
Definition: QuEST.h:232
getRandomStateVector
QVector getRandomStateVector(int numQb)
Returns a random numQb-length L2-normalised state-vector from an undisclosed distribution.
Definition: utilities.cpp:468
ComplexMatrixN::real
qreal ** real
Definition: QuEST.h:189
initDebugState
void initDebugState(Qureg qureg)
Initialises qureg to be in the un-normalised, non-physical state with with -th complex amplitude give...
Definition: QuEST.c:1578
applyPauliSum
void applyPauliSum(Qureg inQureg, enum pauliOpType *allPauliCodes, qreal *termCoeffs, int numSumTerms, Qureg outQureg)
Modifies outQureg to be the result of applying the weighted sum of Pauli products (a Hermitian but no...
Definition: QuEST.c:1048
getRandomQMatrix
QMatrix getRandomQMatrix(int dim)
Returns a dim-by-dim complex matrix, where the real and imaginary value of each element are independe...
Definition: utilities.cpp:379
applyMatrixN
void applyMatrixN(Qureg qureg, int *targs, int numTargs, ComplexMatrixN u)
Apply a general N-by-N matrix, which may be non-unitary, on any number of target qubits.
Definition: QuEST.c:1103
applyMultiVarPhaseFunc
void applyMultiVarPhaseFunc(Qureg qureg, int *qubits, int *numQubitsPerReg, int numRegs, enum bitEncoding encoding, qreal *coeffs, qreal *exponents, int *numTermsPerReg)
Induces a phase change upon each amplitude of qureg, determined by a multi-variable exponential polyn...
Definition: QuEST.c:761
INVERSE_NORM
@ INVERSE_NORM
Definition: QuEST.h:232
getTwosComplement
long long int getTwosComplement(long long int decimal, int numBits)
Returns the two's complement signed encoding of the unsigned number decimal, which must be a number b...
Definition: utilities.cpp:1286
Qureg
Represents a system of qubits.
Definition: QuEST.h:322
ComplexMatrixN::imag
qreal ** imag
Definition: QuEST.h:190
PRODUCT
@ PRODUCT
Definition: QuEST.h:233
applyMultiVarPhaseFuncOverrides
void applyMultiVarPhaseFuncOverrides(Qureg qureg, int *qubits, int *numQubitsPerReg, int numRegs, enum bitEncoding encoding, qreal *coeffs, qreal *exponents, int *numTermsPerReg, long long int *overrideInds, qreal *overridePhases, int numOverrides)
Induces a phase change upon each amplitude of qureg, determined by a multi-variable exponential polyn...
Definition: QuEST.c:778
DiagonalOp::numElemsPerChunk
long long int numElemsPerChunk
The number of the 2^numQubits amplitudes stored on each distributed node.
Definition: QuEST.h:302
createCloneQureg
Qureg createCloneQureg(Qureg qureg, QuESTEnv env)
Create a new Qureg which is an exact clone of the passed qureg, which can be either a state-vector or...
Definition: QuEST.c:64
QMatrix
std::vector< std::vector< qcomp > > QMatrix
A complex square matrix.
Definition: utilities.hpp:49
toComplexMatrix2
ComplexMatrix2 toComplexMatrix2(QMatrix qm)
Returns a ComplexMatrix2 copy of QMatix qm.
Definition: utilities.cpp:1021
applyFullQFT
void applyFullQFT(Qureg qureg)
Applies the quantum Fourier transform (QFT) to the entirety of qureg.
Definition: QuEST.c:876
sublists
Catch::Generators::GeneratorWrapper< int * > sublists(int *list, int len, int sublen)
Returns a Catch2 generator of every length-sublen sublist of length-len list, in increasing lexograph...
Definition: utilities.cpp:1488
getRandomProbabilities
std::vector< qreal > getRandomProbabilities(int numProbs)
Returns a list of random real scalars, each in [0, 1], which sum to unity.
Definition: utilities.cpp:472
SCALED_DISTANCE
@ SCALED_DISTANCE
Definition: QuEST.h:234
DiagonalOp::real
qreal * real
The real values of the 2^numQubits complex elements.
Definition: QuEST.h:308
createQureg
Qureg createQureg(int numQubits, QuESTEnv env)
Creates a state-vector Qureg object representing a set of qubits which will remain in a pure state.
Definition: QuEST.c:36
destroyPauliHamil
void destroyPauliHamil(PauliHamil h)
Destroy a PauliHamil instance, created with either createPauliHamil() or createPauliHamilFromFile().
Definition: QuEST.c:1414
getPureDensityMatrix
QMatrix getPureDensityMatrix(QVector state)
Returns a density matrix initialised into the given pure state.
Definition: utilities.cpp:507
getRandomDensityMatrix
QMatrix getRandomDensityMatrix(int numQb)
Returns a random numQb-by-numQb density matrix, from an undisclosed distribution, in a very mixed sta...
Definition: utilities.cpp:490
applyParamNamedPhaseFunc
void applyParamNamedPhaseFunc(Qureg qureg, int *qubits, int *numQubitsPerReg, int numRegs, enum bitEncoding encoding, enum phaseFunc functionNameCode, qreal *params, int numParams)
Induces a phase change upon each amplitude of qureg, determined by a named, paramaterized (and potent...
Definition: QuEST.c:831
SCALED_INVERSE_SHIFTED_DISTANCE
@ SCALED_INVERSE_SHIFTED_DISTANCE
Definition: QuEST.h:234
getZeroMatrix
QMatrix getZeroMatrix(size_t dim)
Returns a dim-by-dim square complex matrix, initialised to all zeroes.
Definition: utilities.cpp:153
applyMatrix2
void applyMatrix2(Qureg qureg, int targetQubit, ComplexMatrix2 u)
Apply a general 2-by-2 matrix, which may be non-unitary.
Definition: QuEST.c:1084
PREPARE_TEST
#define PREPARE_TEST(quregVec, quregMatr, refVec, refMatr)
Prepares the needed data structures for unit testing some operators.
Definition: test_operators.cpp:10
getRandomQVector
QVector getRandomQVector(int dim)
Returns a dim-length vector with random complex amplitudes in the square joining {-1-i,...
Definition: utilities.cpp:435
SCALED_NORM
@ SCALED_NORM
Definition: QuEST.h:232
SCALED_INVERSE_PRODUCT
@ SCALED_INVERSE_PRODUCT
Definition: QuEST.h:233
applyReferenceMatrix
void applyReferenceMatrix(QVector &state, int *ctrls, int numCtrls, int *targs, int numTargs, QMatrix op)
Modifies the state-vector state to be the result of left-multiplying the multi-target operator matrix...
Definition: utilities.cpp:841
getDFT
QVector getDFT(QVector in)
Returns the discrete fourier transform of vector in.
Definition: utilities.cpp:652
areEqual
bool areEqual(QVector a, QVector b)
Returns true if the absolute value of the difference between every amplitude in vectors a and b is le...
Definition: utilities.cpp:398
createDensityQureg
Qureg createDensityQureg(int numQubits, QuESTEnv env)
Creates a density matrix Qureg object representing a set of qubits which can enter noisy and mixed st...
Definition: QuEST.c:50
applyReferenceOp
void applyReferenceOp(QVector &state, int *ctrls, int numCtrls, int *targs, int numTargs, QMatrix op)
Modifies the state-vector state to be the result of applying the multi-target operator matrix op,...
Definition: utilities.cpp:728
M_PI
#define M_PI
Definition: QuEST_common.c:41
applyTrotterCircuit
void applyTrotterCircuit(Qureg qureg, PauliHamil hamil, qreal time, int order, int reps)
Applies a trotterisation of unitary evolution to qureg.
Definition: QuEST.c:1070
createPauliHamil
PauliHamil createPauliHamil(int numQubits, int numSumTerms)
Dynamically allocates a Hamiltonian expressed as a real-weighted sum of products of Pauli operators.
Definition: QuEST.c:1398
getRandomComplex
qcomp getRandomComplex()
Returns a random complex number within the square closing (-1-i) and (1+i), from a distribution unifo...
Definition: utilities.cpp:431
bitEncoding
bitEncoding
Flags for specifying how the bits in sub-register computational basis states are mapped to indices in...
Definition: QuEST.h:269
createDiagonalOp
DiagonalOp createDiagonalOp(int numQubits, QuESTEnv env)
Creates a DiagonalOp representing a diagonal operator on the full Hilbert space of a Qureg.
Definition: QuEST.c:1518
ComplexMatrix2
Represents a 2x2 matrix of complex numbers.
Definition: QuEST.h:137
SCALED_INVERSE_NORM
@ SCALED_INVERSE_NORM
Definition: QuEST.h:232