1#   Licensed under the Apache License, Version 2.0 (the "License");
2#   you may not use this file except in compliance with the License.
3#   You may obtain a copy of the License at
4#
5#       http://www.apache.org/licenses/LICENSE-2.0
6#
7#   Unless required by applicable law or agreed to in writing, software
8#   distributed under the License is distributed on an "AS IS" BASIS,
9#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
10#   See the License for the specific language governing permissions and
11#   limitations under the License.
12"""Bravyi-Kitaev transform on fermionic operators."""
13
14from openfermion.ops.operators import (FermionOperator, MajoranaOperator,
15                                       QubitOperator)
16from openfermion.ops.representations import InteractionOperator
17from openfermion.utils.operator_utils import count_qubits
18
19
20def bravyi_kitaev(operator, n_qubits=None):
21    """Apply the Bravyi-Kitaev transform.
22
23    Implementation from arXiv:quant-ph/0003137 and
24    "A New Data Structure for Cumulative Frequency Tables" by Peter M. Fenwick.
25
26    Note that this implementation is equivalent to the one described in
27    arXiv:1208.5986, and is different from the one described in
28    arXiv:1701.07072. The one described in arXiv:1701.07072 is implemented
29    in OpenFermion as `bravyi_kitaev_tree`.
30
31    Args:
32        operator (openfermion.ops.FermionOperator):
33            A FermionOperator to transform.
34        n_qubits (int|None):
35            Can force the number of qubits in the resulting operator above the
36            number that appear in the input operator.
37
38    Returns:
39        transformed_operator: An instance of the QubitOperator class.
40
41    Raises:
42        ValueError: Invalid number of qubits specified.
43    """
44    if isinstance(operator, FermionOperator):
45        return _bravyi_kitaev_fermion_operator(operator, n_qubits)
46    if isinstance(operator, MajoranaOperator):
47        return _bravyi_kitaev_majorana_operator(operator, n_qubits)
48    if isinstance(operator, InteractionOperator):
49        return _bravyi_kitaev_interaction_operator(operator, n_qubits)
50    raise TypeError("Couldn't apply the Bravyi-Kitaev Transform to object "
51                    "of type {}.".format(type(operator)))
52
53
54def _update_set(index, n_qubits):
55    """The bits that need to be updated upon flipping the occupancy
56    of a mode."""
57    indices = set()
58
59    # For bit manipulation we need to count from 1 rather than 0
60    index += 1
61    # Ensure index is not a member of the set
62    index += index & -index
63
64    while index <= n_qubits:
65        indices.add(index - 1)
66        # Add least significant one to index
67        # E.g. 00010100 -> 00011000
68        index += index & -index
69    return indices
70
71
72def _occupation_set(index):
73    """The bits whose parity stores the occupation of mode `index`."""
74    indices = set()
75
76    # For bit manipulation we need to count from 1 rather than 0
77    index += 1
78
79    indices.add(index - 1)
80    parent = index & (index - 1)
81    index -= 1
82    while index != parent:
83        indices.add(index - 1)
84        # Remove least significant one from index
85        # E.g. 00010100 -> 00010000
86        index &= index - 1
87    return indices
88
89
90def _parity_set(index):
91    """The bits whose parity stores the parity of the bits 0 .. `index`."""
92    indices = set()
93
94    while index > 0:
95        indices.add(index - 1)
96        # Remove least significant one from index
97        # E.g. 00010100 -> 00010000
98        index &= index - 1
99    return indices
100
101
102def _bravyi_kitaev_majorana_operator(operator, n_qubits):
103    # Compute the number of qubits.
104    N = count_qubits(operator)
105    if n_qubits is None:
106        n_qubits = N
107    if n_qubits < N:
108        raise ValueError('Invalid number of qubits specified.')
109
110    # Compute transformed operator.
111    transformed_terms = (_transform_majorana_term(term=term,
112                                                  coefficient=coeff,
113                                                  n_qubits=n_qubits)
114                         for term, coeff in operator.terms.items())
115    return inline_sum(summands=transformed_terms, seed=QubitOperator())
116
117
118def _transform_majorana_term(term, coefficient, n_qubits):
119    # Build the Bravyi-Kitaev transformed operators.
120    transformed_ops = (_transform_majorana_operator(majorana_index, n_qubits)
121                       for majorana_index in term)
122    return inline_product(factors=transformed_ops,
123                          seed=QubitOperator((), coefficient))
124
125
126def _transform_majorana_operator(majorana_index, n_qubits):
127    q, b = divmod(majorana_index, 2)
128
129    update_set = _update_set(q, n_qubits)
130    update_set.add(q)
131    occupation_set = _occupation_set(q)
132    parity_set = _parity_set(q)
133
134    if b:
135        return QubitOperator([(q, 'Y')] + [(i, 'X') for i in update_set - {q}] +
136                             [(i, 'Z')
137                              for i in (parity_set ^ occupation_set) - {q}])
138    else:
139        return QubitOperator([(i, 'X') for i in update_set] +
140                             [(i, 'Z') for i in parity_set])
141
142
143def _transform_operator_term(term, coefficient, n_qubits):
144    """
145    Args:
146        term (list[tuple[int, int]]):
147            A list of (mode, raising-vs-lowering) ladder operator terms.
148        coefficient (float):
149        n_qubits (int):
150    Returns:
151        QubitOperator:
152    """
153
154    # Build the Bravyi-Kitaev transformed operators.
155    transformed_ladder_ops = (_transform_ladder_operator(
156        ladder_operator, n_qubits) for ladder_operator in term)
157    return inline_product(factors=transformed_ladder_ops,
158                          seed=QubitOperator((), coefficient))
159
160
161def _bravyi_kitaev_fermion_operator(operator, n_qubits):
162    # Compute the number of qubits.
163    N = count_qubits(operator)
164    if n_qubits is None:
165        n_qubits = N
166    if n_qubits < N:
167        raise ValueError('Invalid number of qubits specified.')
168
169    # Compute transformed operator.
170    transformed_terms = (_transform_operator_term(
171        term=term, coefficient=operator.terms[term], n_qubits=n_qubits)
172                         for term in operator.terms)
173    return inline_sum(summands=transformed_terms, seed=QubitOperator())
174
175
176def _transform_ladder_operator(ladder_operator, n_qubits):
177    """
178    Args:
179        ladder_operator (tuple[int, int]): the ladder operator
180        n_qubits (int): the number of qubits
181    Returns:
182        QubitOperator
183    """
184    index, action = ladder_operator
185
186    update_set = _update_set(index, n_qubits)
187    update_set.add(index)
188    occupation_set = _occupation_set(index)
189    parity_set = _parity_set(index)
190
191    # Initialize the transformed majorana operator (a_p^\dagger + a_p) / 2
192    transformed_operator = QubitOperator([(i, 'X') for i in update_set] +
193                                         [(i, 'Z') for i in parity_set], .5)
194    # Get the transformed (a_p^\dagger - a_p) / 2
195    # Below is equivalent to X(update_set) * Z(parity_set ^ occupation_set)
196    transformed_majorana_difference = QubitOperator(
197        [(index, 'Y')] + [(i, 'X') for i in update_set - {index}] +
198        [(i, 'Z') for i in (parity_set ^ occupation_set) - {index}], -.5j)
199
200    # Raising
201    if action == 1:
202        transformed_operator += transformed_majorana_difference
203    # Lowering
204    else:
205        transformed_operator -= transformed_majorana_difference
206
207    return transformed_operator
208
209
210def inline_sum(summands, seed):
211    """Computes a sum, using the __iadd__ operator.
212    Args:
213        seed (T): The starting total. The zero value.
214        summands (iterable[T]): Values to add (with +=) into the total.
215    Returns:
216        T: The result of adding all the factors into the zero value.
217    """
218    for r in summands:
219        seed += r
220    return seed
221
222
223def inline_product(factors, seed):
224    """Computes a product, using the __imul__ operator.
225    Args:
226        seed (T): The starting total. The unit value.
227        factors (iterable[T]): Values to multiply (with *=) into the total.
228    Returns:
229        T: The result of multiplying all the factors into the unit value.
230    """
231    for r in factors:
232        seed *= r
233    return seed
234
235
236def _bravyi_kitaev_interaction_operator(interaction_operator, n_qubits):
237    """Implementation of the Bravyi-Kitaev transformation for OpenFermion
238    Interaction Operators. This implementation is equivalent to that described
239    in arXiv:1208.5986, and has been written to optimize compute time by using
240    algebraic expressions for general products a_i^dagger a_j^dagger as outlined
241    in Table II of Seeley, Richard, Love.
242    """
243
244    one_body = interaction_operator.one_body_tensor
245    two_body = interaction_operator.two_body_tensor
246    constant_term = interaction_operator.constant
247
248    # Compute the number of qubits.
249    N = len(one_body)
250    if n_qubits is None:
251        n_qubits = N
252    if n_qubits < N:
253        raise ValueError('Invalid number of qubits specified.')
254
255    qubit_hamiltonian = QubitOperator()
256    qubit_hamiltonian_op = []
257    qubit_hamiltonian_coef = []
258
259    #  For cases A - F see Table II of Seeley, Richard, Love
260    for i in range(n_qubits):
261        # A. Number operators: n_i
262        if abs(one_body[i, i]) > 0:
263            qubit_hamiltonian += _qubit_operator_creation(
264                *_seeley_richard_love(i, i, one_body[i, i], n_qubits))
265
266        for j in range(i):
267            # Case B: Coulomb and exchange operators
268            if abs(one_body[i, j]) > 0:
269                operators, coef_list = _seeley_richard_love(
270                    i, j, one_body[i, j], n_qubits)
271                qubit_hamiltonian_op.extend(operators)
272                qubit_hamiltonian_coef.extend(coef_list)
273
274                operators, coef_list = _seeley_richard_love(
275                    j, i, one_body[i, j].conj(), n_qubits)
276                qubit_hamiltonian_op.extend(operators)
277                qubit_hamiltonian_coef.extend(coef_list)
278
279            coef = _two_body_coef(two_body, i, j, j, i) / 4
280            if abs(coef) > 0:
281                qubit_hamiltonian_op.append(
282                    tuple((index, "Z") for index in _occupation_set(i)))
283                qubit_hamiltonian_op.append(
284                    tuple((index, "Z") for index in _occupation_set(j)))
285                qubit_hamiltonian_op.append(
286                    tuple((index, "Z") for index in _F_ij_set(i, j)))
287
288                qubit_hamiltonian_coef.append(-coef)
289                qubit_hamiltonian_coef.append(-coef)
290                qubit_hamiltonian_coef.append(coef)
291                constant_term += coef
292
293    # C. Number-excitation operators: n_i a_j^d a_k
294    for i in range(n_qubits):
295        for j in range(n_qubits):
296            for k in range(j):
297                if i not in (j, k):
298                    coef = _two_body_coef(two_body, i, j, k, i)
299
300                    if abs(coef) > 0:
301                        number = _qubit_operator_creation(
302                            *_seeley_richard_love(i, i, 1, n_qubits))
303
304                        excitation_op, excitation_coef = _seeley_richard_love(
305                            j, k, coef, n_qubits)
306                        operators_hc, coef_list_hc = _seeley_richard_love(
307                            k, j, coef.conj(), n_qubits)
308                        excitation_op.extend(operators_hc)
309                        excitation_coef.extend(coef_list_hc)
310                        excitation = _qubit_operator_creation(
311                            excitation_op, excitation_coef)
312
313                        number *= excitation
314                        qubit_hamiltonian += number
315
316    # D. Double-excitation operators: c_i^d c_j^d c_k c_l
317    for i in range(n_qubits):
318        for j in range(i):
319            for k in range(j):
320                for l in range(k):
321                    coef = -_two_body_coef(two_body, i, j, k, l)
322                    if abs(coef) > 0:
323                        qubit_hamiltonian += _hermitian_one_body_product(
324                            i, j, k, l, coef, n_qubits)
325
326                    coef = -_two_body_coef(two_body, i, k, j, l)
327                    if abs(coef) > 0:
328                        qubit_hamiltonian += _hermitian_one_body_product(
329                            i, k, j, l, coef, n_qubits)
330
331                    coef = -_two_body_coef(two_body, i, l, j, k)
332                    if abs(coef) > 0:
333                        qubit_hamiltonian += _hermitian_one_body_product(
334                            i, l, j, k, coef, n_qubits)
335
336    qubit_hamiltonian_op.append(())
337    qubit_hamiltonian_coef.append(constant_term)
338    qubit_hamiltonian += _qubit_operator_creation(qubit_hamiltonian_op,
339                                                  qubit_hamiltonian_coef)
340
341    return qubit_hamiltonian
342
343
344def _two_body_coef(two_body, a, b, c, d):
345    return two_body[a, b, c, d] - two_body[a, b, d, c] + two_body[
346        b, a, d, c] - two_body[b, a, c, d]
347
348
349def _hermitian_one_body_product(a, b, c, d, coef, n_qubits):
350    """ Takes the 4 indices for a two-body operator and constructs the
351    Bravyi-Kitaev form by splitting the two-body into 2 one-body operators,
352    multiplying them together and then re-adding the Hermitian conjugate to
353    give a Hermitian operator. """
354
355    c_dag_c_ac = _qubit_operator_creation(
356        *_seeley_richard_love(a, c, coef, n_qubits))
357    c_dag_c_bd = _qubit_operator_creation(
358        *_seeley_richard_love(b, d, 1, n_qubits))
359    c_dag_c_ac *= c_dag_c_bd
360    hermitian_sum = c_dag_c_ac
361
362    c_dag_c_ca = _qubit_operator_creation(
363        *_seeley_richard_love(c, a, coef.conj(), n_qubits))
364    c_dag_c_db = _qubit_operator_creation(
365        *_seeley_richard_love(d, b, 1, n_qubits))
366    c_dag_c_ca *= c_dag_c_db
367
368    hermitian_sum += c_dag_c_ca
369    return hermitian_sum
370
371
372def _qubit_operator_creation(operators, coefficents):
373    """ Takes a list of tuples for operators/indices, and another for
374    coefficents"""
375
376    qubit_operator = QubitOperator()
377
378    for index in zip(operators, coefficents):
379        qubit_operator += QubitOperator(index[0], index[1])
380
381    return qubit_operator
382
383
384def _seeley_richard_love(i, j, coef, n_qubits):
385    """Algebraic expressions for general products of the form a_i^d a_j term in
386    the Bravyi-Kitaev basis. These expressions vary in form depending on the
387    parity of the indices i and j, as well as onthe overlaps between the parity
388    and update sets of the indices"""
389
390    seeley_richard_love_op = []
391    seeley_richard_love_coef = []
392    coef *= 0.25
393    # Case 0
394    if i == j:  # Simplifies to the number operator
395        seeley_richard_love_op.append(
396            tuple((index, "Z") for index in _occupation_set(i)))
397        seeley_richard_love_coef.append(-coef * 2)
398
399        seeley_richard_love_op.append(())
400        seeley_richard_love_coef.append(coef * 2)
401
402    # Case 1
403    elif i % 2 == 0 and j % 2 == 0:
404        x_pad = tuple((index, "X") for index in _U_diff_a_set(i, j, n_qubits))
405        y_pad = tuple((index, "Y") for index in _alpha_set(i, j, n_qubits))
406        z_pad = tuple(
407            (index, "Z") for index in _P0_ij_diff_a_set(i, j, n_qubits))
408
409        left_pad = x_pad + y_pad + z_pad
410
411        seeley_richard_love_op.append(left_pad + ((j, "Y"), (i, "X")))
412        seeley_richard_love_op.append(left_pad + ((j, "X"), (i, "Y")))
413        seeley_richard_love_op.append(left_pad + ((j, "X"), (i, "X")))
414        seeley_richard_love_op.append(left_pad + ((j, "Y"), (i, "Y")))
415
416        if i < j:
417            seeley_richard_love_coef.append(coef)
418            seeley_richard_love_coef.append(-coef)
419            seeley_richard_love_coef.append(complex(0, -coef))
420            seeley_richard_love_coef.append(complex(0, -coef))
421
422        else:  # whenever i < j introduce phase of -j
423            seeley_richard_love_coef.append(complex(0, -coef))
424            seeley_richard_love_coef.append(complex(0, coef))
425            seeley_richard_love_coef.append(-coef)
426            seeley_richard_love_coef.append(-coef)
427
428    # Case 2
429    elif i % 2 == 1 and j % 2 == 0 and i not in _parity_set(j):
430        x_pad = tuple((index, "X") for index in _U_diff_a_set(i, j, n_qubits))
431        y_pad = tuple((index, "Y") for index in _alpha_set(i, j, n_qubits))
432
433        left_pad = x_pad + y_pad
434
435        right_pad_1 = tuple(
436            (index, "Z")
437            for index in _P0_ij_set(i, j) - _alpha_set(i, j, n_qubits))
438        right_pad_2 = tuple(
439            (index, "Z")
440            for index in _P2_ij_set(i, j) - _alpha_set(i, j, n_qubits))
441
442        seeley_richard_love_op.append(left_pad + ((j, "Y"),
443                                                  (i, "X")) + right_pad_1)
444        seeley_richard_love_op.append(left_pad + ((j, "X"),
445                                                  (i, "X")) + right_pad_1)
446        seeley_richard_love_op.append(left_pad + ((j, "X"),
447                                                  (i, "Y")) + right_pad_2)
448        seeley_richard_love_op.append(left_pad + ((j, "Y"),
449                                                  (i, "Y")) + right_pad_2)
450
451        if i < j:
452            seeley_richard_love_coef.append(coef)
453            seeley_richard_love_coef.append(complex(0, -coef))
454            seeley_richard_love_coef.append(-coef)
455            seeley_richard_love_coef.append(complex(0, -coef))
456
457        else:
458            seeley_richard_love_coef.append(complex(0, -coef))
459            seeley_richard_love_coef.append(-coef)
460            seeley_richard_love_coef.append(complex(0, coef))
461            seeley_richard_love_coef.append(-coef)
462
463    # Case 3
464    elif i % 2 == 1 and j % 2 == 0 and i in _parity_set(j):
465        left_pad = tuple((index, "X") for index in _U_ij_set(i, j, n_qubits))
466        right_pad_1 = tuple((index, "Z") for index in _P0_ij_set(i, j) - {i})
467        right_pad_2 = tuple((index, "Z") for index in _P2_ij_set(i, j) - {i})
468
469        seeley_richard_love_op.append(left_pad + ((j, "Y"),
470                                                  (i, "Y")) + right_pad_1)
471        seeley_richard_love_coef.append(coef)
472
473        seeley_richard_love_op.append(left_pad + ((j, "X"),
474                                                  (i, "Y")) + right_pad_1)
475        seeley_richard_love_coef.append(complex(0, -coef))
476
477        seeley_richard_love_op.append(left_pad + ((j, "X"),
478                                                  (i, "X")) + right_pad_2)
479        seeley_richard_love_coef.append(coef)
480
481        seeley_richard_love_op.append(left_pad + ((j, "Y"),
482                                                  (i, "X")) + right_pad_2)
483        seeley_richard_love_coef.append(complex(0, coef))
484
485    # Case 4
486    elif i % 2 == 0 and j % 2 == 1 and i not in _parity_set(
487            j) and j not in _update_set(i, n_qubits):
488        x_pad = tuple((index, "X") for index in _U_diff_a_set(i, j, n_qubits))
489        y_pad = tuple((index, "Y") for index in _alpha_set(i, j, n_qubits))
490        left_pad = x_pad + y_pad
491
492        right_pad_1 = tuple(
493            (index, "Z")
494            for index in _P0_ij_set(i, j) - _alpha_set(i, j, n_qubits))
495        right_pad_2 = tuple(
496            (index, "Z")
497            for index in _P1_ij_set(i, j) - _alpha_set(i, j, n_qubits))
498
499        seeley_richard_love_op.append(left_pad + ((j, "X"),
500                                                  (i, "Y")) + right_pad_1)
501        seeley_richard_love_op.append(left_pad + ((j, "X"),
502                                                  (i, "X")) + right_pad_1)
503        seeley_richard_love_op.append(left_pad + ((j, "Y"),
504                                                  (i, "X")) + right_pad_2)
505        seeley_richard_love_op.append(left_pad + ((j, "Y"),
506                                                  (i, "Y")) + right_pad_2)
507
508        if i < j:
509            seeley_richard_love_coef.append(-coef)
510            seeley_richard_love_coef.append(complex(0, -coef))
511            seeley_richard_love_coef.append(coef)
512            seeley_richard_love_coef.append(complex(0, -coef))
513        else:
514            seeley_richard_love_coef.append(complex(0, coef))
515            seeley_richard_love_coef.append(-coef)
516            seeley_richard_love_coef.append(complex(0, -coef))
517            seeley_richard_love_coef.append(-coef)
518
519    # Case 5
520    elif i % 2 == 0 and j % 2 == 1 and i not in _parity_set(
521            j) and j in _update_set(i, n_qubits):
522        x_range_1 = _U_ij_set(i, j, n_qubits) - {j}
523        left_pad_1 = tuple((index, "X") for index in x_range_1)
524
525        x_range_2 = x_range_1 - _alpha_set(i, j, n_qubits)
526        left_pad_2 = tuple((index, "X") for index in x_range_2)
527
528        y_pad = tuple((index, "Y") for index in _alpha_set(i, j, n_qubits))
529        z_range_1 = _P0_ij_set(i, j) - _alpha_set(i, j, n_qubits)
530        z_pad = tuple((index, "Z") for index in z_range_1)
531        right_pad_1 = y_pad + z_pad
532
533        z_range_2 = _P1_ij_set(i, j).union({j})
534        right_pad_2 = tuple((index, "Z") for index in z_range_2)
535
536        seeley_richard_love_op.append(left_pad_2 + ((i, "Y"),) + right_pad_1)
537        seeley_richard_love_coef.append(-coef)
538
539        seeley_richard_love_op.append(left_pad_2 + ((i, "X"),) + right_pad_1)
540        seeley_richard_love_coef.append(complex(
541            0, -coef))  # Phase flip of -1 relative to original paper
542
543        seeley_richard_love_op.append(left_pad_1 + ((i, "Y"),) + right_pad_2)
544        seeley_richard_love_coef.append(complex(0, coef))
545
546        seeley_richard_love_op.append(left_pad_1 + ((i, "X"),) + right_pad_2)
547        seeley_richard_love_coef.append(-coef)
548
549    # Case 6
550    elif i % 2 == 0 and j % 2 == 1 and i in _parity_set(j) and j in _update_set(
551            i, n_qubits):
552        left_pad = tuple(
553            (index, "X") for index in _U_ij_set(i, j, n_qubits) - {j})
554        right_pad = tuple((index, "Z") for index in _P1_ij_set(i, j).union({j}))
555
556        seeley_richard_love_op.append(left_pad + ((i, "X"),))
557        seeley_richard_love_coef.append(coef)
558
559        seeley_richard_love_op.append(left_pad + ((i, "Y"),))
560        seeley_richard_love_coef.append(complex(0, -coef))
561
562        seeley_richard_love_op.append(left_pad + ((i, "Y"),) + right_pad)
563        seeley_richard_love_coef.append(complex(0, coef))
564
565        seeley_richard_love_op.append(left_pad + ((i, "X"),) + right_pad)
566        seeley_richard_love_coef.append(-coef)
567
568    # Case 7
569    elif i % 2 == 1 and j % 2 == 1 and i not in _parity_set(
570            j) and j not in _update_set(i, n_qubits):
571        x_pad = tuple((index, "X") for index in _U_diff_a_set(i, j, n_qubits))
572        y_pad = tuple((index, "Y") for index in _alpha_set(i, j, n_qubits))
573        left_pad = x_pad + y_pad
574
575        z_range_1 = _P0_ij_set(i, j) - _alpha_set(i, j, n_qubits)
576        right_pad_1 = tuple((index, "Z") for index in z_range_1)
577
578        z_range_2 = _P1_ij_set(i, j) - _alpha_set(i, j, n_qubits)
579        right_pad_2 = tuple((index, "Z") for index in z_range_2)
580
581        z_range_3 = _P2_ij_set(i, j) - _alpha_set(i, j, n_qubits)
582        right_pad_3 = tuple((index, "Z") for index in z_range_3)
583
584        z_range_4 = _P3_ij_set(i, j) - _alpha_set(i, j, n_qubits)
585        right_pad_4 = tuple((index, "Z") for index in z_range_4)
586
587        if i < j:
588            seeley_richard_love_op.append(left_pad + ((j, "X"),
589                                                      (i, "X")) + right_pad_1)
590            seeley_richard_love_coef.append(complex(0, -coef))
591
592            seeley_richard_love_op.append(left_pad + ((j, "Y"),
593                                                      (i, "X")) + right_pad_2)
594            seeley_richard_love_coef.append(coef)
595
596            seeley_richard_love_op.append(left_pad + ((j, "X"),
597                                                      (i, "Y")) + right_pad_3)
598            seeley_richard_love_coef.append(-coef)
599
600            seeley_richard_love_op.append(left_pad + ((j, "Y"),
601                                                      (i, "Y")) + right_pad_4)
602            seeley_richard_love_coef.append(complex(0, -coef))
603
604        else:
605            seeley_richard_love_op.append(left_pad + ((j, "X"),
606                                                      (i, "X")) + right_pad_1)
607            seeley_richard_love_coef.append(-coef)
608
609            seeley_richard_love_op.append(left_pad + ((j, "Y"),
610                                                      (i, "X")) + right_pad_2)
611            seeley_richard_love_coef.append(complex(0, -coef))
612
613            seeley_richard_love_op.append(left_pad + ((j, "X"),
614                                                      (i, "Y")) + right_pad_3)
615            seeley_richard_love_coef.append(complex(0, coef))
616
617            seeley_richard_love_op.append(left_pad + ((j, "Y"),
618                                                      (i, "Y")) + right_pad_4)
619            seeley_richard_love_coef.append(-coef)
620
621    # Case 8
622    elif i % 2 == 1 and j % 2 == 1 and i in _parity_set(
623            j) and j not in _update_set(i, n_qubits):
624        left_pad = tuple((index, "X") for index in _U_ij_set(i, j, n_qubits))
625
626        z_range_1 = _P0_ij_set(i, j) - {i}
627        right_pad_1 = tuple((index, "Z") for index in z_range_1)
628
629        z_range_2 = _P1_ij_set(i, j) - {i}
630        right_pad_2 = tuple((index, "Z") for index in z_range_2)
631
632        z_range_3 = _P2_ij_set(i, j) - {i}
633        right_pad_3 = tuple((index, "Z") for index in z_range_3)
634
635        z_range_4 = _P3_ij_set(i, j) - {i}
636        right_pad_4 = tuple((index, "Z") for index in z_range_4)
637
638        seeley_richard_love_op.append(left_pad + ((j, "X"),
639                                                  (i, "Y")) + right_pad_1)
640        seeley_richard_love_coef.append(complex(0, -coef))
641
642        seeley_richard_love_op.append(left_pad + ((j, "Y"),
643                                                  (i, "Y")) + right_pad_2)
644        seeley_richard_love_coef.append(coef)
645
646        seeley_richard_love_op.append(left_pad + ((j, "X"),
647                                                  (i, "X")) + right_pad_3)
648        seeley_richard_love_coef.append(coef)
649
650        seeley_richard_love_op.append(left_pad + ((j, "Y"),
651                                                  (i, "X")) + right_pad_4)
652        seeley_richard_love_coef.append(complex(0, coef))
653
654    # Case 9
655    elif i % 2 == 1 and j % 2 == 1 and i not in _parity_set(
656            j) and j in _update_set(i, n_qubits):
657        x_range_1 = _U_ij_set(i, j, n_qubits) - {j}
658        left_pad_3 = tuple((index, "X") for index in x_range_1)
659
660        x_range_2 = x_range_1 - _alpha_set(i, j, n_qubits)
661        left_pad_1 = tuple((index, "X") for index in x_range_2)
662
663        x_range_3 = x_range_1.union({i}) - _alpha_set(i, j, n_qubits)
664        left_pad_2 = tuple((index, "X") for index in x_range_3)
665
666        z_range_1 = _P2_ij_set(i, j) - _alpha_set(i, j, n_qubits)
667        z_pad_1 = tuple((index, "Z") for index in z_range_1)
668        y_range = _alpha_set(i, j, n_qubits)
669        y_pad = tuple((index, "Y") for index in y_range)
670        right_pad_1 = z_pad_1 + y_pad
671
672        z_range_2 = _P0_ij_set(i, j) - _alpha_set(i, j, n_qubits)
673        z_pad_2 = tuple((index, "Z") for index in z_range_2)
674        right_pad_2 = z_pad_2 + y_pad
675
676        z_range_3 = _P1_ij_set(i, j).union({j})
677        right_pad_3 = tuple((index, "Z") for index in z_range_3)
678
679        z_range_4 = _P3_ij_set(i, j).union({j})
680        right_pad_4 = tuple((index, "Z") for index in z_range_4)
681
682        seeley_richard_love_op.append(left_pad_1 + ((i, "Y"),) + right_pad_1)
683        seeley_richard_love_coef.append(
684            -coef)  # phase of -j relative to original paper
685
686        seeley_richard_love_op.append(left_pad_2 + right_pad_2)
687        seeley_richard_love_coef.append(complex(
688            0, -coef))  # phase of -j relative to original paper
689
690        seeley_richard_love_op.append(left_pad_3 + ((i, "X"),) + right_pad_3)
691        seeley_richard_love_coef.append(-coef)
692
693        seeley_richard_love_op.append(left_pad_3 + ((i, "Y"),) + right_pad_4)
694        seeley_richard_love_coef.append(complex(0, coef))
695
696    # Case 10
697    elif i % 2 == 1 and j % 2 == 1 and i in _parity_set(j) and j in _update_set(
698            i, n_qubits):
699        left_pad = tuple(
700            (index, "X") for index in _U_ij_set(i, j, n_qubits) - {j})
701        right_pad_1 = tuple((index, "Z") for index in _P0_ij_set(i, j) - {i})
702        right_pad_2 = tuple((index, "Z") for index in _P2_ij_set(i, j) - {i})
703        right_pad_3 = tuple((index, "Z") for index in _P1_ij_set(i, j))
704        right_pad_4 = tuple((index, "Z") for index in _P3_ij_set(i, j))
705
706        seeley_richard_love_op.append(left_pad + ((i, "Y"),) + right_pad_1)
707        seeley_richard_love_coef.append(complex(0, -coef))
708
709        seeley_richard_love_op.append(left_pad + ((i, "X"),) + right_pad_2)
710        seeley_richard_love_coef.append(coef)
711
712        seeley_richard_love_op.append(left_pad + ((j, "Z"),
713                                                  (i, "X")) + right_pad_3)
714        seeley_richard_love_coef.append(-coef)
715
716        seeley_richard_love_op.append(left_pad + ((j, "Z"),
717                                                  (i, "Y")) + right_pad_4)
718        seeley_richard_love_coef.append(complex(0, coef))
719
720    return seeley_richard_love_op, seeley_richard_love_coef
721
722
723def _remainder_set(index):
724    return _parity_set(index) - _occupation_set(index)
725
726
727def _F_ij_set(i, j):
728    return _occupation_set(i).symmetric_difference(_occupation_set(j))
729
730
731def _P0_ij_set(i, j):
732    """ The symmetric difference of sets P(i) and P(j). """
733    return _parity_set(i).symmetric_difference(_parity_set(j))
734
735
736def _P1_ij_set(i, j):
737    """ The symmetric difference of sets P(i) and R(j). """
738    return _parity_set(i).symmetric_difference(_remainder_set(j))
739
740
741def _P2_ij_set(i, j):
742    """ The symmetric difference of sets R(i) and P(j). """
743    return _remainder_set(i).symmetric_difference(_parity_set(j))
744
745
746def _P3_ij_set(i, j):
747    """ The symmetric difference of sets R(i) and R(j). """
748    return _remainder_set(i).symmetric_difference(_remainder_set(j))
749
750
751def _U_ij_set(i, j, n_qubits):
752    """ The symmetric difference of sets U(i) and U(j)"""
753    return _update_set(i,
754                       n_qubits).symmetric_difference(_update_set(j, n_qubits))
755
756
757def _U_diff_a_set(i, j, n_qubits):
758    """ Calculates the set {member U_ij diff alpha_ij}. """
759    return _U_ij_set(i, j, n_qubits) - _alpha_set(i, j, n_qubits)
760
761
762def _P0_ij_diff_a_set(i, j, n_qubits):
763    """ Calculates the set {member P_ij diff alpha_ij}. """
764    return _P0_ij_set(i, j) - _alpha_set(i, j, n_qubits)
765
766
767def _alpha_set(i, j, n_qubits):
768    return _update_set(i, n_qubits).intersection(_parity_set(j))
769