nlft_qsp

Modules:

Name	Description
approximate	Provides functions to compute Laurent/Chebyshev/Fourier approximations.
file	Module for JSON serialization of various classes in the package.
nlft	Defining the nonlinear Fourier transform, as well as functions to compute the forward NLFT.
numerics	Module dealing with floating point types and error tolerance.
plot	Helper functions to plot polynomials using matplotlib.
poly	Definitions of complex sequences, polynomials and Chebyshev expansions, as well as all operations between them.
qsp	Definitions of quantum signal processing protocols and public API for QSP solvers.
rand	Random generators (mainly used for testing)
solvers	Module containing all the solvers for polynomial completion and inverse nonlinear Fourier transform/QSP synthesis.
util	Utility functions.

Classes:

Name	Description
ChebyshevQSPPhaseFactors	Phase factors for a Chebyshev QSP protocol.
ChebyshevTExpansion	Linear combination of Chebyshev polynomials of the first kind.
GQSPPhaseFactors	Phase factors for a Generalized QSP protocol.
NonLinearFourierSequence	Class representing a finitely supported sequence of complex numbers over \(\mathbb{Z}\).
PhaseFactors	Set of phase factors for a general Quantum Signal Processing protocol.
Polynomial	Represents a general Laurent polynomial of one complex variable.
QSVTPhaseFactors	Phase factors for a QSVT/Reflection QSP protocol.
XQSPPhaseFactors	Phase factors for a XQSP protocol.
YQSPPhaseFactors	Phase factors for a YQSP protocol.

Functions:

Name	Description
chebqsp_approximate	DEPRECATED: use ChebyshevQSPPhaseFactors.approximate() instead.
chebqsp_solve	DEPRECATED: use ChebyshevQSPPhaseFactors.solve() instead.
chebyshev_approximate	Computes the Chebyshev expansion up to \(N\) for a complex-valued function \(f : [-1, 1] \rightarrow \mathbb{C}\).
fourier_approximate	Computes the Fourier series of the given function \(f(z)\), \(z = e^{i\theta}\) being a complex number of unit modulus.
gqsp_solve	DEPRECATED: use GQSPPhaseFactors.solve() instead.
plot_chebyshev	Plots the real part of each object in funcs over the interval \([-1, 1]\).
plot_fourier	Plots the absolute value of each object in funcs over the unit circle, i.e., plugging \(z = e^{it}\) for \(t \in [-\pi, \pi]\).
qsvt_solve	DEPRECATED: use QSVTPhaseFactors.solve() instead.
xqsp_solve	DEPRECATED: use XQSPPhaseFactors.solve() instead.
xqsp_solve_laurent	DEPRECATED: use XQSPPhaseFactors.solve_laurent() instead.
yqsp_solve	DEPRECATED: use YQSPPhaseFactors.solve() instead.
yqsp_solve_laurent	DEPRECATED: use YQSPPhaseFactors.solve_laurent() instead.

ChebyshevQSPPhaseFactors

Bases: XQSPPhaseFactors

Phase factors for a Chebyshev QSP protocol.

\[ e^{i\phi_0 Z} \tilde{x} e^{i\phi_1 Z} \tilde{x} \cdots \tilde{x} e^{i\phi_n Z} = \begin{pmatrix} P(x) & iQ(x) \sqrt{1 - x^2} \\ iQ^*(x) \sqrt{1 - x^2} & P^*(x) \end{pmatrix} \]

where the signal operator \(\tilde{x}\) is

\[ \tilde{x} = \begin{pmatrix} x & i\sqrt{1 - x^2} \\ i\sqrt{1 - x^2} & x \end{pmatrix} \]

Note

This is the ansatz of arXiv:2105.02859 Theorem 9, but the polynomial construction is implemented by conjugating XQSP with a Hadamard gate.

Methods:

Name	Description
approximate	Approximate the given callable object \(f\) (which takes \(x \in [-1, 1]\) and returns a real number)
solve	Returns the set of phase factors for a Chebyshev QSP protocol implementing the polynomial \(P(x)\) (as the real part of the top-left polynomial, see Theorem 9 of arXiv:2105.02859).

Source code in nlft_qsp/qsp.py

@serializable(
    type_tag="@qspx/phase_factors/chebqsp",
    fields={"phi": "phi"}
)
@qsp_variant('cheb', modes=['c'], display_name='Chebyshev QSP')
class ChebyshevQSPPhaseFactors(XQSPPhaseFactors):
    r"""Phase factors for a Chebyshev QSP protocol.

    $$ e^{i\phi_0 Z} \tilde{x} e^{i\phi_1 Z} \tilde{x} \cdots \tilde{x} e^{i\phi_n Z} = \begin{pmatrix} P(x) & iQ(x) \sqrt{1 - x^2} \\ iQ^*(x) \sqrt{1 - x^2} & P^*(x) \end{pmatrix} $$

    where the signal operator $\tilde{x}$ is

    $$ \tilde{x} = \begin{pmatrix} x & i\sqrt{1 - x^2} \\ i\sqrt{1 - x^2} & x \end{pmatrix} $$

    Note:
        This is the ansatz of [arXiv:2105.02859](https://arxiv.org/abs/2105.02859) Theorem 9, but the polynomial
        construction is implemented by conjugating XQSP with a Hadamard gate."""
    def processing_operator(self, k: int): # exp(i phi[k] Z)
        return np.exp(1j*self.phi[k]), 0

    def protocol_conjugation(self, P, Q): # Return the polynomials (P', Q') = H (P, Q) H
        return ((P + P.conjugate()) + (Q - Q.conjugate()))/2, ((P - P.conjugate()) - (Q + Q.conjugate()))/2

    def processing_operator_conjugation(self, a, b):
        return np.real(a) + 1j*np.imag(b), 1j*np.imag(a) - np.real(b)

    def duplicate(self):
        return ChebyshevQSPPhaseFactors(self.phi)

    def iX(self):
        raise ValueError("Multiplying by iX is not possible.")

    def iZ(self):
        return super.iX() # applying iZ is equivalent to applying iX, by the Hadamard conjugation

    @classmethod
    def solve(cls, T: list[complex_type] | Polynomial | ChebyshevTExpansion) -> "ChebyshevQSPPhaseFactors":
        r"""Returns the set of phase factors for a Chebyshev QSP protocol implementing the polynomial $P(x)$ (as the real part of the top-left polynomial, see Theorem 9 of [arXiv:2105.02859](https://arxiv.org/abs/2105.02859)).

        The target polynomial will be $T(x) = \sum_{k = 0}^n c_k T_k(x)$ (if `T` is a `ChebyshevTExpansion`) or $T(x) = \sum_{k = 0}^n c_k x^k$ (if `T` is a `Polynomial`), where $T_k(x)$ are the Chebyshev polynomials of the first kind.

        Args:
            T: a Chebyshev expansion object or the desired Polynomial $P(x)$ (that will be converted to the Chebyshev basis).

        Raises:
            ValueError: If the target polynomial does not have definite parity or is not real.

        Note:
            `T` can also be a list of complex numbers. This will be regarded as coefficients in the Chebyshev basis. Passing the list directly is discouraged and will be removed in future releases."""
        if isinstance(T, list):
            T = ChebyshevTExpansion(T)

        if isinstance(T, Polynomial):
            T = ChebyshevTExpansion.from_polynomial(T)

        if not T.is_real():
            raise ValueError("Only real polynomials are supported.")

        xqsp = XQSPPhaseFactors.solve_laurent(T.to_laurent())
        return ChebyshevQSPPhaseFactors(xqsp.phi)

    @classmethod
    def approximate(cls, f: Callable, deg: int) -> "ChebyshevQSPPhaseFactors":
        r"""Approximate the given callable object $f$ (which takes $x \in [-1, 1]$ and returns a real number)
        and returns the Chebyshev QSP phase factors implementing an approximating polynomial of degree `deg`
        (as the real part of the top-left polynomial, see Theorem 9 of [arXiv:2105.02859](https://arxiv.org/abs/2105.02859)).

        Note: The parity of `deg` should coincide with the parity of $f$, otherwise the Chebyshev approximator might give numerical errors."""
        return cls.solve(chebyshev_approximate(f, deg))

approximate(f: Callable, deg: int) -> ChebyshevQSPPhaseFactors classmethod

Approximate the given callable object \(f\) (which takes \(x \in [-1, 1]\) and returns a real number) and returns the Chebyshev QSP phase factors implementing an approximating polynomial of degree deg (as the real part of the top-left polynomial, see Theorem 9 of arXiv:2105.02859).

Note: The parity of deg should coincide with the parity of \(f\), otherwise the Chebyshev approximator might give numerical errors.

Source code in nlft_qsp/qsp.py

@classmethod
def approximate(cls, f: Callable, deg: int) -> "ChebyshevQSPPhaseFactors":
    r"""Approximate the given callable object $f$ (which takes $x \in [-1, 1]$ and returns a real number)
    and returns the Chebyshev QSP phase factors implementing an approximating polynomial of degree `deg`
    (as the real part of the top-left polynomial, see Theorem 9 of [arXiv:2105.02859](https://arxiv.org/abs/2105.02859)).

    Note: The parity of `deg` should coincide with the parity of $f$, otherwise the Chebyshev approximator might give numerical errors."""
    return cls.solve(chebyshev_approximate(f, deg))

solve(T: list[complex_type] | Polynomial | ChebyshevTExpansion) -> ChebyshevQSPPhaseFactors classmethod

Returns the set of phase factors for a Chebyshev QSP protocol implementing the polynomial \(P(x)\) (as the real part of the top-left polynomial, see Theorem 9 of arXiv:2105.02859).

The target polynomial will be \(T(x) = \sum_{k = 0}^n c_k T_k(x)\) (if T is a ChebyshevTExpansion) or \(T(x) = \sum_{k = 0}^n c_k x^k\) (if T is a Polynomial), where \(T_k(x)\) are the Chebyshev polynomials of the first kind.

Parameters:

Name	Type	Description	Default
T	list[complex_type] | Polynomial | ChebyshevTExpansion	a Chebyshev expansion object or the desired Polynomial \(P(x)\) (that will be converted to the Chebyshev basis).	required

Raises:

Type	Description
ValueError	If the target polynomial does not have definite parity or is not real.

Note

T can also be a list of complex numbers. This will be regarded as coefficients in the Chebyshev basis. Passing the list directly is discouraged and will be removed in future releases.

Source code in nlft_qsp/qsp.py

@classmethod
def solve(cls, T: list[complex_type] | Polynomial | ChebyshevTExpansion) -> "ChebyshevQSPPhaseFactors":
    r"""Returns the set of phase factors for a Chebyshev QSP protocol implementing the polynomial $P(x)$ (as the real part of the top-left polynomial, see Theorem 9 of [arXiv:2105.02859](https://arxiv.org/abs/2105.02859)).

    The target polynomial will be $T(x) = \sum_{k = 0}^n c_k T_k(x)$ (if `T` is a `ChebyshevTExpansion`) or $T(x) = \sum_{k = 0}^n c_k x^k$ (if `T` is a `Polynomial`), where $T_k(x)$ are the Chebyshev polynomials of the first kind.

    Args:
        T: a Chebyshev expansion object or the desired Polynomial $P(x)$ (that will be converted to the Chebyshev basis).

    Raises:
        ValueError: If the target polynomial does not have definite parity or is not real.

    Note:
        `T` can also be a list of complex numbers. This will be regarded as coefficients in the Chebyshev basis. Passing the list directly is discouraged and will be removed in future releases."""
    if isinstance(T, list):
        T = ChebyshevTExpansion(T)

    if isinstance(T, Polynomial):
        T = ChebyshevTExpansion.from_polynomial(T)

    if not T.is_real():
        raise ValueError("Only real polynomials are supported.")

    xqsp = XQSPPhaseFactors.solve_laurent(T.to_laurent())
    return ChebyshevQSPPhaseFactors(xqsp.phi)

ChebyshevTExpansion

Bases: ComplexL0Sequence

Linear combination of Chebyshev polynomials of the first kind.

Parameters:

Name	Type	Description	Default
coeffs	list[complex_type] | Polynomial	Either the coefficients of the linear combination, or the symmetric Laurent polynomial \(P(z)\) which is equal up to the change of variable \(x = \frac{z + z^{-1}}{2}\).	required

Note

The change of variable between Laurent polynomials and Chebyshev expansions is given by the relation \(T_k(\cos \theta) = \cos k\theta = \frac{z^k + z^{-k}}{2}\). Thus the expansion is

\[\sum_{k = 0}^n c_k T_k(x) = \sum_{k = 0}^n c_k \frac{z^k + z^{-k}}{2}\]

Methods:

Name	Description
__call__	Evaluates the Chebyshev expansion at the given number.
__str__	Converts the expansion to a human-readable string representation.
from_laurent_polynomial	Returns the Chebyshev expansion \(T\) satisfying \(T(x) = \frac{P(z) + P^*(z)}{2}\).
from_polynomial	Returns the Chebyshev expansion \(T\) satisfying \(T(x) = P(x)\).
to_laurent	Returns the Laurent polynomial \(P(e^{i\theta}) = T(\cos \theta)\) where \(T(x)\) is represented by self.
to_polynomial	Returns the polynomial \(P\) satisfying \(P(x) = T(x)\).

Source code in nlft_qsp/poly.py

@serializable(
    type_tag="@qspx/chebyshev_t",
    fields={
        "coeffs": "coeffs"
    }
)
class ChebyshevTExpansion(ComplexL0Sequence):
    r"""Linear combination of Chebyshev polynomials of the first kind.

    Args:
        coeffs: Either the coefficients of the linear combination, or the symmetric Laurent polynomial $P(z)$ which is equal up to the change of variable $x = \frac{z + z^{-1}}{2}$.

    Note:
        The change of variable between Laurent polynomials and Chebyshev expansions is given by the relation $T_k(\cos \theta) = \cos k\theta = \frac{z^k + z^{-k}}{2}$. Thus the expansion is

        $$\sum_{k = 0}^n c_k T_k(x) = \sum_{k = 0}^n c_k \frac{z^k + z^{-k}}{2}$$
    """
    def __init__(self, coeffs: list[complex_type] | Polynomial):
        if isinstance(coeffs, list) or isinstance(coeffs, np.ndarray):
            super().__init__(coeffs, support_start=0)
        elif isinstance(coeffs, Polynomial):
            if not coeffs.is_symmetric():
                raise ValueError("The given Laurent polynomial is not symmetric.")

            coeffs = [2*coeffs[k] for k in range(coeffs.support().stop)]
            coeffs[0] /= 2
            super().__init__(coeffs, support_start=0)
        else:
            raise ValueError("Only a coefficient vector or symmetric Laurent polynomials are allowed.")

    def degree(self) -> int:
        return self.coeffs.shape[0] - 1

    def __call__(self, x: float_type) -> complex_type:
        """Evaluates the Chebyshev expansion at the given number.

        Args:
            x (float_type): The point at which to evaluate the expansion.

        Returns:
            complex: The evaluated result.
        """
        theta = np.arccos(x)

        return sum(self[k] * np.cos(k * theta) for k in self.support())

    def __str__(self) -> str:
        """Converts the expansion to a human-readable string representation.

        Returns:
            str: The string representation of the expansion.
        """
        return ' + '.join(f"{c} T_{self.support_start + k}(x)" for k, c in enumerate(self.coeffs))

    def to_laurent(self) -> Polynomial:
        r"""Returns the Laurent polynomial $P(e^{i\theta}) = T(\cos \theta)$ where $T(x)$ is represented by `self`."""
        P = Polynomial(np.concat([np.flipud(self.coeffs), self.coeffs[1:]]), support_start=-self.coeffs.shape[0]+1)
        P[0] *= 2
        return P/2

    @classmethod
    def from_polynomial(cls, P: Polynomial) -> "ChebyshevTExpansion":
        r"""Returns the Chebyshev expansion $T$ satisfying $T(x) = P(x)$."""
        return ChebyshevTExpansion(np.polynomial.chebyshev.poly2cheb(P.coeffs))

    @classmethod
    def from_laurent_polynomial(cls, P: Polynomial):
        r"""Returns the Chebyshev expansion $T$ satisfying $T(x) = \frac{P(z) + P^*(z)}{2}$.

        Raises:
            ValueError: if $P$ is not symmetric.

        Note: $P$ must be symmetric."""
        if not P.is_symmetric():
            raise ValueError("The given Laurent polynomial is not symmetric.")

        coeffs = [2*P[k] for k in range(P.support().stop)]
        coeffs[0] /= 2

        return ChebyshevTExpansion(coeffs)

    def to_polynomial(self) -> Polynomial:
        """Returns the polynomial $P$ satisfying $P(x) = T(x)$."""
        return Polynomial(np.polynomial.chebyshev.cheb2poly(self.coeffs))

__call__(x: float_type) -> complex_type

Evaluates the Chebyshev expansion at the given number.

Parameters:

Name	Type	Description	Default
x	float_type	The point at which to evaluate the expansion.	required

