Source code for qrisp.arithmetic.SBP_arithmetic

"""
\********************************************************************************
* Copyright (c) 2023 the Qrisp authors
*
* This program and the accompanying materials are made available under the
* terms of the Eclipse Public License 2.0 which is available at
* http://www.eclipse.org/legal/epl-2.0.
*
* This Source Code may also be made available under the following Secondary
* Licenses when the conditions for such availability set forth in the Eclipse
* Public License, v. 2.0 are satisfied: GNU General Public License, version 2
* with the GNU Classpath Exception which is
* available at https://www.gnu.org/software/classpath/license.html.
*
* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0
********************************************************************************/
"""


import numpy as np
import sympy as sp

from qrisp.arithmetic.poly_tools import (
    expr_to_list,
    filter_pow,
    get_ordered_symbol_list,
)
from qrisp.core import QuantumArray, QuantumVariable, cp, cx, cz, h, mcx, p, z, rz, crz
from qrisp.misc import gate_wrap
from qrisp.circuit import XGate

# Threshold of rounding used in detecting integer multiples of pi
pi_mult_round_threshold = 11

depth_tracker = {}


# Check if the given expression is a polynomial
def check_for_polynomial(expr):
    for x in sp.preorder_traversal(expr):
        if not isinstance(
            x,
            (
                sp.core.mul.Mul,
                sp.core.add.Add,
                sp.core.symbol.Symbol,
                sp.core.numbers.Integer,
                sp.core.numbers.Float,
                sp.core.power.Pow,
            ),
        ):
            return False

        if isinstance(x, sp.core.power.Pow):
            if not isinstance(x.args[1], sp.core.numbers.Integer):
                return False
    return True


# Efficient implementation of the multicontrolled U_g gate
def multi_controlled_U_g(
    output_qf, control_qb_list, y, phase_tolerant=False, use_gms=False
):
    # Set alias for quantum session
    # qs = output_qf.qs
    qs = control_qb_list[0].qs()

    # For one control qubit, there is no advantage in calculating an ancilla qubit
    if len(control_qb_list) == 1:
        ancilla = control_qb_list
    # Otherwise request an ancilla qubit
    else:
        ancilla = QuantumVariable(1, qs=output_qf.qs, name="sbp_anc*")

    # If an ancilla qb is available
    # we synthesize the result of the boolean multiplication in it

    if len(control_qb_list) != 1:
        if use_gms:
            toffoli_method = "gms"
        else:
            toffoli_method = "gray_pt"
        
        # Apply multi-controlled x gate
        mcx(control_qb_list, ancilla[0], method=toffoli_method)

    # Now we execute the phase gates
    # For this there is 2 possibilies
    # Either use regular phase gates or use GMS gates
    # GMS gates are ion-trap native gates, that allow the entanglement
    # of multiple qubits within a single laser pulse
    
    from qrisp.environments import QuantumEnvironment, control, GMSEnvironment, custom_control
    if use_gms:
        # GMSEnvironment is an environment that allows programming
        # with regular phase gates but converts everything programmed
        # into a GMS gate upon leaving the environment

        env = GMSEnvironment(qs)

    else:
        # This is an environment that simply does nothing
        env = QuantumEnvironment()

    # Enter environment
    env.manual_allocation_management = True
    
    
    @custom_control
    def crz_helper(rot_angle, a, b, ctrl = None, use_gms = False):
        if ctrl is None and not use_gms:
            cx(a, b)
            p(-rot_angle / 2, b)
            cx(a, b)
            p(rot_angle / 2, b)
        elif use_gms:
            crz(rot_angle, ancilla[0], output_qf[i])
        else:
            cx(a, b)
            cp(-rot_angle / 2, ctrl, b)
            cx(a, b)
            cp(rot_angle / 2, ctrl, b)
            
    
    with env:
        phase_accumulator = 0
        cz_counter = 0

        # Execute controlled rotations
        for i in range(len(output_qf)):
            # The -i (instead of i) reverses the gate order implying that
            # we can leave out the swaps of the QFT

            rot_angle = 2 * np.pi * 2 ** (-i - 1) * y
            rot_angle = rot_angle % (2 * np.pi)

            if (
                np.round(rot_angle, pi_mult_round_threshold) != 0
                and np.round((-rot_angle) % (2 * np.pi), pi_mult_round_threshold) != 0
            ):
                # If phase is equal to pi, execute cz gate (costs one CNOT less than CP)
                if False:
                    # if np.round(abs(rot_angle) - np.pi, pi_mult_round_threshold) == 0:
                    # cz(ancilla[0], output_qf.reg[i])
                    
                    z(output_qf.reg[i])
                    cz_counter += 1

                    phase_accumulator += rot_angle / 2
                # Otherwise execute cp gate
                else:
                    crz_helper(rot_angle, ancilla[0], output_qf[i], use_gms = use_gms)
                    phase_accumulator += rot_angle / 2
        
        p(phase_accumulator - np.pi * cz_counter / 2, ancilla[0])

    # Uncompute boolean multiplication
    if len(control_qb_list) != 1:
        if not use_gms:
            toffoli_method += "_inv"

        mcx(control_qb_list, ancilla[0], method=toffoli_method)
        ancilla.delete()

    return phase_accumulator


