The investigation of the various variants of the roots of X and Y.

try_sqrt_matrices_and_corresponding_circuits_of_u_x_r_gates.py

########################################################################################################################
#
# (C) 2020 and later Copyright by Aleksei E. Zhuravlev aka ZA
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
#
########################################################################################################################
#
# The purpose of this program: The investigation of the various variants of the roots of X and Y.
# Inspired by
# https://quantumcomputing.stackexchange.com/questions/6236/how-to-quickly-calculate-the-custom-u3-gate-parameters-theta-phi-and-lamb
#
########################################################################################################################
from qiskit import QuantumCircuit
from qiskit.circuit.library import RXGate, XGate, YGate, U3Gate, U2Gate, RYGate
from qiskit.quantum_info import process_fidelity, Operator
from qiskit.quantum_info.operators.predicates import is_unitary_matrix
from qiskit.quantum_info.synthesis import OneQubitEulerDecomposer
import numpy as np
from math import isclose

sqrtx_matrices = [np.array([[-0.5 + 0.5j, -0.5 - 0.5j], [-0.5 - 0.5j, -0.5 + 0.5j]]),
                  np.array([[-0.5 - 0.5j, -0.5 + 0.5j], [-0.5 + 0.5j, -0.5 - 0.5j]]),
                  np.array([[0.5 + 0.5j, 0.5 - 0.5j], [0.5 - 0.5j, 0.5 + 0.5j]]),  # The most popular!
                  np.array([[0.5 - 0.5j, 0.5 + 0.5j], [0.5 + 0.5j, 0.5 - 0.5j]])]

sqrty_matrices = [np.array([[-0.5 - 0.5j, 0.5 + 0.5j], [-0.5 - 0.5j, -0.5 - 0.5j]]),
                  np.array([[0.5 - 0.5j, 0.5 - 0.5j], [-0.5 + 0.5j, 0.5 - 0.5j]]),
                  np.array([[-0.5 + 0.5j, -0.5 + 0.5j], [0.5 - 0.5j, -0.5 + 0.5j]]),
                  np.array([[0.5 + 0.5j, -0.5 - 0.5j], [0.5 + 0.5j, 0.5 + 0.5j]])]  # The most popular!

# The roots matrices differ slightly (only by signs),
# let's see how similar the Euler angles and phase are obtained from these matrices.

# But first, prepare a set of other input data for checks.

# Other optimized circuits for most popular matrices:

other_circ_sqrtx = QuantumCircuit(1)
other_circ_sqrtx.u2(-np.pi / 4, np.pi / 2, 0)
other_circ_sqrtx.x(0)
other_circ_sqrtx.u3(np.pi, 0, 5 * np.pi / 4, 0)

other_circ_sqrty = QuantumCircuit(1)
other_circ_sqrty.u2(np.pi / 4, 0, 0)
other_circ_sqrty.x(0)
other_circ_sqrty.u3(np.pi, 0, 5 * np.pi / 4, 0)

gates = {'X': (XGate, RXGate, sqrtx_matrices, other_circ_sqrtx),
         'Y': (YGate, RYGate, sqrty_matrices, other_circ_sqrty)}


# Phase shift gate, see:  https://quantumcomputing.stackexchange.com/questions/2477/phase-shift-gate-in-qiskit
# def Ph(qc: QuantumCircuit, gamma: float, qubit):
#     qc.u1(gamma, qubit)
#     qc.x(qubit)
#     qc.u1(gamma, qubit)
#     qc.x(qubit)
#     return qc

# Optimized Phase shift gate, see almost the same with my fix u3(0,...) on u3(pi, ...):
# https://quantumcomputing.stackexchange.com/a/12618/12280
def Ph(qc: QuantumCircuit, gamma: float, qubit):
    qc.x(qubit)
    qc.u3(np.pi, gamma, np.pi + gamma, qubit)
    return qc


def canon(angle: float) -> float:
    if angle < 0 or angle >= np.pi:
        div, mod = divmod(angle, 2 * np.pi)
    else:
        mod = angle
    return round(mod / np.pi, 2)