Returns:

Name	Type	Description
complex	complex_type	The evaluated result.

Source code in nlft_qsp/poly.py

def __call__(self, x: float_type) -> complex_type:
    """Evaluates the Chebyshev expansion at the given number.

    Args:
        x (float_type): The point at which to evaluate the expansion.

    Returns:
        complex: The evaluated result.
    """
    theta = np.arccos(x)

    return sum(self[k] * np.cos(k * theta) for k in self.support())

__str__() -> str

Converts the expansion to a human-readable string representation.

Returns:

Name	Type	Description
str	str	The string representation of the expansion.

Source code in nlft_qsp/poly.py

def __str__(self) -> str:
    """Converts the expansion to a human-readable string representation.

    Returns:
        str: The string representation of the expansion.
    """
    return ' + '.join(f"{c} T_{self.support_start + k}(x)" for k, c in enumerate(self.coeffs))

from_laurent_polynomial(P: Polynomial) classmethod

Returns the Chebyshev expansion \(T\) satisfying \(T(x) = \frac{P(z) + P^*(z)}{2}\).

Raises:

Type	Description
ValueError	if \(P\) is not symmetric.

Note: \(P\) must be symmetric.

Source code in nlft_qsp/poly.py

@classmethod
def from_laurent_polynomial(cls, P: Polynomial):
    r"""Returns the Chebyshev expansion $T$ satisfying $T(x) = \frac{P(z) + P^*(z)}{2}$.

    Raises:
        ValueError: if $P$ is not symmetric.

    Note: $P$ must be symmetric."""
    if not P.is_symmetric():
        raise ValueError("The given Laurent polynomial is not symmetric.")

    coeffs = [2*P[k] for k in range(P.support().stop)]
    coeffs[0] /= 2

    return ChebyshevTExpansion(coeffs)

from_polynomial(P: Polynomial) -> ChebyshevTExpansion classmethod

Returns the Chebyshev expansion \(T\) satisfying \(T(x) = P(x)\).

Source code in nlft_qsp/poly.py

@classmethod
def from_polynomial(cls, P: Polynomial) -> "ChebyshevTExpansion":
    r"""Returns the Chebyshev expansion $T$ satisfying $T(x) = P(x)$."""
    return ChebyshevTExpansion(np.polynomial.chebyshev.poly2cheb(P.coeffs))

to_laurent() -> Polynomial

Returns the Laurent polynomial \(P(e^{i\theta}) = T(\cos \theta)\) where \(T(x)\) is represented by self.

Source code in nlft_qsp/poly.py

def to_laurent(self) -> Polynomial:
    r"""Returns the Laurent polynomial $P(e^{i\theta}) = T(\cos \theta)$ where $T(x)$ is represented by `self`."""
    P = Polynomial(np.concat([np.flipud(self.coeffs), self.coeffs[1:]]), support_start=-self.coeffs.shape[0]+1)
    P[0] *= 2
    return P/2

to_polynomial() -> Polynomial

Returns the polynomial \(P\) satisfying \(P(x) = T(x)\).

Source code in nlft_qsp/poly.py

def to_polynomial(self) -> Polynomial:
    """Returns the polynomial $P$ satisfying $P(x) = T(x)$."""
    return Polynomial(np.polynomial.chebyshev.cheb2poly(self.coeffs))

GQSPPhaseFactors

Bases: PhaseFactors

Phase factors for a Generalized QSP protocol.

\[ e^{i\lambda Z} e^{i\phi_0 X} e^{i\theta_0 Z} W(z) e^{i\phi_1 X} e^{i\theta_1 Z} W(z) \cdots W(z) e^{i\phi_n X} e^{i\theta_n Z} = \begin{pmatrix} P(z) & Q(z) \\ \cdot & \cdot \end{pmatrix} \]

where \(W(z) = \mathrm{diag}(z, 1)\) (mode='analytic') or \(W(z) = \mathrm{diag}(z, z^{-1})\) (mode='laurent').

Note

This class follows the convention of arXiv:2503.03026, Theorem 2, which is different from the original GQSP convention. If the convention of arXiv:2308.01501 Theorem 3 is desired, then one should use the to_mw_gqsp() method.

Methods:

Name	Description
from_nlfs	Computes the GQSP phase factors for a given NLFT sequence.
phase_offset	Returns the phase of the leading coefficient of P, where (P, Q) is the pair of polynomials generated by this set.
solve	Returns the set of phase factors for a Generalized QSP protocol producing the given polynomial. See here for an overview of the QSP variants.
solve_laurent	Returns the set of phase factors for a Generalized QSP protocol producing the given definite-parity polynomial. See here for an overview of the QSP variants.
to_mw_gqsp	Converts the GQSP phase factors to the convention of arXiv:2308.01501 Theorem 3.
to_xqsp	Converts the QSP phase factors into XQSP phase factors.
to_yqsp	Converts the QSP phase factors into YQSP phase factors.

Source code in nlft_qsp/qsp.py

@serializable(
    type_tag="@qspx/phase_factors/gqsp",
    fields={"phi": "phi", "lbd": "lbd", "theta": "theta"}
)
@qsp_variant('g', modes=['a', 'l'], display_name='Generalized QSP')
class GQSPPhaseFactors(PhaseFactors):
    r"""Phase factors for a Generalized QSP protocol.

    $$ e^{i\lambda Z} e^{i\phi_0 X} e^{i\theta_0 Z} W(z) e^{i\phi_1 X} e^{i\theta_1 Z} W(z) \cdots W(z) e^{i\phi_n X} e^{i\theta_n Z} = \begin{pmatrix} P(z) & Q(z) \\ \cdot & \cdot \end{pmatrix} $$

    where $W(z) = \mathrm{diag}(z, 1)$ (`mode='analytic'`) or $W(z) = \mathrm{diag}(z, z^{-1})$ (`mode='laurent'`).

    Note:
        This class follows the convention of [arXiv:2503.03026](https://arxiv.org/abs/2503.03026), Theorem 2, which is different from the original GQSP convention. If the convention of [arXiv:2308.01501](https://arxiv.org/abs/2308.01501) Theorem 3 is desired, then one should use the `to_mw_gqsp()` method."""
    def __init__(self, phi: list[float_type], lbd: float_type=0, theta: list[float_type]=None):
        self.lbd = lbd
        self.phi = list(phi)

        if theta is None:
            theta = [0] * len(phi)

        self.theta = list(theta)

        if len(theta) < len(phi):
            theta += [0] * (len(phi) - len(theta))

        if len(theta) > len(phi):
            phi += [0] * (len(theta) - len(phi))

    def duplicate(self):
        return GQSPPhaseFactors(self.phi, self.lbd, self.theta)

    def processing_operator(self, k: int):
        if k == 0: # exp(i lbd Z) exp(i phi X) exp(i theta Z)
            return np.exp(1j*(self.lbd + self.theta[k]))*np.cos(self.phi[k]), \
                1j*np.exp(1j*(self.lbd - self.theta[k]))*np.sin(self.phi[k])
        else: # exp(i phi X) exp(i theta Z)
            return np.exp(1j*self.theta[k])*np.cos(self.phi[k]), 1j*np.exp(-1j*self.theta[k])*np.sin(self.phi[k])

    def degree(self):
        return len(self.phi) - 1

    def iX(self):
        pf = self.duplicate()
        pf.theta[-1] = -pf.theta[-1]
        pf.phi[-1] += np.pi/2
        return pf

    def iY(self):
        return self.iZ().iX()

    def iZ(self):
        pf = self.duplicate()
        pf.theta[-1] += np.pi/2
        return pf

    def phase_offset(self) -> float_type:
        """Returns the phase of the leading coefficient of P, where (P, Q) is the pair of polynomials generated by this set."""
        return self.lbd + sum(self.theta)

    def to_xqsp(self):
        """Converts the QSP phase factors into XQSP phase factors.

        Raises:
            ValueError: If the phase factors do not lie in the XQSP subalgebra."""
        F = self.to_nlfs()
        if not F.is_imaginary():
            raise ValueError("The phase factors are not reducible to XQSP.")

        return XQSPPhaseFactors.from_nlfs(F)

    def to_yqsp(self):
        """Converts the QSP phase factors into YQSP phase factors.

        Raises:
            ValueError: If the phase factors do not lie in the YQSP subalgebra."""
        F = self.to_nlfs()
        if not F.is_real():
            raise ValueError("The phase factors are not reducible to YQSP.")

        return YQSPPhaseFactors.from_nlfs(F)

    def to_nlfs(self) -> NonLinearFourierSequence:
        n = self.degree()
        alpha = self.phase_offset() # phase of the leading coefficient of P

        psi = [0] * (n+1) # prefactors

        psi[0] = self.lbd - alpha/2
        for k in range(n):
            psi[k+1] = self.theta[k] + psi[k]

        phi = self.phi
        return NonLinearFourierSequence([1j*np.tan(phik)*np.exp(2j*psik) for phik, psik in zip(phi, psi)])

    @classmethod
    def from_nlfs(cls, F: NonLinearFourierSequence, alpha: float_type = 0) -> PhaseFactors:
        r"""Computes the GQSP phase factors for a given NLFT sequence.
        If $\mathrm{NLFT}(F) = (a, b)$, then the returned phase factors will implement $(e^{i \alpha} z^n a, b)$ in the analytic picture.

        Args:
            F (NonLinearFourierSequence): The sequence to be converted to phase factors.
            alpha (float_type): In the pair of polynomials $(P, Q)$ generated by the returned phase factors, $P$ will be multiplied by `e^{i \alpha}`.

        Note:
            The support start of $F$ is ignored, so the support of $b$ is assumed to start at $0$."""
        psi = [phase_prefactor(Fk) for Fk in F.coeffs]
        lbd = psi[0]

        phi = [np.arctan(np.imag(Fk * np.exp(-2j * psik))) for Fk, psik in zip(F, psi)]

        psi += [0] # we add psi[n+1] to compute theta
        theta = [psi[k+1] - psi[k] for k in range(len(phi))]

        lbd += alpha/2
        theta[-1] += alpha/2
        return GQSPPhaseFactors(phi, lbd, theta)

    def to_mw_gqsp(self) -> tuple[list[float_type], list[float_type], float_type]:
        r"""Converts the GQSP phase factors to the convention of [arXiv:2308.01501](https://arxiv.org/abs/2308.01501) Theorem 3.

        In particular, a protocol should be constructed of the form

        $$ R(\theta_0, \phi_0, \lambda) \tilde{w} R(\theta_1, \phi_1, 0) \tilde{w} \cdots \tilde{w} R(\phi_n, \theta_n, 0) = \begin{pmatrix} P(z) & Q(z) \\ \cdot & \cdot \end{pmatrix} $$

        where the matrix $R$ is of the form

        $$ R(\theta, \phi, \lambda) = \begin{pmatrix} e^{i(\lambda + \phi)} \cos \theta & e^{i\phi} \sin \theta \\ e^{i\lambda} \sin \theta & -\cos \theta \end{pmatrix} $$

        Returns:
            The phase factors \vec{\theta}, \vec{\phi} (as lists) and \lambda, in this order."""
        n = self.degree()

        alpha = - self.lbd - sum(self.theta) - (n+1) * np.pi # needed to counterbalance the global phase

        mw_theta = [self.phi[k] for k in range(n+1)]
        mw_phi = [2 * self.lbd - np.pi/2 + alpha] + [2 * self.theta[k-1] - np.pi for k in range(1, n+1)]
        mw_lbd = 2 * self.theta[n] - np.pi/2 + alpha

        return mw_theta, mw_phi, mw_lbd

    @classmethod
    def solve(cls, P: Polynomial, convention='qsp') -> "GQSPPhaseFactors":
        r"""Returns the set of phase factors for a Generalized QSP protocol producing the given polynomial. See [here](https://arxiv.org/abs/2503.03026) for an overview of the QSP variants.
        A complementary $Q$ will be computed with Weiss' algorithm.

        Args:
            convention (str): Whether the phase factors should produce $(P, Q)$ (`'qsp'`), or $(Q, P)$ (`'nlft'`).

        Note:
            The sup norm of $P$ should be bounded by $1 - \eta < 1$.
            The time required by the algorithm to compute the phase factors will scale with $1/\eta$.

            The support_start of $P$ will be ignored."""
        if 1 - P.sup_norm(4*P.effective_degree()) < bd.machine_threshold():
            raise ValueError("The given polynomial cannot be too close to or larger than one on the unit circle.")

        match convention:
            case "qsp":
                P = -1j * Polynomial(P.coeffs, 0) # (Q, -iP) -> (P, iQ)
            case "nlft":
                P = Polynomial(P.coeffs, 0)
            case _:
                raise ValueError("The given mode does not exist. Only modes available are 'qsp', 'nlft'.")

        F = _riemann_hilbert_weiss(P) # NLFT(F) = (Q, P)

        if convention != "qsp":
            return GQSPPhaseFactors.from_nlfs(F)
        return GQSPPhaseFactors.from_nlfs(F).iX()

    @classmethod
    def solve_laurent(cls, P: Polynomial, convention='qsp') -> "GQSPPhaseFactors":
        r"""Returns the set of phase factors for a Generalized QSP protocol producing the given definite-parity polynomial. See [here](https://arxiv.org/abs/2503.03026) for an overview of the QSP variants.
        A complementary $Q$ will be computed with Weiss' algorithm.

        Args:
            convention (str): Whether the phase factors should produce $(P, Q)$ (`'qsp'`), or $(Q, P)$ (`'nlft'`).

        Raises:
            ValueError: If $P$ does not have definite parity.

        Note:
            The sup norm of $P$ should be bounded by $1 - \eta < 1$.
            The time required by the algorithm to compute the phase factors will scale with $1/\eta$.
            """
        if not is_definite_parity(P):
            raise ValueError("Laurent polynomial is not of definite parity.")

        return cls.solve(laurent_to_analytic(P), convention)

from_nlfs(F: NonLinearFourierSequence, alpha: float_type = 0) -> PhaseFactors classmethod

Computes the GQSP phase factors for a given NLFT sequence. If \(\mathrm{NLFT}(F) = (a, b)\), then the returned phase factors will implement \((e^{i \alpha} z^n a, b)\) in the analytic picture.

Parameters:

Name	Type	Description	Default
F	NonLinearFourierSequence	The sequence to be converted to phase factors.	required
alpha	float_type	In the pair of polynomials \((P, Q)\) generated by the returned phase factors, \(P\) will be multiplied by e^{i \alpha}.	0

Note

The support start of \(F\) is ignored, so the support of \(b\) is assumed to start at \(0\).

Source code in nlft_qsp/qsp.py

@classmethod
def from_nlfs(cls, F: NonLinearFourierSequence, alpha: float_type = 0) -> PhaseFactors:
    r"""Computes the GQSP phase factors for a given NLFT sequence.
    If $\mathrm{NLFT}(F) = (a, b)$, then the returned phase factors will implement $(e^{i \alpha} z^n a, b)$ in the analytic picture.

    Args:
        F (NonLinearFourierSequence): The sequence to be converted to phase factors.
        alpha (float_type): In the pair of polynomials $(P, Q)$ generated by the returned phase factors, $P$ will be multiplied by `e^{i \alpha}`.

    Note:
        The support start of $F$ is ignored, so the support of $b$ is assumed to start at $0$."""
    psi = [phase_prefactor(Fk) for Fk in F.coeffs]
    lbd = psi[0]

    phi = [np.arctan(np.imag(Fk * np.exp(-2j * psik))) for Fk, psik in zip(F, psi)]

    psi += [0] # we add psi[n+1] to compute theta
    theta = [psi[k+1] - psi[k] for k in range(len(phi))]

    lbd += alpha/2
    theta[-1] += alpha/2
    return GQSPPhaseFactors(phi, lbd, theta)

phase_offset() -> float_type

Returns the phase of the leading coefficient of P, where (P, Q) is the pair of polynomials generated by this set.

Source code in nlft_qsp/qsp.py

def phase_offset(self) -> float_type:
    """Returns the phase of the leading coefficient of P, where (P, Q) is the pair of polynomials generated by this set."""
    return self.lbd + sum(self.theta)

solve(P: Polynomial, convention='qsp') -> GQSPPhaseFactors classmethod

Returns the set of phase factors for a Generalized QSP protocol producing the given polynomial. See here for an overview of the QSP variants. A complementary \(Q\) will be computed with Weiss' algorithm.

Parameters:

Name	Type	Description	Default
convention	str	Whether the phase factors should produce \((P, Q)\) ('qsp'), or \((Q, P)\) ('nlft').	'qsp'

Note

The sup norm of \(P\) should be bounded by \(1 - \eta < 1\). The time required by the algorithm to compute the phase factors will scale with \(1/\eta\).

The support_start of \(P\) will be ignored.

Source code in nlft_qsp/qsp.py

