Source code for qrisp.alg_primitives.arithmetic.ripple_division

"""
\********************************************************************************
* Copyright (c) 2025 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
********************************************************************************/
"""

from qrisp.core.gate_application_functions import x
from qrisp import gate_wrap


@gate_wrap(is_qfree=True, permeability="args")
def q_int_div(numerator, divisor, adder="cuccaro", n=None, log_output=True):
    # This function performs integer division based on the following algorithm
    # Which is a reformulation of the division algorithm presented in
    # https://en.wikipedia.org/wiki/Division_algorithm#Non-restoring_division

    # Q = 2**(n)-1
    # D *= 2**(n-1)
    # R = int(N)

    # import numpy as np
    # from qrisp.misc import int_as_array

    # for i in range(n-1, -1, -1):

    #     D = D/2
    #     if not np.sign(D)*np.sign(R) != -1:
    #         Q -= 2**(i)

    #     if int_as_array(Q, n)[::-1][i]:
    #         R = R-D
    #     else:
    #         R = R+D

    # Q -= 2**(n-1)
    # Q += 0.5

    # if D < 0:
    #     Q += -0.5
    #     R = R - D
    # else:
    #     Q += -0.5
    #     R = R + D

    # D *= 2

    # The D multiplications/divisions are realized by bitshifts
    # The conditioned Q -= 2**i operation is performed by NOTing the appropriate
    # qubit controlled on the condition

    # The R = R+-D operation is realized by wrapping applying NOT gates on the R
    # qubits before and after the addition controlled on the corresponding Q-qubit
    # This yields the desired behavior because (x' + y)' = x - y
    # where the prime denotes negation

    # The Q -= 2**(n-1) is simply performed by flipping the significance n-1 bit
    # The Q = Q + 0.5 is performed by adding an aditional qubit at the -1 significance
    # of Q and turning it on

    # The final conditional is performed similar to the previous condition
    # Where the variables which are being operated on are wrapped in CNOTs

    qs = numerator.qs

    if numerator.exponent < 0 or divisor.exponent < 0:
        raise Exception(
            "Tried to call integer division with non integer quantum floats"
        )

    # n is the bitsize of the quotient
    # The largest possible number the quotient can take is
    # 2**(numerator_max_sig - divisor_min_sig)
    # We need one additional bit (why?)

    if n is None:
        n = int(numerator.mshape[1] - divisor.exponent) + 1

    # The quotient bits will successively be calculated/
    # added starting at the most significant bit.
    # As the algorithm produces empty bits in the remainder at the same rate as
    # quotient bits are required, we initialize the quotient as a 1 bit variable
    # Since the highest significance ist calculated fist, this bit needs significance n.

    # We initialize this float with a single Qubit which will be freed up
    # once the algorithm produced the first bit
    # We can't initialize variables with 0 qubits
    from qrisp.qtypes.quantum_float import QuantumFloat

    quotient = QuantumFloat(1, n, signed=True)

    # In order to determine the datashape of the remainder, we note that
    # within the algorithm, there are steps where the remainder variable needs to hold
    # larger+finer numbers than the final result for the remainder

    # The crucial step is the addition
    # remainder += divisor

    # Regarding the exponent of the remainder, we know that it has to be
    # less than the exponent of the numerator (so we can properly initialize)
    # But also less than the divisor exponent - 1 in order to properly execute
    # the last iteration which is a D/2 addition
    remainder_exponent = min(numerator.exponent, divisor.exponent - 1)

    # For the maximum significance,
    # we note that the divisor is initially bit shifted by n
    divisor.exp_shift(n - 1)

    # Therefore, the required maximum significance is acquired
    # with the formula for addition max_sig_add = max(max_sig_a, max_sig_b) + 1
    remainder_max_sig = max(numerator.mshape[1], divisor.mshape[1]) + 1

    # Create remainder float
    remainder_size = remainder_max_sig - remainder_exponent
    remainder = QuantumFloat(remainder_size, remainder_exponent, signed=True)

    # We now can initialize the remainder as the value of the numerator
    remainder.init_from(numerator)

    # If the divisor is signed we need to flip the quotient sign bit (why?)
    if divisor.signed:
        qs.cx(divisor[-1], quotient[-1])

    from qrisp.alg_primitives.arithmetic import inpl_add

    # Perform iterations
    for i in range(n - 1, -1, -1):
        if log_output:
            R = list(remainder.get_measurement().keys())[0]
            print("R: ", R)

            D = list(divisor.get_measurement().keys())[0]
            print("D: ", D)

            Q = list(quotient.get_measurement().keys())[0]
            print("Q: ", Q)

            # print("Q: " + bin_rep(int(Q/2.**quotient.exponent),
            # quotient.size-1) + "(" + str(Q) + ")")
            print("---")

        # The sign bit of the remainder is the newly calculated bit
        # of the quotient of this iteration
        # Therefore we need to move this bit from the remainder to the quotient
        # We do this in a very "hacky" fashion

        # Remove from remainder
        remainder_sign_bit = remainder.reg.pop(-1)

        # Add to quotient at position 0
        quotient.reg.insert(0, remainder_sign_bit)

        # This next instruction is a bit involved to understand

        # If we didn't use the technique of transfering the bits from the remainder
        # to the quotient, we would have to call quotient.x() in order to realize
        # the initial Q = 2**n - 1 statement. In this case we would then afterwards
        # CNOT from remainder sign bit it the quotient sign bit from this iteration
        # Since x gate on the quotient qubit and the CNOT commute,
        # we can also perform the x gate after the CNOT

        # Translating this to our case we replaced the CNOT with the bit transfer,
        # we arive at the neccessity of an x at this point
        x(quotient[0])

        # Since the new bit was added at the least significant end,
        # we need to lower the exponent
        quotient.exp_shift(-1)

        # This is to handle the D part of the statement "np.sign(D)*np.sign(R) != 1"
        if divisor.signed:
            qs.cx(divisor[-1], quotient[0])

        # We constructed Q with one initial qubit, which will be freed up after
        # the first iteration. This does not produce a qubit overhead,
        # because the upcomming add function needs more than one ancilla anyway
        if i == n - 1:
            quotient.reduce(quotient[1])

        # Perform the D = D/2 instruction
        divisor.exp_shift(-1)

        # We wrap the remainder in CNOT gates to make use of
        # (R' + D)' = (R - D)
        # where the prime denotes the negation
        for j in range(remainder.size):
            qs.cx(quotient[0], remainder[j])

        # This perform the addition R += D or R -= D
        inpl_add(
            remainder,
            divisor,
            ignore_rounding_error=True,
            ignore_overflow_error=True,
            adder=adder,
        )

        # Instead of executing this layer of CNOT gates,
        # we use the layer from the next iteration and only
        # CNOT the control qubit from this layer
        # for j in range(remainder.size):
        # qs.cx(quotient[0], remainder[j])

        # We still need to make up for the fact that the next layer of CNOT
        # gates doesnt include remainder[-1]

        # qs.cx(quotient[0], remainder[-1])

        # However remainder[-1] is also the control qubit of the next iteration
        # implying our plan of CNOTting the next control qubit cancels the
        # previous line

        # qs.cx(quotient[0], remainder[-1])

        # In any iteration apart from the first, quotient[0] does not contain
        # the desired value because of the sheenanigans of the previous lines.
        # In detail: It contains the desired value (+) the quotient[0] qubit
        # of the previous iteration. We cancel with this command
        if i != n - 1:
            qs.cx(quotient[1], quotient[0])

    # Since there is no more upcoming iteration, we simply perform the CNOT layer
    # without any further tricks.
    for j in range(remainder.size):
        qs.cx(quotient[0], remainder[j])

    # This performs the Q -= 2**(n-1) instruction (the -1 qubit is the sign bit =>
    # the -2 qubit is the significance n-1 bit)
    x(quotient[-2])

    # Handle the case that the divisor is signed
    if divisor.signed:
        # Create Qubit to account for the fact that we will add 0.5 to the quotient,
        # even though it only supports integers
        quotient.extend(1, position=0)
        quotient.exp_shift(-1)
        x(quotient[0])

        # Add/Subtract 0.5 executed via (Q' + 0.5)' = Q - 0.5
        qs.cx(divisor[-1], quotient.reg)
        quotient.incr(-(2**quotient.exponent))
        qs.cx(divisor[-1], quotient.reg)

        # Remove 0.5 significance qubit (contains 0 now)
        quotient.reduce(quotient[0])
        quotient.exp_shift(1)

        # Incase the divisor is not only signed but actually negative, we CNOT the
        # remainder such that the upcoming addition circuit results in a subtraction
        qs.cx(divisor[-1], remainder.reg)

    # Perform R = R +- D
    inpl_add(
        remainder,
        divisor,
        ignore_rounding_error=True,
        ignore_overflow_error=True,
        adder=adder,
    )

    if divisor.signed:
        # Undo the CNOT
        pass
        qs.cx(divisor[-1], remainder.reg)

    # Perform D = 2*D
    divisor.exp_shift(1)

    # Finally, if the numerator was signed, we flip the quotient sign bit
    if numerator.signed:
        qs.cx(numerator[-1], quotient[-1])

    return quotient, remainder


