Coverage for src/algolib/maths/algebra/polynomial.py: 76%
112 statements
« prev ^ index » next coverage.py v7.10.4, created at 2025-08-20 19:37 +0000
1# src/algolib/maths/algebra/polynomial.py
2"""
3Lightweight univariate polynomial with real coefficients.
5Notes
6-----
7- Coefficients are stored in **ascending** degree order: ``p(x) = c0 + c1*x + … + cn*x**n``.
8- The public API is immutable; the internal representation is a tuple of ``float``.
9- Supported operations include addition, subtraction, multiplication (convolution),
10 evaluation (Horner scheme), derivative, and antiderivative (indefinite integral).
12Examples
13--------
14>>> p = Polynomial([1, 2, 3]) # 1 + 2x + 3x^2
15>>> q = Polynomial.identity() # 1
16>>> (p * q).coeffs == p.coeffs
17True
18>>> p.derivative().coeffs # d/dx (1 + 2x + 3x^2) = 2 + 6x
19(2.0, 6.0)
20"""
22from __future__ import annotations
24from dataclasses import dataclass
25from typing import Iterable, Sequence, Tuple, Union
27from algolib.exceptions import InvalidTypeError, InvalidValueError
29Number = Union[int, float]
32def _strip_trailing_zeros(cs: Sequence[Number]) -> Tuple[float, ...]:
33 """Remove trailing zeros to keep a canonical representation.
35 Parameters
36 ----------
37 cs : Sequence[Number]
38 Input coefficients in ascending degree order.
40 Returns
41 -------
42 Tuple[float, ...]
43 Coefficients with all high‑degree zeros removed. At least one
44 coefficient is preserved (the zero polynomial becomes ``(0.0,)``).
45 """
46 if not cs:
47 return (0.0,)
48 out = list(map(float, cs))
49 # remove highest-degree zeros, keep at least one coefficient
50 while len(out) > 1 and out[-1] == 0.0:
51 out.pop()
52 return tuple(out)
54def _horner_kahan(coeffs: Sequence[Number], x: Number) -> Number:
55 """Evaluate ``p(x)`` via Horner's method with Kahan compensation for the
56 accumulating sum.
58 Notes
59 -----
60 The Kahan term stabilizes the *addition* in each fused step ``s*x + c_k``.
61 This is a lightweight compromise that improves accuracy for floating‑point
62 inputs without adding significant overhead.
64 Parameters
65 ----------
66 coeffs : Sequence[Number]
67 Coefficients in ascending degree order.
68 x : Number
69 Evaluation point. Real and complex values are supported.
71 Returns
72 -------
73 Number
74 The evaluated value ``p(x)``.
75 """
76 n = len(coeffs)
77 if n == 1:
78 return float(coeffs[0])
79 s: Number = float(coeffs[-1])
80 c: Number = 0.0
81 for k in range(n - 2, -1, -1):
82 prod = s * x
83 y = float(coeffs[k]) - c
84 t = prod + y
85 c = (t - prod) - y
86 s = t
87 return s
91@dataclass(frozen=True, slots=True)
92class Polynomial:
93 """
94 Univariate polynomial with real coefficients.
96 Parameters
97 ----------
98 coeffs : Sequence[Number]
99 Coefficients in ascending degree order:
100 ``[c0, c1, ..., cn]`` represents :math:`\\sum_{k=0}^n c_k x^k`.
102 Raises
103 ------
104 InvalidTypeError
105 If any coefficient is not a real number.
106 InvalidValueError
107 If ``coeffs`` is empty.
109 Notes
110 -----
111 - Construction enforces a **canonical form** by stripping trailing zeros; the
112 zero polynomial is thus represented as ``(0.0,)`` and has degree 0.
113 - All operations return new ``Polynomial`` instances (immutable API).
114 """
116 coeffs: Tuple[float, ...] # stored canonical (no trailing zeros)
118 # ------------------------------- construction -------------------------------
120 def __init__(self, coeffs: Iterable[Number]) -> None:
121 # Validate and convert to floats, mapping type errors to our domain error.
122 tmp: list[float] = []
123 try:
124 for c in coeffs:
125 if not isinstance(c, (int, float)):
126 raise InvalidTypeError("coeffs must contain real numbers (int or float)")
127 tmp.append(float(c))
128 except InvalidTypeError:
129 # re-raise our own type error untouched
130 raise
131 except Exception as e:
132 # any other conversion error -> InvalidTypeError for a clean API
133 raise InvalidTypeError("coeffs must contain real numbers (int or float)") from e
135 if len(tmp) == 0:
136 raise InvalidValueError("coeffs must be non-empty")
138 # Strip trailing zeros but keep at least one
139 i = len(tmp) - 1
140 while i > 0 and tmp[i] == 0.0:
141 i -= 1
142 object.__setattr__(self, "coeffs", tuple(tmp[: i + 1]))
145 @staticmethod
146 def zeros(deg: int) -> "Polynomial":
147 """
148 Return the (canonical) zero polynomial.
150 Parameters
151 ----------
152 deg : int
153 Requested degree (non‑negative). This value is accepted for API symmetry,
154 but the canonical zero polynomial always has degree 0 after construction.
156 Returns
157 -------
158 Polynomial
159 The zero polynomial.
160 """
161 if deg < 0:
162 raise InvalidValueError("degree must be non-negative.")
163 return Polynomial([0.0] * (deg + 1))
165 @staticmethod
166 def constant(c: Number) -> "Polynomial":
167 """Return a constant polynomial ``p(x) = c``."""
168 return Polynomial([c])
170 @staticmethod
171 def identity() -> "Polynomial":
172 """Return the multiplicative identity polynomial ``1``."""
173 return Polynomial([1.0])
175 # --------------------------------- queries ----------------------------------
177 @property
178 def degree(self) -> int:
179 """Degree of the polynomial (``len(coeffs) - 1``). By convention, ``deg(0) = 0``."""
180 return len(self.coeffs) - 1
182 def __call__(self, x: Number) -> Number:
183 """
184 Evaluate ``p(x)``.
186 Parameters
187 ----------
188 x : float or complex
189 Evaluation point.
191 Returns
192 -------
193 float or complex
194 The value ``p(x)`` computed via Horner's method with lightweight
195 Kahan compensation.
196 """
197 return _horner_kahan(self.coeffs, x)
199 # -------------------------------- calculus ----------------------------------
201 def derivative(self) -> "Polynomial":
202 r"""
203 Return the analytical derivative :math:`p'(x)`.
205 Notes
206 -----
207 If :math:`p(x) = a_0 + a_1 x + \cdots + a_n x^n`, then
208 :math:`p'(x) = a_1 + 2 a_2 x + \cdots + n a_n x^{n-1}`.
209 """
211 if self.degree == 0:
212 return Polynomial([0.0])
213 der = [k * self.coeffs[k] for k in range(1, len(self.coeffs))]
214 return Polynomial(der)
216 def integral(self, c0: Number = 0.0) -> "Polynomial":
217 r"""
218 Return an antiderivative :math:`\int p(x)\,dx` with constant term ``c0``.
220 Parameters
221 ----------
222 c0 : Number, default=0.0
223 Constant of integration (the coefficient of :math:`x^0`).
225 Returns
226 -------
227 Polynomial
228 An antiderivative whose derivative equals ``self``.
229 """
230 out = [float(c0)]
231 out.extend(c / (k + 1.0) for k, c in enumerate(self.coeffs))
232 return Polynomial(out)
234 # --------------------------------- algebra ----------------------------------
236 def __add__(self, other: "Polynomial") -> "Polynomial":
237 if not isinstance(other, Polynomial):
238 raise InvalidTypeError("other must be Polynomial")
239 n = max(self.degree, other.degree) + 1
240 cs = [0.0] * n
241 for i, c in enumerate(self.coeffs):
242 cs[i] += c
243 for i, c in enumerate(other.coeffs):
244 cs[i] += c
245 return Polynomial(cs)
247 def __sub__(self, other: "Polynomial") -> "Polynomial":
248 if not isinstance(other, Polynomial):
249 raise InvalidTypeError("other must be Polynomial")
250 n = max(self.degree, other.degree) + 1
251 cs = [0.0] * n
252 for i, c in enumerate(self.coeffs):
253 cs[i] += c
254 for i, c in enumerate(other.coeffs):
255 cs[i] -= c
256 return Polynomial(cs)
258 def __mul__(self, other: "Polynomial") -> "Polynomial":
259 if not isinstance(other, Polynomial):
260 raise InvalidTypeError("other must be Polynomial")
261 na, nb = self.degree + 1, other.degree + 1
262 out = [0.0] * (na + nb - 1)
263 for i, ai in enumerate(self.coeffs):
264 for j, bj in enumerate(other.coeffs):
265 out[i + j] += ai * bj
266 return Polynomial(out)
268 # --------------------------------- display ----------------------------------
270 def __repr__(self) -> str:
271 """Unambiguous representation, e.g. ``Polynomial(coeffs=(1.0, 2.0))``."""
272 return f"Polynomial(coeffs={self.coeffs})"
275 def __str__(self) -> str:
276 """Human‑readable form such as ``"3x^2 + 2x + 1"`` (or ``"0"`` for the zero polynomial)."""
277 terms = []
278 for k, c in enumerate(self.coeffs):
279 if c == 0.0:
280 continue
281 if k == 0:
282 terms.append(f"{c:.12g}")
283 elif k == 1:
284 terms.append(f"{c:.12g}x")
285 else:
286 terms.append(f"{c:.12g}x^{k}")
287 return " + ".join(reversed(terms)) if terms else "0"