@classmethod
def solve(cls, P: Polynomial, convention='qsp') -> "GQSPPhaseFactors":
    r"""Returns the set of phase factors for a Generalized QSP protocol producing the given polynomial. See [here](https://arxiv.org/abs/2503.03026) for an overview of the QSP variants.
    A complementary $Q$ will be computed with Weiss' algorithm.

    Args:
        convention (str): Whether the phase factors should produce $(P, Q)$ (`'qsp'`), or $(Q, P)$ (`'nlft'`).

    Note:
        The sup norm of $P$ should be bounded by $1 - \eta < 1$.
        The time required by the algorithm to compute the phase factors will scale with $1/\eta$.

        The support_start of $P$ will be ignored."""
    if 1 - P.sup_norm(4*P.effective_degree()) < bd.machine_threshold():
        raise ValueError("The given polynomial cannot be too close to or larger than one on the unit circle.")

    match convention:
        case "qsp":
            P = -1j * Polynomial(P.coeffs, 0) # (Q, -iP) -> (P, iQ)
        case "nlft":
            P = Polynomial(P.coeffs, 0)
        case _:
            raise ValueError("The given mode does not exist. Only modes available are 'qsp', 'nlft'.")

    F = _riemann_hilbert_weiss(P) # NLFT(F) = (Q, P)

    if convention != "qsp":
        return GQSPPhaseFactors.from_nlfs(F)
    return GQSPPhaseFactors.from_nlfs(F).iX()

solve_laurent(P: Polynomial, convention='qsp') -> GQSPPhaseFactors classmethod

Returns the set of phase factors for a Generalized QSP protocol producing the given definite-parity polynomial. See here for an overview of the QSP variants. A complementary \(Q\) will be computed with Weiss' algorithm.

Parameters:

Name	Type	Description	Default
convention	str	Whether the phase factors should produce \((P, Q)\) ('qsp'), or \((Q, P)\) ('nlft').	'qsp'

Raises:

Type	Description
ValueError	If \(P\) does not have definite parity.

Note

The sup norm of \(P\) should be bounded by \(1 - \eta < 1\). The time required by the algorithm to compute the phase factors will scale with \(1/\eta\).

Source code in nlft_qsp/qsp.py

@classmethod
def solve_laurent(cls, P: Polynomial, convention='qsp') -> "GQSPPhaseFactors":
    r"""Returns the set of phase factors for a Generalized QSP protocol producing the given definite-parity polynomial. See [here](https://arxiv.org/abs/2503.03026) for an overview of the QSP variants.
    A complementary $Q$ will be computed with Weiss' algorithm.

    Args:
        convention (str): Whether the phase factors should produce $(P, Q)$ (`'qsp'`), or $(Q, P)$ (`'nlft'`).

    Raises:
        ValueError: If $P$ does not have definite parity.

    Note:
        The sup norm of $P$ should be bounded by $1 - \eta < 1$.
        The time required by the algorithm to compute the phase factors will scale with $1/\eta$.
        """
    if not is_definite_parity(P):
        raise ValueError("Laurent polynomial is not of definite parity.")

    return cls.solve(laurent_to_analytic(P), convention)

to_mw_gqsp() -> tuple[list[float_type], list[float_type], float_type]

Converts the GQSP phase factors to the convention of arXiv:2308.01501 Theorem 3.

In particular, a protocol should be constructed of the form

\[ R(\theta_0, \phi_0, \lambda) \tilde{w} R(\theta_1, \phi_1, 0) \tilde{w} \cdots \tilde{w} R(\phi_n, \theta_n, 0) = \begin{pmatrix} P(z) & Q(z) \\ \cdot & \cdot \end{pmatrix} \]

where the matrix \(R\) is of the form

\[ R(\theta, \phi, \lambda) = \begin{pmatrix} e^{i(\lambda + \phi)} \cos \theta & e^{i\phi} \sin \theta \\ e^{i\lambda} \sin \theta & -\cos \theta \end{pmatrix} \]

Returns:

Type	Description
tuple[list[float_type], list[float_type], float_type]	The phase factors \vec{\theta}, \vec{\phi} (as lists) and \lambda, in this order.

Source code in nlft_qsp/qsp.py

def to_mw_gqsp(self) -> tuple[list[float_type], list[float_type], float_type]:
    r"""Converts the GQSP phase factors to the convention of [arXiv:2308.01501](https://arxiv.org/abs/2308.01501) Theorem 3.

    In particular, a protocol should be constructed of the form

    $$ R(\theta_0, \phi_0, \lambda) \tilde{w} R(\theta_1, \phi_1, 0) \tilde{w} \cdots \tilde{w} R(\phi_n, \theta_n, 0) = \begin{pmatrix} P(z) & Q(z) \\ \cdot & \cdot \end{pmatrix} $$

    where the matrix $R$ is of the form

    $$ R(\theta, \phi, \lambda) = \begin{pmatrix} e^{i(\lambda + \phi)} \cos \theta & e^{i\phi} \sin \theta \\ e^{i\lambda} \sin \theta & -\cos \theta \end{pmatrix} $$

    Returns:
        The phase factors \vec{\theta}, \vec{\phi} (as lists) and \lambda, in this order."""
    n = self.degree()

    alpha = - self.lbd - sum(self.theta) - (n+1) * np.pi # needed to counterbalance the global phase

    mw_theta = [self.phi[k] for k in range(n+1)]
    mw_phi = [2 * self.lbd - np.pi/2 + alpha] + [2 * self.theta[k-1] - np.pi for k in range(1, n+1)]
    mw_lbd = 2 * self.theta[n] - np.pi/2 + alpha

    return mw_theta, mw_phi, mw_lbd

to_xqsp()

Converts the QSP phase factors into XQSP phase factors.

Raises:

Type	Description
ValueError	If the phase factors do not lie in the XQSP subalgebra.

Source code in nlft_qsp/qsp.py

def to_xqsp(self):
    """Converts the QSP phase factors into XQSP phase factors.

    Raises:
        ValueError: If the phase factors do not lie in the XQSP subalgebra."""
    F = self.to_nlfs()
    if not F.is_imaginary():
        raise ValueError("The phase factors are not reducible to XQSP.")

    return XQSPPhaseFactors.from_nlfs(F)

to_yqsp()

Converts the QSP phase factors into YQSP phase factors.

Raises:

Type	Description
ValueError	If the phase factors do not lie in the YQSP subalgebra.

Source code in nlft_qsp/qsp.py

def to_yqsp(self):
    """Converts the QSP phase factors into YQSP phase factors.

    Raises:
        ValueError: If the phase factors do not lie in the YQSP subalgebra."""
    F = self.to_nlfs()
    if not F.is_real():
        raise ValueError("The phase factors are not reducible to YQSP.")

    return YQSPPhaseFactors.from_nlfs(F)

NonLinearFourierSequence

Bases: ComplexL0Sequence

Class representing a finitely supported sequence of complex numbers over \(\mathbb{Z}\). The class provides methods to compute the nonlinear Fourier transform (NLFT) associated with the sequence.

Methods:

Name	Description
__init__	Initializes a nonlinear Fourier sequence with a given list of complex values and support starting index.
transform	Computes the nonlinear Fourier transform \((a(z), b(z))\) over \(SU(2)\) associated with this sequence.
transform_bounds	Computes the nonlinear Fourier transform over \(SU(2)\) for the subsequence within the specified range.

Source code in nlft_qsp/nlft.py

class NonLinearFourierSequence(ComplexL0Sequence):
    r"""Class representing a finitely supported sequence of complex numbers over $\mathbb{Z}$.
    The class provides methods to compute the nonlinear Fourier transform (NLFT) associated with the sequence.
    """

    def __init__(self, coeffs: list[complex_type] | np.ndarray=[], support_start: int=0):
        r"""Initializes a nonlinear Fourier sequence with a given list of complex values and support starting index.

        Args:
            coeffs (list[complex_type] | np.ndarray): A list of complex numbers representing the sequence. The list includes both the lower 
                                      and upper bounds of the sequence.
            support_start (int): The index of the first element of the sequence in $\mathbb{Z}$. The support of the sequence will 
                                 be in the range [`support_start`, `support_start + coeffs.shape[0]`].
        """
        super().__init__(coeffs, support_start)

    def transform_bounds(self, inf, sup) -> tuple[Polynomial, Polynomial]:
        """
        Computes the nonlinear Fourier transform over $SU(2)$ for the subsequence within the specified range.

        Args:
            inf (int): The lower bound (included) index of the sequence for the transformation.
            sup (int): The upper bound (excluded) index of the sequence for the transformation.

        Returns:
            The $SU(2)$-NLFT of the subsequence in [`inf`, `sup`).

        Note:
            This is used only internally in order to compute the polynomials through a divide-and-conquer strategy.
            If only interested in the final polynomials, please refer to `transform()`.
        """
        if sup - inf <= 0:
            return Polynomial([1]), Polynomial([0])

        if sup - inf <= 1:
            F = self[inf]
            den = np.sqrt(1 + F * np.conj(F))
            return Polynomial([1/den]), Polynomial([F/den], inf)  # (1/den, F/den z^inf)

        mid = (sup + inf) // 2
        a1, b1 = self.transform_bounds(inf, mid)
        a2, b2 = self.transform_bounds(mid, sup)

        return a1 * a2 - b1 * b2.conjugate(), a1 * b2 + b1 * a2.conjugate()

    def transform(self) -> tuple[Polynomial, Polynomial]:
        r"""Computes the nonlinear Fourier transform $(a(z), b(z))$ over $SU(2)$ associated with this sequence.

        $$\mathrm{NLFT}(F) = (a, b)$$

        See [here](https://arxiv.org/abs/2503.03026) for a definition of the nonlinear Fourier transform.

        Returns:
            The $SU(2)$-NLFT of the sequence.
        """
        n = self.coeffs.shape[0] - 1
        a, b = self.transform_bounds(self.support_start, self.support_start + n + 1)

        if n < 0:
            return a, b

        return a.truncate(-n, 0), b.truncate(self.support_start, self.support_start + n)

__init__(coeffs: list[complex_type] | np.ndarray = [], support_start: int = 0)

Initializes a nonlinear Fourier sequence with a given list of complex values and support starting index.

Parameters:

Name	Type	Description	Default
coeffs	list[complex_type] | ndarray	A list of complex numbers representing the sequence. The list includes both the lower and upper bounds of the sequence.	[]
support_start	int	The index of the first element of the sequence in \(\mathbb{Z}\). The support of the sequence will be in the range [support_start, support_start + coeffs.shape[0]].	0

Source code in nlft_qsp/nlft.py

def __init__(self, coeffs: list[complex_type] | np.ndarray=[], support_start: int=0):
    r"""Initializes a nonlinear Fourier sequence with a given list of complex values and support starting index.

    Args:
        coeffs (list[complex_type] | np.ndarray): A list of complex numbers representing the sequence. The list includes both the lower 
                                  and upper bounds of the sequence.
        support_start (int): The index of the first element of the sequence in $\mathbb{Z}$. The support of the sequence will 
                             be in the range [`support_start`, `support_start + coeffs.shape[0]`].
    """
    super().__init__(coeffs, support_start)

transform() -> tuple[Polynomial, Polynomial]

Computes the nonlinear Fourier transform \((a(z), b(z))\) over \(SU(2)\) associated with this sequence.

\[\mathrm{NLFT}(F) = (a, b)\]

See here for a definition of the nonlinear Fourier transform.

Returns:

Type	Description
tuple[Polynomial, Polynomial]	The \(SU(2)\)-NLFT of the sequence.

Source code in nlft_qsp/nlft.py

def transform(self) -> tuple[Polynomial, Polynomial]:
    r"""Computes the nonlinear Fourier transform $(a(z), b(z))$ over $SU(2)$ associated with this sequence.

    $$\mathrm{NLFT}(F) = (a, b)$$

    See [here](https://arxiv.org/abs/2503.03026) for a definition of the nonlinear Fourier transform.

    Returns:
        The $SU(2)$-NLFT of the sequence.
    """
    n = self.coeffs.shape[0] - 1
    a, b = self.transform_bounds(self.support_start, self.support_start + n + 1)

    if n < 0:
        return a, b

    return a.truncate(-n, 0), b.truncate(self.support_start, self.support_start + n)

transform_bounds(inf, sup) -> tuple[Polynomial, Polynomial]

Computes the nonlinear Fourier transform over \(SU(2)\) for the subsequence within the specified range.

Parameters:

Name	Type	Description	Default
inf	int	The lower bound (included) index of the sequence for the transformation.	required
sup	int	The upper bound (excluded) index of the sequence for the transformation.	required

Returns:

Type	Description
tuple[Polynomial, Polynomial]	The \(SU(2)\)-NLFT of the subsequence in [inf, sup).

Note

This is used only internally in order to compute the polynomials through a divide-and-conquer strategy. If only interested in the final polynomials, please refer to transform().

Source code in nlft_qsp/nlft.py

def transform_bounds(self, inf, sup) -> tuple[Polynomial, Polynomial]:
    """
    Computes the nonlinear Fourier transform over $SU(2)$ for the subsequence within the specified range.

    Args:
        inf (int): The lower bound (included) index of the sequence for the transformation.
        sup (int): The upper bound (excluded) index of the sequence for the transformation.

    Returns:
        The $SU(2)$-NLFT of the subsequence in [`inf`, `sup`).

    Note:
        This is used only internally in order to compute the polynomials through a divide-and-conquer strategy.
        If only interested in the final polynomials, please refer to `transform()`.
    """
    if sup - inf <= 0:
        return Polynomial([1]), Polynomial([0])

    if sup - inf <= 1:
        F = self[inf]
        den = np.sqrt(1 + F * np.conj(F))
        return Polynomial([1/den]), Polynomial([F/den], inf)  # (1/den, F/den z^inf)

    mid = (sup + inf) // 2
    a1, b1 = self.transform_bounds(inf, mid)
    a2, b2 = self.transform_bounds(mid, sup)

    return a1 * a2 - b1 * b2.conjugate(), a1 * b2 + b1 * a2.conjugate()

PhaseFactors

Set of phase factors for a general Quantum Signal Processing protocol. It also provides methods to construct polynomials generated by QSP protocols. Each subclass of this class represents a different QSP ansatz.

Methods:

Name	Description
degree	Returns the degree of the polynomials generated by the QSP protocol.
iX	Returns a new QSP protocol, obtained by multiplying the given QSP protocol by \(iX\) on the right, where \(X\) is the Pauli matrix.
iY	Returns a new QSP protocol, obtained by multiplying the given QSP protocol by \(iY\) on the right, where \(Y\) is the Pauli matrix.
iZ	Returns a new QSP protocol, obtained by multiplying the given QSP protocol by \(iZ\) on the right, where \(Z\) is the Pauli matrix.
polynomials	Returns the pair of polynomials \((P, Q)\) generated by the given set of phase factors.
polynomials_bounds	Returns the pair of polynomials \((P, Q) = A_{inf} \tilde{v} A_{inf+1} \tilde{v} ... \tilde{v} A_{sup-1} \tilde{v} A_{sup}\),
processing_operator	Returns the \(k\)-th signal processing operator according to the given QSP variant.
processing_operator_conjugation	This applies a conjugation to each signal processing operator returned by processing_operator(), e.g., with the Hadamard gate.
protocol_conjugation	Given \((P, Q)\), this applies a final conjugation to the whole protocol, e.g., with the Hadamard gate.
signal_operator	Returns the pair of polynomials given by multiplying \((P_1, Q_1) W(z) (P_2, Q_2)\), where W(z) is the signal operator.
to_nlfs	Returns the nonlinear Fourier sequence generating \((z^{-n} P, Q)\),

Source code in nlft_qsp/qsp.py

