Skip to content

Polynomials

qsp_proc.polynomials

Polynomial types and approximations for QSP (Chebyshev, Laurent, Jacobi-Anger, etc.)

ChebyshevPoly

Polynomial represented in first-kind Chebyshev basis.

For coefficients c_k this class represents P(x) = sum_k c_k T_k(x).

Source code in src\qsp_proc\polynomials\chebyshev.py 
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
class ChebyshevPoly:
    """Polynomial represented in first-kind Chebyshev basis.

    For coefficients ``c_k`` this class represents
    ``P(x) = sum_k c_k T_k(x)``.
    """

    def __init__(self, coeffs: Sequence[complex] | np.ndarray | Chebyshev) -> None:
        """Initialize the polynomial from coefficients or a Chebyshev object."""
        if isinstance(coeffs, Chebyshev):
            self._poly = coeffs
        else:
            self._poly = Chebyshev(np.asarray(coeffs, dtype=np.complex128))

    @property
    def coeffs(self) -> np.ndarray:
        """Return a copy of Chebyshev coefficients."""
        return np.asarray(self._poly.coef, dtype=np.complex128).copy()

    @property
    def degree(self) -> int:
        """Return polynomial degree."""
        return int(self._poly.degree())

    def __call__(self, x: float | np.ndarray) -> complex | np.ndarray:
        """Evaluate the polynomial at ``x``."""
        values = self._poly(x)
        if np.ndim(values) == 0:
            return complex(values)
        return np.asarray(values)

    def parity(self, tol: float = 1e-14) -> int:
        """Return parity by Chebyshev index support.

        Args:
            tol: Numerical zero threshold for coefficients.

        Returns:
            ``0`` if even, ``1`` if odd, and ``-1`` if mixed.
        """
        nz = np.flatnonzero(np.abs(self.coeffs) > tol)
        if nz.size == 0:
            return 0
        if np.all((nz % 2) == 0):
            return 0
        if np.all((nz % 2) == 1):
            return 1
        return -1

    def decompose_eq17(self) -> tuple[ChebyshevPoly, ChebyshevPoly, ChebyshevPoly, ChebyshevPoly]:
        """Decompose ``P`` according to Eq. (17) into four real components.

        Eq. (17) in the paper writes
        ``P(x) = f_R,E(x) + f_R,O(x) + i f_I,E(x) + i f_I,O(x)``.

        Returns:
            Tuple ``(f_R,E, f_R,O, f_I,E, f_I,O)`` represented in Chebyshev basis.
        """
        coeffs = self.coeffs
        real_coeffs = np.real(coeffs)
        imag_coeffs = np.imag(coeffs)
        idx = np.arange(len(coeffs), dtype=np.int64)
        even_mask = (idx % 2) == 0
        odd_mask = ~even_mask

        fr_even = np.zeros_like(real_coeffs, dtype=np.complex128)
        fr_odd = np.zeros_like(real_coeffs, dtype=np.complex128)
        fi_even = np.zeros_like(imag_coeffs, dtype=np.complex128)
        fi_odd = np.zeros_like(imag_coeffs, dtype=np.complex128)

        fr_even[even_mask] = real_coeffs[even_mask]
        fr_odd[odd_mask] = real_coeffs[odd_mask]
        fi_even[even_mask] = imag_coeffs[even_mask]
        fi_odd[odd_mask] = imag_coeffs[odd_mask]

        return (
            ChebyshevPoly(fr_even),
            ChebyshevPoly(fr_odd),
            ChebyshevPoly(fi_even),
            ChebyshevPoly(fi_odd),
        )

    @classmethod
    def from_eq17_components(
        cls,
        f_re_even: ChebyshevPoly,
        f_re_odd: ChebyshevPoly,
        f_im_even: ChebyshevPoly,
        f_im_odd: ChebyshevPoly,
    ) -> ChebyshevPoly:
        """Reconstruct ``P`` from Eq. (17) components."""
        degree = max(f_re_even.degree, f_re_odd.degree, f_im_even.degree, f_im_odd.degree)
        out = np.zeros(degree + 1, dtype=np.complex128)
        parts = (
            (f_re_even.coeffs, 1.0 + 0.0j),
            (f_re_odd.coeffs, 1.0 + 0.0j),
            (f_im_even.coeffs, 1.0j),
            (f_im_odd.coeffs, 1.0j),
        )
        for coeffs, scale in parts:
            out[: len(coeffs)] += scale * coeffs
        return cls(out)

    def linf_norm_on_unit_circle(self, n_samples: int = 4096) -> float:
        """Estimate ``||P||_inf`` on ``U(1)`` using Eq. (16) coordinate map.

        Eq. (16) uses ``z = exp(i*theta)`` and ``x = cos(theta)``, which here is
        evaluated as ``x = (z + z^{-1}) / 2`` for sampled unit-circle points.
        """
        return sampled_linf_on_unit_circle(
            evaluator=lambda z: np.asarray(self._poly((z + 1.0 / z) / 2.0)),
            num_samples=n_samples,
        )

    def subnormalize(self, bound: float = 0.5, n_samples: int = 4096) -> ChebyshevPoly:
        """Return a scaled polynomial with sampled norm at most ``bound``."""
        if bound < 0:
            raise ValueError("bound must be non-negative.")
        current = self.linf_norm_on_unit_circle(n_samples=n_samples)
        if current == 0.0 or current <= bound:
            return ChebyshevPoly(self._poly)
        return ChebyshevPoly(self._poly * (bound / current))

    def to_laurent(self) -> LaurentPoly:
        """Convert first-kind Chebyshev series to Laurent form.

        Uses Appendix A identities with ``x = (z + z^{-1})/2``:
        ``T_0(x)=1`` and ``T_n(x)=(z^n + z^{-n})/2`` for ``n>=1``.
        """
        from .laurent import LaurentPoly

        c = self.coeffs
        degree = len(c) - 1
        dense = np.zeros(2 * degree + 1, dtype=np.complex128)
        center = degree
        dense[center] = c[0]
        for n in range(1, degree + 1):
            half = 0.5 * c[n]
            dense[center - n] += half
            dense[center + n] += half
        return LaurentPoly(coefficients=dense, min_degree=-degree)

coeffs property

Return a copy of Chebyshev coefficients.

degree property

Return polynomial degree.

__init__(coeffs)

Initialize the polynomial from coefficients or a Chebyshev object.

Source code in src\qsp_proc\polynomials\chebyshev.py 
24
25
26
27
28
29
def __init__(self, coeffs: Sequence[complex] | np.ndarray | Chebyshev) -> None:
    """Initialize the polynomial from coefficients or a Chebyshev object."""
    if isinstance(coeffs, Chebyshev):
        self._poly = coeffs
    else:
        self._poly = Chebyshev(np.asarray(coeffs, dtype=np.complex128))

__call__(x)

Evaluate the polynomial at x.

Source code in src\qsp_proc\polynomials\chebyshev.py 
41
42
43
44
45
46
def __call__(self, x: float | np.ndarray) -> complex | np.ndarray:
    """Evaluate the polynomial at ``x``."""
    values = self._poly(x)
    if np.ndim(values) == 0:
        return complex(values)
    return np.asarray(values)

parity(tol=1e-14)

Return parity by Chebyshev index support.

Parameters:

Name Type Description Default
tol float

Numerical zero threshold for coefficients.

 1e-14

Returns:

Type Description
int

0 if even, 1 if odd, and -1 if mixed.
Source code in src\qsp_proc\polynomials\chebyshev.py 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
def parity(self, tol: float = 1e-14) -> int:
    """Return parity by Chebyshev index support.

    Args:
        tol: Numerical zero threshold for coefficients.

    Returns:
        ``0`` if even, ``1`` if odd, and ``-1`` if mixed.
    """
    nz = np.flatnonzero(np.abs(self.coeffs) > tol)
    if nz.size == 0:
        return 0
    if np.all((nz % 2) == 0):
        return 0
    if np.all((nz % 2) == 1):
        return 1
    return -1

decompose_eq17()

Decompose P according to Eq. (17) into four real components.

Eq. (17) in the paper writes P(x) = f_R,E(x) + f_R,O(x) + i f_I,E(x) + i f_I,O(x).

Returns:

Type Description
tuple[ChebyshevPoly, ChebyshevPoly, ChebyshevPoly, ChebyshevPoly]

Tuple (f_R,E, f_R,O, f_I,E, f_I,O) represented in Chebyshev basis.
Source code in src\qsp_proc\polynomials\chebyshev.py 
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def decompose_eq17(self) -> tuple[ChebyshevPoly, ChebyshevPoly, ChebyshevPoly, ChebyshevPoly]:
    """Decompose ``P`` according to Eq. (17) into four real components.

    Eq. (17) in the paper writes
    ``P(x) = f_R,E(x) + f_R,O(x) + i f_I,E(x) + i f_I,O(x)``.

    Returns:
        Tuple ``(f_R,E, f_R,O, f_I,E, f_I,O)`` represented in Chebyshev basis.
    """
    coeffs = self.coeffs
    real_coeffs = np.real(coeffs)
    imag_coeffs = np.imag(coeffs)
    idx = np.arange(len(coeffs), dtype=np.int64)
    even_mask = (idx % 2) == 0
    odd_mask = ~even_mask

    fr_even = np.zeros_like(real_coeffs, dtype=np.complex128)
    fr_odd = np.zeros_like(real_coeffs, dtype=np.complex128)
    fi_even = np.zeros_like(imag_coeffs, dtype=np.complex128)
    fi_odd = np.zeros_like(imag_coeffs, dtype=np.complex128)

    fr_even[even_mask] = real_coeffs[even_mask]
    fr_odd[odd_mask] = real_coeffs[odd_mask]
    fi_even[even_mask] = imag_coeffs[even_mask]
    fi_odd[odd_mask] = imag_coeffs[odd_mask]

    return (
        ChebyshevPoly(fr_even),
        ChebyshevPoly(fr_odd),
        ChebyshevPoly(fi_even),
        ChebyshevPoly(fi_odd),
    )

from_eq17_components(f_re_even, f_re_odd, f_im_even, f_im_odd) classmethod

Reconstruct P from Eq. (17) components.

Source code in src\qsp_proc\polynomials\chebyshev.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
@classmethod
def from_eq17_components(
    cls,
    f_re_even: ChebyshevPoly,
    f_re_odd: ChebyshevPoly,
    f_im_even: ChebyshevPoly,
    f_im_odd: ChebyshevPoly,
) -> ChebyshevPoly:
    """Reconstruct ``P`` from Eq. (17) components."""
    degree = max(f_re_even.degree, f_re_odd.degree, f_im_even.degree, f_im_odd.degree)
    out = np.zeros(degree + 1, dtype=np.complex128)
    parts = (
        (f_re_even.coeffs, 1.0 + 0.0j),
        (f_re_odd.coeffs, 1.0 + 0.0j),
        (f_im_even.coeffs, 1.0j),
        (f_im_odd.coeffs, 1.0j),
    )
    for coeffs, scale in parts:
        out[: len(coeffs)] += scale * coeffs
    return cls(out)

