1 /* gates.c: Basic gates for quantum register manipulation
2 
3    Copyright 2003, 2004 Bjoern Butscher, Hendrik Weimer
4 
5    This file is part of libquantum
6 
7    libquantum is free software; you can redistribute it and/or modify
8    it under the terms of the GNU General Public License as published
9    by the Free Software Foundation; either version 3 of the License,
10    or (at your option) any later version.
11 
12    libquantum is distributed in the hope that it will be useful, but
13    WITHOUT ANY WARRANTY; without even the implied warranty of
14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15    General Public License for more details.
16 
17    You should have received a copy of the GNU General Public License
18    along with libquantum; if not, write to the Free Software
19    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
20    MA 02110-1301, USA
21 
22 */
23 
24 #include <stdlib.h>
25 #include <stdio.h>
26 #include <math.h>
27 #include <stdarg.h>
28 
29 #include "matrix.h"
30 #include "defs.h"
31 #include "complex.h"
32 #include "qureg.h"
33 #include "decoherence.h"
34 #include "qec.h"
35 #include "objcode.h"
36 #include "error.h"
37 
38 /* Apply a controlled-not gate */
39 
40 void
quantum_cnot(int control,int target,quantum_reg * reg)41 quantum_cnot(int control, int target, quantum_reg *reg)
42 {
43   int i;
44   int qec;
45 
46   quantum_qec_get_status(&qec, NULL);
47 
48   if(qec)
49     quantum_cnot_ft(control, target, reg);
50   else
51     {
52       if(quantum_objcode_put(CNOT, control, target))
53 	return;
54 
55 #ifdef _OPENMP
56 #pragma omp parallel for
57 #endif
58       for(i=0; i<reg->size; i++)
59 	{
60 	  /* Flip the target bit of a basis state if the control bit is set */
61 
62 	  if((reg->state[i] & ((MAX_UNSIGNED) 1 << control)))
63 	    reg->state[i] ^= ((MAX_UNSIGNED) 1 << target);
64 	}
65       quantum_decohere(reg);
66     }
67 }
68 
69 /* Apply a toffoli (or controlled-controlled-not) gate */
70 
71 void
quantum_toffoli(int control1,int control2,int target,quantum_reg * reg)72 quantum_toffoli(int control1, int control2, int target, quantum_reg *reg)
73 {
74   int i;
75   int qec;
76 
77   quantum_qec_get_status(&qec, NULL);
78 
79   if(qec)
80     quantum_toffoli_ft(control1, control2, target, reg);
81   else
82     {
83       if(quantum_objcode_put(TOFFOLI, control1, control2, target))
84 	return;
85 
86 #ifdef _OPENMP
87 #pragma omp parallel for
88 #endif
89       for(i=0; i<reg->size; i++)
90 	{
91 	  /* Flip the target bit of a basis state if both control bits are
92 	     set */
93 
94 	  if(reg->state[i] & ((MAX_UNSIGNED) 1 << control1))
95 	    {
96 	      if(reg->state[i] & ((MAX_UNSIGNED) 1 << control2))
97 		{
98 		  reg->state[i] ^= ((MAX_UNSIGNED) 1 << target);
99 		}
100 	    }
101 	}
102       quantum_decohere(reg);
103     }
104 }
105 
106 /* Apply an unbounded toffoli gate. This gate is not considered
107 elementary and is not available on all physical realizations of a
108 quantum computer. Be sure to pass the function the correct number of
109 controlling qubits. The target is given in the last argument. */
110 
111 void
quantum_unbounded_toffoli(int controlling,quantum_reg * reg,...)112 quantum_unbounded_toffoli(int controlling, quantum_reg *reg, ...)
113 {
114   va_list bits;
115   int target;
116   int *controls;
117   int i, j;
118 
119   controls = malloc(controlling * sizeof(int));
120 
121   if(!controls)
122     quantum_error(QUANTUM_ENOMEM);
123 
124   quantum_memman(controlling * sizeof(int));
125 
126   va_start(bits, reg);
127 
128   for(i=0; i<controlling; i++)
129     controls[i] = va_arg(bits, int);
130 
131   target = va_arg(bits, int);
132 
133   va_end(bits);
134 
135 #ifdef _OPENMP
136 #pragma omp parallel for private (j)
137 #endif
138   for(i=0; i<reg->size; i++)
139     {
140       for(j=0; (j < controlling) &&
141 	    (reg->state[i] & (MAX_UNSIGNED) 1 << controls[j]); j++);
142 
143       if(j == controlling) /* all control bits are set */
144 	reg->state[i] ^= ((MAX_UNSIGNED) 1 << target);
145     }
146 
147   free(controls);
148   quantum_memman(-controlling * sizeof(int));
149 
150   quantum_decohere(reg);
151 
152 }
153 
154 
155 /* Apply a sigma_x (or not) gate */
156 
157 void
quantum_sigma_x(int target,quantum_reg * reg)158 quantum_sigma_x(int target, quantum_reg *reg)
159 {
160   int i;
161   int qec;
162 
163   quantum_qec_get_status(&qec, NULL);
164 
165   if(qec)
166     quantum_sigma_x_ft(target, reg);
167   else
168     {
169       if(quantum_objcode_put(SIGMA_X, target))
170 	return;
171 
172 #ifdef _OPENMP
173 #pragma omp parallel for
174 #endif
175       for(i=0; i<reg->size; i++)
176 	{
177 	  /* Flip the target bit of each basis state */
178 
179 	  reg->state[i] ^= ((MAX_UNSIGNED) 1 << target);
180 	}
181       quantum_decohere(reg);
182     }
183 }
184 
185 /* Apply a sigma_y gate */
186 
187 void
quantum_sigma_y(int target,quantum_reg * reg)188 quantum_sigma_y(int target, quantum_reg *reg)
189 {
190   int i;
191 
192   if(quantum_objcode_put(SIGMA_Y, target))
193     return;
194 
195 #ifdef _OPENMP
196 #pragma omp parallel for
197 #endif
198   for(i=0; i<reg->size;i++)
199     {
200       /* Flip the target bit of each basis state and multiply with
201 	 +/- i */
202 
203       reg->state[i] ^= ((MAX_UNSIGNED) 1 << target);
204 
205       if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
206 	reg->amplitude[i] *= IMAGINARY;
207       else
208 	reg->amplitude[i] *= -IMAGINARY;
209     }
210 
211   quantum_decohere(reg);
212 }
213 
214 /* Apply a sigma_y gate */
215 
216 void
quantum_sigma_z(int target,quantum_reg * reg)217 quantum_sigma_z(int target, quantum_reg *reg)
218 {
219   int i;
220 
221   if(quantum_objcode_put(SIGMA_Z, target))
222     return;
223 
224 #ifdef _OPENMP
225 #pragma omp parallel for
226 #endif
227   for(i=0; i<reg->size; i++)
228     {
229       /* Multiply with -1 if the target bit is set */
230 
231       if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
232 	reg->amplitude[i] *= -1;
233     }
234   quantum_decohere(reg);
235 }
236 
237 /* Swap the first WIDTH bits of the quantum register. This is done
238    classically by renaming the bits, unless QEC is enabled. */
239 
240 void
quantum_swaptheleads(int width,quantum_reg * reg)241 quantum_swaptheleads(int width, quantum_reg *reg)
242 {
243   int i, j;
244   int pat1, pat2;
245   int qec;
246   MAX_UNSIGNED l;
247 
248   quantum_qec_get_status(&qec, NULL);
249 
250   if(qec)
251     {
252       for(i=0; i<width; i++)
253 	{
254 	  quantum_cnot(i, width+i, reg);
255 	  quantum_cnot(width+i, i, reg);
256 	  quantum_cnot(i, width+i, reg);
257 	}
258     }
259   else
260     {
261       for(i=0; i<reg->size; i++)
262 	{
263 
264 	  if(quantum_objcode_put(SWAPLEADS, width))
265 	    return;
266 
267 	  /* calculate left bit pattern */
268 
269 	  pat1 = reg->state[i] % ((MAX_UNSIGNED) 1 << width);
270 
271 	  /*calculate right but pattern */
272 
273 	  pat2 = 0;
274 
275 	  for(j=0; j<width; j++)
276 	    pat2 += reg->state[i] & ((MAX_UNSIGNED) 1 << (width + j));
277 
278 	  /* construct the new basis state */
279 
280 	  l = reg->state[i] - (pat1 + pat2);
281 	  l += (pat1 << width);
282 	  l += (pat2 >> width);
283 	  reg->state[i] = l;
284 	}
285     }
286 }
287 
288 /* Swap WIDTH bits starting at WIDTH and 2*WIDTH+2 controlled by
289    CONTROL */
290 
291 void
quantum_swaptheleads_omuln_controlled(int control,int width,quantum_reg * reg)292 quantum_swaptheleads_omuln_controlled(int control, int width, quantum_reg *reg)
293 {
294    int i;
295 
296   for(i=0; i<width; i++)
297     {
298       quantum_toffoli(control, width+i, 2*width+i+2, reg);
299       quantum_toffoli(control, 2*width+i+2, width+i, reg);
300       quantum_toffoli(control, width+i, 2*width+i+2, reg);
301     }
302 }
303 
304 /* Apply the 2x2 matrix M to the target bit. M should be unitary. */
305 
306 void
quantum_gate1(int target,quantum_matrix m,quantum_reg * reg)307 quantum_gate1(int target, quantum_matrix m, quantum_reg *reg)
308 {
309   int i, j, k, iset;
310   int addsize=0, decsize=0;
311   COMPLEX_FLOAT t, tnot=0;
312   float limit;
313   char *done;
314 
315   if((m.cols != 2) || (m.rows != 2))
316     quantum_error(QUANTUM_EMSIZE);
317 
318   if(reg->hashw)
319     {
320       quantum_reconstruct_hash(reg);
321 
322       /* calculate the number of basis states to be added */
323 
324       for(i=0; i<reg->size; i++)
325 	{
326 	  /* determine whether XORed basis state already exists */
327 
328 	  if(quantum_get_state(reg->state[i]
329 			       ^ ((MAX_UNSIGNED) 1 << target), *reg) == -1)
330 	    addsize++;
331 	}
332 
333       /* allocate memory for the new basis states */
334 
335       reg->state = realloc(reg->state,
336 			       (reg->size + addsize) * sizeof(MAX_UNSIGNED));
337       reg->amplitude = realloc(reg->amplitude,
338 			       (reg->size + addsize) * sizeof(COMPLEX_FLOAT));
339 
340       if(reg->size && !(reg->state && reg->amplitude))
341 	quantum_error(QUANTUM_ENOMEM);
342 
343       quantum_memman(addsize*(sizeof(COMPLEX_FLOAT) + sizeof(MAX_UNSIGNED)));
344 
345       for(i=0; i<addsize; i++)
346 	{
347 	  reg->state[i+reg->size] = 0;
348 	  reg->amplitude[i+reg->size] = 0;
349 	}
350 
351     }
352 
353   done = calloc(reg->size + addsize, sizeof(char));
354 
355   if(!done)
356     quantum_error(QUANTUM_ENOMEM);
357 
358   quantum_memman(reg->size + addsize * sizeof(char));
359 
360   k = reg->size;
361 
362   limit = (1.0 / ((MAX_UNSIGNED) 1 << reg->width)) * epsilon;
363 
364   /* perform the actual matrix multiplication */
365 
366   for(i=0; i<reg->size; i++)
367     {
368       if(!done[i])
369 	{
370 	  /* determine if the target of the basis state is set */
371 
372 	  iset = reg->state[i] & ((MAX_UNSIGNED) 1 << target);
373 
374 	  tnot = 0;
375 	  j = quantum_get_state(reg->state[i]
376 				^ ((MAX_UNSIGNED) 1<<target), *reg);
377 	  if(j >= 0)
378 	    tnot = reg->amplitude[j];
379 
380 	  t = reg->amplitude[i];
381 
382 	  if(iset)
383 	    reg->amplitude[i] = m.t[2] * tnot + m.t[3] * t;
384 
385 	  else
386 	    reg->amplitude[i] = m.t[0] * t + m.t[1] * tnot;
387 
388 	  if(j >= 0)
389 	    {
390 	      if(iset)
391 		reg->amplitude[j] = m.t[0] * tnot + m.t[1] * t;
392 
393 	      else
394 		reg->amplitude[j] = m.t[2] * t + m.t[3] * tnot;
395 	    }
396 
397 
398 	  else /* new basis state will be created */
399 	    {
400 
401 	      if((m.t[1] == 0) && (iset))
402 		break;
403 	      if((m.t[2] == 0) && !(iset))
404 		 break;
405 
406 	      reg->state[k] = reg->state[i]
407 		^ ((MAX_UNSIGNED) 1 << target);
408 
409 	      if(iset)
410 		reg->amplitude[k] = m.t[1] * t;
411 
412 	      else
413 		reg->amplitude[k] = m.t[2] * t;
414 
415 	      k++;
416 	    }
417 
418 	  if(j >= 0)
419 	    done[j] = 1;
420 
421 	}
422     }
423 
424   reg->size += addsize;
425 
426   free(done);
427   quantum_memman(-reg->size * sizeof(char));
428 
429   /* remove basis states with extremely small amplitude */
430 
431   if(reg->hashw)
432     {
433       for(i=0, j=0; i<reg->size; i++)
434 	{
435 	  if(quantum_prob_inline(reg->amplitude[i]) < limit)
436 	    {
437 	      j++;
438 	      decsize++;
439 	    }
440 
441 	  else if(j)
442 	    {
443 	      reg->state[i-j] = reg->state[i];
444 	      reg->amplitude[i-j] = reg->amplitude[i];
445 	    }
446 	}
447 
448       if(decsize)
449 	{
450 	  reg->size -= decsize;
451 	  reg->amplitude = realloc(reg->amplitude,
452 				   reg->size * sizeof(COMPLEX_FLOAT));
453 	  reg->state = realloc(reg->state,
454 			       reg->size * sizeof(MAX_UNSIGNED));
455 
456 
457 	  if(reg->size && !(reg->state && reg->amplitude))
458 	    quantum_error(QUANTUM_ENOMEM);
459 
460 	  quantum_memman(-decsize * (sizeof(MAX_UNSIGNED)
461 				     + sizeof(COMPLEX_FLOAT)));
462 	}
463     }
464 
465   if(reg->size > (1 << (reg->hashw-1)))
466     fprintf(stderr, "Warning: inefficient hash table (size %i vs hash %i)\n",
467 	    reg->size, 1<<reg->hashw);
468 
469 
470   quantum_decohere(reg);
471 }
472 
473 /* Apply the 4x4 matrix M to the bits TARGET1 and TARGET2. M should be
474    unitary.
475 
476    Warning: code is mostly untested.*/
477 
478 void
quantum_gate2(int target1,int target2,quantum_matrix m,quantum_reg * reg)479 quantum_gate2(int target1, int target2, quantum_matrix m, quantum_reg *reg)
480 {
481   int i, j, k, l;
482   int addsize=0, decsize=0;
483   COMPLEX_FLOAT psi_sub[4];
484   int base[4];
485   int bits[2];
486   float limit;
487   char *done;
488 
489   if((m.cols != 4) || (m.rows != 4))
490     quantum_error(QUANTUM_EMSIZE);
491 
492   /* Build hash table */
493 
494   for(i=0; i<(1 << reg->hashw); i++)
495     reg->hash[i] = 0;
496 
497   for(i=0; i<reg->size; i++)
498     quantum_add_hash(reg->state[i], i, reg);
499 
500   /* calculate the number of basis states to be added */
501 
502   for(i=0; i<reg->size; i++)
503     {
504       if(quantum_get_state(reg->state[i] ^ ((MAX_UNSIGNED) 1 << target1),
505 			   *reg) == -1)
506 	addsize++;
507       if(quantum_get_state(reg->state[i] ^ ((MAX_UNSIGNED) 1 << target2),
508 			   *reg) == -1)
509 	addsize++;
510     }
511 
512   /* allocate memory for the new basis states */
513 
514   reg->state = realloc(reg->state,
515 		       (reg->size + addsize) * sizeof(MAX_UNSIGNED));
516   reg->amplitude = realloc(reg->amplitude,
517 			   (reg->size + addsize) * sizeof(COMPLEX_FLOAT));
518 
519   if(reg->size && !(reg->state && reg->amplitude))
520     quantum_error(QUANTUM_ENOMEM);
521 
522   quantum_memman(addsize*(sizeof(COMPLEX_FLOAT) + sizeof(MAX_UNSIGNED)));
523 
524   for(i=0; i<addsize; i++)
525     {
526       reg->state[i+reg->size] = 0;
527       reg->amplitude[i+reg->size] = 0;
528     }
529 
530   done = calloc(reg->size + addsize, sizeof(char));
531 
532   if(!done)
533     quantum_error(QUANTUM_EMSIZE);
534 
535   quantum_memman(reg->size + addsize * sizeof(char));
536 
537   l = reg->size;
538 
539   limit = (1.0 / ((MAX_UNSIGNED) 1 << reg->width)) / 1000000;
540 
541   bits[0] = target1;
542   bits[1] = target2;
543 
544   /* perform the actual matrix multiplication */
545 
546   for(i=0; i<reg->size; i++)
547     {
548       if(!done[i])
549 	{
550 	  j = quantum_bitmask(reg->state[i], 2, bits);
551 	  base[j] = i;
552 	  base[j ^ 1] = quantum_get_state(reg->state[i]
553 					  ^ ((MAX_UNSIGNED) 1 << target2),
554 					  *reg);
555 	  base[j ^ 2] = quantum_get_state(reg->state[i]
556 					  ^ ((MAX_UNSIGNED) 1 << target1),
557 					  *reg);
558 	  base[j ^ 3] = quantum_get_state(reg->state[i]
559 					  ^ ((MAX_UNSIGNED) 1 << target1)
560 					  ^ ((MAX_UNSIGNED) 1 << target2),
561 					  *reg);
562 
563 	  for(j=0; j<4; j++)
564 	    {
565 	      if(base[j] == -1)
566 		{
567 		  base[j] = l;
568 		  //		  reg->node[l].state = reg->state[i]
569 		  l++;
570 		}
571 	      psi_sub[j] = reg->amplitude[base[j]];
572 	    }
573 
574 	  for(j=0; j<4; j++)
575 	    {
576 	      reg->amplitude[base[j]] = 0;
577 	      for(k=0; k<4; k++)
578 		reg->amplitude[base[j]] += M(m, k, j) * psi_sub[k];
579 
580 	      done[base[j]] = 1;
581 	    }
582 
583 	}
584     }
585 
586   reg->size += addsize;
587 
588   free(done);
589 
590   quantum_memman(-reg->size * sizeof(char));
591 
592   /* remove basis states with extremely small amplitude */
593 
594   for(i=0, j=0; i<reg->size; i++)
595     {
596       if(quantum_prob_inline(reg->amplitude[i]) < limit)
597 	{
598 	  j++;
599 	  decsize++;
600 	}
601 
602       else if(j)
603 	{
604 	  reg->state[i-j] = reg->state[i];
605 	  reg->amplitude[i-j] = reg->amplitude[i];
606 	}
607     }
608 
609   if(decsize)
610     {
611       reg->size -= decsize;
612       reg->amplitude = realloc(reg->amplitude,
613 			       reg->size * sizeof(COMPLEX_FLOAT));
614       reg->state = realloc(reg->state,
615 			   reg->size * sizeof(MAX_UNSIGNED));
616 
617 
618       if(reg->size && !(reg->state && reg->amplitude))
619 	quantum_error(QUANTUM_ENOMEM);
620 
621       quantum_memman(-decsize * (sizeof(MAX_UNSIGNED)
622 				 + sizeof(COMPLEX_FLOAT)));
623 
624     }
625 
626   quantum_decohere(reg);
627 }
628 
629 /* Apply a hadamard gate */
630 
631 void
quantum_hadamard(int target,quantum_reg * reg)632 quantum_hadamard(int target, quantum_reg *reg)
633 {
634   quantum_matrix m;
635 
636   if(quantum_objcode_put(HADAMARD, target))
637     return;
638 
639   m = quantum_new_matrix(2, 2);
640 
641   m.t[0] = sqrt(1.0/2);  m.t[1] = sqrt(1.0/2);
642   m.t[2] = sqrt(1.0/2);  m.t[3] = -sqrt(1.0/2);
643 
644   quantum_gate1(target, m, reg);
645 
646   quantum_delete_matrix(&m);
647 
648 }
649 
650 /* Apply a walsh-hadamard transform */
651 
652 void
quantum_walsh(int width,quantum_reg * reg)653 quantum_walsh(int width, quantum_reg *reg)
654 {
655   int i;
656 
657   for(i=0; i<width; i++)
658     quantum_hadamard(i, reg);
659 
660 }
661 
662 /* Apply a rotation about the x-axis by the angle GAMMA */
663 
664 void
quantum_r_x(int target,float gamma,quantum_reg * reg)665 quantum_r_x(int target, float gamma, quantum_reg *reg)
666 {
667   quantum_matrix m;
668 
669   if(quantum_objcode_put(ROT_X, target, (double) gamma))
670     return;
671 
672   m = quantum_new_matrix(2, 2);
673 
674   m.t[0] = cos(gamma / 2);              m.t[1] = -IMAGINARY * sin(gamma / 2);
675   m.t[2] = -IMAGINARY * sin(gamma / 2); m.t[3] = cos(gamma / 2);
676 
677   quantum_gate1(target, m, reg);
678 
679   quantum_delete_matrix(&m);
680 
681 }
682 
683 /* Apply a rotation about the y-axis by the angle GAMMA */
684 
685 void
quantum_r_y(int target,float gamma,quantum_reg * reg)686 quantum_r_y(int target, float gamma, quantum_reg *reg)
687 {
688   quantum_matrix m;
689 
690   if(quantum_objcode_put(ROT_Y, target, (double) gamma))
691     return;
692 
693   m = quantum_new_matrix(2, 2);
694 
695   m.t[0] = cos(gamma / 2);  m.t[1] = -sin(gamma / 2);
696   m.t[2] = sin(gamma / 2);  m.t[3] = cos(gamma / 2);
697 
698   quantum_gate1(target, m, reg);
699 
700   quantum_delete_matrix(&m);
701 
702 }
703 
704 /* Apply a rotation about the z-axis by the angle GAMMA */
705 
706 void
quantum_r_z(int target,float gamma,quantum_reg * reg)707 quantum_r_z(int target, float gamma, quantum_reg *reg)
708 {
709   int i;
710   COMPLEX_FLOAT z;
711 
712   if(quantum_objcode_put(ROT_Z, target, (double) gamma))
713     return;
714 
715   z = quantum_cexp(gamma/2);
716 
717   for(i=0; i<reg->size; i++)
718     {
719       if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
720 	reg->amplitude[i] *= z;
721       else
722 	reg->amplitude[i] /= z;
723     }
724 
725   quantum_decohere(reg);
726 }
727 
728 /* Scale the phase of qubit */
729 
730 void
quantum_phase_scale(int target,float gamma,quantum_reg * reg)731 quantum_phase_scale(int target, float gamma, quantum_reg *reg)
732 {
733   int i;
734   COMPLEX_FLOAT z;
735 
736   if(quantum_objcode_put(PHASE_SCALE, target, (double) gamma))
737     return;
738 
739   z = quantum_cexp(gamma);
740 
741 #ifdef _OPENMP
742 #pragma omp parallel for
743 #endif
744   for(i=0; i<reg->size; i++)
745     {
746       reg->amplitude[i] *= z;
747     }
748 
749   quantum_decohere(reg);
750 }
751 
752 
753 /* Apply a phase kick by the angle GAMMA */
754 
755 void
quantum_phase_kick(int target,float gamma,quantum_reg * reg)756 quantum_phase_kick(int target, float gamma, quantum_reg *reg)
757 {
758   int i;
759   COMPLEX_FLOAT z;
760 
761   if(quantum_objcode_put(PHASE_KICK, target, (double) gamma))
762     return;
763 
764   z = quantum_cexp(gamma);
765 
766 #ifdef _OPENMP
767 #pragma omp parallel for
768 #endif
769   for(i=0; i<reg->size; i++)
770     {
771       if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
772 	reg->amplitude[i] *= z;
773     }
774 
775   quantum_decohere(reg);
776 }
777 
778 /* Apply a conditional phase shift by PI / 2^(CONTROL - TARGET) */
779 
780 void
quantum_cond_phase(int control,int target,quantum_reg * reg)781 quantum_cond_phase(int control, int target, quantum_reg *reg)
782 {
783   int i;
784   COMPLEX_FLOAT z;
785 
786   if(quantum_objcode_put(COND_PHASE, control, target))
787     return;
788 
789   z = quantum_cexp(pi / ((MAX_UNSIGNED) 1 << (control - target)));
790 
791 #ifdef _OPENMP
792 #pragma omp parallel for
793 #endif
794   for(i=0; i<reg->size; i++)
795     {
796       if(reg->state[i] & ((MAX_UNSIGNED) 1 << control))
797 	{
798 	  if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
799 	    reg->amplitude[i] *= z;
800 	}
801     }
802 
803   quantum_decohere(reg);
804 }
805 
806 
807 void
quantum_cond_phase_inv(int control,int target,quantum_reg * reg)808 quantum_cond_phase_inv(int control, int target, quantum_reg *reg)
809 {
810   int i;
811   COMPLEX_FLOAT z;
812 
813   z = quantum_cexp(-pi / ((MAX_UNSIGNED) 1 << (control - target)));
814 
815 #ifdef _OPENMP
816 #pragma omp parallel for
817 #endif
818   for(i=0; i<reg->size; i++)
819     {
820       if(reg->state[i] & ((MAX_UNSIGNED) 1 << control))
821 	{
822 	  if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
823 	    reg->amplitude[i] *= z;
824 	}
825     }
826 
827   quantum_decohere(reg);
828 }
829 
830 
831 void
quantum_cond_phase_kick(int control,int target,float gamma,quantum_reg * reg)832 quantum_cond_phase_kick(int control, int target, float gamma, quantum_reg *reg)
833 {
834   int i;
835   COMPLEX_FLOAT z;
836 
837   if(quantum_objcode_put(COND_PHASE, control, target, (double) gamma))
838     return;
839 
840   z = quantum_cexp(gamma);
841 
842 #ifdef _OPENMP
843 #pragma omp parallel for
844 #endif
845   for(i=0; i<reg->size; i++)
846     {
847       if(reg->state[i] & ((MAX_UNSIGNED) 1 << control))
848 	{
849 	  if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
850 	    reg->amplitude[i] *= z;
851 	}
852      }
853   quantum_decohere(reg);
854 }
855 
856 void
quantum_cond_phase_shift(int control,int target,float gamma,quantum_reg * reg)857 quantum_cond_phase_shift(int control, int target, float gamma, quantum_reg *reg)
858 {
859   int i;
860   COMPLEX_FLOAT z;
861 
862   if(quantum_objcode_put(COND_PHASE, control, target, (double) gamma))
863     return;
864 
865   z = quantum_cexp(gamma/2);
866 
867 #ifdef _OPENMP
868 #pragma omp parallel for
869 #endif
870   for(i=0; i<reg->size; i++)
871     {
872       if(reg->state[i] & ((MAX_UNSIGNED) 1 << control))
873 	{
874 	  if(reg->state[i] & ((MAX_UNSIGNED) 1 << target))
875 	    reg->amplitude[i] *= z;
876 	  else
877 	    reg->amplitude[i] /= z;
878 	}
879      }
880   quantum_decohere(reg);
881 }
882 
883 
884 
885 
886 /* Increase the gate counter by INC steps or reset it if INC < 0. The
887    current value of the counter is returned. */
888 
889 int
quantum_gate_counter(int inc)890 quantum_gate_counter(int inc)
891 {
892   static int counter = 0;
893 
894   if(inc > 0)
895     counter += inc;
896   else if(inc < 0)
897     counter = 0;
898 
899   return counter;
900 }
901