class PhaseFactors:
    """Set of phase factors for a general Quantum Signal Processing protocol.
    It also provides methods to construct polynomials generated by QSP protocols.
    Each subclass of this class represents a different QSP ansatz."""
    def duplicate(self):
        raise NotImplementedError()

    def processing_operator(self, k: int):
        r"""Returns the $k$-th signal processing operator according to the given QSP variant."""
        raise NotImplementedError()

    def signal_operator(self, P1: Polynomial, Q1: Polynomial, P2: Polynomial, Q2: Polynomial) -> tuple[Polynomial, Polynomial]:
        """Returns the pair of polynomials given by multiplying $(P_1, Q_1) W(z) (P_2, Q_2)$, where `W(z)` is the signal operator.

        Note:
            This might not reflect the signal operator as expressed in the original papers. Some basis transformations are implicitly made for computational efficiency."""
        # W(z) = diag(z, z^(-1))
        zP1 = P1.shift(1)  # z*P1
        zQ1 = Q1.shift(-1) # z^(-1)*Q1
        return zP1 * P2 - zQ1 * Q2.conjugate(), zP1 * Q2 + zQ1 * P2.conjugate()

    def protocol_conjugation(self, P, Q):
        r"""Given $(P, Q)$, this applies a final conjugation to the whole protocol, e.g., with the Hadamard gate."""
        return P, Q

    def processing_operator_conjugation(self, a, b):
        """This applies a conjugation to each signal processing operator returned by `processing_operator()`, e.g., with the Hadamard gate.
        This is done mainly to make `processing_operator()` return the operators as expressed by the original ansatze in the papers, while keeping computational efficiency and numerical stability."""
        return a, b

    def degree(self):
        """Returns the degree of the polynomials generated by the QSP protocol."""
        raise NotImplementedError()

    def iX(self):
        r"""Returns a new QSP protocol, obtained by multiplying the given QSP protocol by $iX$ on the right, where $X$ is the Pauli matrix.
        This will make the generated polynomials undergo the transformation $(P, Q) \rightarrow (iQ, iP)$.
        This method is useful to bring to swap the places of the two polynomials, to switch between NLFT and QSP conventions.

        Raises:
            ValueError: if multiplying by $iX$ does not preserve the subalgebra of the phase factors."""
        raise ValueError("Multiplying by iX is not possible.")

    def iY(self):
        r"""Returns a new QSP protocol, obtained by multiplying the given QSP protocol by $iY$ on the right, where $Y$ is the Pauli matrix.
        This will make the generated polynomials undergo the transformation $(P, Q) \rightarrow (-Q, P)$.
        This method is useful to bring to swap the places of the two polynomials, to switch between NLFT and QSP conventions.

        Raises:
            ValueError: if multiplying by $iY$ does not preserve the subalgebra of the phase factors."""
        raise ValueError("Multiplying by iY is not possible.")

    def iZ(self):
        r"""Returns a new QSP protocol, obtained by multiplying the given QSP protocol by $iZ$ on the right, where $Z$ is the Pauli matrix.
        This will make the generated polynomials undergo the transformation $(P, Q) \rightarrow (iP, -iQ)$.
        This method is useful to bring to swap the places of the two polynomials, to switch between NLFT and QSP conventions.

        Raises:
            ValueError: if multiplying by $iZ$ does not preserve the subalgebra of the phase factors."""
        raise ValueError("Multiplying by iZ is not possible.")

    def polynomials_bounds(self, inf: int, sup: int) -> tuple[Polynomial, Polynomial]:
        r"""Returns the pair of polynomials $(P, Q) = A_{inf} \tilde{v} A_{inf+1} \tilde{v} ... \tilde{v} A_{sup-1} \tilde{v} A_{sup}$,
        where $\tilde{v} = \mathrm{diag}(z, z^{-1})$ is the signal operator.

        Note:
            This assumes the Laurent picture ($\tilde{v} = \mathrm{diag}(z, z^{-1})$). For the analytic picture $\tilde{w} = \mathrm{diag}(z, 1)$ use `polynomials()`."""
        if sup - inf < 0:
            return Polynomial([1]), Polynomial([0])
        if sup - inf <= 0:
            p, q = self.processing_operator(inf)
            p, q = self.processing_operator_conjugation(p, q)
            return Polynomial([p]), Polynomial([q])

        mid = (sup + inf) // 2
        P1, Q1 = self.polynomials_bounds(inf, mid)
        P2, Q2 = self.polynomials_bounds(mid+1, sup)

        return self.signal_operator(P1, Q1, P2, Q2) # (P1, Q1) W(z) (P2, Q2)

    def polynomials(self, inf: int = 0, sup: int = -1, mode: str = 'analytic') -> tuple[Polynomial, Polynomial]:
        """Returns the pair of polynomials $(P, Q)$ generated by the given set of phase factors.
        The polynomials are computed with a divide-and-conquer strategy (see [here](https://arxiv.org/abs/2410.06409)).

        Args:
            mode (str): Either `'analytic'` or `'laurent'`, indicating whether an analytic or a Laurent QSP protocol should be composed.
        """
        if sup < 0:
            sup = self.degree()

        Pl, Ql = self.polynomials_bounds(inf, sup)
        Pl, Ql = self.protocol_conjugation(Pl, Ql)
        match mode:
            case 'analytic': # convert from Laurent to analytic picture
                return laurent_to_analytic(Pl), laurent_to_analytic(Ql)
            case 'laurent':
                return Pl, Ql
            case _:
                raise ValueError("mode can only be 'analytic' or 'laurent'.")

    def to_nlfs(self) -> NonLinearFourierSequence:
        """Returns the nonlinear Fourier sequence generating $(z^{-n} P, Q)$,
        where $(P, Q)$ is the pair of polynomial generated by the given set of GQSP phase factors.

        Note: if the phase factors are not canonical, then the phase of the leading coefficient of $P$ is adjusted
        so that it becomes real and positive, and $(z^{-n} P, Q)$ is in the image of the NLFT."""
        raise NotImplementedError()

degree()

Returns the degree of the polynomials generated by the QSP protocol.

Source code in nlft_qsp/qsp.py

def degree(self):
    """Returns the degree of the polynomials generated by the QSP protocol."""
    raise NotImplementedError()

iX()

Returns a new QSP protocol, obtained by multiplying the given QSP protocol by \(iX\) on the right, where \(X\) is the Pauli matrix. This will make the generated polynomials undergo the transformation \((P, Q) \rightarrow (iQ, iP)\). This method is useful to bring to swap the places of the two polynomials, to switch between NLFT and QSP conventions.

Raises:

Type	Description
ValueError	if multiplying by \(iX\) does not preserve the subalgebra of the phase factors.

Source code in nlft_qsp/qsp.py

def iX(self):
    r"""Returns a new QSP protocol, obtained by multiplying the given QSP protocol by $iX$ on the right, where $X$ is the Pauli matrix.
    This will make the generated polynomials undergo the transformation $(P, Q) \rightarrow (iQ, iP)$.
    This method is useful to bring to swap the places of the two polynomials, to switch between NLFT and QSP conventions.

    Raises:
        ValueError: if multiplying by $iX$ does not preserve the subalgebra of the phase factors."""
    raise ValueError("Multiplying by iX is not possible.")

iY()

Returns a new QSP protocol, obtained by multiplying the given QSP protocol by \(iY\) on the right, where \(Y\) is the Pauli matrix. This will make the generated polynomials undergo the transformation \((P, Q) \rightarrow (-Q, P)\). This method is useful to bring to swap the places of the two polynomials, to switch between NLFT and QSP conventions.

Raises:

Type	Description
ValueError	if multiplying by \(iY\) does not preserve the subalgebra of the phase factors.

Source code in nlft_qsp/qsp.py

def iY(self):
    r"""Returns a new QSP protocol, obtained by multiplying the given QSP protocol by $iY$ on the right, where $Y$ is the Pauli matrix.
    This will make the generated polynomials undergo the transformation $(P, Q) \rightarrow (-Q, P)$.
    This method is useful to bring to swap the places of the two polynomials, to switch between NLFT and QSP conventions.

    Raises:
        ValueError: if multiplying by $iY$ does not preserve the subalgebra of the phase factors."""
    raise ValueError("Multiplying by iY is not possible.")

iZ()

Returns a new QSP protocol, obtained by multiplying the given QSP protocol by \(iZ\) on the right, where \(Z\) is the Pauli matrix. This will make the generated polynomials undergo the transformation \((P, Q) \rightarrow (iP, -iQ)\). This method is useful to bring to swap the places of the two polynomials, to switch between NLFT and QSP conventions.

Raises:

Type	Description
ValueError	if multiplying by \(iZ\) does not preserve the subalgebra of the phase factors.

Source code in nlft_qsp/qsp.py

def iZ(self):
    r"""Returns a new QSP protocol, obtained by multiplying the given QSP protocol by $iZ$ on the right, where $Z$ is the Pauli matrix.
    This will make the generated polynomials undergo the transformation $(P, Q) \rightarrow (iP, -iQ)$.
    This method is useful to bring to swap the places of the two polynomials, to switch between NLFT and QSP conventions.

    Raises:
        ValueError: if multiplying by $iZ$ does not preserve the subalgebra of the phase factors."""
    raise ValueError("Multiplying by iZ is not possible.")

polynomials(inf: int = 0, sup: int = -1, mode: str = 'analytic') -> tuple[Polynomial, Polynomial]

Returns the pair of polynomials \((P, Q)\) generated by the given set of phase factors. The polynomials are computed with a divide-and-conquer strategy (see here).

Parameters:

Name	Type	Description	Default
mode	str	Either 'analytic' or 'laurent', indicating whether an analytic or a Laurent QSP protocol should be composed.	'analytic'

Source code in nlft_qsp/qsp.py

def polynomials(self, inf: int = 0, sup: int = -1, mode: str = 'analytic') -> tuple[Polynomial, Polynomial]:
    """Returns the pair of polynomials $(P, Q)$ generated by the given set of phase factors.
    The polynomials are computed with a divide-and-conquer strategy (see [here](https://arxiv.org/abs/2410.06409)).

    Args:
        mode (str): Either `'analytic'` or `'laurent'`, indicating whether an analytic or a Laurent QSP protocol should be composed.
    """
    if sup < 0:
        sup = self.degree()

    Pl, Ql = self.polynomials_bounds(inf, sup)
    Pl, Ql = self.protocol_conjugation(Pl, Ql)
    match mode:
        case 'analytic': # convert from Laurent to analytic picture
            return laurent_to_analytic(Pl), laurent_to_analytic(Ql)
        case 'laurent':
            return Pl, Ql
        case _:
            raise ValueError("mode can only be 'analytic' or 'laurent'.")

polynomials_bounds(inf: int, sup: int) -> tuple[Polynomial, Polynomial]

Returns the pair of polynomials \((P, Q) = A_{inf} \tilde{v} A_{inf+1} \tilde{v} ... \tilde{v} A_{sup-1} \tilde{v} A_{sup}\), where \(\tilde{v} = \mathrm{diag}(z, z^{-1})\) is the signal operator.

Note

This assumes the Laurent picture (\(\tilde{v} = \mathrm{diag}(z, z^{-1})\)). For the analytic picture \(\tilde{w} = \mathrm{diag}(z, 1)\) use polynomials().

Source code in nlft_qsp/qsp.py

def polynomials_bounds(self, inf: int, sup: int) -> tuple[Polynomial, Polynomial]:
    r"""Returns the pair of polynomials $(P, Q) = A_{inf} \tilde{v} A_{inf+1} \tilde{v} ... \tilde{v} A_{sup-1} \tilde{v} A_{sup}$,
    where $\tilde{v} = \mathrm{diag}(z, z^{-1})$ is the signal operator.

    Note:
        This assumes the Laurent picture ($\tilde{v} = \mathrm{diag}(z, z^{-1})$). For the analytic picture $\tilde{w} = \mathrm{diag}(z, 1)$ use `polynomials()`."""
    if sup - inf < 0:
        return Polynomial([1]), Polynomial([0])
    if sup - inf <= 0:
        p, q = self.processing_operator(inf)
        p, q = self.processing_operator_conjugation(p, q)
        return Polynomial([p]), Polynomial([q])

    mid = (sup + inf) // 2
    P1, Q1 = self.polynomials_bounds(inf, mid)
    P2, Q2 = self.polynomials_bounds(mid+1, sup)

    return self.signal_operator(P1, Q1, P2, Q2) # (P1, Q1) W(z) (P2, Q2)

processing_operator(k: int)

Returns the \(k\)-th signal processing operator according to the given QSP variant.

Source code in nlft_qsp/qsp.py

def processing_operator(self, k: int):
    r"""Returns the $k$-th signal processing operator according to the given QSP variant."""
    raise NotImplementedError()

processing_operator_conjugation(a, b)

This applies a conjugation to each signal processing operator returned by processing_operator(), e.g., with the Hadamard gate. This is done mainly to make processing_operator() return the operators as expressed by the original ansatze in the papers, while keeping computational efficiency and numerical stability.

Source code in nlft_qsp/qsp.py

def processing_operator_conjugation(self, a, b):
    """This applies a conjugation to each signal processing operator returned by `processing_operator()`, e.g., with the Hadamard gate.
    This is done mainly to make `processing_operator()` return the operators as expressed by the original ansatze in the papers, while keeping computational efficiency and numerical stability."""
    return a, b

protocol_conjugation(P, Q)

Given \((P, Q)\), this applies a final conjugation to the whole protocol, e.g., with the Hadamard gate.

Source code in nlft_qsp/qsp.py

def protocol_conjugation(self, P, Q):
    r"""Given $(P, Q)$, this applies a final conjugation to the whole protocol, e.g., with the Hadamard gate."""
    return P, Q

signal_operator(P1: Polynomial, Q1: Polynomial, P2: Polynomial, Q2: Polynomial) -> tuple[Polynomial, Polynomial]

Returns the pair of polynomials given by multiplying \((P_1, Q_1) W(z) (P_2, Q_2)\), where W(z) is the signal operator.

Note

This might not reflect the signal operator as expressed in the original papers. Some basis transformations are implicitly made for computational efficiency.

Source code in nlft_qsp/qsp.py

def signal_operator(self, P1: Polynomial, Q1: Polynomial, P2: Polynomial, Q2: Polynomial) -> tuple[Polynomial, Polynomial]:
    """Returns the pair of polynomials given by multiplying $(P_1, Q_1) W(z) (P_2, Q_2)$, where `W(z)` is the signal operator.

    Note:
        This might not reflect the signal operator as expressed in the original papers. Some basis transformations are implicitly made for computational efficiency."""
    # W(z) = diag(z, z^(-1))
    zP1 = P1.shift(1)  # z*P1
    zQ1 = Q1.shift(-1) # z^(-1)*Q1
    return zP1 * P2 - zQ1 * Q2.conjugate(), zP1 * Q2 + zQ1 * P2.conjugate()

to_nlfs() -> NonLinearFourierSequence

Returns the nonlinear Fourier sequence generating \((z^{-n} P, Q)\), where \((P, Q)\) is the pair of polynomial generated by the given set of GQSP phase factors.

Note: if the phase factors are not canonical, then the phase of the leading coefficient of \(P\) is adjusted so that it becomes real and positive, and \((z^{-n} P, Q)\) is in the image of the NLFT.

Source code in nlft_qsp/qsp.py

def to_nlfs(self) -> NonLinearFourierSequence:
    """Returns the nonlinear Fourier sequence generating $(z^{-n} P, Q)$,
    where $(P, Q)$ is the pair of polynomial generated by the given set of GQSP phase factors.

    Note: if the phase factors are not canonical, then the phase of the leading coefficient of $P$ is adjusted
    so that it becomes real and positive, and $(z^{-n} P, Q)$ is in the image of the NLFT."""
    raise NotImplementedError()

Polynomial

Bases: ComplexL0Sequence

Represents a general Laurent polynomial of one complex variable.

Attributes:

Name	Type	Description
coeffs	list[complex_type]	List of complex coefficients.
shape	tuple	The shape of the coefficients. This is always () for scalar polynomials, but is used for matrix polynomials.
support_start	int	Minimum degree that appears in the polynomial.

Methods:

Name	Description
__call__	Evaluates the polynomial using Horner's method.
__init__	Initializes a Polynomial instance.
__str__	Converts the polynomial to a human-readable string representation.
analytic_part	Discards all the negative degrees, keeping only the non-negative ones.
anti_analytic_part	Discards all the positive degrees, keeping only the non-positive ones.
block_matrix	Returns a matrix polynomial containing the given polynomials/constants as blocks.
conjugate	Returns the conjugate polynomial on the unit circle. If \(p(z) = \sum_k p_k z^k\), then its conjugate is defined as \(p^*(z) = \sum_k p_k^* z^{-k}\). This is also known as the Schwarz reflection across \(\mathbb{T}\).
diagonal_block_matrix	Returns a matrix polynomial containing the given polynomials/constants as blocks along the diagonal.
duplicate	Creates a duplicate of the current polynomial.
effective_degree	Returns the size of the support of the polynomial minus 1 (max degree - min degree).
eval_at_roots_of_unity	Evaluates the polynomial at the \(N\)-th roots of unity using the inverse fast Fourier transform.
hilbert_transform	Returns the polynomial \(Q\) such that \(P + Q\) yields an analytic polynomial (\(P\) being self).
only_negative_degrees	DEPRECATED: use anti_analytic_part() instead.
only_positive_degrees	DEPRECATED: use analytic_part() instead.
schwarz_transform	Returns the anti-analytic polynomial whose real part gives the current polynomial.
sharp	Same as conjugate(), but support_start is left unchanged.
shift	Creates a new polynomial equal to the current one, multiplied by z^k.
sup_norm	Estimates the supremum norm of the polynomial over the unit circle
truncate	Keeps only the coefficients in \([m, n]\), discarding the others.

Source code in nlft_qsp/poly.py