# Wrapper for the integer division
# Supports division for quantum float of arbitrary exponent.
# And allows to give a precision threshold (prec)
# The resulting quotient variable Q_res has the property |Q_res - Q_real| < 2**prec

[docs]
def q_divmod(numerator, divisor, adder="thapliyal", prec=0):
    """
    Performs division up to arbitrary precision. Returns the quotient and the remainder.

    Parameters
    ----------
    numerator : QuantumFloat
        The QuantumFloat to divide.
    divisor : QuantumFloat
        The QuantumFloat to divide by.
    adder : str, optional
        The type of adder to use. Available are "thapliyal" and "cuccarro".
        The default is "thapliyal".
    prec : int, optional
        The precision of the division. If the precision is set to $k$,
        the approximated quotient $q_{apr}$ and the true quotient $q_{true}$
        satisfy $|q_{apr} - q_{true}|<2^{-k}$. The default is 0.

    Returns
    -------
    quotient : QuantumFloat
        The approximated quotient.
    remainder : QuantumFloat
        The remainder which satisfying $q*d + r = n$.

    Examples
    --------

    We calculate 10/8 with varying precision:

    >>> from qrisp import QuantumFloat, q_divmod, multi_measurement
    >>> num = QuantumFloat(4)
    >>> div = QuantumFloat(4)
    >>> num[:] = 10
    >>> div[:] = 8
    >>> quotient, remainder = q_divmod(num, div, prec = 1)
    >>> multi_measurement([quotient, remainder])
    {(1.0, 2.0): 1.0}

    Now with higher precision

    >>> num = QuantumFloat(4)
    >>> div = QuantumFloat(4)
    >>> num[:] = 10
    >>> div[:] = 8
    >>> quotient, remainder = q_divmod(num, div, prec = 3)
    >>> multi_measurement([quotient, remainder])
    {(1.25, 0.0): 1.0}

    """

    # The idea is to bit shift numerator and divisor by s such that,
    # they are both integers. To increase the precision we shift N by -prec
    # (we'll see shortly why this increases the precision).
    # We write

    # N_tilde = 2**(s-prec) * N
    # D_tilde = 2**s * D

    # We then perform integer division giving Q_tilde, R_tilde such that

    # N_tilde/D_tilde = Q_tilde + R_tilde/D_tilde
    # where |R_tilde| < |D_tilde|

    # We now set

    # Q = Q_tilde 2**(prec)
    # and insert, which yields

    # N_tilde/D_tilde = 2**(-prec) * N/D  = 2**(-prec) * Q + R_tilde/D_tilde

    # =>  N/D = Q + 2**prec*R_tilde/D_tilde

    # So the error is
    # |2**prec*R_tilde/D_tilde| < 2**prec

    # Because R_tilde < D_tilde

    prec = -prec
    s = -min(numerator.exponent, divisor.exponent)

    num_exp_shift = s - prec
    div_exp_shift = s

    numerator.exp_shift(num_exp_shift)
    divisor.exp_shift(div_exp_shift)

    # print("N_tilde: ", numerator.get_measurement())
    # print("D_tilde: ", divisor.get_measurement())

    quotient, remainder = q_int_div(numerator, divisor, adder=adder, log_output=False)

    # print("Q_tilde: ", quotient.get_measurement())
    # print("R_tilde: ", remainder.get_measurement())

    quotient.exp_shift(prec)
    remainder.exp_shift(-num_exp_shift)

    divisor.exp_shift(-div_exp_shift)
    numerator.exp_shift(-num_exp_shift)

    return quotient, remainder