# This function returns the U_g gate (up to qiskit endian conventions)
# This is used for applying U_g gates which have no control-knobs
def U_g(y, qv):
    for i in range(len(qv)):
        rot_angle = np.pi * 2 ** (-i) * y
        rot_angle = rot_angle % (2 * np.pi)
        if (
            np.round(rot_angle, pi_mult_round_threshold) != 0
            and np.round((-rot_angle) % (2 * np.pi), pi_mult_round_threshold) != 0
        ):
            p(rot_angle, qv[i])


# This function encodes a semi-boolean polynomial into a circuit i.e.
# For the polynomial p(x_0, x_1, x_2) = 4*x_0*x_1 + 2*x_0*x_2
# the effect of this function is:
# U_f |x_0, x_1, x_2>|0> = |x_0, x_1, x_2>|p(x_1,x_2,x_3)>
def sb_polynomial_encoder(
    input_qf_list, output_qf, poly, inplace_mult=1, use_gms=False, init_op="auto"
):
    # As the polynomial has only boolean variables,
    # powers can be ignored since x**k = x for x in GF(2)
    poly = filter_pow(poly.expand()).expand() / 2.0**output_qf.exponent

    # Acquire list of symbols present in the polynomial
    symbol_list = []

    for qf in input_qf_list:
        temp_var_list = list(sp.symbols(qf.name + "_" + "0:" + str(qf.size)))
        symbol_list += temp_var_list

    n = len(symbol_list)

    if n != sum([var.size for var in input_qf_list]):
        raise Exception(
            "Input variables do not the required amount of qubits to encode polynomial"
        )

    # Acquire monomials in list form
    monomial_list = expr_to_list(poly)

    from qrisp import QFT, conjugate, QuantumEnvironment
    

    if inplace_mult == 1:
        env = conjugate(QFT)(
            output_qf,
            inv=False,
            exec_swap=False,
            inplace_mult=inplace_mult,
            use_gms=use_gms,
        )
    else:
        QFT(
            output_qf,
            inv=False,
            exec_swap=False,
            inplace_mult=inplace_mult,
            use_gms=use_gms,
        )
        env = QuantumEnvironment()

    with env:
        # The list of qubits contained in the variables of input_var_list
        input_qubits = sum([list(var.reg) for var in input_qf_list], [])
    
        control_qubit_list = []
        y_list = []
    
        # Iterate through the monomials
        for monom in monomial_list:
            # Prepare the two variables coeff (which is the coefficient of the monomial)
            # And the list of variables which appear in the monomial
    
            # For this, go through the cases which can appear
    
            # This describes the case where there is only a single term in the monomial
            # Either a constant or a variable
            if len(monom) == 1:
                if isinstance(monom[0], sp.core.symbol.Symbol):
                    coeff = 1
                    variables = list(monom)
                else:
                    coeff = float(monom[0])
                    variables = []
    
            # This describes the case where there is multiple terms in the monomial
            elif not isinstance(monom[0], sp.core.symbol.Symbol):
                coeff = monom[0]
                variables = list(monom[1:])
            else:
                coeff = 1
                variables = list(monom)
    
            # Check if the coefficient is an integer (up to float errors)
            if abs(int(np.round(float(coeff))) - coeff) > 1e-14:
                pass
                # raise Exception("Tried to encode sb-polynomial
                # with non-integer coefficient")
    
            # Append coefficient to y_list
            y_list.append(int(np.round(float(coeff))))
    
            # Prepare the qubits on which the U_g should be controlled
            control_qubit_numbers = [symbol_list.index(var) for var in variables]
    
            control_qubits = [input_qubits[nr] for nr in control_qubit_numbers]
    
            control_qubits = list(set(control_qubits))
    
            control_qubits.sort(key=lambda x: x.identifier)
    
            control_qubit_list.append(control_qubits)
    
        # Now we apply the multi controlled U_g gate
        # Here the order in which they are applied makes a huge difference
        # In order to determine, which U_g to apply next, we evaluate a cost function
        # (which we determined through trial and error)
        # and choose the U_g with the lowest cost
    
        # Note that for this feature to yield an improvement, the quantum session requires
        # multiple free ancilla qubits to work on
        def delay_cost(depth_array):
            return max(depth_array)
    
        def find_best_monomial(control_qubit_list, depth_dic):
            delay_cost_list = []
    
            for i in range(len(control_qubit_list)):
                depth_array = np.array([depth_dic[qb] for qb in control_qubit_list[i]])
    
                delay_cost_list.append(delay_cost(depth_array))
    
            return np.argmin(delay_cost_list)
    
        # Iterate through the list of U_g gates
        while control_qubit_list:
            # TO-DO fix depth calculation inside environment
            # Update depth_dic (contains the depth of each qubit)
            # depth_dic = output_qf.qs.get_depth_dic()
    
            # Determine best U_g
            # monomial_index = find_best_monomial(control_qubit_list, depth_dic)
            monomial_index = 0
            # Find control qubits and their coefficient
            control_qubits = control_qubit_list.pop(monomial_index)
            y = y_list.pop(monomial_index)
    
            # Apply (controlled) U_g
            if len(control_qubits):
                multi_controlled_U_g(output_qf, control_qubits, y, use_gms=use_gms)
            else:
                U_g(y, output_qf)

    if inplace_mult != 1:
        # Apply QFT
        QFT(output_qf, inv=True, exec_swap=False, use_gms=use_gms)