@serializable(
    type_tag="@qspx/polynomial",
    fields={
        "support_start": "support_start",
        "coeffs": "coeffs"
    }
)
class Polynomial(ComplexL0Sequence):
    """Represents a general Laurent polynomial of one complex variable.

    Attributes:
        coeffs (list[complex_type]): List of complex coefficients.
        shape (tuple): The shape of the coefficients. This is always `()` for scalar polynomials, but is used for matrix polynomials.
        support_start (int): Minimum degree that appears in the polynomial.
    """

    def __init__(self, coeffs: list[complex_type] | np.ndarray = None, support_start: int = 0, shape: tuple[int] = None):
        """Initializes a Polynomial instance.

        Args:
            coeffs: List of complex numbers as coefficients.
            shape (optional): Shape of the coefficient array.
            support_start (optional): Minimum degree in the polynomial. Defaults to 0.

        Note:
            Please provide exactly one of coeffs or shape.
        """
        super().__init__(coeffs, support_start, shape)

    def duplicate(self) -> "Polynomial":
        """Creates a duplicate of the current polynomial.

        Returns:
            Polynomial: A new Polynomial instance with the same coefficients and support.
        """
        return type(self)(self.coeffs, self.support_start)

    def shift(self, k: int):
        """Creates a new polynomial equal to the current one, multiplied by `z^k`."""
        return type(self)(self.coeffs, self.support_start + k)

    def effective_degree(self) -> int:
        """Returns the size of the support of the polynomial minus 1 (max degree - min degree).

        Note:
            This does not check for leading or trailing zeros in the coefficient array.

        Returns:
            int: The effective degree of the polynomial.
        """
        return self.coeffs.shape[0] - 1

    def conjugate(self) -> "Polynomial":
        r"""Returns the conjugate polynomial on the unit circle. If $p(z) = \sum_k p_k z^k$, then its conjugate is defined as $p^*(z) = \sum_k p_k^* z^{-k}$. This is also known as the Schwarz reflection across $\mathbb{T}$.

        Note:
            For a matrix polynomial, each coefficient is conjugate-transposed.
        """
        if self.shape == ():
            conj_coeffs = np.flip(np.conj(self.coeffs), axis=0)
        else:
            axes = [0] + list(range(1, self.coeffs.ndim))[::-1] # transpose = reverse the order of the axes except the first one, which is the index of the coefficient
            conj_coeffs = np.flip(np.transpose(np.conj(self.coeffs), axes=axes), axis=0)

        return Polynomial(conj_coeffs, -(self.support_start + self.coeffs.shape[0] - 1))

    def sharp(self) -> "Polynomial":
        r"""Same as `conjugate()`, but `support_start` is left unchanged.

        Returns:
            Polynomial: The sharp-conjugate polynomial.
        """
        p = self.conjugate()
        p.support_start += self.effective_degree() + 1
        return p

    def schwarz_transform(self) -> "Polynomial":
        r"""Returns the anti-analytic polynomial whose real part gives the current polynomial.

        Note:
            This is equivalent to adding $i \mathcal{H}[p]`, where $\mathcal{H}[p]$ is the Hilbert transform of $p$.
        """
        schwarz_coeffs = []
        for k in self.support():
            if k < 0:
                schwarz_coeffs.append(2*self[k])
            elif k == 0:
                schwarz_coeffs.append(self[k])

        return Polynomial(np.array(schwarz_coeffs, dtype=complex_type), self.support_start)

    def hilbert_transform(self) -> "Polynomial":
        r"""Returns the polynomial $Q$ such that $P + Q$ yields an analytic polynomial ($P$ being `self`).

        Note:
            This is actually $i \mathcal{H}[P]$, i.e., the Hilbert transform as returned is already multiplied by $i$.
        """
        hilbert_coeffs = []
        for k in self.support():
            if k < 0:
                hilbert_coeffs.append(-self[k])
            elif k > 0:
                hilbert_coeffs.append(self[k])
            else:
                hilbert_coeffs.append(0)

        return Polynomial(np.array(hilbert_coeffs, dtype=complex_type), self.support_start)

    def __mul__(self, other):
        if isinstance(other, Number):
            return Polynomial(self.coeffs * other, self.support_start)
        elif isinstance(other, np.ndarray):
            if other.shape != self.shape:
                raise ValueError(f"The shape of the given matrix {other.shape} does not match the expected shape {self.shape}.")

            return Polynomial(self.coeffs * other, self.support_start)
        elif not isinstance(other, Polynomial):
            raise TypeError("Polynomial multiplication admits only other polynomials or constants with compatible shape.")

        # Pad so they end up with the same length
        target_len = self.coeffs.shape[0] + other.coeffs.shape[0] - 1
        pad_a = [(0, target_len - self.coeffs.shape[0])] + [(0, 0)] * (self.coeffs.ndim - 1)
        pad_b = [(0, target_len - other.coeffs.shape[0])] + [(0, 0)] * (other.coeffs.ndim - 1)

        coeffs_a = np.fft.fft(np.pad(self.coeffs, pad_width=pad_a), axis=0)
        coeffs_b = np.fft.fft(np.pad(other.coeffs, pad_width=pad_b), axis=0)

        # Multiply in the Fourier domain
        coeffs_c = [a * b for a, b in zip(coeffs_a, coeffs_b)]

        # Inverse FFT to get the result
        new_coeffs = np.fft.ifft(coeffs_c, axis=0)
        support_start = self.support_start + other.support_start  # Lowest degree of the new poly

        return Polynomial(new_coeffs, support_start)

    def __rmul__(self, other):
        return self * other

    def __truediv__(self, other):
        if isinstance(other, Number):
            return Polynomial(self.coeffs / other, self.support_start)

        raise TypeError("Polynomial division is only possible with scalars.")

    def __matmul__(self, other):
        if isinstance(other, np.ndarray):
            if self.coeffs.shape[-1] != other.shape[-2]:
                raise ValueError(f"Incompatible shapes: {self.shape} vs {other.shape}.")

            return Polynomial(self.coeffs @ other, self.support_start)
        elif not isinstance(other, Polynomial):
            raise TypeError("Polynomial multiplication admits only other polynomials or constants with compatible shape.")

        if self.coeffs.shape[-1] != other.shape[-2]:
            raise ValueError(f"Incompatible shapes: {self.shape} vs {other.shape}.")

        # Pad so they end up with the same length
        target_len = self.coeffs.shape[0] + other.coeffs.shape[0] - 1
        pad_a = [(0, target_len - self.coeffs.shape[0])] + [(0, 0)] * (self.coeffs.ndim - 1)
        pad_b = [(0, target_len - other.coeffs.shape[0])] + [(0, 0)] * (other.coeffs.ndim - 1)

        coeffs_a = np.fft.fft(np.pad(self.coeffs, pad_width=pad_a), axis=0)
        coeffs_b = np.fft.fft(np.pad(other.coeffs, pad_width=pad_b), axis=0)

        # Multiply in the Fourier domain
        coeffs_c = [a @ b for a, b in zip(coeffs_a, coeffs_b)]

        # Inverse FFT to get the result
        new_coeffs = np.fft.ifft(coeffs_c, axis=0)
        support_start = self.support_start + other.support_start  # Lowest degree of the new poly

        return Polynomial(new_coeffs, support_start)

    def __rmatmul__(self, other):
        if isinstance(other, np.ndarray):
            if other.shape[-1] != self.coeffs.shape[-2]:
                raise ValueError(f"The shape of the given matrix {other.shape} does not match the expected shape {self.shape}.")

            return Polynomial(other @ self.coeffs, self.support_start)

        raise TypeError("Polynomial multiplication admits only other polynomials or constants with compatible shape.")

    def __call__(self, z) -> complex_type | np.ndarray:
        """Evaluates the polynomial using Horner's method.

        Args:
            z (complex): The point at which to evaluate the polynomial.

        Returns:
            complex: The evaluated result.
        """
        res = self.coeffs[-1]
        for k in reversed(range(self.coeffs.shape[0] - 1)):
            res = res * z + self.coeffs[k]
        return res * (z ** self.support_start)

    def eval_at_roots_of_unity(self, N: int) -> list[complex_type]: # TODO remove power-of-two assumption
        r"""Evaluates the polynomial at the $N$-th roots of unity using the inverse fast Fourier transform.

        Args:
            N (int): A power of two specifying the number of roots. If $N$ is not a power of two, then the next power of two is implicitly taken.

        Returns:
            list[complex]: List of evaluations at the N-th roots of unity.
            The $k$-th element of the list will be $P(e^{2\pi i k/N})$, $P$ being `self`.
        """
        N = next_power_of_two(N)
        M = next_power_of_two(max(N, self.coeffs.shape[0]))

        pad_width = [(0, M - self.coeffs.shape[0])] + [(0, 0)] * (self.coeffs.ndim - 1)
        coeffs = np.pad(self.coeffs, pad_width=pad_width)
        coeffs = np.roll(coeffs, self.support_start, axis=0)
        # This has the effect of having everything multiplied by z^s

        evals = np.fft.ifft(coeffs, norm='forward', axis=0) # M evaluations at the M-th roots of unity
        return evals[::M//N]

    def sup_norm(self, N=1024) -> float_type:
        """Estimates the supremum norm of the polynomial over the unit circle

        Args:
            N (int): the number of samples to compute the maximum from. If $N$ is not a power of two, then the next power of two is taken.

        Returns:
            float_type: An estimate for the supremum norm of the polynomial over the unit circle.
        """
        if self.shape == ():
            return np.max([np.abs(sample) for sample in self.eval_at_roots_of_unity(N)])

        return np.max([np.linalg.norm(sample, ord=2) for sample in self.eval_at_roots_of_unity(N)])

    def truncate(self, m: int, n: int) -> "Polynomial":
        """Keeps only the coefficients in $[m, n]$, discarding the others.

        Returns:
            Polynomial: A new, truncated polynomial.
        """
        if self.shape == ():
            return Polynomial([self[k] for k in range(m, n+1)], m)

        return Polynomial(np.array([self[k] for k in range(m, n+1)], dtype=complex_type), m)

    def only_positive_degrees(self) -> "Polynomial":
        """DEPRECATED: use `analytic_part()` instead."""
        return self.analytic_part()

    def analytic_part(self) -> "Polynomial":
        """Discards all the negative degrees, keeping only the non-negative ones.

        Returns:
            Polynomial: A new polynomial containing only the positive-degree coefficients."""
        return self.truncate(0, self.support_start + self.coeffs.shape[0] - 1)

    def only_negative_degrees(self) -> "Polynomial":
        """DEPRECATED: use `anti_analytic_part()` instead."""
        return self.anti_analytic_part()

    def anti_analytic_part(self) -> "Polynomial":
        """Discards all the positive degrees, keeping only the non-positive ones.

        Returns:
            Polynomial: A new polynomial containing only the non-positive-degree coefficients."""
        return self.truncate(self.support_start, 0)

    def __str__(self) -> str:
        """Converts the polynomial to a human-readable string representation.

        Returns:
            str: The string representation of the polynomial.
        """
        return ' + '.join(f"{c} z^{self.support_start + k}" for k, c in enumerate(self.coeffs))

    @classmethod
    def block_matrix(cls, blocks: list[list["Polynomial"]]):
        """Returns a matrix polynomial containing the given polynomials/constants as blocks.

        Raises: ValueError if the block shapes do not match."""
        joint_support = range(min(min(p.support().start for p in l) for l in blocks),
                              max(max(p.support().stop for p in l) for l in blocks))

        coeffs = []
        for k in joint_support:
            coeffs.append(np.block([[p[k] for p in l] for l in blocks]))

        return Polynomial(np.array(coeffs, dtype=complex_type), support_start=joint_support.start)

    @classmethod
    def diagonal_block_matrix(cls, blocks: list["Polynomial"]):
        """Returns a matrix polynomial containing the given polynomials/constants as blocks along the diagonal.

        Raises: ValueError if the block shapes do not match."""
        joint_support = range(min(p.support().start for p in blocks), max(p.support().stop for p in blocks))

        coeffs = []
        for k in joint_support:
            coeffs.append(sp.linalg.block_diag(*[p[k] for p in blocks]))

        return Polynomial(np.array(coeffs, dtype=complex_type), support_start=joint_support.start)

__call__(z) -> complex_type | np.ndarray

Evaluates the polynomial using Horner's method.

Parameters:

Name	Type	Description	Default
z	complex	The point at which to evaluate the polynomial.	required

Returns:

Name	Type	Description
complex	complex_type | ndarray	The evaluated result.

Source code in nlft_qsp/poly.py

def __call__(self, z) -> complex_type | np.ndarray:
    """Evaluates the polynomial using Horner's method.

    Args:
        z (complex): The point at which to evaluate the polynomial.

    Returns:
        complex: The evaluated result.
    """
    res = self.coeffs[-1]
    for k in reversed(range(self.coeffs.shape[0] - 1)):
        res = res * z + self.coeffs[k]
    return res * (z ** self.support_start)

__init__(coeffs: list[complex_type] | np.ndarray = None, support_start: int = 0, shape: tuple[int] = None)

Initializes a Polynomial instance.

Parameters:

Name	Type	Description	Default
coeffs	list[complex_type] | ndarray	List of complex numbers as coefficients.	None
shape	optional	Shape of the coefficient array.	None
support_start	optional	Minimum degree in the polynomial. Defaults to 0.	0

Note

Please provide exactly one of coeffs or shape.

Source code in nlft_qsp/poly.py

def __init__(self, coeffs: list[complex_type] | np.ndarray = None, support_start: int = 0, shape: tuple[int] = None):
    """Initializes a Polynomial instance.

    Args:
        coeffs: List of complex numbers as coefficients.
        shape (optional): Shape of the coefficient array.
        support_start (optional): Minimum degree in the polynomial. Defaults to 0.

    Note:
        Please provide exactly one of coeffs or shape.
    """
    super().__init__(coeffs, support_start, shape)

__str__() -> str

Converts the polynomial to a human-readable string representation.

Returns:

Name	Type	Description
str	str	The string representation of the polynomial.

Source code in nlft_qsp/poly.py

def __str__(self) -> str:
    """Converts the polynomial to a human-readable string representation.

    Returns:
        str: The string representation of the polynomial.
    """
    return ' + '.join(f"{c} z^{self.support_start + k}" for k, c in enumerate(self.coeffs))

analytic_part() -> Polynomial

Discards all the negative degrees, keeping only the non-negative ones.

Returns:

Name	Type	Description
Polynomial	Polynomial	A new polynomial containing only the positive-degree coefficients.

Source code in nlft_qsp/poly.py

def analytic_part(self) -> "Polynomial":
    """Discards all the negative degrees, keeping only the non-negative ones.

    Returns:
        Polynomial: A new polynomial containing only the positive-degree coefficients."""
    return self.truncate(0, self.support_start + self.coeffs.shape[0] - 1)

anti_analytic_part() -> Polynomial

Discards all the positive degrees, keeping only the non-positive ones.

Returns:

Name	Type	Description
Polynomial	Polynomial	A new polynomial containing only the non-positive-degree coefficients.

Source code in nlft_qsp/poly.py

def anti_analytic_part(self) -> "Polynomial":
    """Discards all the positive degrees, keeping only the non-positive ones.

    Returns:
        Polynomial: A new polynomial containing only the non-positive-degree coefficients."""
    return self.truncate(self.support_start, 0)

block_matrix(blocks: list[list[Polynomial]]) classmethod

Returns a matrix polynomial containing the given polynomials/constants as blocks.

Raises: ValueError if the block shapes do not match.

Source code in nlft_qsp/poly.py

@classmethod
def block_matrix(cls, blocks: list[list["Polynomial"]]):
    """Returns a matrix polynomial containing the given polynomials/constants as blocks.

    Raises: ValueError if the block shapes do not match."""
    joint_support = range(min(min(p.support().start for p in l) for l in blocks),
                          max(max(p.support().stop for p in l) for l in blocks))

    coeffs = []
    for k in joint_support:
        coeffs.append(np.block([[p[k] for p in l] for l in blocks]))

    return Polynomial(np.array(coeffs, dtype=complex_type), support_start=joint_support.start)

conjugate() -> Polynomial

Returns the conjugate polynomial on the unit circle. If \(p(z) = \sum_k p_k z^k\), then its conjugate is defined as \(p^*(z) = \sum_k p_k^* z^{-k}\). This is also known as the Schwarz reflection across \(\mathbb{T}\).

Note

For a matrix polynomial, each coefficient is conjugate-transposed.

Source code in nlft_qsp/poly.py

def conjugate(self) -> "Polynomial":
    r"""Returns the conjugate polynomial on the unit circle. If $p(z) = \sum_k p_k z^k$, then its conjugate is defined as $p^*(z) = \sum_k p_k^* z^{-k}$. This is also known as the Schwarz reflection across $\mathbb{T}$.

    Note:
        For a matrix polynomial, each coefficient is conjugate-transposed.
    """
    if self.shape == ():
        conj_coeffs = np.flip(np.conj(self.coeffs), axis=0)
    else:
        axes = [0] + list(range(1, self.coeffs.ndim))[::-1] # transpose = reverse the order of the axes except the first one, which is the index of the coefficient
        conj_coeffs = np.flip(np.transpose(np.conj(self.coeffs), axes=axes), axis=0)

    return Polynomial(conj_coeffs, -(self.support_start + self.coeffs.shape[0] - 1))

diagonal_block_matrix(blocks: list[Polynomial]) classmethod

Returns a matrix polynomial containing the given polynomials/constants as blocks along the diagonal.

Raises: ValueError if the block shapes do not match.

Source code in nlft_qsp/poly.py

@classmethod
def diagonal_block_matrix(cls, blocks: list["Polynomial"]):
    """Returns a matrix polynomial containing the given polynomials/constants as blocks along the diagonal.

    Raises: ValueError if the block shapes do not match."""
    joint_support = range(min(p.support().start for p in blocks), max(p.support().stop for p in blocks))

    coeffs = []
    for k in joint_support:
        coeffs.append(sp.linalg.block_diag(*[p[k] for p in blocks]))

    return Polynomial(np.array(coeffs, dtype=complex_type), support_start=joint_support.start)

duplicate() -> Polynomial

Creates a duplicate of the current polynomial.

Returns:

Name	Type	Description
Polynomial	Polynomial	A new Polynomial instance with the same coefficients and support.

Source code in nlft_qsp/poly.py

def duplicate(self) -> "Polynomial":
    """Creates a duplicate of the current polynomial.

    Returns:
        Polynomial: A new Polynomial instance with the same coefficients and support.
    """
    return type(self)(self.coeffs, self.support_start)

effective_degree() -> int

Returns the size of the support of the polynomial minus 1 (max degree - min degree).

Note

This does not check for leading or trailing zeros in the coefficient array.

Returns:

Name	Type	Description
int	int	The effective degree of the polynomial.

Source code in nlft_qsp/poly.py

def effective_degree(self) -> int:
    """Returns the size of the support of the polynomial minus 1 (max degree - min degree).

    Note:
        This does not check for leading or trailing zeros in the coefficient array.

    Returns:
        int: The effective degree of the polynomial.
    """
    return self.coeffs.shape[0] - 1

eval_at_roots_of_unity(N: int) -> list[complex_type]

Evaluates the polynomial at the \(N\)-th roots of unity using the inverse fast Fourier transform.

Parameters:

Name	Type	Description	Default
N	int	A power of two specifying the number of roots. If \(N\) is not a power of two, then the next power of two is implicitly taken.	required

Returns:

Type	Description
list[complex_type]	list[complex]: List of evaluations at the N-th roots of unity.
list[complex_type]	The \(k\)-th element of the list will be \(P(e^{2\pi i k/N})\), \(P\) being self.

Source code in nlft_qsp/poly.py

def eval_at_roots_of_unity(self, N: int) -> list[complex_type]: # TODO remove power-of-two assumption
    r"""Evaluates the polynomial at the $N$-th roots of unity using the inverse fast Fourier transform.

    Args:
        N (int): A power of two specifying the number of roots. If $N$ is not a power of two, then the next power of two is implicitly taken.

    Returns:
        list[complex]: List of evaluations at the N-th roots of unity.
        The $k$-th element of the list will be $P(e^{2\pi i k/N})$, $P$ being `self`.
    """
    N = next_power_of_two(N)
    M = next_power_of_two(max(N, self.coeffs.shape[0]))

    pad_width = [(0, M - self.coeffs.shape[0])] + [(0, 0)] * (self.coeffs.ndim - 1)
    coeffs = np.pad(self.coeffs, pad_width=pad_width)
    coeffs = np.roll(coeffs, self.support_start, axis=0)
    # This has the effect of having everything multiplied by z^s

    evals = np.fft.ifft(coeffs, norm='forward', axis=0) # M evaluations at the M-th roots of unity
    return evals[::M//N]

hilbert_transform() -> Polynomial

Returns the polynomial \(Q\) such that \(P + Q\) yields an analytic polynomial (\(P\) being self).

Note

This is actually \(i \mathcal{H}[P]\), i.e., the Hilbert transform as returned is already multiplied by \(i\).

Source code in nlft_qsp/poly.py

def hilbert_transform(self) -> "Polynomial":
    r"""Returns the polynomial $Q$ such that $P + Q$ yields an analytic polynomial ($P$ being `self`).

    Note:
        This is actually $i \mathcal{H}[P]$, i.e., the Hilbert transform as returned is already multiplied by $i$.
    """
    hilbert_coeffs = []
    for k in self.support():
        if k < 0:
            hilbert_coeffs.append(-self[k])
        elif k > 0:
            hilbert_coeffs.append(self[k])
        else:
            hilbert_coeffs.append(0)

    return Polynomial(np.array(hilbert_coeffs, dtype=complex_type), self.support_start)

only_negative_degrees() -> Polynomial

DEPRECATED: use anti_analytic_part() instead.

Source code in nlft_qsp/poly.py

def only_negative_degrees(self) -> "Polynomial":
    """DEPRECATED: use `anti_analytic_part()` instead."""
    return self.anti_analytic_part()

only_positive_degrees() -> Polynomial

DEPRECATED: use analytic_part() instead.

Source code in nlft_qsp/poly.py

def only_positive_degrees(self) -> "Polynomial":
    """DEPRECATED: use `analytic_part()` instead."""
    return self.analytic_part()

schwarz_transform() -> Polynomial

Returns the anti-analytic polynomial whose real part gives the current polynomial.

Note

This is equivalent to adding $i \mathcal{H}[p]`, where \(\mathcal{H}[p]\) is the Hilbert transform of \(p\).

Source code in nlft_qsp/poly.py

def schwarz_transform(self) -> "Polynomial":
    r"""Returns the anti-analytic polynomial whose real part gives the current polynomial.

    Note:
        This is equivalent to adding $i \mathcal{H}[p]`, where $\mathcal{H}[p]$ is the Hilbert transform of $p$.
    """
    schwarz_coeffs = []
    for k in self.support():
        if k < 0:
            schwarz_coeffs.append(2*self[k])
        elif k == 0:
            schwarz_coeffs.append(self[k])

    return Polynomial(np.array(schwarz_coeffs, dtype=complex_type), self.support_start)

sharp() -> Polynomial

Same as conjugate(), but support_start is left unchanged.

Returns:

Name	Type	Description
Polynomial	Polynomial	The sharp-conjugate polynomial.

Source code in nlft_qsp/poly.py

def sharp(self) -> "Polynomial":
    r"""Same as `conjugate()`, but `support_start` is left unchanged.

    Returns:
        Polynomial: The sharp-conjugate polynomial.
    """
    p = self.conjugate()
    p.support_start += self.effective_degree() + 1
    return p

shift(k: int)

Creates a new polynomial equal to the current one, multiplied by z^k.

Source code in nlft_qsp/poly.py

def shift(self, k: int):
    """Creates a new polynomial equal to the current one, multiplied by `z^k`."""
    return type(self)(self.coeffs, self.support_start + k)

sup_norm(N=1024) -> float_type

Estimates the supremum norm of the polynomial over the unit circle

Parameters:

Name	Type	Description	Default
N	int	the number of samples to compute the maximum from. If \(N\) is not a power of two, then the next power of two is taken.	1024

Returns:

Name	Type	Description
float_type	float_type	An estimate for the supremum norm of the polynomial over the unit circle.

Source code in nlft_qsp/poly.py

def sup_norm(self, N=1024) -> float_type:
    """Estimates the supremum norm of the polynomial over the unit circle

    Args:
        N (int): the number of samples to compute the maximum from. If $N$ is not a power of two, then the next power of two is taken.

    Returns:
        float_type: An estimate for the supremum norm of the polynomial over the unit circle.
    """
    if self.shape == ():
        return np.max([np.abs(sample) for sample in self.eval_at_roots_of_unity(N)])

    return np.max([np.linalg.norm(sample, ord=2) for sample in self.eval_at_roots_of_unity(N)])

truncate(m: int, n: int) -> Polynomial

Keeps only the coefficients in \([m, n]\), discarding the others.

Returns:

Name	Type	Description
Polynomial	Polynomial	A new, truncated polynomial.

Source code in nlft_qsp/poly.py

def truncate(self, m: int, n: int) -> "Polynomial":
    """Keeps only the coefficients in $[m, n]$, discarding the others.

    Returns:
        Polynomial: A new, truncated polynomial.
    """
    if self.shape == ():
        return Polynomial([self[k] for k in range(m, n+1)], m)

    return Polynomial(np.array([self[k] for k in range(m, n+1)], dtype=complex_type), m)

QSVTPhaseFactors

Bases: ChebyshevQSPPhaseFactors

Phase factors for a QSVT/Reflection QSP protocol.

\[ e^{i\phi_0 Z} \tilde{r} e^{i\phi_1 Z} \tilde{r} \cdots \tilde{r} e^{i\phi_n Z} = \begin{pmatrix} P(x) & \cdot \\ \cdot & \cdot \end{pmatrix} \]

with the Hermitian signal operator \(\tilde{r}\):

\[ \tilde{r} = \begin{pmatrix} x & \sqrt{1 - x^2} \\ \sqrt{1 - x^2} & -x \end{pmatrix} \]

Note

This is the ansatz of Corollary 8 arXiv:1806.01838, but the polynomial construction is implemented by implicitly adjusting the phase factors from Chebyshev QSP, see arXiv:2105.02859, (A5).

Methods:

Name	Description
approximate	Approximate the given callable object \(f\) (which takes \(x \in [-1, 1]\) and returns a real number)
from_chebqsp	Returns a set of QSVT phase factors constructing the same left chebyshev expansion as the given QSP protocol.
solve	Returns the set of phase factors for a QSVT protocol implementing the polynomial \(P(x)\) (as the real part of the top-left polynomial, see Theorem 9 of arXiv:2105.02859).
to_chebqsp	Returns a set of Chebyshev QSP phase factors constructing the same left chebyshev expansion as the given QSVT protocol.

Source code in nlft_qsp/qsp.py

@serializable(
    type_tag="@qspx/phase_factors/qsvt",
    fields={"phi": "phi"}
)
@qsp_variant('qsvt', modes=['c'], display_name='QSVT')
class QSVTPhaseFactors(ChebyshevQSPPhaseFactors):
    r"""Phase factors for a QSVT/Reflection QSP protocol.

    $$ e^{i\phi_0 Z} \tilde{r} e^{i\phi_1 Z} \tilde{r} \cdots \tilde{r} e^{i\phi_n Z} = \begin{pmatrix} P(x) & \cdot \\ \cdot & \cdot \end{pmatrix} $$

    with the Hermitian signal operator $\tilde{r}$:

    $$ \tilde{r} = \begin{pmatrix} x & \sqrt{1 - x^2} \\ \sqrt{1 - x^2} & -x \end{pmatrix} $$

    Note:
        This is the ansatz of Corollary 8 [arXiv:1806.01838](https://arxiv.org/abs/1806.01838), but the polynomial
        construction is implemented by implicitly adjusting the phase factors from Chebyshev QSP, see [arXiv:2105.02859](https://arxiv.org/abs/2105.02859), (A5)."""
    def protocol_conjugation(self, P, Q):
        d = self.degree()
        return ((P + P.conjugate()) + (Q - Q.conjugate()))*((-1j)**(d+1)/2), \
               ((P - P.conjugate()) - (Q + Q.conjugate()))*((-1j)**d/2)

    def processing_operator_conjugation(self, a, b):
        a *= 1j # R(x) = -I*exp(I*pi*Z/4)*W(x)*exp(I*pi*Z/4)
        return np.real(a) + 1j*np.imag(b), 1j*np.imag(a) - np.real(b)

    def duplicate(self):
        return QSVTPhaseFactors(self.phi)

    @classmethod
    def from_chebqsp(cls, pf: ChebyshevQSPPhaseFactors):
        """Returns a set of QSVT phase factors constructing the same left chebyshev expansion as the given QSP protocol."""
        d = pf.degree()
        phi = [0] * (d+1)

        phi[0] = pf.phi[0] + (2*d - 1)*np.pi/4
        for k in range(1, d):
            phi[k] = pf.phi[k] - np.pi/2
        phi[d] = pf.phi[d] - np.pi/4

        return QSVTPhaseFactors(phi)

    def to_chebqsp(self) -> ChebyshevQSPPhaseFactors:
        """Returns a set of Chebyshev QSP phase factors constructing the same left chebyshev expansion as the given QSVT protocol."""
        d = self.degree()
        phi = [0] * (d+1)

        phi[0] = self.phi[0] - (2*d - 1)*np.pi/4
        for k in range(1, d):
            phi[k] = self.phi[k] + np.pi/2
        phi[d] = self.phi[d] + np.pi/4

        return ChebyshevQSPPhaseFactors(phi)

    @classmethod
    def solve(cls, T: list[complex_type] | Polynomial | ChebyshevTExpansion) -> "QSVTPhaseFactors":
        r"""Returns the set of phase factors for a QSVT protocol implementing the polynomial $P(x)$ (as the real part of the top-left polynomial, see Theorem 9 of [arXiv:2105.02859](https://arxiv.org/abs/2105.02859)).

        The target polynomial will be $T(x) = \sum_{k = 0}^n c_k T_k(x)$ (if `T` is a `ChebyshevTExpansion`) or $T(x) = \sum_{k = 0}^n c_k x^k$ (if `T` is a `Polynomial`), where $T_k(x)$ are the Chebyshev polynomials of the first kind.

        Args:
            T: a Chebyshev expansion object or the desired Polynomial $P(x)$ (that will be converted to the Chebyshev basis).

        Raises:
            ValueError: If the target polynomial does not have definite parity or is not real.

        Note:
            `T` can also be a list of complex numbers. This will be regarded as coefficients in the Chebyshev basis. Passing the list directly is discouraged and will be removed in future releases."""
        return QSVTPhaseFactors.from_chebqsp(ChebyshevQSPPhaseFactors.solve(T))

    @classmethod
    def approximate(cls, f: Callable, deg: int) -> "QSVTPhaseFactors":
        r"""Approximate the given callable object $f$ (which takes $x \in [-1, 1]$ and returns a real number)
            and returns the QSVT phase factors implementing an approximating polynomial of degree `deg`
            (as the real part of the top-left polynomial, see Theorem 9 of [arXiv:2105.02859](https://arxiv.org/abs/2105.02859)).

            Note: The parity of `deg` should coincide with the parity of $f$, otherwise the Chebyshev approximator might give numerical errors."""
        return cls.solve(chebyshev_approximate(f, deg))

approximate(f: Callable, deg: int) -> QSVTPhaseFactors classmethod

Approximate the given callable object \(f\) (which takes \(x \in [-1, 1]\) and returns a real number) and returns the QSVT phase factors implementing an approximating polynomial of degree deg (as the real part of the top-left polynomial, see Theorem 9 of arXiv:2105.02859).

Note: The parity of deg should coincide with the parity of \(f\), otherwise the Chebyshev approximator might give numerical errors.

Source code in nlft_qsp/qsp.py

@classmethod
def approximate(cls, f: Callable, deg: int) -> "QSVTPhaseFactors":
    r"""Approximate the given callable object $f$ (which takes $x \in [-1, 1]$ and returns a real number)
        and returns the QSVT phase factors implementing an approximating polynomial of degree `deg`
        (as the real part of the top-left polynomial, see Theorem 9 of [arXiv:2105.02859](https://arxiv.org/abs/2105.02859)).

        Note: The parity of `deg` should coincide with the parity of $f$, otherwise the Chebyshev approximator might give numerical errors."""
    return cls.solve(chebyshev_approximate(f, deg))

from_chebqsp(pf: ChebyshevQSPPhaseFactors) classmethod

Returns a set of QSVT phase factors constructing the same left chebyshev expansion as the given QSP protocol.

Source code in nlft_qsp/qsp.py

@classmethod
def from_chebqsp(cls, pf: ChebyshevQSPPhaseFactors):
    """Returns a set of QSVT phase factors constructing the same left chebyshev expansion as the given QSP protocol."""
    d = pf.degree()
    phi = [0] * (d+1)

    phi[0] = pf.phi[0] + (2*d - 1)*np.pi/4
    for k in range(1, d):
        phi[k] = pf.phi[k] - np.pi/2
    phi[d] = pf.phi[d] - np.pi/4

    return QSVTPhaseFactors(phi)

solve(T: list[complex_type] | Polynomial | ChebyshevTExpansion) -> QSVTPhaseFactors classmethod

Returns the set of phase factors for a QSVT protocol implementing the polynomial \(P(x)\) (as the real part of the top-left polynomial, see Theorem 9 of arXiv:2105.02859).

The target polynomial will be \(T(x) = \sum_{k = 0}^n c_k T_k(x)\) (if T is a ChebyshevTExpansion) or \(T(x) = \sum_{k = 0}^n c_k x^k\) (if T is a Polynomial), where \(T_k(x)\) are the Chebyshev polynomials of the first kind.

Parameters:

Name	Type	Description	Default
T	list[complex_type] | Polynomial | ChebyshevTExpansion	a Chebyshev expansion object or the desired Polynomial \(P(x)\) (that will be converted to the Chebyshev basis).	required

Raises:

Type	Description
ValueError	If the target polynomial does not have definite parity or is not real.

Note

T can also be a list of complex numbers. This will be regarded as coefficients in the Chebyshev basis. Passing the list directly is discouraged and will be removed in future releases.

Source code in nlft_qsp/qsp.py

@classmethod
def solve(cls, T: list[complex_type] | Polynomial | ChebyshevTExpansion) -> "QSVTPhaseFactors":
    r"""Returns the set of phase factors for a QSVT protocol implementing the polynomial $P(x)$ (as the real part of the top-left polynomial, see Theorem 9 of [arXiv:2105.02859](https://arxiv.org/abs/2105.02859)).

    The target polynomial will be $T(x) = \sum_{k = 0}^n c_k T_k(x)$ (if `T` is a `ChebyshevTExpansion`) or $T(x) = \sum_{k = 0}^n c_k x^k$ (if `T` is a `Polynomial`), where $T_k(x)$ are the Chebyshev polynomials of the first kind.

    Args:
        T: a Chebyshev expansion object or the desired Polynomial $P(x)$ (that will be converted to the Chebyshev basis).

    Raises:
        ValueError: If the target polynomial does not have definite parity or is not real.

    Note:
        `T` can also be a list of complex numbers. This will be regarded as coefficients in the Chebyshev basis. Passing the list directly is discouraged and will be removed in future releases."""
    return QSVTPhaseFactors.from_chebqsp(ChebyshevQSPPhaseFactors.solve(T))

to_chebqsp() -> ChebyshevQSPPhaseFactors

Returns a set of Chebyshev QSP phase factors constructing the same left chebyshev expansion as the given QSVT protocol.

Source code in nlft_qsp/qsp.py

def to_chebqsp(self) -> ChebyshevQSPPhaseFactors:
    """Returns a set of Chebyshev QSP phase factors constructing the same left chebyshev expansion as the given QSVT protocol."""
    d = self.degree()
    phi = [0] * (d+1)

    phi[0] = self.phi[0] - (2*d - 1)*np.pi/4
    for k in range(1, d):
        phi[k] = self.phi[k] + np.pi/2
    phi[d] = self.phi[d] + np.pi/4

    return ChebyshevQSPPhaseFactors(phi)

XQSPPhaseFactors

Bases: PhaseFactors

Phase factors for a XQSP protocol.

\[ e^{i\phi_0 X} W(z) e^{i\phi_1 X} W(z) \cdots W(z) e^{i\phi_n X} = \begin{pmatrix} P(z) & Q(z) \\ \cdot & \cdot \end{pmatrix} \]

where \(W(z) = \mathrm{diag}(z, 1)\) (mode='analytic') or \(W(z) = \mathrm{diag}(z, z^{-1})\) (mode='laurent').

Methods:

Name	Description
from_nlfs	Computes the XQSP phase factors for a given imaginary NLFT sequence.
solve	Returns the set of phase factors for a XQSP protocol producing the given polynomial. See here for an overview of the QSP variants.
solve_laurent	Returns the set of phase factors for a XQSP protocol producing the given definite-parity Laurent polynomial. See here for an overview of the QSP variants.

Source code in nlft_qsp/qsp.py

@serializable(
    type_tag="@qspx/phase_factors/xqsp",
    fields={"phi": "phi"}
)
@qsp_variant('x', modes=['l', 'a'], display_name='XQSP')
class XQSPPhaseFactors(PhaseFactors):
    r"""Phase factors for a XQSP protocol.

    $$ e^{i\phi_0 X} W(z) e^{i\phi_1 X} W(z) \cdots W(z) e^{i\phi_n X} = \begin{pmatrix} P(z) & Q(z) \\ \cdot & \cdot \end{pmatrix} $$

    where $W(z) = \mathrm{diag}(z, 1)$ (`mode='analytic'`) or $W(z) = \mathrm{diag}(z, z^{-1})$ (`mode='laurent'`).
    """
    def __init__(self, phi: list[float_type]):
        self.phi = list(phi)

    def duplicate(self):
        return XQSPPhaseFactors(self.phi)

    def processing_operator(self, k: int): # exp(i phi[k] X)
        return np.cos(self.phi[k]), 1j*np.sin(self.phi[k])

    def degree(self):
        return len(self.phi) - 1

    def iX(self):
        pf = self.duplicate()
        pf.phi[-1] += np.pi/2
        return pf

    def to_nlfs(self) -> NonLinearFourierSequence:
        return NonLinearFourierSequence([1j*np.tan(phik) for phik in self.phi])

    @classmethod
    def from_nlfs(cls, F: NonLinearFourierSequence) -> PhaseFactors:
        r"""Computes the XQSP phase factors for a given imaginary NLFT sequence.
        If $\mathrm{NLFT}(F) = (a, b)$, then the returned phase factors will implement $(z^n a, b)$ in the analytic picture.

        Args:
            F (NonLinearFourierSequence): The imaginary sequence to be converted to phase factors.

        Raises:
            ValueError: if $F$ is not imaginary.

        Note:
            The support start of $F$ is ignored, so the support of $b$ is assumed to start at $0$."""
        if not F.is_imaginary():
            raise ValueError("The nonlinear Fourier sequence must be imaginary in order to be turned into a XQSP protocol.")

        return XQSPPhaseFactors([np.arctan(np.imag(Fk)) for Fk in F.coeffs])

    @classmethod
    def solve(cls, P: Polynomial, convention='qsp') -> "XQSPPhaseFactors":
        r"""Returns the set of phase factors for a XQSP protocol producing the given polynomial. See [here](https://arxiv.org/abs/2503.03026) for an overview of the QSP variants.
            A complementary $Q$ will be computed with Weiss' algorithm.

            Args:
                convention (str): Whether the phase factors should produce $(P, Q)$ (`'qsp'`), or $(Q, P)$ (`'nlft'`).

            Raises:
                ValueError: If $P$ does not lie in the XQSP subalgebra.

            Note:
                The sup norm of $P$ should be bounded by $1 - \eta < 1$.
                The time required by the algorithm to compute the phase factors will scale with $1/\eta$.

                The support_start of $P$ will be ignored. This is a solver for analytic QSP. In order to obtain phase factors for Laurent XQSP, use `XQSPPhaseFactors.solve_laurent()`."""
        if 1 - P.sup_norm(4*P.effective_degree()) < bd.machine_threshold():
            raise ValueError("The given polynomial cannot be too close to or larger than one on the unit circle.")

        match convention:
            case "qsp":
                P = -1j * Polynomial(P.coeffs, 0) # (Q, -iP) -> (P, iQ)
            case "nlft":
                P = Polynomial(P.coeffs, 0)
            case _:
                raise ValueError("The given mode does not exist. Only modes available are 'qsp', 'nlft'.")

        F = _riemann_hilbert_weiss(P) # NLFT(F) = (Q, P)

        if convention != "qsp":
            return XQSPPhaseFactors.from_nlfs(F)
        return XQSPPhaseFactors.from_nlfs(F).iX()

    @classmethod
    def solve_laurent(cls, P: Polynomial, convention='qsp') -> "XQSPPhaseFactors":
        r"""Returns the set of phase factors for a XQSP protocol producing the given definite-parity Laurent polynomial. See [here](https://arxiv.org/abs/2503.03026) for an overview of the QSP variants.
        A complementary $Q$ will be computed with Weiss' algorithm.

        Args:
            convention (str): Whether the phase factors should produce $(P, Q)$ (`'qsp'`), or $(Q, P)$ (`'nlft'`).

        Raises:
            ValueError: If $P$ does not lie in the X-constrained subalgebra or $P$ has not definite-parity.

        Note:
            The sup norm of $P$ should be bounded by $1 - \eta < 1$.
            The time required by the algorithm to compute the phase factors will scale with $1/\eta$.

            The support_start of $P$ will be ignored. In order to obtain phase factors for Laurent XQSP, first convert the polynomial into analytic form."""
        if not is_definite_parity(P):
            raise ValueError("Laurent polynomial is not of definite parity.")

        return cls.solve(laurent_to_analytic(P), convention)

from_nlfs(F: NonLinearFourierSequence) -> PhaseFactors classmethod

Computes the XQSP phase factors for a given imaginary NLFT sequence. If \(\mathrm{NLFT}(F) = (a, b)\), then the returned phase factors will implement \((z^n a, b)\) in the analytic picture.

Parameters:

Name	Type	Description	Default
F	NonLinearFourierSequence	The imaginary sequence to be converted to phase factors.	required

Raises:

Type	Description
ValueError	if \(F\) is not imaginary.

Note

The support start of \(F\) is ignored, so the support of \(b\) is assumed to start at \(0\).

Source code in nlft_qsp/qsp.py

@classmethod
def from_nlfs(cls, F: NonLinearFourierSequence) -> PhaseFactors:
    r"""Computes the XQSP phase factors for a given imaginary NLFT sequence.
    If $\mathrm{NLFT}(F) = (a, b)$, then the returned phase factors will implement $(z^n a, b)$ in the analytic picture.

    Args:
        F (NonLinearFourierSequence): The imaginary sequence to be converted to phase factors.

    Raises:
        ValueError: if $F$ is not imaginary.

    Note:
        The support start of $F$ is ignored, so the support of $b$ is assumed to start at $0$."""
    if not F.is_imaginary():
        raise ValueError("The nonlinear Fourier sequence must be imaginary in order to be turned into a XQSP protocol.")

    return XQSPPhaseFactors([np.arctan(np.imag(Fk)) for Fk in F.coeffs])

solve(P: Polynomial, convention='qsp') -> XQSPPhaseFactors classmethod

Returns the set of phase factors for a XQSP protocol producing the given polynomial. See here for an overview of the QSP variants. A complementary \(Q\) will be computed with Weiss' algorithm.

Parameters:

Name	Type	Description	Default
convention	str	Whether the phase factors should produce \((P, Q)\) ('qsp'), or \((Q, P)\) ('nlft').	'qsp'

Raises:

Type	Description
ValueError	If \(P\) does not lie in the XQSP subalgebra.

Note

The sup norm of \(P\) should be bounded by \(1 - \eta < 1\). The time required by the algorithm to compute the phase factors will scale with \(1/\eta\).

The support_start of \(P\) will be ignored. This is a solver for analytic QSP. In order to obtain phase factors for Laurent XQSP, use XQSPPhaseFactors.solve_laurent().

Source code in nlft_qsp/qsp.py

@classmethod
def solve(cls, P: Polynomial, convention='qsp') -> "XQSPPhaseFactors":
    r"""Returns the set of phase factors for a XQSP protocol producing the given polynomial. See [here](https://arxiv.org/abs/2503.03026) for an overview of the QSP variants.
        A complementary $Q$ will be computed with Weiss' algorithm.

        Args:
            convention (str): Whether the phase factors should produce $(P, Q)$ (`'qsp'`), or $(Q, P)$ (`'nlft'`).

        Raises:
            ValueError: If $P$ does not lie in the XQSP subalgebra.

        Note:
            The sup norm of $P$ should be bounded by $1 - \eta < 1$.
            The time required by the algorithm to compute the phase factors will scale with $1/\eta$.

            The support_start of $P$ will be ignored. This is a solver for analytic QSP. In order to obtain phase factors for Laurent XQSP, use `XQSPPhaseFactors.solve_laurent()`."""
    if 1 - P.sup_norm(4*P.effective_degree()) < bd.machine_threshold():
        raise ValueError("The given polynomial cannot be too close to or larger than one on the unit circle.")

    match convention:
        case "qsp":
            P = -1j * Polynomial(P.coeffs, 0) # (Q, -iP) -> (P, iQ)
        case "nlft":
            P = Polynomial(P.coeffs, 0)
        case _:
            raise ValueError("The given mode does not exist. Only modes available are 'qsp', 'nlft'.")

    F = _riemann_hilbert_weiss(P) # NLFT(F) = (Q, P)

    if convention != "qsp":
        return XQSPPhaseFactors.from_nlfs(F)
    return XQSPPhaseFactors.from_nlfs(F).iX()

solve_laurent(P: Polynomial, convention='qsp') -> XQSPPhaseFactors classmethod

Returns the set of phase factors for a XQSP protocol producing the given definite-parity Laurent polynomial. See here for an overview of the QSP variants. A complementary \(Q\) will be computed with Weiss' algorithm.

Parameters:

Name	Type	Description	Default
convention	str	Whether the phase factors should produce \((P, Q)\) ('qsp'), or \((Q, P)\) ('nlft').	'qsp'

Raises:

Type	Description
ValueError	If \(P\) does not lie in the X-constrained subalgebra or \(P\) has not definite-parity.

Note

The sup norm of \(P\) should be bounded by \(1 - \eta < 1\). The time required by the algorithm to compute the phase factors will scale with \(1/\eta\).

The support_start of \(P\) will be ignored. In order to obtain phase factors for Laurent XQSP, first convert the polynomial into analytic form.

Source code in nlft_qsp/qsp.py

@classmethod
def solve_laurent(cls, P: Polynomial, convention='qsp') -> "XQSPPhaseFactors":
    r"""Returns the set of phase factors for a XQSP protocol producing the given definite-parity Laurent polynomial. See [here](https://arxiv.org/abs/2503.03026) for an overview of the QSP variants.
    A complementary $Q$ will be computed with Weiss' algorithm.

    Args:
        convention (str): Whether the phase factors should produce $(P, Q)$ (`'qsp'`), or $(Q, P)$ (`'nlft'`).

    Raises:
        ValueError: If $P$ does not lie in the X-constrained subalgebra or $P$ has not definite-parity.

    Note:
        The sup norm of $P$ should be bounded by $1 - \eta < 1$.
        The time required by the algorithm to compute the phase factors will scale with $1/\eta$.

        The support_start of $P$ will be ignored. In order to obtain phase factors for Laurent XQSP, first convert the polynomial into analytic form."""
    if not is_definite_parity(P):
        raise ValueError("Laurent polynomial is not of definite parity.")

    return cls.solve(laurent_to_analytic(P), convention)

YQSPPhaseFactors

Bases: PhaseFactors

Phase factors for a YQSP protocol.

\[ e^{i\phi_0 Y} W(z) e^{i\phi_1 Y} W(z) \cdots W(z) e^{i\phi_n Y} = \begin{pmatrix} P(z) & Q(z) \\ \cdot & \cdot \end{pmatrix} \]

where \(W(z) = \mathrm{diag}(z, 1)\) (mode='analytic') or \(W(z) = \mathrm{diag}(z, z^{-1})\) (mode='laurent').

Methods:

Name	Description
from_nlfs	Computes the YQSP phase factors for a given real NLFT sequence.
solve	Returns the set of phase factors for a YQSP protocol producing the given polynomial. See here for an overview of the QSP variants.
solve_laurent	Returns the set of phase factors for a YQSP protocol producing the given definite-parity polynomial. See here for an overview of the QSP variants.

Source code in nlft_qsp/qsp.py

@serializable(
    type_tag="@qspx/phase_factors/yqsp",
    fields={"phi": "phi"}
)
@qsp_variant('y', modes=['l', 'a'], display_name='YQSP')
class YQSPPhaseFactors(PhaseFactors):
    r"""Phase factors for a YQSP protocol.

        $$ e^{i\phi_0 Y} W(z) e^{i\phi_1 Y} W(z) \cdots W(z) e^{i\phi_n Y} = \begin{pmatrix} P(z) & Q(z) \\ \cdot & \cdot \end{pmatrix} $$

        where $W(z) = \mathrm{diag}(z, 1)$ (`mode='analytic'`) or $W(z) = \mathrm{diag}(z, z^{-1})$ (`mode='laurent'`).
        """
    def __init__(self, phi: list[float_type]):
        self.phi = list(phi)

    def duplicate(self):
        return YQSPPhaseFactors(self.phi)

    def processing_operator(self, k: int): # exp(i phi[k] Y)
        return np.cos(self.phi[k]), np.sin(self.phi[k])

    def degree(self):
        return len(self.phi) - 1

    def iY(self):
        pf = self.duplicate()
        pf.phi[-1] += np.pi/2
        return pf

    def to_nlfs(self) -> NonLinearFourierSequence:
        return NonLinearFourierSequence([np.tan(phik) for phik in self.phi])

    @classmethod
    def from_nlfs(cls, F: NonLinearFourierSequence) -> PhaseFactors:
        r"""Computes the YQSP phase factors for a given real NLFT sequence.
        If $\mathrm{NLFT}(F) = (a, b)$, then the returned phase factors will implement $(z^n a, b)$ in the analytic picture.

        Args:
            F (NonLinearFourierSequence): The real sequence to be converted to phase factors.

        Raises:
            ValueError: if $F$ is not real.

        Note:
            The support start of $F$ is ignored, so the support of $b$ is assumed to start at $0$."""
        if not F.is_real():
            raise ValueError("The Non-Linear Fourier sequence must be real in order to be turned into a YQSP protocol.")

        return YQSPPhaseFactors([np.arctan(np.real(Fk)) for Fk in F.coeffs])

    @classmethod
    def solve(cls, P: Polynomial, convention='qsp') -> "YQSPPhaseFactors":
        r"""Returns the set of phase factors for a YQSP protocol producing the given polynomial. See [here](https://arxiv.org/abs/2503.03026) for an overview of the QSP variants.
        A complementary $Q$ will be computed with Weiss' algorithm.

        Args:
            convention (str): Whether the phase factors should produce $(P, Q)$ (`'qsp'`), or $(Q, P)$ (`'nlft'`).

        Raises:
            ValueError: If $P$ does not lie in the YQSP subalgebra.

        Note:
            The sup norm of $P$ should be bounded by $1 - \eta < 1$.
            The time required by the algorithm to compute the phase factors will scale with $1/\eta$.

            The support_start of $P$ will be ignored. This is a solver for analytic QSP. In order to obtain phase factors for Laurent YQSP, use `YQSPPhaseFactors.solve_laurent()`."""
        if 1 - P.sup_norm(4*P.effective_degree()) < bd.machine_threshold():
            raise ValueError("The given polynomial cannot be too close to or larger than one on the unit circle.")

        match convention:
            case "qsp":
                P = -Polynomial(P.coeffs, 0) # (Q, -P) -> (P, Q)
            case "nlft":
                P = Polynomial(P.coeffs, 0)
            case _:
                raise ValueError("The given mode does not exist. Only modes available are 'qsp', 'nlft'.")

        F = _riemann_hilbert_weiss(P) # NLFT(F) = (Q, P)

        if convention != "qsp":
            return YQSPPhaseFactors.from_nlfs(F)
        return YQSPPhaseFactors.from_nlfs(F).iY()

    @classmethod
    def solve_laurent(cls, P: Polynomial, convention='qsp') -> "YQSPPhaseFactors":
        r"""Returns the set of phase factors for a YQSP protocol producing the given definite-parity polynomial. See [here](https://arxiv.org/abs/2503.03026) for an overview of the QSP variants.
        A complementary $Q$ will be computed with Weiss' algorithm.

        Args:
            convention (str): Whether the phase factors should produce $(P, Q)$ (`'qsp'`), or $(Q, P)$ (`'nlft'`).

        Raises:
            ValueError: If $P$ does not lie in the YQSP subalgebra or $P$ has not definite-parity.

        Note:
            The sup norm of $P$ should be bounded by $1 - \eta < 1$.
            The time required by the algorithm to compute the phase factors will scale with $1/\eta$.
            """
        if not is_definite_parity(P):
            raise ValueError("Laurent polynomial is not of definite parity.")

        return cls.solve(laurent_to_analytic(P), convention)

from_nlfs(F: NonLinearFourierSequence) -> PhaseFactors classmethod

Computes the YQSP phase factors for a given real NLFT sequence. If \(\mathrm{NLFT}(F) = (a, b)\), then the returned phase factors will implement \((z^n a, b)\) in the analytic picture.

Parameters:

Name	Type	Description	Default
F	NonLinearFourierSequence	The real sequence to be converted to phase factors.	required

Raises:

Type	Description
ValueError	if \(F\) is not real.

Note

The support start of \(F\) is ignored, so the support of \(b\) is assumed to start at \(0\).

Source code in nlft_qsp/qsp.py

@classmethod
def from_nlfs(cls, F: NonLinearFourierSequence) -> PhaseFactors:
    r"""Computes the YQSP phase factors for a given real NLFT sequence.
    If $\mathrm{NLFT}(F) = (a, b)$, then the returned phase factors will implement $(z^n a, b)$ in the analytic picture.

    Args:
        F (NonLinearFourierSequence): The real sequence to be converted to phase factors.

    Raises:
        ValueError: if $F$ is not real.

    Note:
        The support start of $F$ is ignored, so the support of $b$ is assumed to start at $0$."""
    if not F.is_real():
        raise ValueError("The Non-Linear Fourier sequence must be real in order to be turned into a YQSP protocol.")

    return YQSPPhaseFactors([np.arctan(np.real(Fk)) for Fk in F.coeffs])

solve(P: Polynomial, convention='qsp') -> YQSPPhaseFactors classmethod

Returns the set of phase factors for a YQSP protocol producing the given polynomial. See here for an overview of the QSP variants. A complementary \(Q\) will be computed with Weiss' algorithm.

Parameters:

Name	Type	Description	Default
convention	str	Whether the phase factors should produce \((P, Q)\) ('qsp'), or \((Q, P)\) ('nlft').	'qsp'

Raises:

Type	Description
ValueError	If \(P\) does not lie in the YQSP subalgebra.

Note

The sup norm of \(P\) should be bounded by \(1 - \eta < 1\). The time required by the algorithm to compute the phase factors will scale with \(1/\eta\).

The support_start of \(P\) will be ignored. This is a solver for analytic QSP. In order to obtain phase factors for Laurent YQSP, use YQSPPhaseFactors.solve_laurent().

Source code in nlft_qsp/qsp.py

@classmethod
def solve(cls, P: Polynomial, convention='qsp') -> "YQSPPhaseFactors":
    r"""Returns the set of phase factors for a YQSP protocol producing the given polynomial. See [here](https://arxiv.org/abs/2503.03026) for an overview of the QSP variants.
    A complementary $Q$ will be computed with Weiss' algorithm.

    Args:
        convention (str): Whether the phase factors should produce $(P, Q)$ (`'qsp'`), or $(Q, P)$ (`'nlft'`).

    Raises:
        ValueError: If $P$ does not lie in the YQSP subalgebra.

    Note:
        The sup norm of $P$ should be bounded by $1 - \eta < 1$.
        The time required by the algorithm to compute the phase factors will scale with $1/\eta$.

        The support_start of $P$ will be ignored. This is a solver for analytic QSP. In order to obtain phase factors for Laurent YQSP, use `YQSPPhaseFactors.solve_laurent()`."""
    if 1 - P.sup_norm(4*P.effective_degree()) < bd.machine_threshold():
        raise ValueError("The given polynomial cannot be too close to or larger than one on the unit circle.")

    match convention:
        case "qsp":
            P = -Polynomial(P.coeffs, 0) # (Q, -P) -> (P, Q)
        case "nlft":
            P = Polynomial(P.coeffs, 0)
        case _:
            raise ValueError("The given mode does not exist. Only modes available are 'qsp', 'nlft'.")

    F = _riemann_hilbert_weiss(P) # NLFT(F) = (Q, P)

    if convention != "qsp":
        return YQSPPhaseFactors.from_nlfs(F)
    return YQSPPhaseFactors.from_nlfs(F).iY()

solve_laurent(P: Polynomial, convention='qsp') -> YQSPPhaseFactors classmethod

Returns the set of phase factors for a YQSP protocol producing the given definite-parity polynomial. See here for an overview of the QSP variants. A complementary \(Q\) will be computed with Weiss' algorithm.

Parameters:

Name	Type	Description	Default
convention	str	Whether the phase factors should produce \((P, Q)\) ('qsp'), or \((Q, P)\) ('nlft').	'qsp'

Raises:

Type	Description
ValueError	If \(P\) does not lie in the YQSP subalgebra or \(P\) has not definite-parity.

Note

The sup norm of \(P\) should be bounded by \(1 - \eta < 1\). The time required by the algorithm to compute the phase factors will scale with \(1/\eta\).

Source code in nlft_qsp/qsp.py

@classmethod
def solve_laurent(cls, P: Polynomial, convention='qsp') -> "YQSPPhaseFactors":
    r"""Returns the set of phase factors for a YQSP protocol producing the given definite-parity polynomial. See [here](https://arxiv.org/abs/2503.03026) for an overview of the QSP variants.
    A complementary $Q$ will be computed with Weiss' algorithm.

    Args:
        convention (str): Whether the phase factors should produce $(P, Q)$ (`'qsp'`), or $(Q, P)$ (`'nlft'`).

    Raises:
        ValueError: If $P$ does not lie in the YQSP subalgebra or $P$ has not definite-parity.

    Note:
        The sup norm of $P$ should be bounded by $1 - \eta < 1$.
        The time required by the algorithm to compute the phase factors will scale with $1/\eta$.
        """
    if not is_definite_parity(P):
        raise ValueError("Laurent polynomial is not of definite parity.")

    return cls.solve(laurent_to_analytic(P), convention)

chebqsp_approximate(f, deg: int) -> ChebyshevQSPPhaseFactors

DEPRECATED: use ChebyshevQSPPhaseFactors.approximate() instead.

Source code in nlft_qsp/qsp.py

def chebqsp_approximate(f, deg: int) -> ChebyshevQSPPhaseFactors:
    r"""DEPRECATED: use `ChebyshevQSPPhaseFactors.approximate()` instead."""
    return ChebyshevQSPPhaseFactors.approximate(f, deg)

chebqsp_solve(T: list[complex_type] | ChebyshevTExpansion) -> ChebyshevQSPPhaseFactors

DEPRECATED: use ChebyshevQSPPhaseFactors.solve() instead.

Source code in nlft_qsp/qsp.py

def chebqsp_solve(T: list[complex_type] | ChebyshevTExpansion) -> ChebyshevQSPPhaseFactors:
    r"""DEPRECATED: use `ChebyshevQSPPhaseFactors.solve()` instead."""
    return ChebyshevQSPPhaseFactors.solve(T)

chebyshev_approximate(f: Callable, N: int) -> ChebyshevTExpansion

Computes the Chebyshev expansion up to \(N\) for a complex-valued function \(f : [-1, 1] \rightarrow \mathbb{C}\).

Parameters:

Name	Type	Description	Default
f	Callable	complex-valued function f(x)	required
N	int	degree of Chebyshev approximation	required

Source code in nlft_qsp/approximate.py

def chebyshev_approximate(f: Callable, N: int) -> ChebyshevTExpansion:
    r"""Computes the Chebyshev expansion up to $N$ for a complex-valued function $f : [-1, 1] \rightarrow \mathbb{C}$.

    Args:
        f (Callable): complex-valued function f(x)
        N (int): degree of Chebyshev approximation
    """
    x = np.cos(np.pi * np.arange(N + 1) / N)
    fx = [f(xk) for xk in x]

    fx_m = fx + fx[1:-1][::-1]

    F = np.fft.fft(fx_m)
    F = F[:N + 1] / N
    F[0]  /= 2
    F[-1] /= 2

    return ChebyshevTExpansion(F)

fourier_approximate(f: Callable, N: int) -> Polynomial

Computes the Fourier series of the given function \(f(z)\), \(z = e^{i\theta}\) being a complex number of unit modulus.

Parameters:

Name	Type	Description	Default
f	Callable	a function taking a complex number \(z\) and returning a complex number.	required
N	int	The degree of approximation.	required

Returns:

Name	Type	Description
Polynomial	Polynomial	A Laurent polynomial of degrees in \(\{-N, N+1, \ldots, N-1\}\) approximating \(f\)

Source code in nlft_qsp/approximate.py

def fourier_approximate(f: Callable, N: int) -> Polynomial:
    r"""Computes the Fourier series of the given function $f(z)$, $z = e^{i\theta}$ being a complex number of unit modulus.

    Args:
        f (Callable): a function taking a complex number $z$ and returning a complex number.
        N (int): The degree of approximation.

    Returns:
        Polynomial: A Laurent polynomial of degrees in $\{-N, N+1, \ldots, N-1\}$ approximating $f$
    """
    M = next_power_of_two(2*N+1)

    zk = np.exp(2j * np.pi * np.arange(M) / M)
    fk = np.array([f(z) for z in zk], dtype=complex)

    pk = (np.fft.fft(fk) / M).tolist()

    return Polynomial(pk[M-N:] + pk[:N+1], support_start=-N)

gqsp_solve(P: Polynomial, mode='qsp') -> GQSPPhaseFactors

DEPRECATED: use GQSPPhaseFactors.solve() instead.

Source code in nlft_qsp/qsp.py

def gqsp_solve(P: Polynomial, mode='qsp') -> GQSPPhaseFactors:
    r"""DEPRECATED: use `GQSPPhaseFactors.solve()` instead."""
    return GQSPPhaseFactors.solve(P, mode)

plot_chebyshev(funcs: dict, num_points: int = 1000)

Plots the real part of each object in funcs over the interval \([-1, 1]\). These can be Python functions, Polynomial objects, ChebyshevTExpansion objects, or any callable object.

Parameters:

Name	Type	Description	Default
funcs	dict	a dictionary where each key is the name appearing in the legend of the corresponding function plot.	required
num_points	int	number of sampling points.	1000

Source code in nlft_qsp/plot.py

def plot_chebyshev(funcs: dict, num_points: int=1000):
    r"""
    Plots the real part of each object in funcs over the interval $[-1, 1]$.
    These can be Python functions, `Polynomial` objects, `ChebyshevTExpansion` objects, or
    any callable object.

    Args:
        funcs (dict): a dictionary where each key is the name appearing in the legend of the corresponding function plot.
        num_points (int): number of sampling points.
    """
    plt.figure(figsize=(6, 3))

    for name, f in funcs.items():
        try:
            x_vals = [-1 + 2*k/num_points for k in range(num_points+1)]
            y_vals = [np.real(f(x)) for x in x_vals]
            plt.plot(x_vals, y_vals, label=name)
        except Exception as e:
            print(f"Error evaluating function {name}: {e}")

    plt.xlabel("x")
    plt.legend()
    plt.grid(True)
    plt.show()

plot_fourier(funcs: dict, num_points: int = 1000)

Plots the absolute value of each object in funcs over the unit circle, i.e., plugging \(z = e^{it}\) for \(t \in [-\pi, \pi]\). These can be Python functions, Polynomial objects, or any callable object.

Parameters:

Name	Type	Description	Default
funcs	dict	a dictionary where each key is the name appearing in the legend of the corresponding function plot.	required
num_points	int	number of sampling points.	1000

Source code in nlft_qsp/plot.py

def plot_fourier(funcs: dict, num_points: int=1000):
    r"""
    Plots the absolute value of each object in funcs over the unit circle, i.e., plugging $z = e^{it}$ for $t \in [-\pi, \pi]$.
    These can be Python functions, `Polynomial` objects, or any callable object.

    Args:
        funcs (dict): a dictionary where each key is the name appearing in the legend of the corresponding function plot.
        num_points (int): number of sampling points.
    """
    plt.figure(figsize=(6, 3))

    for name, f in funcs.items():
        try:
            x_vals = [-np.pi + 2*np.pi*k/num_points for k in range(num_points+1)]
            y_vals = [np.abs(f(np.exp(1j * x))) for x in x_vals]
            plt.plot(x_vals, y_vals, label=name)
        except Exception as e:
            print(f"Error evaluating function {name}: {e}")

    plt.xlabel("x")
    plt.xlim(left=-np.pi, right=np.pi)
    plt.legend()
    plt.grid(True)
    plt.show()

qsvt_solve(T: list[complex_type] | ChebyshevTExpansion) -> QSVTPhaseFactors

DEPRECATED: use QSVTPhaseFactors.solve() instead.

Source code in nlft_qsp/qsp.py

def qsvt_solve(T: list[complex_type] | ChebyshevTExpansion) -> QSVTPhaseFactors:
    r"""DEPRECATED: use `QSVTPhaseFactors.solve()` instead."""
    return QSVTPhaseFactors.solve(T)

xqsp_solve(P: Polynomial, mode='qsp') -> XQSPPhaseFactors

DEPRECATED: use XQSPPhaseFactors.solve() instead.

Source code in nlft_qsp/qsp.py

def xqsp_solve(P: Polynomial, mode='qsp') -> XQSPPhaseFactors:
    r"""DEPRECATED: use `XQSPPhaseFactors.solve()` instead."""
    return XQSPPhaseFactors.solve(P, mode)

xqsp_solve_laurent(P: Polynomial, mode='qsp') -> XQSPPhaseFactors

DEPRECATED: use XQSPPhaseFactors.solve_laurent() instead.

Source code in nlft_qsp/qsp.py

def xqsp_solve_laurent(P: Polynomial, mode='qsp') -> XQSPPhaseFactors:
    r"""DEPRECATED: use `XQSPPhaseFactors.solve_laurent()` instead."""
    return XQSPPhaseFactors.solve_laurent(P, mode)

yqsp_solve(P: Polynomial, mode='qsp') -> YQSPPhaseFactors

DEPRECATED: use YQSPPhaseFactors.solve() instead.

Source code in nlft_qsp/qsp.py

def yqsp_solve(P: Polynomial, mode='qsp') -> YQSPPhaseFactors:
    r"""DEPRECATED: use `YQSPPhaseFactors.solve()` instead."""
    return YQSPPhaseFactors.solve(P, mode)

yqsp_solve_laurent(P: Polynomial, mode='qsp') -> YQSPPhaseFactors

DEPRECATED: use YQSPPhaseFactors.solve_laurent() instead.

Source code in nlft_qsp/qsp.py

def yqsp_solve_laurent(P: Polynomial, mode='qsp') -> YQSPPhaseFactors:
    r"""DEPRECATED: use `YQSPPhaseFactors.solve_laurent()` instead."""
    return YQSPPhaseFactors.solve_laurent(P, mode)