linf_norm_on_unit_circle(n_samples=4096)

Estimate ||P||_inf on U(1) using Eq. (16) coordinate map.

Eq. (16) uses z = exp(i*theta) and x = cos(theta), which here is evaluated as x = (z + z^{-1}) / 2 for sampled unit-circle points.

Source code in src\qsp_proc\polynomials\chebyshev.py 
120
121
122
123
124
125
126
127
128
129
def linf_norm_on_unit_circle(self, n_samples: int = 4096) -> float:
    """Estimate ``||P||_inf`` on ``U(1)`` using Eq. (16) coordinate map.

    Eq. (16) uses ``z = exp(i*theta)`` and ``x = cos(theta)``, which here is
    evaluated as ``x = (z + z^{-1}) / 2`` for sampled unit-circle points.
    """
    return sampled_linf_on_unit_circle(
        evaluator=lambda z: np.asarray(self._poly((z + 1.0 / z) / 2.0)),
        num_samples=n_samples,
    )

subnormalize(bound=0.5, n_samples=4096)

Return a scaled polynomial with sampled norm at most bound.

Source code in src\qsp_proc\polynomials\chebyshev.py 
131
132
133
134
135
136
137
138
def subnormalize(self, bound: float = 0.5, n_samples: int = 4096) -> ChebyshevPoly:
    """Return a scaled polynomial with sampled norm at most ``bound``."""
    if bound < 0:
        raise ValueError("bound must be non-negative.")
    current = self.linf_norm_on_unit_circle(n_samples=n_samples)
    if current == 0.0 or current <= bound:
        return ChebyshevPoly(self._poly)
    return ChebyshevPoly(self._poly * (bound / current))

to_laurent()

Convert first-kind Chebyshev series to Laurent form.

Uses Appendix A identities with x = (z + z^{-1})/2: T_0(x)=1 and T_n(x)=(z^n + z^{-n})/2 for n>=1.

Source code in src\qsp_proc\polynomials\chebyshev.py 
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
def to_laurent(self) -> LaurentPoly:
    """Convert first-kind Chebyshev series to Laurent form.

    Uses Appendix A identities with ``x = (z + z^{-1})/2``:
    ``T_0(x)=1`` and ``T_n(x)=(z^n + z^{-n})/2`` for ``n>=1``.
    """
    from .laurent import LaurentPoly

    c = self.coeffs
    degree = len(c) - 1
    dense = np.zeros(2 * degree + 1, dtype=np.complex128)
    center = degree
    dense[center] = c[0]
    for n in range(1, degree + 1):
        half = 0.5 * c[n]
        dense[center - n] += half
        dense[center + n] += half
    return LaurentPoly(coefficients=dense, min_degree=-degree)

LaurentPoly

Finite Laurent polynomial P(z)=sum_k a_k z^k with dense storage.