[docs]
def q_div(numerator, divisor, prec=None):
    """
    Performs division up to arbitrary precision and uncomputes the remainder.

    Parameters
    ----------
    numerator : QuantumFloat
        The QuantumFloat to divide.
    divisor : QuantumFloat
        The QuantumFloat to divide by.
    prec : int, optional
        The precision of the division. If the precision is set to $k$,
        the approximated quotient $q_{apr}$ and the true quotient $q_{true}$
        satisfy $|q_{apr} - q_{true}|<2^{-k}$.
        By default, a suited precision will be determined from the other inputs.

    Returns
    -------
    quotient : QuantumFloat
        The result of the division.

    Examples
    --------

    We calculate 10/8:

    >>> from qrisp import QuantumFloat, q_div
    >>> num = QuantumFloat(4)
    >>> div = QuantumFloat(4)
    >>> num[:] = 10
    >>> div[:] = 8
    >>> quotient = q_div(num, div, prec = 2)
    >>> print(quotient)
    {1.25: 1.0}

    """

    from qrisp import U_g_inpl_adder, h, hybrid_mult

    if prec is None:
        prec = divisor.size - numerator.exponent

    quotient, remainder = q_divmod(numerator, divisor, prec=prec)

    hybrid_mult(quotient, divisor, remainder, init_op="qft", terminal_op=None)

    U_g_inpl_adder(remainder, numerator, mult_factor=-1)

    h(remainder)
    # QFT(remainder, inv = True, exec_swap = False)
    # remainder += -1

    remainder.delete()

    return quotient




[docs]
def qf_inversion(qf, prec=None):
    """
    Calculates the multiplicative inverse of a QuantumFloat.

    Parameters
    ----------
    qf : QuantumFloat
        The QuantumFloat to invert.
    prec : int, optional
        The precision of the inversion. If the precision is set to $k$,
        the approximated inverse $q_{apr}$ and the true inverse $q_{true}$ satisfy
        $|q_{res} - q_{true}|<2^{-k}$.
        By default, a suited precision will be determined from the other input.

    Returns
    -------
    result : QuantumFloat
        A QuantumFloat containing the inverse.

    Examples
    --------

    We calculate the inverse of 0.75

    >>> from qrisp import QuantumFloat, qf_inversion
    >>> qf = QuantumFloat(2, -2)
    >>> qf[:] = 0.75
    >>> qf_inv = qf_inversion(qf, prec = 4)
    >>> print(qf_inv)
    {1.3125: 1.0}
    >>> 0.75**-1
    1.3333333333333333

    """

    if prec is None:
        prec = qf.size

    from qrisp import QuantumFloat, x

    numerator = QuantumFloat(1, signed=False)
    numerator.encode(1)

    result = q_div(numerator, qf, prec)

    x(numerator)
    numerator.delete()

    return result