# Multiplies two integers
# The general idea is to represent the integers as values of two polynomials
# p(x) = x_0 + 2*x_1 + 4*x_2 ...
# p(y) = y_0 + 2*y_1 + 4*y_2 ...
# these polynomials get multiplied and which forms a bigger polynomial
# which can be encoded using the polynomial encoder
[docs]def sbp_mult(factor_1_qf, factor_2_qf, output_qf=None): """ Performs multiplication based on the evaluation of `semi-boolean polynomials <https://ieeexplore.ieee.org/document/9815035>`_. Parameters ---------- factor_1_qf : QuantumFloat The first factor to multiply. factor_2_qf : QuantumFloat The second factor to multiply. output_qf : QuantumFloat, optional The QuantumFloat to store the result in. By default, a suited new QuantumFloat is created. Returns ------- output_qf : QuantumFloat A QuantumFloat containing the result of the multiplication. Examples -------- We multiply two QuantumFloats: :: from qrisp import QuantumFloat, sbp_mult qf_0 = QuantumFloat(3) qf_1 = QuantumFloat(3) qf_0[:] = 3 qf_1[:] = 4 qf_res = sbp_mult(qf_0, qf_1) print(qf_res) :: #Yields: {12: 1.0} """ if output_qf is None: from qrisp.qtypes.quantum_float import create_output_qf output_qf = create_output_qf([factor_1_qf, factor_2_qf], op="mul") # Multiply the polynmials mult_poly = factor_1_qf.sb_poly(output_qf.msize) * factor_2_qf.sb_poly( output_qf.msize ) # Apply sb encoder sb_polynomial_encoder([factor_1_qf, factor_2_qf], output_qf, mult_poly) return output_qf
[docs]def sbp_add(summand_1_qf, summand_2_qf, output_qf=None): """ Performs addition based on the evaluation of `semi-boolean polynomials <https://ieeexplore.ieee.org/document/9815035>`_. Parameters ---------- summand_1_qf : QuantumFloat The first summand to add. summand_2_qf : QuantumFloat The second summand to add. output_qf : QuantumFloat, optional The QuantumFloat to store the result in. By default, a suited new QuantumFloat is created. Returns ------- output_qf : QuantumFloat A QuantumFloat containing the result of the addition. Examples -------- We add two QuantumFloats: :: from qrisp import QuantumFloat, sbp_add qf_0 = QuantumFloat(3) qf_1 = QuantumFloat(3) qf_0[:] = 3 qf_1[:] = 4 qf_res = sbp_add(qf_0, qf_1) print(qf_res) :: # Yields: {7: 1.0} """ if output_qf is None: from qrisp.qtypes.quantum_float import create_output_qf output_qf = create_output_qf([summand_1_qf, summand_2_qf], op="add") sum_poly = summand_1_qf.sb_poly(output_qf.msize) + summand_2_qf.sb_poly( output_qf.msize ) # Apply sb encoder sb_polynomial_encoder([summand_1_qf, summand_2_qf], output_qf, sum_poly) return output_qf
[docs]def sbp_sub(summand_1_qf, summand_2_qf, output_qf=None): """ Performs subtraction based on the evaluation of `semi-boolean polynomials <https://ieeexplore.ieee.org/document/9815035>`_. Parameters ---------- summand_1_qf : QuantumFloat The QuantumFloat to subtract from. summand_2_qf : QuantumFloat The QuantumFloat to subtract. output_qf : QuantumFloat, optional The QuantumFloat to store the result in. By default, a suited new QuantumFloat is created. Returns ------- output_qf : QuantumFloat A QuantumFloat containing the result of the subtraction. Examples -------- We add two QuantumFloats: :: from qrisp import QuantumFloat, sbp_sub qf_0 = QuantumFloat(3) qf_1 = QuantumFloat(3) qf_0[:] = 3 qf_1[:] = 4 qf_res = sbp_sub(qf_0, qf_1) print(qf_res) :: # Yields: {-1: 1.0} """ if output_qf is None: from qrisp.qtypes.quantum_float import create_output_qf output_qf = create_output_qf([summand_1_qf, summand_2_qf], op="sub") dif_poly = summand_1_qf.sb_poly(output_qf.msize) - summand_2_qf.sb_poly( output_qf.msize ) # Apply sb encoder sb_polynomial_encoder([summand_1_qf, summand_2_qf], output_qf, dif_poly) return output_qf
# Encode the polynomial given in poly in output var # depending on the input variables qv_list @gate_wrap(is_qfree=True, permeability=[0]) def polynomial_encoder(qf_list, output_qf, poly, encoding_dic=None, inplace_mult=1): """ Evaluates a (multivariate) sympy polynomial on a list of QuantumFloats using `semi-boolean polynomials <https://ieeexplore.ieee.org/document/9815035>`_. Parameters ---------- qf_list : list[QuantumFloat] The list of QuantumFloats to evaluate the polynomial on. output_qf : QuantumFloat The QuantumFloat to evaluate into. poly : sympy expression The polynomial to evaluate. encoding_dic : dict, optional A dictionary which has the QuantumFloats of qf_list as keys and the associated sympy symbols as values. By default, the symbols of the polynomial will be ordered alphabetically and then matched to the order in qf_list. inplace_mult : int, optional This integer allow to perform an inplace multiplication on output_qf, before the polynomial is evaluated. Note that due to reversibility only odd numbers are supported. Raises ------ Exception Provided QuantumFloat list does not include the appropriate amount of elements to encode given polynomial". Returns ------- None. Examples -------- We evaluate the polynomial $x^2 + 2y^2$ on two QuantumFloats: :: from sympy import Symbol x = Symbol("x") y = Symbol("y") poly = x**2 + 2*y**2 from qrisp import QuantumFloat, polynomial_encoder x_qf = QuantumFloat(3) y_qf = QuantumFloat(3) x_qf[:] = 3 y_qf[:] = 2 res_qf = QuantumFloat(7) encoding_dic = {x_qf : x, y_qf : y} polynomial_encoder([x_qf, y_qf], res_qf, poly, encoding_dic) print(res_qf) :: # Yields: {17.0: 1.0} """ if isinstance(qf_list, QuantumArray): qf_list = list(qf_list.flatten().qv_array) if encoding_dic is not None: symbol_list = [encoding_dic[qv.name] for qv in qf_list] else: symbol_list = get_ordered_symbol_list(poly) if len(symbol_list) != len(qf_list): raise Exception( "Provided QuantumFloat list does not include the appropriate amount" "of elements to encode given polynomial" ) if not output_qf.signed: for qf in qf_list: if qf.signed: raise Exception( "When encoding into an unsigned quantum float" "provide only unsigned inputs" ) sb_poly_list = [qf.sb_poly(output_qf.size) for qf in qf_list] repl_dic = {symbol_list[i]: sb_poly_list[i] for i in range(len(qf_list))} # Substitute SB-polynomials sb_polynomial = poly.subs(repl_dic).expand() # Apply sb encoder sb_polynomial_encoder(qf_list, output_qf, sb_polynomial, inplace_mult=inplace_mult) def split(a, n): k, m = divmod(len(a), n) return (a[i * k + min(i, m) : (i + 1) * k + min(i + 1, m)] for i in range(n)) # Performs inplace addition on a QFT-ed variable # by applying successive U_g gates def U_g_inpl_adder(modified_var, summand, mult_factor=1): applied_phases = [] # Extract the coefficients of the SB-polynomial of x from sympy import Poly summand_coeffs = list(Poly(summand.sb_poly(modified_var.msize)).as_dict().values()) summand_coeffs = [float(coeff) for coeff in summand_coeffs][::-1] # summand_coeffs.sort() summand_coeffs = summand_coeffs + (summand.size - len(summand_coeffs)) * [0] for j in range(summand.size): if summand_coeffs[j] == 0: continue applied_phases.append( multi_controlled_U_g( modified_var, [summand[j]], mult_factor * summand_coeffs[j] * 2**-modified_var.exponent, ) ) return applied_phases # This algorithm is based on the equation X*U_g(y)*X = U_g(-y)*G with G as global phase # This is used to perform conditional addition/subtraction required for multiplication # using algorithm 3 from https://arxiv.org/abs/2112.10537 # s = x << (n + 1) # s− = x # for i in range(y.size): # if y[i]: # s+= (x << i) # else: # s−= (x << i) # return s >> 1 # We will use the function U_g_inply_adder for the inplace additions # while performing conditional additions using CNOT gates. # The approach is therefore a hybrid of SBP-method and traditional logic approaches
[docs]def hybrid_mult( x, y, output_qf=None, init_op="h", terminal_op="qft", phase_tolerant=False, cl_factor=1, ): """ An advanced algorithm for multiplication which has better depth, gate-count and compile time than :meth:`sbp_mult <qrisp.sbp_mult>`. It does not support squaring a single QuantumFloat though. This algorithm also operates on the Fourier transform. Because of this, between successive multiplications targeting the same QuantumFloat it is not neccessary to Fourier-Transform. This advantage is expressed in the parameters init_op and terminal_op. These can be set to either 'h', 'qft' or None to leave out self canceling Fourier-transforms. Parameters ---------- x : QuantumFloat The first factor to multiply. y : QuantumFloat The second factor to multiply. output_qf : QuantumFloat, optional The QuantumFloat to store the result in. By default a suited QuantumFloat is created. init_op : str, optional The operation to bring output_qf into it's Fourier-transform. The default is 'h'. terminal_op : str, optional The operation to bring output_qf back from it's Fourier-transform. The default is "qft". phase_tolerant : bool, optional If set to True, differing results introduce differing extra phases. This can be usefull to save resources incase this functions will get uncomputed. The default is False. cl_factor : float, optional Allows to multiply the result by a classical factor without any extra gates. The default is 1. Returns ------- output_qf : QuantumFloat The QuantumFloat containing the result. Examples -------- We multiply two QuantumFloat with eachother and an additional classical factor :: from qrisp import QuantumFloat, hybrid_mult qf_0 = QuantumFloat(3) qf_1 = QuantumFloat(3) qf_0[:] = 3 qf_1[:] = 4 qf_res = hybrid_mult(qf_0, qf_1, cl_factor = 2) print(qf_res) :: # Yields: {24: 1.0} """ from qrisp import QFT, cx, h, merge, z # The two factors take asymetrical roles in this algorithm # This implies that there is most likely a prefered choice # for the roles depending on the size of the factors # Several trials showed that these roles work best if not x.size > y.size: x, y = y, x # We shift the exponent of both factors such that they are integers # and shift the result in the end, by the sum of both exponents # This allows convenient treatment of non-integer inputs # while only constructing the algorithm for integers x_exp = int(x.exponent) y_exp = int(y.exponent) x.exp_shift(-x_exp) y.exp_shift(-y_exp) # If no output_qf is given, create one from qrisp.qtypes.quantum_float import create_output_qf if isinstance(output_qf, type(None)): output_qf = create_output_qf([y, x], op="mul") else: output_qf.exp_shift(-(x_exp + y_exp)) output_qf.extend(1, position=0) # output_qf.exp_shift(-1) # Since the result is right shifted at the end, this implies that the zero-th qubit # of the output will not contain any information (otherwise this would imply, # that integer multiplications can yield non-integer results). # However the result still needs to be able to # display every possible result. This is why it needs "extra working-space" # We therefore extend it by one qubit # Merge the sessions of all involved variables merge([output_qf, x, y]) # Perform initial operation # (check the general documentation of SBP arithmetic for more details) if init_op == "h": h(output_qf) elif init_op == "qft": QFT(output_qf, exec_swap=False) else: h(output_qf[0]) # We treat the case that y is signed by applying a negation on output_qf that is # conditioned on the sign qubit of y (see command [2]). This will be reversed # at the end of the algorithm. This implies that every addition that happens # in between is actually a subtraction. # So this acts as a sign reversal of the result. # We then need to make sure we multiply with the absolute of y # We do this by applying a bitwise negation on the mantissa of y conditioned on # the sign qubit (command [4]). This bitwise negation acts as # NOT y = -y + 1 # This implies we need to correct the additional +1. # This additional +1 is equivalent to an extra +x in the result # Therefore we need to subtract this +x from output_qf. # Note that this subtraction also has to be performed # conditioned on the sign qubit of y. # Instead of adding an additional control qubit, we again follow the strategy # of turning a conditional subtraction into a conditional subtraction/addition if y.signed: # This performs the first part of the conditional subtraction/addition # for j in range(x.size): # multi_controlled_U_g(output_qf, [x[j]], -x_coeffs[j]) U_g_inpl_adder(output_qf, x, -1 * cl_factor) # command [1] # Negate output_qf conditioned on the sign qubit of y cx(y[-1], output_qf) # command [2] # This would perform the second part of the subtraction/addition # However this command is merged into command [5] for more performance gains. # If commented in, this command reverses command [1] # if output_qf has not been negated IF the sign qubit # has not negated output_qf. Otherwise the subtraction of x is completed. # for j in range(x.size): # multi_controlled_U_g(output_qf, [x[j]], x_coeffs[j]) #command [3] # Negate mantissa of y cx(y[-1], y[:-1]) # command [4] # This performs the initial two lines of the initially mentioned # multiplication algorithm: # s = x << (n + 1) # s− = x # Note that the boolean y.signed adds the phase that would have been added # in command [3] # however not requiring another round of U_g gates applied_phases = U_g_inpl_adder( output_qf, x, cl_factor * (2 ** (y.msize) - 1 + y.signed) ) # We now come to the loop of the multiplication algorithm # for i in range(y.size): # if y[i]: # s+= (x << i) # else: # s−= (x << i) # Negate output_qf in order to perform the conditional addition/subtraction if not phase_tolerant and terminal_op != "qft": hybrid_mult_anc = QuantumVariable(1) if y.signed: cx(y[-1], hybrid_mult_anc[0]) for k in range(len(applied_phases)): # cp(-applied_phases[k] * 2, hybrid_mult_anc[0], x[k]) cp(-applied_phases[k] * 2, x[k], hybrid_mult_anc[0]) cx(y[0], output_qf) for i in range(y.msize): # Perform in-place addition # (note that if y[i] is active an addition has to be performed). # Therefore, we need to perform a subtraction here, # because if y[i] is active output_qf has been negated applied_phases = U_g_inpl_adder(output_qf, x, -(2**i) * cl_factor) if not phase_tolerant and terminal_op != "qft": if i != 0: cx(y[i - 1], hybrid_mult_anc[0]) cx(y[i], hybrid_mult_anc[0]) for k in range(len(applied_phases)): # cp(-applied_phases[k] * 2, hybrid_mult_anc[0], x[k]) cp(-applied_phases[k] * 2, x[k], hybrid_mult_anc[0]) # sbp+= -2*applied_phases[k]*Symbol("y_" + str(i))*Symbol("x_" + str(k)) # This command is equivalent to # cx(y[i], output_qf) # cx(y[i+1], output_qf) # ie. performs the output_qf negation but requires less cnot gates if i != y.msize - 1: cx(y[i], y[i + 1]) cx(y[i + 1], output_qf) cx(y[i], y[i + 1]) else: if not y.signed: # This command performs the final negation # if y is signed, we can perform the negation # together with the negation conditioned on the sign qubit # of y (command [6]) cx(y[i], output_qf) # command [5] if not phase_tolerant and terminal_op != "qft": cx(y[i], hybrid_mult_anc[0]) if y.signed: cx(y[-1], hybrid_mult_anc[0]) # from sympy import Symbol, simplify # temp_x = sum([2**i*Symbol("x_" + str(i)) for i in range(3)], 0) # temp_y = sum([2**i*Symbol("y_" + str(i)) for i in range(3)], 0) # sbp = sbp/(2*np.pi) # sbp = simplify((sbp + (temp_x*temp_y).expand()/2**7)) # print(sbp) hybrid_mult_anc.delete() if y.signed: # This makes sure the negation from command [5] is performed cx(y[i], y[-1]) cx(y[-1], output_qf) # command[6] cx(y[i], y[-1]) # Reverse command [4] cx(y[-1], y[:-1]) # Free up the qubit which we identified as containing no information h(output_qf[0]) output_qf.reduce(output_qf[0]) # Perform terminal qft if terminal_op == "qft": QFT(output_qf, inv=True, exec_swap=False) if not phase_tolerant and terminal_op == "qft": # This is a phase tolerant correction using only single qubit gates. # We found it by looking at the sbp of the correction using cp gates and # identified the sbp of x*y in it. for i in range(output_qf.size): p(-(2 * np.pi) * 2**i / 2 ** (output_qf.size + 1), output_qf[i]) if output_qf.signed: z(output_qf[-1]) output_qf.exp_shift((x_exp + y_exp)) # Perform exponent shifts x.exp_shift(x_exp) y.exp_shift(y_exp) return output_qf
# Wrapper for choosing the best multiplication algorithm def q_mult(factor_1, factor_2, target=None, method="auto"): if method == "auto": if factor_1.reg == factor_2.reg: return q_mult(factor_1, factor_2, target, method="sbp") else: return q_mult(factor_1, factor_2, target, method="hybrid") elif method == "sbp": from qrisp.qtypes.quantum_float import create_output_qf if target is None: target = create_output_qf([factor_1, factor_2], op="mul") sbp_mult(factor_1, factor_2, target) return target elif method == "hybrid": from qrisp.arithmetic import hybrid_mult return hybrid_mult(factor_1, factor_2) def QFT_inpl_mult(qv, inplace_mult=1): from qrisp.misc import is_inv qv = list(qv) qv = qv[::-1] n = len(qv) if not is_inv(inplace_mult, n): raise Exception( "Tried to perform non-invertible inplace multiplication during Fourier-Transform" ) # Perform QFT with inplace multiplication for i in range(n): if i != n - 1: h(qv[i]) if i == n - 1: break for k in range(n - i - 1): if k + i + 1 != n - 1: cp(inplace_mult * 2 * np.pi / 2 ** (k + 2), qv[k + i + 1], qv[i]) else: # The -1 here cancels some cp gates of the inverse qft cp((inplace_mult - 1) * 2 * np.pi / 2 ** (k + 2), qv[k + i + 1], qv[i]) # Perform reversed QFT without inplace multiplication and without canceled steps for i in range(n)[::-1]: for k in range(n - i - 1)[::-1]: if k + i + 1 != n - 1: cp(-2 * np.pi / 2 ** (k + 2), qv[k + i + 1], qv[i]) if i != n - 1: h(qv[i]) return qv
[docs]def inpl_mult(qf, mult_int, treat_overflow=True): """ Performs inplace multiplication of a :ref:`QuantumFloat` with a classical integer. To prevent overflow errors, this function automatically adjusts the mantissa size. If you want to prevent this behavior, set ``treat_overflow = False``. Parameters ---------- qf : QuantumFloat The QuantumFloat to inplace-multiply. mult_int : int The integer to perform the multiplication with. treat_overflow : bool If set to ``False``, the mantissa will not be extended to prevent overflow errors. The default is ``True``. Examples -------- We create a QuantumFloat, bring it to superposition and perform an inplace multiplication. :: from qrisp import QuantumFloat, h, inpl_mult a = QuantumFloat(5, signed = True) h(a[0]) h(a[-1]) print(a) :: # Yields: {0: 0.25, 1: 0.25, -32: 0.25, -31: 0.25} :: inpl_mult(a, -5) print(a) :: # Yields: {0: 0.25, 155: 0.25, 160: 0.25, -5: 0.25} """ if mult_int < 0 and not qf.signed: raise Exception( "Tried to inplace-multiply unsigned QuantumFloat with negative factor" ) bit_shift = 0 if int(mult_int) != mult_int: c = abs(mult_int) for i in range(32): if int(2**i*c) == 2**i*c: break else: raise Exception("Tried to inplace multiply with number of to much precision") bit_shift = -i mult_int = 2**i*mult_int else: while not mult_int % 2: bit_shift += 1 mult_int = mult_int // 2 if mult_int != 1: if treat_overflow: extension_size = int(np.ceil(np.log2(abs(mult_int)))) qf.extend(extension_size, position=qf.size - 1) if qf.signed: cx(qf[-1], qf[-1 - extension_size : -1]) from qrisp.arithmetic.SBP_arithmetic import QFT_inpl_mult QFT_inpl_mult(qf, inplace_mult=mult_int) quantum_bit_shift(qf, bit_shift, treat_overflow) if treat_overflow and bit_shift<0 and qf.signed: cx(qf[-1], qf[bit_shift-1:-1])
def quantum_bit_shift(qf, bit_shift, treat_overflow = True): from qrisp import cyclic_shift, control, QuantumFloat if isinstance(bit_shift, QuantumFloat): if bit_shift.signed or qf.signed: raise Exception("Quantum-quantum bitshifting is currently only supported for unsigned arguments") for i in range(*bit_shift.mshape): with control(bit_shift.significant(i)): quantum_bit_shift(qf, 2**i) return if treat_overflow: if bit_shift > 0: if qf.signed: qf.extend(bit_shift, position=qf.size-1) else: qf.extend(bit_shift, position=qf.size) else: qf.extend(abs(bit_shift), position=0) qf.exp_shift(bit_shift) if qf.signed: cyclic_shift(qf[:-1], bit_shift) else: cyclic_shift(qf, bit_shift)