for gname, (gate, rgate, sqrtg_matrices, other_circ) in gates.items():
    print('=' * 101)
    print(f'Take {gname} gate:')
    for sqrtg_matrix in sqrtg_matrices:
        print('=' * 101)
        print(f'Take the matrix of the square root of {gname}:')
        print(sqrtg_matrix)
        # make sure the matrix is correct:
        assert is_unitary_matrix(sqrtg_matrix)
        assert Operator(sqrtg_matrix).power(2) == Operator(gate())
        # get Euler angles and phase:
        theta, phi, lam, gamma = OneQubitEulerDecomposer('U3').angles_and_phase(sqrtg_matrix)
        print('Euler decomposition results (in rad):')
        print(f'theta={theta}, phi={phi}, lambda={lam}, gamma={gamma}')
        # print(f'theta={theta / np.pi}, phi={round(phi / np.pi, 2)}, lambda={lam / np.pi}, gamma={gamma / np.pi} (in
        # pi)')
        print('or canonical (0 <= angle < 2pi) values (in pi):')
        print(f'theta={canon(theta)}, phi={canon(phi)}, lambda={canon(lam)}, gamma={canon(gamma)}')
        assert 0 <= theta <= np.pi
        # build circuit:
        circuit = QuantumCircuit(1)
        circuit.u3(theta, phi, lam, 0)
        Ph(circuit, gamma, 0)
        print('The circuit with phase shift:')
        print(circuit)
        sgn = 'SQRT' + gname
        print(f'Compliance checks (process_fidelity(U3(angles), {sgn})=1, circuit={sgn}_matrix, circuit**2={gname}):')
        print(check1 := isclose(process_fidelity(Operator(U3Gate(theta, phi, lam)), Operator(sqrtg_matrix)), 1))
        print(check2 := Operator(circuit) == Operator(sqrtg_matrix))
        print(check3 := Operator(circuit).power(2) == Operator(gate()))
        if check1 and check2 and check3:
            print(f'Apparently it is possible to use this circuit as the root of {gname},',
                  'especially if the phase is important!')
        else:
            print(f'This circuit is NOT recommended to use as a square root of {gname}, if the phase is important!')

        if np.allclose(sqrtg_matrix, gate().power(1/2).to_matrix()):
            # sqrtg_matrices[2]):
            print('--------------------------------------------------------------------------------------------')
            print(f'Special checks for the most popular matrix of the square root of {gname} and its angles and phase.')
            print('--------------------------------------------------------------------------------------------')
            print('Other optimized circuit (with phase shift):')
            print(other_circ)
            print(f'Compliance checks (circuit={sgn}_matrix, circuit**2={gname}):')
            print(check4 := Operator(other_circ) == Operator(sqrtg_matrix))
            print(check5 := Operator(other_circ).power(2) == Operator(gate()))
            if check4 and check5:
                print(f'Apparently it is possible to use this circuit as the root of {gname},',
                      'especially if the phase is important!')
            else:
                print(f'This circuit is NOT recommended to use as a square root of {gname}, if the phase is important!')

            print('-' * 101)
            print('As an insufficiently good circuit take the circuit of one U3-gate only ')
            print('(with same angles but without phase shift).')
            print('At first quite obvious, but still gates equivalency check',
                  f'(U3(theta, phi, lam)=U2(phi, lam)=R{gname}(pi/2):')
            print(Operator(U3Gate(theta, phi, lam)) == Operator(U2Gate(phi, lam)) == Operator(rgate(np.pi / 2)))
            print('And so the circuit:')
            circ_without_phase = QuantumCircuit(1)
            # circ_without_phase.u3(np.pi / 2, - np.pi / 2, np.pi / 2, 0)
            circ_without_phase.u3(theta, phi, lam, 0)
            print(circ_without_phase)
            print(f'Compliance checks (circuit={sgn}_matrix, circuit**2={gname}):')
            print(check6 := Operator(circ_without_phase) == Operator(sqrtg_matrix))
            print(check7 := Operator(circ_without_phase).power(2) == Operator(gate()))
            if check6 and check7:
                print(f'Apparently it is possible to use this circuit as the root of {gname},',
                      'especially if the phase is important!')
            else:
                print(f'This circuit is NOT recommended to use as a square root of {gname}, if the phase is important!')