Source code in src\qsp_proc\polynomials\laurent.py 
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
class LaurentPoly:
    """Finite Laurent polynomial ``P(z)=sum_k a_k z^k`` with dense storage."""

    def __init__(
        self,
        coefficients: Sequence[complex] | np.ndarray | Mapping[int, complex],
        min_degree: int | None = None,
        *,
        trim_tol: float = 0.0,
    ) -> None:
        """Create a Laurent polynomial from dense or sparse coefficients."""
        if isinstance(coefficients, Mapping):
            if not coefficients:
                self._coeffs = np.array([0.0 + 0.0j], dtype=np.complex128)
                self._min_degree = 0
                return
            min_k = int(min(coefficients))
            max_k = int(max(coefficients))
            dense = np.zeros(max_k - min_k + 1, dtype=np.complex128)
            for exp, value in coefficients.items():
                dense[int(exp) - min_k] = complex(value)
            raw_coeffs = dense
            raw_min = min_k
        else:
            dense = np.asarray(coefficients, dtype=np.complex128)
            if dense.ndim != 1 or dense.size == 0:
                raise ValueError("Dense coefficients must be a non-empty 1D sequence.")
            raw_coeffs = dense.copy()
            raw_min = -((dense.size - 1) // 2) if min_degree is None else int(min_degree)

        self._coeffs, self._min_degree = trim_dense_support(raw_coeffs, raw_min, tol=trim_tol)

    @property
    def coeffs(self) -> np.ndarray:
        """Return dense coefficients as a copy."""
        return self._coeffs.copy()

    @property
    def coefficients(self) -> np.ndarray:
        """Alias for :attr:`coeffs`."""
        return self.coeffs

    @property
    def min_degree(self) -> int:
        """Return the minimum represented exponent."""
        return self._min_degree

    @property
    def max_degree(self) -> int:
        """Return the maximum represented exponent."""
        return self._min_degree + len(self._coeffs) - 1

    @property
    def degree(self) -> int:
        """Return half-width degree ``max(|min_degree|, |max_degree|)``."""
        return max(abs(self.min_degree), abs(self.max_degree))

    def coefficient(self, exponent: int) -> complex:
        """Return coefficient of ``z**exponent`` (or ``0`` outside support)."""
        if exponent < self.min_degree or exponent > self.max_degree:
            return 0.0 + 0.0j
        return complex(self._coeffs[exponent - self.min_degree])

    def __call__(self, z: complex | np.ndarray) -> complex | np.ndarray:
        """Evaluate the Laurent polynomial at scalar or array ``z``."""
        exponents = np.arange(self.min_degree, self.max_degree + 1, dtype=np.int64)
        z_arr = np.asarray(z, dtype=np.complex128)
        values = np.sum(self._coeffs * z_arr[..., None] ** exponents, axis=-1)
        if np.ndim(z_arr) == 0:
            return complex(values)
        return np.asarray(values)

    def evaluate_on_unit_circle(self, theta: float | np.ndarray) -> complex | np.ndarray:
        """Evaluate on ``U(1)`` via ``z = exp(i*theta)``."""
        theta_arr = np.asarray(theta, dtype=np.float64)
        return self(np.exp(1j * theta_arr))

    def linf_norm_on_unit_circle(self, num_samples: int = 4096) -> float:
        """Estimate sampled ``||P||_inf`` on ``U(1)``."""
        return sampled_linf_on_unit_circle(evaluator=self.__call__, num_samples=num_samples)

    def with_bounded_linf_on_unit_circle(
        self, bound: float, num_samples: int = 4096
    ) -> LaurentPoly:
        """Return scaled polynomial with sampled unit-circle norm at most ``bound``."""
        if bound < 0:
            raise ValueError("bound must be non-negative.")
        current = self.linf_norm_on_unit_circle(num_samples=num_samples)
        if current == 0.0 or current <= bound:
            return LaurentPoly(self._coeffs, self._min_degree)
        return LaurentPoly((bound / current) * self._coeffs, self._min_degree)

    def is_reciprocal(self, tol: float = 1e-12) -> bool:
        """Check Appendix A reciprocal condition ``F(z)=F(z^{-1})``."""
        if self.min_degree != -self.max_degree:
            return False
        return bool(np.allclose(self._coeffs, self._coeffs[::-1], atol=tol, rtol=0.0))

    def is_anti_reciprocal(self, tol: float = 1e-12) -> bool:
        """Check Appendix A anti-reciprocal condition ``F(z)=-F(z^{-1})``."""
        if self.min_degree != -self.max_degree:
            return False
        return bool(np.allclose(self._coeffs, -self._coeffs[::-1], atol=tol, rtol=0.0))

    def to_chebyshev(self, tol: float = 1e-12) -> ChebyshevPoly:
        """Convert reciprocal Laurent polynomial to first-kind Chebyshev series."""
        from .chebyshev import ChebyshevPoly

        first_kind, second_kind = self.to_appendix_a_components()
        if second_kind.size > 0 and np.max(np.abs(second_kind)) > tol:
            raise ValueError(
                "Cannot convert to Chebyshev first-kind only: anti-reciprocal component is non-zero."
            )
        return first_kind

    def to_appendix_a_components(self) -> tuple[ChebyshevPoly, np.ndarray]:
        """Return Appendix A Chebyshev decomposition of a Laurent polynomial.

        Any Laurent series can be written as
        ``sum a_n T_n(x) + sqrt(1-x^2) * sum b_n U_n(x)``
        under ``x=(z+z^{-1})/2``.

        Returns:
            A pair ``(a_poly, b_coeffs)`` where:
            - ``a_poly`` is first-kind Chebyshev polynomial ``sum a_n T_n``.
            - ``b_coeffs[n]`` multiplies ``U_n`` in the second-kind series.
        """
        from .chebyshev import ChebyshevPoly

        d = self.degree
        first = np.zeros(d + 1, dtype=np.complex128)
        second = np.zeros(d, dtype=np.complex128)

        first[0] = self.coefficient(0)
        for n in range(1, d + 1):
            cp = self.coefficient(n)
            cm = self.coefficient(-n)
            first[n] = cp + cm
            second[n - 1] = 1j * (cp - cm)

        return ChebyshevPoly(first), second

    @classmethod
    def from_appendix_a_components(
        cls, first_kind: ChebyshevPoly, second_kind_coeffs: Sequence[complex] | np.ndarray
    ) -> LaurentPoly:
        """Build Laurent polynomial from Appendix A Chebyshev components.

        The input represents
        ``sum a_n T_n(x) + sqrt(1-x^2) * sum b_n U_n(x)``.
        """
        a = first_kind.coeffs
        b = np.asarray(second_kind_coeffs, dtype=np.complex128)
        degree = max(len(a) - 1, len(b))
        center = degree
        dense = np.zeros(2 * degree + 1, dtype=np.complex128)

        dense[center] += a[0]
        for n in range(1, len(a)):
            half = 0.5 * a[n]
            dense[center + n] += half
            dense[center - n] += half
        for n in range(1, len(b) + 1):
            coeff = b[n - 1]
            dense[center + n] += -0.5j * coeff
            dense[center - n] += 0.5j * coeff
        return cls(dense, min_degree=-degree)

    @classmethod
    def from_chebyshev(cls, poly: ChebyshevPoly) -> LaurentPoly:
        """Construct Laurent polynomial from a first-kind Chebyshev series."""
        return poly.to_laurent()

    def __mul__(self, other: complex | float | int) -> LaurentPoly:
        """Scalar multiplication ``P * c``."""
        return LaurentPoly(self._coeffs * complex(other), self._min_degree)

    def __rmul__(self, other: complex | float | int) -> LaurentPoly:
        """Scalar multiplication ``c * P``."""
        return self.__mul__(other)

    def __repr__(self) -> str:
        """Return unambiguous string representation."""
        return f"LaurentPoly(coefficients={self._coeffs!r}, min_degree={self._min_degree})"

coeffs property

Return dense coefficients as a copy.

coefficients property

Alias for :attr:coeffs.

min_degree property

Return the minimum represented exponent.

max_degree property

Return the maximum represented exponent.

degree property

Return half-width degree max(|min_degree|, |max_degree|).

__init__(coefficients, min_degree=None, *, trim_tol=0.0)

Create a Laurent polynomial from dense or sparse coefficients.

Source code in src\qsp_proc\polynomials\laurent.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
def __init__(
    self,
    coefficients: Sequence[complex] | np.ndarray | Mapping[int, complex],
    min_degree: int | None = None,
    *,
    trim_tol: float = 0.0,
) -> None:
    """Create a Laurent polynomial from dense or sparse coefficients."""
    if isinstance(coefficients, Mapping):
        if not coefficients:
            self._coeffs = np.array([0.0 + 0.0j], dtype=np.complex128)
            self._min_degree = 0
            return
        min_k = int(min(coefficients))
        max_k = int(max(coefficients))
        dense = np.zeros(max_k - min_k + 1, dtype=np.complex128)
        for exp, value in coefficients.items():
            dense[int(exp) - min_k] = complex(value)
        raw_coeffs = dense
        raw_min = min_k
    else:
        dense = np.asarray(coefficients, dtype=np.complex128)
        if dense.ndim != 1 or dense.size == 0:
            raise ValueError("Dense coefficients must be a non-empty 1D sequence.")
        raw_coeffs = dense.copy()
        raw_min = -((dense.size - 1) // 2) if min_degree is None else int(min_degree)

    self._coeffs, self._min_degree = trim_dense_support(raw_coeffs, raw_min, tol=trim_tol)

coefficient(exponent)

Return coefficient of z**exponent (or 0 outside support).

Source code in src\qsp_proc\polynomials\laurent.py 
73
74
75
76
77
def coefficient(self, exponent: int) -> complex:
    """Return coefficient of ``z**exponent`` (or ``0`` outside support)."""
    if exponent < self.min_degree or exponent > self.max_degree:
        return 0.0 + 0.0j
    return complex(self._coeffs[exponent - self.min_degree])

__call__(z)

Evaluate the Laurent polynomial at scalar or array z.

Source code in src\qsp_proc\polynomials\laurent.py 
79
80
81
82
83
84
85
86
def __call__(self, z: complex | np.ndarray) -> complex | np.ndarray:
    """Evaluate the Laurent polynomial at scalar or array ``z``."""
    exponents = np.arange(self.min_degree, self.max_degree + 1, dtype=np.int64)
    z_arr = np.asarray(z, dtype=np.complex128)
    values = np.sum(self._coeffs * z_arr[..., None] ** exponents, axis=-1)
    if np.ndim(z_arr) == 0:
        return complex(values)
    return np.asarray(values)

evaluate_on_unit_circle(theta)

Evaluate on U(1) via z = exp(i*theta).

Source code in src\qsp_proc\polynomials\laurent.py 
88
89
90
91
def evaluate_on_unit_circle(self, theta: float | np.ndarray) -> complex | np.ndarray:
    """Evaluate on ``U(1)`` via ``z = exp(i*theta)``."""
    theta_arr = np.asarray(theta, dtype=np.float64)
    return self(np.exp(1j * theta_arr))

linf_norm_on_unit_circle(num_samples=4096)

Estimate sampled ||P||_inf on U(1).

Source code in src\qsp_proc\polynomials\laurent.py 
93
94
95
def linf_norm_on_unit_circle(self, num_samples: int = 4096) -> float:
    """Estimate sampled ``||P||_inf`` on ``U(1)``."""
    return sampled_linf_on_unit_circle(evaluator=self.__call__, num_samples=num_samples)

with_bounded_linf_on_unit_circle(bound, num_samples=4096)

Return scaled polynomial with sampled unit-circle norm at most bound.

Source code in src\qsp_proc\polynomials\laurent.py 
 97
 98
 99
100
101
102
103
104
105
106
def with_bounded_linf_on_unit_circle(
    self, bound: float, num_samples: int = 4096
) -> LaurentPoly:
    """Return scaled polynomial with sampled unit-circle norm at most ``bound``."""
    if bound < 0:
        raise ValueError("bound must be non-negative.")
    current = self.linf_norm_on_unit_circle(num_samples=num_samples)
    if current == 0.0 or current <= bound:
        return LaurentPoly(self._coeffs, self._min_degree)
    return LaurentPoly((bound / current) * self._coeffs, self._min_degree)

is_reciprocal(tol=1e-12)

Check Appendix A reciprocal condition F(z)=F(z^{-1}).

Source code in src\qsp_proc\polynomials\laurent.py 
108
109
110
111
112
def is_reciprocal(self, tol: float = 1e-12) -> bool:
    """Check Appendix A reciprocal condition ``F(z)=F(z^{-1})``."""
    if self.min_degree != -self.max_degree:
        return False
    return bool(np.allclose(self._coeffs, self._coeffs[::-1], atol=tol, rtol=0.0))

is_anti_reciprocal(tol=1e-12)

Check Appendix A anti-reciprocal condition F(z)=-F(z^{-1}).

Source code in src\qsp_proc\polynomials\laurent.py 
114
115
116
117
118
def is_anti_reciprocal(self, tol: float = 1e-12) -> bool:
    """Check Appendix A anti-reciprocal condition ``F(z)=-F(z^{-1})``."""
    if self.min_degree != -self.max_degree:
        return False
    return bool(np.allclose(self._coeffs, -self._coeffs[::-1], atol=tol, rtol=0.0))

to_chebyshev(tol=1e-12)

Convert reciprocal Laurent polynomial to first-kind Chebyshev series.

Source code in src\qsp_proc\polynomials\laurent.py 
120
121
122
123
124
125
126
127
128
129
def to_chebyshev(self, tol: float = 1e-12) -> ChebyshevPoly:
    """Convert reciprocal Laurent polynomial to first-kind Chebyshev series."""
    from .chebyshev import ChebyshevPoly

    first_kind, second_kind = self.to_appendix_a_components()
    if second_kind.size > 0 and np.max(np.abs(second_kind)) > tol:
        raise ValueError(
            "Cannot convert to Chebyshev first-kind only: anti-reciprocal component is non-zero."
        )
    return first_kind

to_appendix_a_components()

Return Appendix A Chebyshev decomposition of a Laurent polynomial.

Any Laurent series can be written as sum a_n T_n(x) + sqrt(1-x^2) * sum b_n U_n(x) under x=(z+z^{-1})/2.

Returns:

Type Description
ChebyshevPoly

A pair (a_poly, b_coeffs) where:
ndarray
  • a_poly is first-kind Chebyshev polynomial sum a_n T_n.
tuple[ChebyshevPoly, ndarray]
  • b_coeffs[n] multiplies U_n in the second-kind series.
Source code in src\qsp_proc\polynomials\laurent.py 
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
def to_appendix_a_components(self) -> tuple[ChebyshevPoly, np.ndarray]:
    """Return Appendix A Chebyshev decomposition of a Laurent polynomial.

    Any Laurent series can be written as
    ``sum a_n T_n(x) + sqrt(1-x^2) * sum b_n U_n(x)``
    under ``x=(z+z^{-1})/2``.

    Returns:
        A pair ``(a_poly, b_coeffs)`` where:
        - ``a_poly`` is first-kind Chebyshev polynomial ``sum a_n T_n``.
        - ``b_coeffs[n]`` multiplies ``U_n`` in the second-kind series.
    """
    from .chebyshev import ChebyshevPoly

    d = self.degree
    first = np.zeros(d + 1, dtype=np.complex128)
    second = np.zeros(d, dtype=np.complex128)

    first[0] = self.coefficient(0)
    for n in range(1, d + 1):
        cp = self.coefficient(n)
        cm = self.coefficient(-n)
        first[n] = cp + cm
        second[n - 1] = 1j * (cp - cm)

    return ChebyshevPoly(first), second

from_appendix_a_components(first_kind, second_kind_coeffs) classmethod

Build Laurent polynomial from Appendix A Chebyshev components.

The input represents sum a_n T_n(x) + sqrt(1-x^2) * sum b_n U_n(x).

Source code in src\qsp_proc\polynomials\laurent.py 
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
@classmethod
def from_appendix_a_components(
    cls, first_kind: ChebyshevPoly, second_kind_coeffs: Sequence[complex] | np.ndarray
) -> LaurentPoly:
    """Build Laurent polynomial from Appendix A Chebyshev components.

    The input represents
    ``sum a_n T_n(x) + sqrt(1-x^2) * sum b_n U_n(x)``.
    """
    a = first_kind.coeffs
    b = np.asarray(second_kind_coeffs, dtype=np.complex128)
    degree = max(len(a) - 1, len(b))
    center = degree
    dense = np.zeros(2 * degree + 1, dtype=np.complex128)

    dense[center] += a[0]
    for n in range(1, len(a)):
        half = 0.5 * a[n]
        dense[center + n] += half
        dense[center - n] += half
    for n in range(1, len(b) + 1):
        coeff = b[n - 1]
        dense[center + n] += -0.5j * coeff
        dense[center - n] += 0.5j * coeff
    return cls(dense, min_degree=-degree)

from_chebyshev(poly) classmethod

Construct Laurent polynomial from a first-kind Chebyshev series.

Source code in src\qsp_proc\polynomials\laurent.py 
184
185
186
187
@classmethod
def from_chebyshev(cls, poly: ChebyshevPoly) -> LaurentPoly:
    """Construct Laurent polynomial from a first-kind Chebyshev series."""
    return poly.to_laurent()

__mul__(other)

Scalar multiplication P * c.

Source code in src\qsp_proc\polynomials\laurent.py 
189
190
191
def __mul__(self, other: complex | float | int) -> LaurentPoly:
    """Scalar multiplication ``P * c``."""
    return LaurentPoly(self._coeffs * complex(other), self._min_degree)

__rmul__(other)

Scalar multiplication c * P.

Source code in src\qsp_proc\polynomials\laurent.py 
193
194
195
def __rmul__(self, other: complex | float | int) -> LaurentPoly:
    """Scalar multiplication ``c * P``."""
    return self.__mul__(other)

__repr__()

Return unambiguous string representation.

Source code in src\qsp_proc\polynomials\laurent.py 
197
198
199
def __repr__(self) -> str:
    """Return unambiguous string representation."""
    return f"LaurentPoly(coefficients={self._coeffs!r}, min_degree={self._min_degree})"

hamiltonian_sim_degree_bound(tau, epsilon)

Return the Appendix B.2 piecewise truncation bound r_tilde(tau, epsilon).

This implements Eq. (73) from the paper for tau >= 0 and 0 < epsilon < 1. For negative tau, the bound is computed with |tau|.

Source code in src\qsp_proc\polynomials\approximations.py 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
def hamiltonian_sim_degree_bound(tau: float, epsilon: float) -> int:
    """Return the Appendix B.2 piecewise truncation bound ``r_tilde(tau, epsilon)``.

    This implements Eq. (73) from the paper for ``tau >= 0`` and ``0 < epsilon < 1``.
    For negative ``tau``, the bound is computed with ``|tau|``.
    """
    if not np.isfinite(tau):
        raise ValueError("tau must be finite.")
    if not (0.0 < epsilon < 1.0):
        raise ValueError("epsilon must satisfy 0 < epsilon < 1.")

    abs_tau = float(abs(tau))
    if abs_tau == 0.0:
        return 0

    log_term = math.log(1.0 / epsilon)
    threshold = log_term / math.e

    if abs_tau >= threshold:
        return int(math.ceil(math.e * abs_tau + log_term))

    denominator = math.log(math.e + (log_term / abs_tau))
    return int(math.ceil((4.0 * log_term) / denominator))

jacobi_anger_laurent(tau, epsilon, *, max_degree=None, subnormalization=math.sqrt(0.5))

Build truncated Jacobi-Anger Laurent target for Hamiltonian simulation.

The target corresponds to exp(i * tau * x) with x = cos(theta) and z = exp(i * theta), truncated to degree r where r is either max_degree or the Appendix B.2 Eq. (73) bound.

Coefficients are generated from the Fourier/Jacobi-Anger identity exp(i * tau * cos(theta)) = sum_k i^k J_k(tau) z^k for k in [-r, r].

Parameters:

Name Type Description Default
tau float

Simulation time parameter.

 required
epsilon float

Approximation precision for selecting truncation degree.

 required
max_degree int | None

Explicit truncation degree override.

 None
subnormalization float

Optional multiplicative scaling (paper uses sqrt(1/2)).

 sqrt(0.5)

Returns:

Type Description
LaurentPoly

Truncated Laurent polynomial target.
Source code in src\qsp_proc\polynomials\approximations.py 
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
def jacobi_anger_laurent(
    tau: float,
    epsilon: float,
    *,
    max_degree: int | None = None,
    subnormalization: float = math.sqrt(0.5),
) -> LaurentPoly:
    """Build truncated Jacobi-Anger Laurent target for Hamiltonian simulation.

    The target corresponds to ``exp(i * tau * x)`` with ``x = cos(theta)`` and
    ``z = exp(i * theta)``, truncated to degree ``r`` where ``r`` is either
    ``max_degree`` or the Appendix B.2 Eq. (73) bound.

    Coefficients are generated from the Fourier/Jacobi-Anger identity
    ``exp(i * tau * cos(theta)) = sum_k i^k J_k(tau) z^k`` for ``k in [-r, r]``.

    Args:
        tau: Simulation time parameter.
        epsilon: Approximation precision for selecting truncation degree.
        max_degree: Explicit truncation degree override.
        subnormalization: Optional multiplicative scaling (paper uses ``sqrt(1/2)``).

    Returns:
        Truncated Laurent polynomial target.
    """
    if not np.isfinite(tau):
        raise ValueError("tau must be finite.")
    if not (0.0 < epsilon < 1.0):
        raise ValueError("epsilon must satisfy 0 < epsilon < 1.")
    if max_degree is not None and max_degree < 0:
        raise ValueError("max_degree must be non-negative when provided.")
    if subnormalization <= 0.0:
        raise ValueError("subnormalization must be positive.")

    try:
        from scipy.special import jv
    except Exception as exc:  # pragma: no cover
        raise ImportError("jacobi_anger_laurent requires scipy.special.jv.") from exc

    degree = hamiltonian_sim_degree_bound(tau, epsilon) if max_degree is None else int(max_degree)
    coeffs: dict[int, complex] = {}
    for k in range(-degree, degree + 1):
        coeffs[k] = complex((1j**k) * jv(k, tau))

    dense = LaurentPoly(coeffs)
    if subnormalization != 1.0:
        return LaurentPoly(dense.coeffs * subnormalization, dense.min_degree)
    return dense

random_complex_target_laurent(degree, *, nz=None, decay_rate=1.5, norm_bound=0.5, rng=None)

Generate Appendix B.1 style sparse random complex target.

The target is sampled in theta-form: P(theta) = sum_j a_j cos(j theta) + i sum_j b_j sin(j theta), then converted to Laurent coefficients.

Parameters:

Name Type Description Default
degree int

Maximum Fourier/Chebyshev order.

 required
nz int | None

Number of non-zero coefficients across {a_j} and {b_j}. Defaults to max(5, degree // 15) as in Appendix B.1.

 None
decay_rate float

Magnitude decay exponent for higher-order terms.

 1.5
norm_bound float

Post-normalization bound for sampled ||P||_inf on U(1).

 0.5
rng Generator | None

Optional NumPy random generator.

 None

Returns:

Type Description
LaurentPoly

Sparse random Laurent polynomial with sampled norm at most norm_bound.
Source code in src\qsp_proc\polynomials\approximations.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def random_complex_target_laurent(
    degree: int,
    *,
    nz: int | None = None,
    decay_rate: float = 1.5,
    norm_bound: float = 0.5,
    rng: np.random.Generator | None = None,
) -> LaurentPoly:
    """Generate Appendix B.1 style sparse random complex target.

    The target is sampled in theta-form:
    ``P(theta) = sum_j a_j cos(j theta) + i sum_j b_j sin(j theta)``,
    then converted to Laurent coefficients.

    Args:
        degree: Maximum Fourier/Chebyshev order.
        nz: Number of non-zero coefficients across ``{a_j}`` and ``{b_j}``.
            Defaults to ``max(5, degree // 15)`` as in Appendix B.1.
        decay_rate: Magnitude decay exponent for higher-order terms.
        norm_bound: Post-normalization bound for sampled ``||P||_inf`` on ``U(1)``.
        rng: Optional NumPy random generator.

    Returns:
        Sparse random Laurent polynomial with sampled norm at most ``norm_bound``.
    """
    if degree < 0:
        raise ValueError("degree must be non-negative.")
    if decay_rate <= 0.0:
        raise ValueError("decay_rate must be positive.")
    if norm_bound < 0.0:
        raise ValueError("norm_bound must be non-negative.")

    rng = np.random.default_rng() if rng is None else rng
    default_nz = max(5, degree // 15) if degree > 0 else 1
    total_candidates = (degree + 1) + degree  # a_0..a_degree and b_1..b_degree
    nnz = default_nz if nz is None else nz
    if nnz < 0 or nnz > total_candidates:
        raise ValueError(f"nz must satisfy 0 <= nz <= {total_candidates}.")

    a = np.zeros(degree + 1, dtype=np.float64)
    b = np.zeros(degree + 1, dtype=np.float64)
    candidates = [("a", j) for j in range(degree + 1)] + [("b", j) for j in range(1, degree + 1)]
    if nnz > 0:
        chosen = rng.choice(len(candidates), size=nnz, replace=False)
        for idx in chosen:
            label, order = candidates[int(idx)]
            scale = float((order + 1) ** (-decay_rate))
            value = float(rng.uniform(-1.0, 1.0) * scale)
            if label == "a":
                a[order] = value
            else:
                b[order] = value

    coeffs: dict[int, complex] = {0: complex(a[0])}
    for j in range(1, degree + 1):
        c_pos = 0.5 * (a[j] + b[j])
        c_neg = 0.5 * (a[j] - b[j])
        if c_pos != 0.0:
            coeffs[j] = complex(c_pos)
        if c_neg != 0.0:
            coeffs[-j] = complex(c_neg)

    return LaurentPoly(coeffs).with_bounded_linf_on_unit_circle(bound=norm_bound)

sampled_linf_on_unit_circle(evaluator, num_samples)

Estimate ||f||_inf on U(1) via uniform sampling.

Parameters:

Name Type Description Default
evaluator Callable[[ndarray], ndarray]

Vectorized callable that accepts sampled z values.

 required
num_samples int

Number of points used for sampling.

 required

Returns:

Type Description
float

Estimated sup norm on sampled points.
Source code in src\qsp_proc\polynomials\utils.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def sampled_linf_on_unit_circle(
    evaluator: Callable[[np.ndarray], np.ndarray], num_samples: int
) -> float:
    """Estimate ``||f||_inf`` on ``U(1)`` via uniform sampling.

    Args:
        evaluator: Vectorized callable that accepts sampled ``z`` values.
        num_samples: Number of points used for sampling.

    Returns:
        Estimated sup norm on sampled points.
    """
    z = unit_circle_points(num_samples)
    values = evaluator(z)
    return float(np.max(np.abs(values)))

trim_dense_support(coeffs, min_degree, *, tol)

Trim near-zero boundary entries in a dense Laurent support.

Parameters:

Name Type Description Default
coeffs ndarray

Dense coefficient array where coeffs[i] maps to exponent min_degree + i.

 required
min_degree int

Smallest exponent represented by coeffs.

 required
tol float

Trimming tolerance.

 required

Returns:

Type Description
tuple[ndarray, int]

A pair (trimmed_coeffs, trimmed_min_degree).
Source code in src\qsp_proc\polynomials\utils.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
def trim_dense_support(
    coeffs: np.ndarray, min_degree: int, *, tol: float
) -> tuple[np.ndarray, int]:
    """Trim near-zero boundary entries in a dense Laurent support.

    Args:
        coeffs: Dense coefficient array where ``coeffs[i]`` maps to exponent
            ``min_degree + i``.
        min_degree: Smallest exponent represented by ``coeffs``.
        tol: Trimming tolerance.

    Returns:
        A pair ``(trimmed_coeffs, trimmed_min_degree)``.
    """
    if coeffs.size == 0:
        return np.array([0.0 + 0.0j], dtype=np.complex128), 0

    mask = np.abs(coeffs) > tol
    if not np.any(mask):
        return np.array([0.0 + 0.0j], dtype=np.complex128), 0

    left = int(np.argmax(mask))
    right = int(len(mask) - np.argmax(mask[::-1]))
    return coeffs[left:right].copy(), min_degree + left

unit_circle_points(num_samples)

Return uniformly sampled points on the unit circle.

Parameters:

Name Type Description Default
num_samples int

Number of points to sample on U(1).

 required

Returns:

Type Description
ndarray

Complex array z_k = exp(i * theta_k) with theta_k uniformly
ndarray

spaced in [0, 2*pi).

Raises:

Type Description
ValueError

If num_samples is not positive.
Source code in src\qsp_proc\polynomials\utils.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def unit_circle_points(num_samples: int) -> np.ndarray:
    """Return uniformly sampled points on the unit circle.

    Args:
        num_samples: Number of points to sample on ``U(1)``.

    Returns:
        Complex array ``z_k = exp(i * theta_k)`` with ``theta_k`` uniformly
        spaced in ``[0, 2*pi)``.

    Raises:
        ValueError: If ``num_samples`` is not positive.
    """
    if num_samples <= 0:
        raise ValueError("num_samples must be positive.")
    theta = np.linspace(0.0, 2.0 * np.pi, num_samples, endpoint=False, dtype=np.float64)
    return np.exp(1j * theta)

qsp_proc.polynomials.chebyshev

Chebyshev polynomial utilities for QSP

ChebyshevPoly

Polynomial represented in first-kind Chebyshev basis.

For coefficients c_k this class represents P(x) = sum_k c_k T_k(x).

Source code in src\qsp_proc\polynomials\chebyshev.py 
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
class ChebyshevPoly:
    """Polynomial represented in first-kind Chebyshev basis.

    For coefficients ``c_k`` this class represents
    ``P(x) = sum_k c_k T_k(x)``.
    """

    def __init__(self, coeffs: Sequence[complex] | np.ndarray | Chebyshev) -> None:
        """Initialize the polynomial from coefficients or a Chebyshev object."""
        if isinstance(coeffs, Chebyshev):
            self._poly = coeffs
        else:
            self._poly = Chebyshev(np.asarray(coeffs, dtype=np.complex128))

    @property
    def coeffs(self) -> np.ndarray:
        """Return a copy of Chebyshev coefficients."""
        return np.asarray(self._poly.coef, dtype=np.complex128).copy()

    @property
    def degree(self) -> int:
        """Return polynomial degree."""
        return int(self._poly.degree())

    def __call__(self, x: float | np.ndarray) -> complex | np.ndarray:
        """Evaluate the polynomial at ``x``."""
        values = self._poly(x)
        if np.ndim(values) == 0:
            return complex(values)
        return np.asarray(values)

    def parity(self, tol: float = 1e-14) -> int:
        """Return parity by Chebyshev index support.

        Args:
            tol: Numerical zero threshold for coefficients.

        Returns:
            ``0`` if even, ``1`` if odd, and ``-1`` if mixed.
        """
        nz = np.flatnonzero(np.abs(self.coeffs) > tol)
        if nz.size == 0:
            return 0
        if np.all((nz % 2) == 0):
            return 0
        if np.all((nz % 2) == 1):
            return 1
        return -1

    def decompose_eq17(self) -> tuple[ChebyshevPoly, ChebyshevPoly, ChebyshevPoly, ChebyshevPoly]:
        """Decompose ``P`` according to Eq. (17) into four real components.

        Eq. (17) in the paper writes
        ``P(x) = f_R,E(x) + f_R,O(x) + i f_I,E(x) + i f_I,O(x)``.

        Returns:
            Tuple ``(f_R,E, f_R,O, f_I,E, f_I,O)`` represented in Chebyshev basis.
        """
        coeffs = self.coeffs
        real_coeffs = np.real(coeffs)
        imag_coeffs = np.imag(coeffs)
        idx = np.arange(len(coeffs), dtype=np.int64)
        even_mask = (idx % 2) == 0
        odd_mask = ~even_mask

        fr_even = np.zeros_like(real_coeffs, dtype=np.complex128)
        fr_odd = np.zeros_like(real_coeffs, dtype=np.complex128)
        fi_even = np.zeros_like(imag_coeffs, dtype=np.complex128)
        fi_odd = np.zeros_like(imag_coeffs, dtype=np.complex128)

        fr_even[even_mask] = real_coeffs[even_mask]
        fr_odd[odd_mask] = real_coeffs[odd_mask]
        fi_even[even_mask] = imag_coeffs[even_mask]
        fi_odd[odd_mask] = imag_coeffs[odd_mask]

        return (
            ChebyshevPoly(fr_even),
            ChebyshevPoly(fr_odd),
            ChebyshevPoly(fi_even),
            ChebyshevPoly(fi_odd),
        )

    @classmethod
    def from_eq17_components(
        cls,
        f_re_even: ChebyshevPoly,
        f_re_odd: ChebyshevPoly,
        f_im_even: ChebyshevPoly,
        f_im_odd: ChebyshevPoly,
    ) -> ChebyshevPoly:
        """Reconstruct ``P`` from Eq. (17) components."""
        degree = max(f_re_even.degree, f_re_odd.degree, f_im_even.degree, f_im_odd.degree)
        out = np.zeros(degree + 1, dtype=np.complex128)
        parts = (
            (f_re_even.coeffs, 1.0 + 0.0j),
            (f_re_odd.coeffs, 1.0 + 0.0j),
            (f_im_even.coeffs, 1.0j),
            (f_im_odd.coeffs, 1.0j),
        )
        for coeffs, scale in parts:
            out[: len(coeffs)] += scale * coeffs
        return cls(out)

    def linf_norm_on_unit_circle(self, n_samples: int = 4096) -> float:
        """Estimate ``||P||_inf`` on ``U(1)`` using Eq. (16) coordinate map.

        Eq. (16) uses ``z = exp(i*theta)`` and ``x = cos(theta)``, which here is
        evaluated as ``x = (z + z^{-1}) / 2`` for sampled unit-circle points.
        """
        return sampled_linf_on_unit_circle(
            evaluator=lambda z: np.asarray(self._poly((z + 1.0 / z) / 2.0)),
            num_samples=n_samples,
        )

    def subnormalize(self, bound: float = 0.5, n_samples: int = 4096) -> ChebyshevPoly:
        """Return a scaled polynomial with sampled norm at most ``bound``."""
        if bound < 0:
            raise ValueError("bound must be non-negative.")
        current = self.linf_norm_on_unit_circle(n_samples=n_samples)
        if current == 0.0 or current <= bound:
            return ChebyshevPoly(self._poly)
        return ChebyshevPoly(self._poly * (bound / current))

    def to_laurent(self) -> LaurentPoly:
        """Convert first-kind Chebyshev series to Laurent form.

        Uses Appendix A identities with ``x = (z + z^{-1})/2``:
        ``T_0(x)=1`` and ``T_n(x)=(z^n + z^{-n})/2`` for ``n>=1``.
        """
        from .laurent import LaurentPoly

        c = self.coeffs
        degree = len(c) - 1
        dense = np.zeros(2 * degree + 1, dtype=np.complex128)
        center = degree
        dense[center] = c[0]
        for n in range(1, degree + 1):
            half = 0.5 * c[n]
            dense[center - n] += half
            dense[center + n] += half
        return LaurentPoly(coefficients=dense, min_degree=-degree)

coeffs property

Return a copy of Chebyshev coefficients.

degree property

Return polynomial degree.

__init__(coeffs)

Initialize the polynomial from coefficients or a Chebyshev object.

Source code in src\qsp_proc\polynomials\chebyshev.py 
24
25
26
27
28
29
def __init__(self, coeffs: Sequence[complex] | np.ndarray | Chebyshev) -> None:
    """Initialize the polynomial from coefficients or a Chebyshev object."""
    if isinstance(coeffs, Chebyshev):
        self._poly = coeffs
    else:
        self._poly = Chebyshev(np.asarray(coeffs, dtype=np.complex128))

__call__(x)

Evaluate the polynomial at x.

Source code in src\qsp_proc\polynomials\chebyshev.py 
41
42
43
44
45
46
def __call__(self, x: float | np.ndarray) -> complex | np.ndarray:
    """Evaluate the polynomial at ``x``."""
    values = self._poly(x)
    if np.ndim(values) == 0:
        return complex(values)
    return np.asarray(values)

parity(tol=1e-14)

Return parity by Chebyshev index support.

Parameters:

Name Type Description Default
tol float

Numerical zero threshold for coefficients.

 1e-14

Returns:

Type Description
int

0 if even, 1 if odd, and -1 if mixed.
Source code in src\qsp_proc\polynomials\chebyshev.py 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
def parity(self, tol: float = 1e-14) -> int:
    """Return parity by Chebyshev index support.

    Args:
        tol: Numerical zero threshold for coefficients.

    Returns:
        ``0`` if even, ``1`` if odd, and ``-1`` if mixed.
    """
    nz = np.flatnonzero(np.abs(self.coeffs) > tol)
    if nz.size == 0:
        return 0
    if np.all((nz % 2) == 0):
        return 0
    if np.all((nz % 2) == 1):
        return 1
    return -1

decompose_eq17()

Decompose P according to Eq. (17) into four real components.

Eq. (17) in the paper writes P(x) = f_R,E(x) + f_R,O(x) + i f_I,E(x) + i f_I,O(x).

Returns:

Type Description
tuple[ChebyshevPoly, ChebyshevPoly, ChebyshevPoly, ChebyshevPoly]

Tuple (f_R,E, f_R,O, f_I,E, f_I,O) represented in Chebyshev basis.
Source code in src\qsp_proc\polynomials\chebyshev.py 
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def decompose_eq17(self) -> tuple[ChebyshevPoly, ChebyshevPoly, ChebyshevPoly, ChebyshevPoly]:
    """Decompose ``P`` according to Eq. (17) into four real components.

    Eq. (17) in the paper writes
    ``P(x) = f_R,E(x) + f_R,O(x) + i f_I,E(x) + i f_I,O(x)``.

    Returns:
        Tuple ``(f_R,E, f_R,O, f_I,E, f_I,O)`` represented in Chebyshev basis.
    """
    coeffs = self.coeffs
    real_coeffs = np.real(coeffs)
    imag_coeffs = np.imag(coeffs)
    idx = np.arange(len(coeffs), dtype=np.int64)
    even_mask = (idx % 2) == 0
    odd_mask = ~even_mask

    fr_even = np.zeros_like(real_coeffs, dtype=np.complex128)
    fr_odd = np.zeros_like(real_coeffs, dtype=np.complex128)
    fi_even = np.zeros_like(imag_coeffs, dtype=np.complex128)
    fi_odd = np.zeros_like(imag_coeffs, dtype=np.complex128)

    fr_even[even_mask] = real_coeffs[even_mask]
    fr_odd[odd_mask] = real_coeffs[odd_mask]
    fi_even[even_mask] = imag_coeffs[even_mask]
    fi_odd[odd_mask] = imag_coeffs[odd_mask]

    return (
        ChebyshevPoly(fr_even),
        ChebyshevPoly(fr_odd),
        ChebyshevPoly(fi_even),
        ChebyshevPoly(fi_odd),
    )

from_eq17_components(f_re_even, f_re_odd, f_im_even, f_im_odd) classmethod

Reconstruct P from Eq. (17) components.

Source code in src\qsp_proc\polynomials\chebyshev.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
@classmethod
def from_eq17_components(
    cls,
    f_re_even: ChebyshevPoly,
    f_re_odd: ChebyshevPoly,
    f_im_even: ChebyshevPoly,
    f_im_odd: ChebyshevPoly,
) -> ChebyshevPoly:
    """Reconstruct ``P`` from Eq. (17) components."""
    degree = max(f_re_even.degree, f_re_odd.degree, f_im_even.degree, f_im_odd.degree)
    out = np.zeros(degree + 1, dtype=np.complex128)
    parts = (
        (f_re_even.coeffs, 1.0 + 0.0j),
        (f_re_odd.coeffs, 1.0 + 0.0j),
        (f_im_even.coeffs, 1.0j),
        (f_im_odd.coeffs, 1.0j),
    )
    for coeffs, scale in parts:
        out[: len(coeffs)] += scale * coeffs
    return cls(out)

linf_norm_on_unit_circle(n_samples=4096)

Estimate ||P||_inf on U(1) using Eq. (16) coordinate map.

Eq. (16) uses z = exp(i*theta) and x = cos(theta), which here is evaluated as x = (z + z^{-1}) / 2 for sampled unit-circle points.

Source code in src\qsp_proc\polynomials\chebyshev.py 
120
121
122
123
124
125
126
127
128
129
def linf_norm_on_unit_circle(self, n_samples: int = 4096) -> float:
    """Estimate ``||P||_inf`` on ``U(1)`` using Eq. (16) coordinate map.

    Eq. (16) uses ``z = exp(i*theta)`` and ``x = cos(theta)``, which here is
    evaluated as ``x = (z + z^{-1}) / 2`` for sampled unit-circle points.
    """
    return sampled_linf_on_unit_circle(
        evaluator=lambda z: np.asarray(self._poly((z + 1.0 / z) / 2.0)),
        num_samples=n_samples,
    )

subnormalize(bound=0.5, n_samples=4096)

Return a scaled polynomial with sampled norm at most bound.

Source code in src\qsp_proc\polynomials\chebyshev.py 
131
132
133
134
135
136
137
138
def subnormalize(self, bound: float = 0.5, n_samples: int = 4096) -> ChebyshevPoly:
    """Return a scaled polynomial with sampled norm at most ``bound``."""
    if bound < 0:
        raise ValueError("bound must be non-negative.")
    current = self.linf_norm_on_unit_circle(n_samples=n_samples)
    if current == 0.0 or current <= bound:
        return ChebyshevPoly(self._poly)
    return ChebyshevPoly(self._poly * (bound / current))

to_laurent()

Convert first-kind Chebyshev series to Laurent form.

Uses Appendix A identities with x = (z + z^{-1})/2: T_0(x)=1 and T_n(x)=(z^n + z^{-n})/2 for n>=1.

Source code in src\qsp_proc\polynomials\chebyshev.py 
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
def to_laurent(self) -> LaurentPoly:
    """Convert first-kind Chebyshev series to Laurent form.

    Uses Appendix A identities with ``x = (z + z^{-1})/2``:
    ``T_0(x)=1`` and ``T_n(x)=(z^n + z^{-n})/2`` for ``n>=1``.
    """
    from .laurent import LaurentPoly

    c = self.coeffs
    degree = len(c) - 1
    dense = np.zeros(2 * degree + 1, dtype=np.complex128)
    center = degree
    dense[center] = c[0]
    for n in range(1, degree + 1):
        half = 0.5 * c[n]
        dense[center - n] += half
        dense[center + n] += half
    return LaurentPoly(coefficients=dense, min_degree=-degree)

qsp_proc.polynomials.laurent

Laurent polynomial utilities.

LaurentPoly

Finite Laurent polynomial P(z)=sum_k a_k z^k with dense storage.

Source code in src\qsp_proc\polynomials\laurent.py 
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
class LaurentPoly:
    """Finite Laurent polynomial ``P(z)=sum_k a_k z^k`` with dense storage."""

    def __init__(
        self,
        coefficients: Sequence[complex] | np.ndarray | Mapping[int, complex],
        min_degree: int | None = None,
        *,
        trim_tol: float = 0.0,
    ) -> None:
        """Create a Laurent polynomial from dense or sparse coefficients."""
        if isinstance(coefficients, Mapping):
            if not coefficients:
                self._coeffs = np.array([0.0 + 0.0j], dtype=np.complex128)
                self._min_degree = 0
                return
            min_k = int(min(coefficients))
            max_k = int(max(coefficients))
            dense = np.zeros(max_k - min_k + 1, dtype=np.complex128)
            for exp, value in coefficients.items():
                dense[int(exp) - min_k] = complex(value)
            raw_coeffs = dense
            raw_min = min_k
        else:
            dense = np.asarray(coefficients, dtype=np.complex128)
            if dense.ndim != 1 or dense.size == 0:
                raise ValueError("Dense coefficients must be a non-empty 1D sequence.")
            raw_coeffs = dense.copy()
            raw_min = -((dense.size - 1) // 2) if min_degree is None else int(min_degree)

        self._coeffs, self._min_degree = trim_dense_support(raw_coeffs, raw_min, tol=trim_tol)

    @property
    def coeffs(self) -> np.ndarray:
        """Return dense coefficients as a copy."""
        return self._coeffs.copy()

    @property
    def coefficients(self) -> np.ndarray:
        """Alias for :attr:`coeffs`."""
        return self.coeffs

    @property
    def min_degree(self) -> int:
        """Return the minimum represented exponent."""
        return self._min_degree

    @property
    def max_degree(self) -> int:
        """Return the maximum represented exponent."""
        return self._min_degree + len(self._coeffs) - 1

    @property
    def degree(self) -> int:
        """Return half-width degree ``max(|min_degree|, |max_degree|)``."""
        return max(abs(self.min_degree), abs(self.max_degree))

    def coefficient(self, exponent: int) -> complex:
        """Return coefficient of ``z**exponent`` (or ``0`` outside support)."""
        if exponent < self.min_degree or exponent > self.max_degree:
            return 0.0 + 0.0j
        return complex(self._coeffs[exponent - self.min_degree])

    def __call__(self, z: complex | np.ndarray) -> complex | np.ndarray:
        """Evaluate the Laurent polynomial at scalar or array ``z``."""
        exponents = np.arange(self.min_degree, self.max_degree + 1, dtype=np.int64)
        z_arr = np.asarray(z, dtype=np.complex128)
        values = np.sum(self._coeffs * z_arr[..., None] ** exponents, axis=-1)
        if np.ndim(z_arr) == 0:
            return complex(values)
        return np.asarray(values)

    def evaluate_on_unit_circle(self, theta: float | np.ndarray) -> complex | np.ndarray:
        """Evaluate on ``U(1)`` via ``z = exp(i*theta)``."""
        theta_arr = np.asarray(theta, dtype=np.float64)
        return self(np.exp(1j * theta_arr))

    def linf_norm_on_unit_circle(self, num_samples: int = 4096) -> float:
        """Estimate sampled ``||P||_inf`` on ``U(1)``."""
        return sampled_linf_on_unit_circle(evaluator=self.__call__, num_samples=num_samples)

    def with_bounded_linf_on_unit_circle(
        self, bound: float, num_samples: int = 4096
    ) -> LaurentPoly:
        """Return scaled polynomial with sampled unit-circle norm at most ``bound``."""
        if bound < 0:
            raise ValueError("bound must be non-negative.")
        current = self.linf_norm_on_unit_circle(num_samples=num_samples)
        if current == 0.0 or current <= bound:
            return LaurentPoly(self._coeffs, self._min_degree)
        return LaurentPoly((bound / current) * self._coeffs, self._min_degree)

    def is_reciprocal(self, tol: float = 1e-12) -> bool:
        """Check Appendix A reciprocal condition ``F(z)=F(z^{-1})``."""
        if self.min_degree != -self.max_degree:
            return False
        return bool(np.allclose(self._coeffs, self._coeffs[::-1], atol=tol, rtol=0.0))

    def is_anti_reciprocal(self, tol: float = 1e-12) -> bool:
        """Check Appendix A anti-reciprocal condition ``F(z)=-F(z^{-1})``."""
        if self.min_degree != -self.max_degree:
            return False
        return bool(np.allclose(self._coeffs, -self._coeffs[::-1], atol=tol, rtol=0.0))

    def to_chebyshev(self, tol: float = 1e-12) -> ChebyshevPoly:
        """Convert reciprocal Laurent polynomial to first-kind Chebyshev series."""
        from .chebyshev import ChebyshevPoly

        first_kind, second_kind = self.to_appendix_a_components()
        if second_kind.size > 0 and np.max(np.abs(second_kind)) > tol:
            raise ValueError(
                "Cannot convert to Chebyshev first-kind only: anti-reciprocal component is non-zero."
            )
        return first_kind

    def to_appendix_a_components(self) -> tuple[ChebyshevPoly, np.ndarray]:
        """Return Appendix A Chebyshev decomposition of a Laurent polynomial.

        Any Laurent series can be written as
        ``sum a_n T_n(x) + sqrt(1-x^2) * sum b_n U_n(x)``
        under ``x=(z+z^{-1})/2``.

        Returns:
            A pair ``(a_poly, b_coeffs)`` where:
            - ``a_poly`` is first-kind Chebyshev polynomial ``sum a_n T_n``.
            - ``b_coeffs[n]`` multiplies ``U_n`` in the second-kind series.
        """
        from .chebyshev import ChebyshevPoly

        d = self.degree
        first = np.zeros(d + 1, dtype=np.complex128)
        second = np.zeros(d, dtype=np.complex128)

        first[0] = self.coefficient(0)
        for n in range(1, d + 1):
            cp = self.coefficient(n)
            cm = self.coefficient(-n)
            first[n] = cp + cm
            second[n - 1] = 1j * (cp - cm)

        return ChebyshevPoly(first), second

    @classmethod
    def from_appendix_a_components(
        cls, first_kind: ChebyshevPoly, second_kind_coeffs: Sequence[complex] | np.ndarray
    ) -> LaurentPoly:
        """Build Laurent polynomial from Appendix A Chebyshev components.

        The input represents
        ``sum a_n T_n(x) + sqrt(1-x^2) * sum b_n U_n(x)``.
        """
        a = first_kind.coeffs
        b = np.asarray(second_kind_coeffs, dtype=np.complex128)
        degree = max(len(a) - 1, len(b))
        center = degree
        dense = np.zeros(2 * degree + 1, dtype=np.complex128)

        dense[center] += a[0]
        for n in range(1, len(a)):
            half = 0.5 * a[n]
            dense[center + n] += half
            dense[center - n] += half
        for n in range(1, len(b) + 1):
            coeff = b[n - 1]
            dense[center + n] += -0.5j * coeff
            dense[center - n] += 0.5j * coeff
        return cls(dense, min_degree=-degree)

    @classmethod
    def from_chebyshev(cls, poly: ChebyshevPoly) -> LaurentPoly:
        """Construct Laurent polynomial from a first-kind Chebyshev series."""
        return poly.to_laurent()

    def __mul__(self, other: complex | float | int) -> LaurentPoly:
        """Scalar multiplication ``P * c``."""
        return LaurentPoly(self._coeffs * complex(other), self._min_degree)

    def __rmul__(self, other: complex | float | int) -> LaurentPoly:
        """Scalar multiplication ``c * P``."""
        return self.__mul__(other)

    def __repr__(self) -> str:
        """Return unambiguous string representation."""
        return f"LaurentPoly(coefficients={self._coeffs!r}, min_degree={self._min_degree})"

coeffs property

Return dense coefficients as a copy.

coefficients property

Alias for :attr:coeffs.

min_degree property

Return the minimum represented exponent.

max_degree property

Return the maximum represented exponent.

degree property

Return half-width degree max(|min_degree|, |max_degree|).

__init__(coefficients, min_degree=None, *, trim_tol=0.0)

Create a Laurent polynomial from dense or sparse coefficients.

Source code in src\qsp_proc\polynomials\laurent.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
def __init__(
    self,
    coefficients: Sequence[complex] | np.ndarray | Mapping[int, complex],
    min_degree: int | None = None,
    *,
    trim_tol: float = 0.0,
) -> None:
    """Create a Laurent polynomial from dense or sparse coefficients."""
    if isinstance(coefficients, Mapping):
        if not coefficients:
            self._coeffs = np.array([0.0 + 0.0j], dtype=np.complex128)
            self._min_degree = 0
            return
        min_k = int(min(coefficients))
        max_k = int(max(coefficients))
        dense = np.zeros(max_k - min_k + 1, dtype=np.complex128)
        for exp, value in coefficients.items():
            dense[int(exp) - min_k] = complex(value)
        raw_coeffs = dense
        raw_min = min_k
    else:
        dense = np.asarray(coefficients, dtype=np.complex128)
        if dense.ndim != 1 or dense.size == 0:
            raise ValueError("Dense coefficients must be a non-empty 1D sequence.")
        raw_coeffs = dense.copy()
        raw_min = -((dense.size - 1) // 2) if min_degree is None else int(min_degree)

    self._coeffs, self._min_degree = trim_dense_support(raw_coeffs, raw_min, tol=trim_tol)

coefficient(exponent)

Return coefficient of z**exponent (or 0 outside support).

Source code in src\qsp_proc\polynomials\laurent.py 
73
74
75
76
77
def coefficient(self, exponent: int) -> complex:
    """Return coefficient of ``z**exponent`` (or ``0`` outside support)."""
    if exponent < self.min_degree or exponent > self.max_degree:
        return 0.0 + 0.0j
    return complex(self._coeffs[exponent - self.min_degree])

__call__(z)

Evaluate the Laurent polynomial at scalar or array z.

Source code in src\qsp_proc\polynomials\laurent.py 
79
80
81
82
83
84
85
86
def __call__(self, z: complex | np.ndarray) -> complex | np.ndarray:
    """Evaluate the Laurent polynomial at scalar or array ``z``."""
    exponents = np.arange(self.min_degree, self.max_degree + 1, dtype=np.int64)
    z_arr = np.asarray(z, dtype=np.complex128)
    values = np.sum(self._coeffs * z_arr[..., None] ** exponents, axis=-1)
    if np.ndim(z_arr) == 0:
        return complex(values)
    return np.asarray(values)

evaluate_on_unit_circle(theta)

Evaluate on U(1) via z = exp(i*theta).

Source code in src\qsp_proc\polynomials\laurent.py 
88
89
90
91
def evaluate_on_unit_circle(self, theta: float | np.ndarray) -> complex | np.ndarray:
    """Evaluate on ``U(1)`` via ``z = exp(i*theta)``."""
    theta_arr = np.asarray(theta, dtype=np.float64)
    return self(np.exp(1j * theta_arr))

linf_norm_on_unit_circle(num_samples=4096)

Estimate sampled ||P||_inf on U(1).

Source code in src\qsp_proc\polynomials\laurent.py 
93
94
95
def linf_norm_on_unit_circle(self, num_samples: int = 4096) -> float:
    """Estimate sampled ``||P||_inf`` on ``U(1)``."""
    return sampled_linf_on_unit_circle(evaluator=self.__call__, num_samples=num_samples)

with_bounded_linf_on_unit_circle(bound, num_samples=4096)

Return scaled polynomial with sampled unit-circle norm at most bound.

Source code in src\qsp_proc\polynomials\laurent.py 
 97
 98
 99
100
101
102
103
104
105
106
def with_bounded_linf_on_unit_circle(
    self, bound: float, num_samples: int = 4096
) -> LaurentPoly:
    """Return scaled polynomial with sampled unit-circle norm at most ``bound``."""
    if bound < 0:
        raise ValueError("bound must be non-negative.")
    current = self.linf_norm_on_unit_circle(num_samples=num_samples)
    if current == 0.0 or current <= bound:
        return LaurentPoly(self._coeffs, self._min_degree)
    return LaurentPoly((bound / current) * self._coeffs, self._min_degree)

is_reciprocal(tol=1e-12)

Check Appendix A reciprocal condition F(z)=F(z^{-1}).

Source code in src\qsp_proc\polynomials\laurent.py 
108
109
110
111
112
def is_reciprocal(self, tol: float = 1e-12) -> bool:
    """Check Appendix A reciprocal condition ``F(z)=F(z^{-1})``."""
    if self.min_degree != -self.max_degree:
        return False
    return bool(np.allclose(self._coeffs, self._coeffs[::-1], atol=tol, rtol=0.0))

is_anti_reciprocal(tol=1e-12)

Check Appendix A anti-reciprocal condition F(z)=-F(z^{-1}).

Source code in src\qsp_proc\polynomials\laurent.py 
114
115
116
117
118
def is_anti_reciprocal(self, tol: float = 1e-12) -> bool:
    """Check Appendix A anti-reciprocal condition ``F(z)=-F(z^{-1})``."""
    if self.min_degree != -self.max_degree:
        return False
    return bool(np.allclose(self._coeffs, -self._coeffs[::-1], atol=tol, rtol=0.0))

to_chebyshev(tol=1e-12)

Convert reciprocal Laurent polynomial to first-kind Chebyshev series.

Source code in src\qsp_proc\polynomials\laurent.py 
120
121
122
123
124
125
126
127
128
129
def to_chebyshev(self, tol: float = 1e-12) -> ChebyshevPoly:
    """Convert reciprocal Laurent polynomial to first-kind Chebyshev series."""
    from .chebyshev import ChebyshevPoly

    first_kind, second_kind = self.to_appendix_a_components()
    if second_kind.size > 0 and np.max(np.abs(second_kind)) > tol:
        raise ValueError(
            "Cannot convert to Chebyshev first-kind only: anti-reciprocal component is non-zero."
        )
    return first_kind

to_appendix_a_components()

Return Appendix A Chebyshev decomposition of a Laurent polynomial.

Any Laurent series can be written as sum a_n T_n(x) + sqrt(1-x^2) * sum b_n U_n(x) under x=(z+z^{-1})/2.

Returns:

Type Description
ChebyshevPoly

A pair (a_poly, b_coeffs) where:
ndarray
  • a_poly is first-kind Chebyshev polynomial sum a_n T_n.
tuple[ChebyshevPoly, ndarray]
  • b_coeffs[n] multiplies U_n in the second-kind series.
Source code in src\qsp_proc\polynomials\laurent.py 
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
def to_appendix_a_components(self) -> tuple[ChebyshevPoly, np.ndarray]:
    """Return Appendix A Chebyshev decomposition of a Laurent polynomial.

    Any Laurent series can be written as
    ``sum a_n T_n(x) + sqrt(1-x^2) * sum b_n U_n(x)``
    under ``x=(z+z^{-1})/2``.

    Returns:
        A pair ``(a_poly, b_coeffs)`` where:
        - ``a_poly`` is first-kind Chebyshev polynomial ``sum a_n T_n``.
        - ``b_coeffs[n]`` multiplies ``U_n`` in the second-kind series.
    """
    from .chebyshev import ChebyshevPoly

    d = self.degree
    first = np.zeros(d + 1, dtype=np.complex128)
    second = np.zeros(d, dtype=np.complex128)

    first[0] = self.coefficient(0)
    for n in range(1, d + 1):
        cp = self.coefficient(n)
        cm = self.coefficient(-n)
        first[n] = cp + cm
        second[n - 1] = 1j * (cp - cm)

    return ChebyshevPoly(first), second

from_appendix_a_components(first_kind, second_kind_coeffs) classmethod

Build Laurent polynomial from Appendix A Chebyshev components.

The input represents sum a_n T_n(x) + sqrt(1-x^2) * sum b_n U_n(x).

Source code in src\qsp_proc\polynomials\laurent.py 
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
@classmethod
def from_appendix_a_components(
    cls, first_kind: ChebyshevPoly, second_kind_coeffs: Sequence[complex] | np.ndarray
) -> LaurentPoly:
    """Build Laurent polynomial from Appendix A Chebyshev components.

    The input represents
    ``sum a_n T_n(x) + sqrt(1-x^2) * sum b_n U_n(x)``.
    """
    a = first_kind.coeffs
    b = np.asarray(second_kind_coeffs, dtype=np.complex128)
    degree = max(len(a) - 1, len(b))
    center = degree
    dense = np.zeros(2 * degree + 1, dtype=np.complex128)

    dense[center] += a[0]
    for n in range(1, len(a)):
        half = 0.5 * a[n]
        dense[center + n] += half
        dense[center - n] += half
    for n in range(1, len(b) + 1):
        coeff = b[n - 1]
        dense[center + n] += -0.5j * coeff
        dense[center - n] += 0.5j * coeff
    return cls(dense, min_degree=-degree)

from_chebyshev(poly) classmethod

Construct Laurent polynomial from a first-kind Chebyshev series.

Source code in src\qsp_proc\polynomials\laurent.py 
184
185
186
187
@classmethod
def from_chebyshev(cls, poly: ChebyshevPoly) -> LaurentPoly:
    """Construct Laurent polynomial from a first-kind Chebyshev series."""
    return poly.to_laurent()

__mul__(other)

Scalar multiplication P * c.

Source code in src\qsp_proc\polynomials\laurent.py 
189
190
191
def __mul__(self, other: complex | float | int) -> LaurentPoly:
    """Scalar multiplication ``P * c``."""
    return LaurentPoly(self._coeffs * complex(other), self._min_degree)

__rmul__(other)

Scalar multiplication c * P.

Source code in src\qsp_proc\polynomials\laurent.py 
193
194
195
def __rmul__(self, other: complex | float | int) -> LaurentPoly:
    """Scalar multiplication ``c * P``."""
    return self.__mul__(other)

__repr__()

Return unambiguous string representation.

Source code in src\qsp_proc\polynomials\laurent.py 
197
198
199
def __repr__(self) -> str:
    """Return unambiguous string representation."""
    return f"LaurentPoly(coefficients={self._coeffs!r}, min_degree={self._min_degree})"

qsp_proc.polynomials.approximations

Target polynomial constructors

hamiltonian_sim_degree_bound(tau, epsilon)

Return the Appendix B.2 piecewise truncation bound r_tilde(tau, epsilon).

This implements Eq. (73) from the paper for tau >= 0 and 0 < epsilon < 1. For negative tau, the bound is computed with |tau|.

Source code in src\qsp_proc\polynomials\approximations.py 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
def hamiltonian_sim_degree_bound(tau: float, epsilon: float) -> int:
    """Return the Appendix B.2 piecewise truncation bound ``r_tilde(tau, epsilon)``.

    This implements Eq. (73) from the paper for ``tau >= 0`` and ``0 < epsilon < 1``.
    For negative ``tau``, the bound is computed with ``|tau|``.
    """
    if not np.isfinite(tau):
        raise ValueError("tau must be finite.")
    if not (0.0 < epsilon < 1.0):
        raise ValueError("epsilon must satisfy 0 < epsilon < 1.")

    abs_tau = float(abs(tau))
    if abs_tau == 0.0:
        return 0

    log_term = math.log(1.0 / epsilon)
    threshold = log_term / math.e

    if abs_tau >= threshold:
        return int(math.ceil(math.e * abs_tau + log_term))

    denominator = math.log(math.e + (log_term / abs_tau))
    return int(math.ceil((4.0 * log_term) / denominator))

jacobi_anger_laurent(tau, epsilon, *, max_degree=None, subnormalization=math.sqrt(0.5))

Build truncated Jacobi-Anger Laurent target for Hamiltonian simulation.

The target corresponds to exp(i * tau * x) with x = cos(theta) and z = exp(i * theta), truncated to degree r where r is either max_degree or the Appendix B.2 Eq. (73) bound.

Coefficients are generated from the Fourier/Jacobi-Anger identity exp(i * tau * cos(theta)) = sum_k i^k J_k(tau) z^k for k in [-r, r].

Parameters:

Name Type Description Default
tau float

Simulation time parameter.

 required
epsilon float

Approximation precision for selecting truncation degree.

 required
max_degree int | None

Explicit truncation degree override.

 None
subnormalization float

Optional multiplicative scaling (paper uses sqrt(1/2)).

 sqrt(0.5)

Returns:

Type Description
LaurentPoly

Truncated Laurent polynomial target.
Source code in src\qsp_proc\polynomials\approximations.py 
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
def jacobi_anger_laurent(
    tau: float,
    epsilon: float,
    *,
    max_degree: int | None = None,
    subnormalization: float = math.sqrt(0.5),
) -> LaurentPoly:
    """Build truncated Jacobi-Anger Laurent target for Hamiltonian simulation.

    The target corresponds to ``exp(i * tau * x)`` with ``x = cos(theta)`` and
    ``z = exp(i * theta)``, truncated to degree ``r`` where ``r`` is either
    ``max_degree`` or the Appendix B.2 Eq. (73) bound.

    Coefficients are generated from the Fourier/Jacobi-Anger identity
    ``exp(i * tau * cos(theta)) = sum_k i^k J_k(tau) z^k`` for ``k in [-r, r]``.

    Args:
        tau: Simulation time parameter.
        epsilon: Approximation precision for selecting truncation degree.
        max_degree: Explicit truncation degree override.
        subnormalization: Optional multiplicative scaling (paper uses ``sqrt(1/2)``).

    Returns:
        Truncated Laurent polynomial target.
    """
    if not np.isfinite(tau):
        raise ValueError("tau must be finite.")
    if not (0.0 < epsilon < 1.0):
        raise ValueError("epsilon must satisfy 0 < epsilon < 1.")
    if max_degree is not None and max_degree < 0:
        raise ValueError("max_degree must be non-negative when provided.")
    if subnormalization <= 0.0:
        raise ValueError("subnormalization must be positive.")

    try:
        from scipy.special import jv
    except Exception as exc:  # pragma: no cover
        raise ImportError("jacobi_anger_laurent requires scipy.special.jv.") from exc

    degree = hamiltonian_sim_degree_bound(tau, epsilon) if max_degree is None else int(max_degree)
    coeffs: dict[int, complex] = {}
    for k in range(-degree, degree + 1):
        coeffs[k] = complex((1j**k) * jv(k, tau))

    dense = LaurentPoly(coeffs)
    if subnormalization != 1.0:
        return LaurentPoly(dense.coeffs * subnormalization, dense.min_degree)
    return dense

random_complex_target_laurent(degree, *, nz=None, decay_rate=1.5, norm_bound=0.5, rng=None)

Generate Appendix B.1 style sparse random complex target.

The target is sampled in theta-form: P(theta) = sum_j a_j cos(j theta) + i sum_j b_j sin(j theta), then converted to Laurent coefficients.

Parameters:

Name Type Description Default
degree int

Maximum Fourier/Chebyshev order.

 required
nz int | None

Number of non-zero coefficients across {a_j} and {b_j}. Defaults to max(5, degree // 15) as in Appendix B.1.

 None
decay_rate float

Magnitude decay exponent for higher-order terms.

 1.5
norm_bound float

Post-normalization bound for sampled ||P||_inf on U(1).

 0.5
rng Generator | None

Optional NumPy random generator.

 None

Returns:

Type Description
LaurentPoly

Sparse random Laurent polynomial with sampled norm at most norm_bound.
Source code in src\qsp_proc\polynomials\approximations.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def random_complex_target_laurent(
    degree: int,
    *,
    nz: int | None = None,
    decay_rate: float = 1.5,
    norm_bound: float = 0.5,
    rng: np.random.Generator | None = None,
) -> LaurentPoly:
    """Generate Appendix B.1 style sparse random complex target.

    The target is sampled in theta-form:
    ``P(theta) = sum_j a_j cos(j theta) + i sum_j b_j sin(j theta)``,
    then converted to Laurent coefficients.

    Args:
        degree: Maximum Fourier/Chebyshev order.
        nz: Number of non-zero coefficients across ``{a_j}`` and ``{b_j}``.
            Defaults to ``max(5, degree // 15)`` as in Appendix B.1.
        decay_rate: Magnitude decay exponent for higher-order terms.
        norm_bound: Post-normalization bound for sampled ``||P||_inf`` on ``U(1)``.
        rng: Optional NumPy random generator.

    Returns:
        Sparse random Laurent polynomial with sampled norm at most ``norm_bound``.
    """
    if degree < 0:
        raise ValueError("degree must be non-negative.")
    if decay_rate <= 0.0:
        raise ValueError("decay_rate must be positive.")
    if norm_bound < 0.0:
        raise ValueError("norm_bound must be non-negative.")

    rng = np.random.default_rng() if rng is None else rng
    default_nz = max(5, degree // 15) if degree > 0 else 1
    total_candidates = (degree + 1) + degree  # a_0..a_degree and b_1..b_degree
    nnz = default_nz if nz is None else nz
    if nnz < 0 or nnz > total_candidates:
        raise ValueError(f"nz must satisfy 0 <= nz <= {total_candidates}.")

    a = np.zeros(degree + 1, dtype=np.float64)
    b = np.zeros(degree + 1, dtype=np.float64)
    candidates = [("a", j) for j in range(degree + 1)] + [("b", j) for j in range(1, degree + 1)]
    if nnz > 0:
        chosen = rng.choice(len(candidates), size=nnz, replace=False)
        for idx in chosen:
            label, order = candidates[int(idx)]
            scale = float((order + 1) ** (-decay_rate))
            value = float(rng.uniform(-1.0, 1.0) * scale)
            if label == "a":
                a[order] = value
            else:
                b[order] = value

    coeffs: dict[int, complex] = {0: complex(a[0])}
    for j in range(1, degree + 1):
        c_pos = 0.5 * (a[j] + b[j])
        c_neg = 0.5 * (a[j] - b[j])
        if c_pos != 0.0:
            coeffs[j] = complex(c_pos)
        if c_neg != 0.0:
            coeffs[-j] = complex(c_neg)

    return LaurentPoly(coeffs).with_bounded_linf_on_unit_circle(bound=norm_bound)

qsp_proc.polynomials.utils

Shared helpers for polynomial operations

unit_circle_points(num_samples)

Return uniformly sampled points on the unit circle.

Parameters:

Name Type Description Default
num_samples int

Number of points to sample on U(1).

 required

Returns:

Type Description
ndarray

Complex array z_k = exp(i * theta_k) with theta_k uniformly
ndarray

spaced in [0, 2*pi).

Raises:

Type Description
ValueError

If num_samples is not positive.
Source code in src\qsp_proc\polynomials\utils.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def unit_circle_points(num_samples: int) -> np.ndarray:
    """Return uniformly sampled points on the unit circle.

    Args:
        num_samples: Number of points to sample on ``U(1)``.

    Returns:
        Complex array ``z_k = exp(i * theta_k)`` with ``theta_k`` uniformly
        spaced in ``[0, 2*pi)``.

    Raises:
        ValueError: If ``num_samples`` is not positive.
    """
    if num_samples <= 0:
        raise ValueError("num_samples must be positive.")
    theta = np.linspace(0.0, 2.0 * np.pi, num_samples, endpoint=False, dtype=np.float64)
    return np.exp(1j * theta)

sampled_linf_on_unit_circle(evaluator, num_samples)

Estimate ||f||_inf on U(1) via uniform sampling.

Parameters:

Name Type Description Default
evaluator Callable[[ndarray], ndarray]

Vectorized callable that accepts sampled z values.

 required
num_samples int

Number of points used for sampling.

 required

Returns:

Type Description
float

Estimated sup norm on sampled points.
Source code in src\qsp_proc\polynomials\utils.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def sampled_linf_on_unit_circle(
    evaluator: Callable[[np.ndarray], np.ndarray], num_samples: int
) -> float:
    """Estimate ``||f||_inf`` on ``U(1)`` via uniform sampling.

    Args:
        evaluator: Vectorized callable that accepts sampled ``z`` values.
        num_samples: Number of points used for sampling.

    Returns:
        Estimated sup norm on sampled points.
    """
    z = unit_circle_points(num_samples)
    values = evaluator(z)
    return float(np.max(np.abs(values)))

trim_dense_support(coeffs, min_degree, *, tol)

Trim near-zero boundary entries in a dense Laurent support.

Parameters:

Name Type Description Default
coeffs ndarray

Dense coefficient array where coeffs[i] maps to exponent min_degree + i.

 required
min_degree int

Smallest exponent represented by coeffs.

 required
tol float

Trimming tolerance.

 required

Returns:

Type Description
tuple[ndarray, int]

A pair (trimmed_coeffs, trimmed_min_degree).
Source code in src\qsp_proc\polynomials\utils.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
def trim_dense_support(
    coeffs: np.ndarray, min_degree: int, *, tol: float
) -> tuple[np.ndarray, int]:
    """Trim near-zero boundary entries in a dense Laurent support.

    Args:
        coeffs: Dense coefficient array where ``coeffs[i]`` maps to exponent
            ``min_degree + i``.
        min_degree: Smallest exponent represented by ``coeffs``.
        tol: Trimming tolerance.

    Returns:
        A pair ``(trimmed_coeffs, trimmed_min_degree)``.
    """
    if coeffs.size == 0:
        return np.array([0.0 + 0.0j], dtype=np.complex128), 0

    mask = np.abs(coeffs) > tol
    if not np.any(mask):
        return np.array([0.0 + 0.0j], dtype=np.complex128), 0

    left = int(np.argmax(mask))
    right = int(len(mask) - np.argmax(mask[::-1]))
    return coeffs[left:right].copy(), min_degree + left