Coverage for src/chebpy/algorithms.py: 97%

1"""Numerical algorithms for Chebyshev approximation and manipulation. 

2 

3This module provides core numerical algorithms used throughout the ChebPy package, 

4including rootfinding, barycentric interpolation, Chebyshev coefficient manipulation, 

5and adaptive approximation techniques. 

6 

7The algorithms implemented here are based on established numerical methods for 

8working with Chebyshev polynomials and approximations, many of which are described 

9in Trefethen's "Approximation Theory and Approximation Practice". 

10""" 

11 

12import warnings 

13from collections.abc import Callable 

14from typing import Any, cast 

15 

16import numpy as np 

17from numpy.fft import fft, ifft 

18 

19from .decorators import preandpostprocess 

20from .settings import _preferences as prefs 

21from .utilities import Interval, infnorm 

22 

23# supress numpy division and multiply warnings 

24np.seterr(divide="ignore", invalid="ignore") 

25 

26# constants 

27SPLITPOINT = -0.004849834917525 

28 

29 

30# local helpers 

31def find(x: np.ndarray) -> np.ndarray: 

32 """Find the indices of non-zero elements in an array. 

33 

34 A simple wrapper around numpy.where that returns only the indices. 

35 

36 Args: 

37 x (array-like): Input array. 

38 

39 Returns: 

40 numpy.ndarray: Indices of non-zero elements in the input array. 

41 """ 

42 return np.where(x)[0] 

43 

44 

45def rootsunit(ak: np.ndarray, htol: float | None = None) -> np.ndarray: 

46 """Compute the roots of a function on [-1,1] using Chebyshev coefficients. 

47 

48 This function finds the real roots of a function on the interval [-1,1] 

49 using the coefficients in its Chebyshev series representation. For large 

50 degree polynomials, it uses a recursive subdivision approach. 

51 

52 Args: 

53 ak (numpy.ndarray): Coefficients of the Chebyshev series. 

54 htol (float, optional): Tolerance for determining which roots to keep. 

55 Defaults to 100 * machine epsilon. 

56 

57 Returns: 

58 numpy.ndarray: Array of roots in the interval [-1,1], sorted in 

59 ascending order. 

60 

61 References: 

62 I. J. Good, "The colleague matrix, a Chebyshev analogue of the 

63 companion matrix", Quarterly Journal of Mathematics 12 (1961). 

64 

65 J. A. Boyd, "Computing zeros on a real interval through 

66 Chebyshev expansion and polynomial rootfinding", SIAM Journal on 

67 Numerical Analysis 40 (2002). 

68 

69 L. N. Trefethen, Approximation Theory and Approximation 

70 Practice, SIAM, 2013, chapter 18. 

71 """ 

72 htol = htol if htol is not None else 1e2 * prefs.eps 

73 n = standard_chop(ak, tol=htol) 

74 ak = ak[:n] 

75 

76 # if n > 50, we split and recurse 

77 if n > 50: 

78 chebpts = chebpts2(ak.size) 

79 lmap = Interval(-1, SPLITPOINT) 

80 rmap = Interval(SPLITPOINT, 1) 

81 lpts = lmap(chebpts) 

82 rpts = rmap(chebpts) 

83 lval = clenshaw(lpts, ak) 

84 rval = clenshaw(rpts, ak) 

85 lcfs = vals2coeffs2(lval) 

86 rcfs = vals2coeffs2(rval) 

87 lrts = rootsunit(lcfs, 2 * htol) 

88 rrts = rootsunit(rcfs, 2 * htol) 

89 return np.append(lmap(lrts), rmap(rrts)) 

90 

91 # trivial base case 

92 if n <= 1: 

93 return np.array([]) 

94 

95 # nontrivial base case: either compute directly or solve 

96 # a Colleague Matrix eigenvalue problem 

97 if n == 2: 

98 rts = np.array([-ak[0] / ak[1]]) 

99 elif n <= 50: 

100 v = 0.5 * np.ones(n - 2) 

101 colleague_matrix = np.diag(v, -1) + np.diag(v, 1) 

102 colleague_matrix[0, 1] = 1 

103 coeffs_matrix = np.zeros(colleague_matrix.shape, dtype=ak.dtype) 

104 coeffs_matrix[-1, :] = ak[:-1] 

105 eigenvalue_matrix = colleague_matrix - 0.5 * 1.0 / ak[-1] * coeffs_matrix 

106 rts = np.linalg.eigvals(eigenvalue_matrix) 

107 

108 # discard values with large imaginary part and treat the remaining 

109 # ones as real; then sort and retain only the roots inside [-1,1] 

110 mask = abs(np.imag(rts)) < htol 

111 rts = np.real(rts[mask]) 

112 rts = rts[abs(rts) <= 1.0 + htol] 

113 rts = np.sort(rts) 

114 if rts.size >= 2: 

115 rts[0] = max([rts[0], -1]) 

116 rts[-1] = min([rts[-1], 1]) 

117 return rts 

118 

119 

120@preandpostprocess 

121def bary(xx: np.ndarray, fk: np.ndarray, xk: np.ndarray, vk: np.ndarray) -> np.ndarray: 

122 """Evaluate a function using the barycentric interpolation formula. 

123 

124 This function implements the barycentric interpolation formula for evaluating 

125 a function at arbitrary points given its values at a set of nodes. It uses 

126 an efficient algorithm that switches between two implementations based on 

127 the number of evaluation points. 

128 

129 Args: 

130 xx (numpy.ndarray): Array of evaluation points. 

131 fk (numpy.ndarray): Array of function values at the interpolation nodes xk. 

132 xk (numpy.ndarray): Array of interpolation nodes. 

133 vk (numpy.ndarray): Barycentric weights corresponding to the interpolation nodes xk. 

134 

135 Returns: 

136 numpy.ndarray: Function values at the evaluation points xx. 

137 

138 References: 

139 J.P. Berrut, L.N. Trefethen, Barycentric Lagrange Interpolation, SIAM 

140 Review (2004) 

141 """ 

142 # either iterate over the evaluation points, or ... 

143 dtype = np.result_type(xx, fk, xk, vk, 1.0) 

144 if xx.size < 4 * xk.size: 

145 out = np.zeros(xx.size, dtype=dtype) 

146 for i in range(xx.size): 

147 tt = vk / (xx[i] - xk) 

148 out[i] = np.dot(tt, fk) / tt.sum() 

149 

150 # ... iterate over the barycenters 

151 else: 

152 numer = np.zeros(xx.size, dtype=dtype) 

153 denom = np.zeros(xx.size) 

154 for j in range(xk.size): 

155 temp = vk[j] / (xx - xk[j]) 

156 numer = numer + temp * fk[j] 

157 denom = denom + temp 

158 out = numer / denom 

159 

160 # replace NaNs 

161 for k in find(np.isnan(out)): 

162 idx = find(xx[k] == xk) 

163 if idx.size > 0: 

164 out[k] = fk[idx[0]] 

165 

166 return out 

167 

168 

169def fh_barywts(x: np.ndarray, d: int, maxind: int | None = None) -> np.ndarray: 

170 """Compute Floater-Hormann barycentric weights. 

171 

172 Args: 

173 x: Interpolation nodes. 

174 d: Floater-Hormann blending degree. 

175 maxind: Optional number of leading weights to compute. 

176 

177 Returns: 

178 Barycentric weights for the Floater-Hormann rational interpolant. 

179 """ 

180 x = np.asarray(x, dtype=float) 

181 n = x.size - 1 

182 if n < 0: 

183 return np.array([]) 

184 if d < 0 or d > n: 

185 msg = f"d must satisfy 0 <= d <= {n}; got {d}" 

186 raise ValueError(msg) 

187 num_weights = min(n + 1, n + 1 if maxind is None else maxind) 

188 w = np.zeros(num_weights) 

189 for k in range(num_weights): 

190 for i in range(max(0, k - d), min(k, n - d) + 1): 

191 s = 1.0 

192 for j in range(i, min(n, i + d) + 1): 

193 if j != k: 

194 s /= x[k] - x[j] 

195 w[k] += (-1.0) ** i * s 

196 return w 

197 

198 

199def _funqui_max_degree(n: int) -> int: 

200 """Return Chebfun's sample-count dependent maximum FH degree.""" 

201 if n < 60: 

202 return min(n, 35) 

203 if n < 100: 

204 return 30 

205 if n < 1000: 

206 return 25 

207 if n < 5000: 

208 return 20 

209 return 15 

210 

211 

212def _fh_weights(x: np.ndarray, d: int) -> np.ndarray: 

213 """Compute FH weights using Chebfun's symmetric equispaced-data shortcut.""" 

214 n = x.size - 1 

215 if d <= (n + 1) / 2: 

216 wl = np.abs(fh_barywts(x, d, d + 1)) 

217 wm = np.full(max(n - 1 - 2 * d, 0), wl[-1]) 

218 w = np.concatenate((wl, wm)) 

219 w = w[: int(np.ceil((n + 1) / 2))] 

220 w = np.concatenate((w, w[::-1])) if n % 2 else np.concatenate((w[:-1], w[::-1])) 

221 w[::2] *= -1 

222 return w 

223 return fh_barywts(x, d) 

224 

225 

226def _funqui_degree(x: np.ndarray, values: np.ndarray) -> int: 

227 """Choose the Floater-Hormann degree used for equispaced data.""" 

228 n = values.size - 1 

229 maxd = _funqui_max_degree(n) 

230 if n <= 2: 

231 return min(4, n) 

232 

233 rm_index = np.unique(np.array([1, n - 2])) 

234 xrm = np.delete(x, rm_index) 

235 valsrm = np.delete(values, rm_index) 

236 

237 if np.linalg.norm(values[rm_index], ord=np.inf) < 2 * np.finfo(float).eps * np.linalg.norm(values, ord=np.inf): 

238 return min(4, n) 

239 

240 errs = [] 

241 for d in range(min(n - 2, maxd) + 1): 

242 if d <= (n - 5) / 2: 

243 wl = np.abs(fh_barywts(xrm, d, d + 2)) 

244 wr = np.abs(fh_barywts(xrm[::-1], d, d + 2))[::-1] 

245 wm = np.full(max(n - 5 - 2 * d, 0), wl[-1]) 

246 w = np.concatenate((wl, wm, wr)) 

247 w[::2] *= -1 

248 else: 

249 w = fh_barywts(xrm, d) 

250 held_out = bary(x[rm_index], valsrm, xrm, w) 

251 errs.append(float(np.max(np.abs(held_out - values[rm_index])))) 

252 if errs[-1] > 1000 * min(errs): 

253 break 

254 return int(np.argmin(errs)) 

255 

256 

257def funqui(values: np.ndarray, domain: np.ndarray) -> Callable[..., Any]: 

258 """Return a Floater-Hormann interpolant for equispaced endpoint data. 

259 

260 Args: 

261 values: One-dimensional sample values on an equispaced grid. 

262 domain: Two finite endpoints defining the sample interval. 

263 

264 Returns: 

265 Callable rational interpolant through the supplied data. 

266 """ 

267 values = np.asarray(values) 

268 domain = np.asarray(domain, dtype=float) 

269 x = np.linspace(domain[0], domain[1], values.size) 

270 d = _funqui_degree(x, values) 

271 w = _fh_weights(x, d) 

272 return lambda zz: bary(zz, values, x, w) 

273 

274 

275@preandpostprocess 

276def clenshaw(xx: np.ndarray, ak: np.ndarray) -> np.ndarray: 

277 """Evaluate a Chebyshev series using Clenshaw's algorithm. 

278 

279 This function implements Clenshaw's algorithm for the evaluation of a 

280 first-kind Chebyshev series expansion at an array of points. 

281 

282 Args: 

283 xx (numpy.ndarray): Array of points at which to evaluate the series. 

284 ak (numpy.ndarray): Coefficients of the Chebyshev series. 

285 

286 Returns: 

287 numpy.ndarray: Values of the Chebyshev series at the points xx. 

288 

289 References: 

290 C. W. Clenshaw, "A note on the summation of Chebyshev series", 

291 Mathematics of Computation, Vol. 9, No. 51, 1955, pp. 118-120. 

292 

293 Examples: 

294 >>> import numpy as np 

295 >>> coeffs = np.array([1.0]) 

296 >>> x = np.array([0.0]) 

297 >>> result = clenshaw(x, coeffs) 

298 >>> bool(abs(float(result[0]) - 1.0) < 1e-10) 

299 True 

300 """ 

301 bk1 = 0 * xx 

302 bk2 = 0 * xx 

303 xx = 2 * xx 

304 idx = range(ak.size) 

305 for k in idx[ak.size : 1 : -2]: 

306 bk2 = ak[k] + xx * bk1 - bk2 

307 bk1 = ak[k - 1] + xx * bk2 - bk1 

308 if np.mod(ak.size - 1, 2) == 1: 

309 bk1, bk2 = ak[1] + xx * bk1 - bk2, bk1 

310 out: np.ndarray = ak[0] + 0.5 * xx * bk1 - bk2 

311 return out 

312 

313 

314def standard_chop(coeffs: np.ndarray, tol: float | None = None) -> int: 

315 """Determine where to truncate a Chebyshev series based on coefficient decay. 

316 

317 This function determines an appropriate cutoff point for a Chebyshev series 

318 by analyzing the decay of its coefficients. It implements the algorithm 

319 described by Aurentz and Trefethen. 

320 

321 Args: 

322 coeffs (numpy.ndarray): Coefficients of the Chebyshev series. 

323 tol (float, optional): Tolerance for determining the cutoff point. 

324 Defaults to machine epsilon from preferences. 

325 

326 Returns: 

327 int: Index at which to truncate the series. 

328 

329 References: 

330 J. Aurentz and L.N. Trefethen, "Chopping a Chebyshev series" (2015) 

331 (http://arxiv.org/pdf/1512.01803v1.pdf) 

332 """ 

333 # check magnitude of tol: 

334 tol = tol if tol is not None else prefs.eps 

335 if tol >= 1: # pragma: no cover 

336 cutoff = 1 

337 return cutoff 

338 

339 # ensure length at least 17: 

340 n = coeffs.size 

341 cutoff = n 

342 if n < 17: 

343 return cutoff 

344 

345 # Step 1: Convert coeffs input to a new monotonically nonincreasing 

346 # vector (envelope) normalized to begin with the value 1. 

347 b = np.flipud(np.abs(coeffs)) 

348 m = np.flipud(np.maximum.accumulate(b)) 

349 if m[0] == 0.0: 

350 cutoff = 1 

351 return cutoff 

352 envelope = m / m[0] 

353 

354 # Step 2: Scan envelope for a value plateauPoint, the first point J-1, 

355 # if any, that is followed by a plateau. Uses 1-based j to match the 

356 # MATLAB reference implementation; envelope is indexed with [j-1]. 

357 for j in range(2, n + 1): 

358 j2 = round(1.25 * j + 5) 

359 if j2 > n: 

360 # there is no plateau: exit 

361 return cutoff 

362 e1 = envelope[j - 1] 

363 e2 = envelope[int(j2) - 1] 

364 r = 3 * (1 - np.log(e1) / np.log(tol)) 

365 plateau = (e1 == 0.0) | (e2 / e1 > r) 

366 if plateau: 

367 # a plateau has been found: go to Step 3 

368 plateau_point = j - 1 

369 break 

370 

371 # Step 3: Fix cutoff at a point where envelope, plus a linear function 

372 # included to bias the result towards the left end, is minimal. 

373 if envelope[plateau_point - 1] == 0.0: # pragma: no cover 

374 cutoff = plateau_point 

375 else: 

376 j3 = int(np.sum(envelope >= tol ** (7.0 / 6.0))) 

377 if j3 < j2: 

378 j2 = j3 + 1 

379 envelope[int(j2) - 1] = tol ** (7.0 / 6.0) 

380 cc = np.log10(envelope[: int(j2)]) 

381 cc = cc + np.linspace(0, (-1.0 / 3.0) * np.log10(tol), int(j2)) 

382 d = np.argmin(cc) 

383 cutoff = max(int(d), 1) 

384 return cutoff 

385 

386 

387def adaptive(cls: Any, fun: Callable[..., Any], hscale: float = 1, maxpow2: int | None = None) -> np.ndarray: 

388 """Adaptively determine the number of points needed to represent a function. 

389 

390 This function implements an adaptive algorithm to determine the appropriate 

391 number of points needed to represent a function to a specified tolerance. 

392 It cycles over powers of two, evaluating the function at Chebyshev points 

393 and checking if the resulting coefficients can be truncated. 

394 

395 Args: 

396 cls: The class that provides the _chebpts and _vals2coeffs methods. 

397 fun (callable): The function to be approximated. 

398 hscale (float, optional): Scale factor for the tolerance. Defaults to 1. 

399 maxpow2 (int, optional): Maximum power of 2 to try. If None, uses the 

400 value from preferences. 

401 

402 Returns: 

403 numpy.ndarray: Coefficients of the Chebyshev series representing the function. 

404 

405 Warns: 

406 UserWarning: If the constructor does not converge within the maximum 

407 number of iterations. 

408 """ 

409 minpow2 = 4 # 17 points 

410 maxpow2 = maxpow2 if maxpow2 is not None else prefs.maxpow2 

411 tol = prefs.eps * max(hscale, 1) 

412 coeffs: np.ndarray = np.array([]) 

413 for k in range(minpow2, max(minpow2, maxpow2) + 1): 

414 n = 2**k + 1 

415 points = cls._chebpts(n) 

416 values = fun(points) 

417 coeffs = cls._vals2coeffs(values) 

418 # If function values are at or below tolerance the function is 

419 # indistinguishable from zero (cf. classicCheck.m vscale==0 guard). 

420 vscale = np.max(np.abs(values)) 

421 if vscale <= tol: 

422 coeffs = np.array([0.0]) 

423 break 

424 chplen = standard_chop(coeffs, tol=tol) 

425 if chplen < coeffs.size: 

426 coeffs = coeffs[:chplen] 

427 break 

428 if k == maxpow2: 

429 warnings.warn(f"The {cls.__name__} constructor did not converge: using {n} points", stacklevel=2) 

430 break 

431 return coeffs 

432 

433 

434def coeffmult(fc: np.ndarray, gc: np.ndarray) -> np.ndarray: 

435 """Multiply two Chebyshev series in coefficient space. 

436 

437 This function performs multiplication of two Chebyshev series represented by 

438 their coefficients. It uses FFT-based convolution for efficiency. 

439 

440 Args: 

441 fc (numpy.ndarray): Coefficients of the first Chebyshev series. 

442 gc (numpy.ndarray): Coefficients of the second Chebyshev series. 

443 

444 Returns: 

445 numpy.ndarray: Coefficients of the product series. 

446 

447 Note: 

448 The input series must have the same length. 

449 """ 

450 fc_extended = np.append(2.0 * fc[:1], (fc[1:], fc[:0:-1])) 

451 gc_extended = np.append(2.0 * gc[:1], (gc[1:], gc[:0:-1])) 

452 ak = ifft(fft(fc_extended) * fft(gc_extended)) 

453 ak = np.append(ak[:1], ak[1:] + ak[:0:-1]) * 0.25 

454 ak = ak[: fc.size] 

455 inputcfs = np.append(fc, gc) 

456 out = np.real(ak) if np.isreal(inputcfs).all() else ak 

457 return out 

458 

459 

460def barywts2(n: int) -> np.ndarray: 

461 """Compute barycentric weights for Chebyshev points of the second kind. 

462 

463 This function calculates the barycentric weights used in the barycentric 

464 interpolation formula for Chebyshev points of the second kind. 

465 

466 Args: 

467 n (int): Number of points (n+1 weights will be computed). 

468 

469 Returns: 

470 numpy.ndarray: Array of barycentric weights. 

471 

472 Note: 

473 For Chebyshev points of the second kind, the weights have a simple 

474 explicit formula with alternating signs. 

475 """ 

476 if n == 0: 

477 wts = np.array([]) 

478 elif n == 1: 

479 wts = np.array([1]) 

480 else: 

481 wts = np.append(np.ones(n - 1), 0.5) 

482 wts[n - 2 :: -2] = -1 

483 wts[0] = 0.5 * wts[0] 

484 return wts 

485 

486 

487def chebpts2(n: int) -> np.ndarray: 

488 """Compute Chebyshev points of the second kind. 

489 

490 This function calculates the n Chebyshev points of the second kind in the 

491 interval [-1, 1], which are the extrema of the Chebyshev polynomial T_{n-1} 

492 together with the endpoints ±1. 

493 

494 Args: 

495 n (int): Number of points to compute. 

496 

497 Returns: 

498 numpy.ndarray: Array of n Chebyshev points of the second kind. 

499 

500 Note: 

501 The points are ordered from left to right on the interval [-1, 1]. 

502 """ 

503 if n == 1: 

504 pts = np.array([0.0]) 

505 else: 

506 nn = np.arange(n) 

507 pts = np.cos(nn[::-1] * np.pi / (n - 1)) 

508 return pts 

509 

510 

511def vals2coeffs2(vals: np.ndarray) -> np.ndarray: 

512 """Convert function values to Chebyshev coefficients. 

513 

514 This function maps function values at Chebyshev points of the second kind 

515 to coefficients of the corresponding first-kind Chebyshev polynomial expansion. 

516 It uses an FFT-based algorithm for efficiency. 

517 

518 Args: 

519 vals (numpy.ndarray): Function values at Chebyshev points of the second kind. 

520 

521 Returns: 

522 numpy.ndarray: Coefficients of the first-kind Chebyshev polynomial expansion. 

523 

524 Note: 

525 This transformation is the discrete cosine transform of type I (DCT-I), 

526 which is implemented here using FFT for efficiency. 

527 """ 

528 n = vals.size 

529 if n <= 1: 

530 coeffs = vals 

531 return coeffs 

532 tmp = np.append(vals[::-1], vals[1:-1]) 

533 if np.isreal(vals).all(): 

534 coeffs = ifft(tmp) 

535 coeffs = np.real(coeffs) 

536 elif np.isreal(1j * vals).all(): # pragma: no cover 

537 coeffs = ifft(np.imag(tmp)) 

538 coeffs = 1j * np.real(coeffs) 

539 else: 

540 coeffs = ifft(tmp) 

541 coeffs = coeffs[:n] 

542 coeffs[1 : n - 1] = 2 * coeffs[1 : n - 1] 

543 return coeffs 

544 

545 

546def coeffs2vals2(coeffs: np.ndarray) -> np.ndarray: 

547 """Convert Chebyshev coefficients to function values. 

548 

549 This function maps coefficients of a first-kind Chebyshev polynomial expansion 

550 to function values at Chebyshev points of the second kind. It uses an FFT-based 

551 algorithm for efficiency. 

552 

553 Args: 

554 coeffs (numpy.ndarray): Coefficients of the first-kind Chebyshev polynomial expansion. 

555 

556 Returns: 

557 numpy.ndarray: Function values at Chebyshev points of the second kind. 

558 

559 Note: 

560 This transformation is the inverse discrete cosine transform of type I (IDCT-I), 

561 which is implemented here using FFT for efficiency. It is the inverse of vals2coeffs2. 

562 """ 

563 n = coeffs.size 

564 if n <= 1: 

565 vals = coeffs 

566 return vals 

567 coeffs = coeffs.copy() 

568 coeffs[1 : n - 1] = 0.5 * coeffs[1 : n - 1] 

569 tmp = np.append(coeffs, coeffs[n - 2 : 0 : -1]) 

570 if np.isreal(coeffs).all(): 

571 vals = fft(tmp) 

572 vals = np.real(vals) 

573 elif np.isreal(1j * coeffs).all(): # pragma: no cover 

574 vals = fft(np.imag(tmp)) 

575 vals = 1j * np.real(vals) 

576 else: 

577 vals = fft(tmp) 

578 vals = vals[n - 1 :: -1] 

579 return vals 

580 

581 

582def cheb2leg(c: np.ndarray) -> np.ndarray: 

583 """Convert Chebyshev coefficients to Legendre coefficients. 

584 

585 Converts the vector ``c`` of Chebyshev coefficients to a vector of Legendre 

586 coefficients such that:: 

587 

588 c[0]*T_0 + c[1]*T_1 + ... = l[0]*P_0 + l[1]*P_1 + ... 

589 

590 Uses a stable O(n²) three-term recurrence derived from the Chebyshev 

591 recurrence ``T_n = 2x T_{n-1} - T_{n-2}``. 

592 

593 Args: 

594 c (array-like): Chebyshev coefficients. 

595 

596 Returns: 

597 numpy.ndarray: Legendre coefficients of the same polynomial. 

598 """ 

599 c = np.asarray(c, dtype=float) 

600 n = c.size 

601 if n <= 1: 

602 return c.copy() 

603 

604 # Build Legendre coefficients via the recurrence: 

605 # M[j, col] = coeff of P_j in T_col 

606 # Recurrence: M[j,col] = 2j/(2j-1)*M[j-1,col-1] 

607 # + 2(j+1)/(2j+3)*M[j+1,col-1] 

608 # - M[j,col-2] 

609 # Initial columns: M[:,0] = [1,0,...], M[:,1] = [0,1,0,...] 

610 leg_coeffs = np.zeros(n) 

611 

612 prev_prev = np.zeros(n) 

613 prev_prev[0] = 1.0 # T_0 = P_0 

614 leg_coeffs += c[0] * prev_prev 

615 

616 prev = np.zeros(n) 

617 prev[1] = 1.0 # T_1 = P_1 

618 leg_coeffs += c[1] * prev 

619 

620 j = np.arange(n) 

621 for col in range(2, n): 

622 curr = np.zeros(n) 

623 # 2j/(2j-1) * prev[j-1] (for j >= 1) 

624 curr[1:] += 2.0 * j[1:] / (2.0 * j[1:] - 1.0) * prev[:-1] 

625 # 2(j+1)/(2j+3) * prev[j+1] (for j+1 <= n-1) 

626 curr[:-1] += 2.0 * (j[:-1] + 1.0) / (2.0 * j[:-1] + 3.0) * prev[1:] 

627 curr -= prev_prev 

628 leg_coeffs += c[col] * curr 

629 prev_prev = prev 

630 prev = curr 

631 

632 return leg_coeffs 

633 

634 

635def leg2cheb(c: np.ndarray) -> np.ndarray: 

636 """Convert Legendre coefficients to Chebyshev coefficients. 

637 

638 Converts the vector ``c`` of Legendre coefficients to a vector of Chebyshev 

639 coefficients such that:: 

640 

641 c[0]*P_0 + c[1]*P_1 + ... = l[0]*T_0 + l[1]*T_1 + ... 

642 

643 Uses a stable O(n²) three-term recurrence derived from the Legendre 

644 recurrence ``(n+1) P_{n+1} = (2n+1) x P_n - n P_{n-1}``. 

645 

646 Args: 

647 c (array-like): Legendre coefficients. 

648 

649 Returns: 

650 numpy.ndarray: Chebyshev coefficients of the same polynomial. 

651 """ 

652 c = np.asarray(c, dtype=float) 

653 n = c.size 

654 if n == 0: 

655 return np.zeros(0) 

656 if n == 1: 

657 return np.array([c[0]]) 

658 

659 # Build Chebyshev coefficients via the Legendre recurrence. 

660 # The Chebyshev representation of P_j is computed column by column. 

661 # Multiplication by x in Chebyshev basis: 

662 # (x*f)[0] = f[1]/2 

663 # (x*f)[1] = f[0] + f[2]/2 

664 # (x*f)[k] = (f[k-1] + f[k+1])/2 for k >= 2 

665 result = np.zeros(n) 

666 

667 prev_prev = np.zeros(n) 

668 prev_prev[0] = 1.0 # P_0 = T_0 

669 result += c[0] * prev_prev 

670 

671 prev = np.zeros(n) 

672 prev[1] = 1.0 # P_1 = T_1 

673 result += c[1] * prev 

674 

675 for j in range(2, n): 

676 # x * prev in Chebyshev basis 

677 xprev = np.zeros(n) 

678 xprev[1] += prev[0] # from x*T_0 = T_1 

679 xprev[: n - 1] += prev[1:] / 2.0 # T_{k-1} from x*T_k for k>=1 

680 xprev[2:] += prev[1 : n - 1] / 2.0 # T_{k+1} from x*T_k for k>=1 

681 

682 curr = ((2 * j - 1) * xprev - (j - 1) * prev_prev) / j 

683 result += c[j] * curr 

684 prev_prev = prev 

685 prev = curr 

686 

687 return result 

688 

689 

690def _conv_legendre(a: np.ndarray, b: np.ndarray) -> tuple[np.ndarray, np.ndarray]: 

691 """Convolve two Legendre series using the Hale-Townsend algorithm. 

692 

693 Computes the convolution of two functions expressed as Legendre series on 

694 [-1, 1]. The result is a piecewise polynomial on [-2, 2], split into a 

695 left piece on [-2, 0] and a right piece on [0, 2]. The Legendre 

696 coefficients of each piece (with respect to the linear map of the piece 

697 to [-1, 1]) are returned. 

698 

699 The algorithm is based on: 

700 N. Hale and A. Townsend, "An algorithm for the convolution of Legendre 

701 series", SIAM J. Sci. Comput., 36(3), A1207-A1220, 2014. 

702 

703 Args: 

704 a (array-like): Legendre coefficients of the first function on [-1, 1]. 

705 b (array-like): Legendre coefficients of the second function on [-1, 1]. 

706 

707 Returns: 

708 tuple: (gamma_left, gamma_right) where each element is a 1-D array of 

709 Legendre coefficients for the left [-2, 0] and right [0, 2] pieces 

710 respectively. 

711 """ 

712 a = np.asarray(a, dtype=float).ravel() 

713 b = np.asarray(b, dtype=float).ravel() 

714 

715 # Ensure a has the higher (or equal) degree 

716 if len(b) > len(a): 

717 a, b = b, a 

718 

719 na, nb = len(a), len(b) 

720 mn = na + nb 

721 

722 # Pad a to length mn 

723 alpha = np.zeros(mn) 

724 alpha[:na] = a 

725 

726 # Build the tridiagonal S matrix (mn x mn), the Legendre cumulative-integral 

727 # operator (S f)(x) = ∫_{-1}^{x} f(t) dt. Entries (using column index n): 

728 # S[0, 0] = 1 

729 # S[n+1, n] = 1/(2n+1) for n >= 0 (sub-diagonal) 

730 # S[n-1, n] = -1/(2n+1) for n >= 1 (super-diagonal) 

731 k = np.arange(mn) 

732 main = np.zeros(mn) 

733 main[0] = 1.0 

734 sub = 1.0 / (2.0 * k[:-1] + 1.0) # [1, 1/3, 1/5, ...], length mn-1 

735 supra = -1.0 / (2.0 * k[1:] + 1.0) # [-1/3, -1/5, -1/7, ...], length mn-1 

736 

737 def _s_apply(v: np.ndarray) -> np.ndarray: 

738 """Apply the S matrix to vector v.""" 

739 res = main * v 

740 res[1:] += sub * v[:-1] 

741 res[:-1] += supra * v[1:] 

742 return cast(np.ndarray, res) 

743 

744 def _rec(alpha_arg: np.ndarray, beta: np.ndarray, sgn: float, s00: float) -> np.ndarray: 

745 """Compute Legendre coefficients of the convolution on one piece. 

746 

747 Uses the recurrence from Theorem 4.1 of Hale & Townsend (2014). 

748 """ 

749 n_beta = len(beta) 

750 # Save / restore main[0] for S 

751 save_main0 = main[0] 

752 main[0] = s00 

753 

754 # scl[k] = (-1)^k / (2k-1) for k=1,...,n_beta (1-indexed) 

755 scl = np.ones(n_beta) / (2.0 * np.arange(1, n_beta + 1) - 1.0) 

756 scl[1::2] = -scl[1::2] 

757 

758 # First column 

759 v_new = _s_apply(alpha_arg) 

760 v = v_new.copy() 

761 gamma = beta[0] * v_new.copy() 

762 beta_scl = scl * beta 

763 beta_scl[0] = 0.0 

764 gamma[0] += float(v_new[:n_beta].dot(beta_scl)) 

765 

766 if n_beta > 1: 

767 # Second column 

768 v_new = _s_apply(v) + sgn * v 

769 v_old = v.copy() 

770 v = v_new.copy() 

771 v_new[0] = 0.0 

772 gamma += beta[1] * v_new 

773 beta_scl = -beta_scl * (2.0 - 0.5) / (2.0 - 1.5) 

774 beta_scl[1] = 0.0 

775 gamma[1] += float(v_new[:n_beta].dot(beta_scl)) 

776 

777 # Remaining columns 

778 for nn in range(3, n_beta + 1): 

779 v_new = (2 * nn - 3) * _s_apply(v) + v_old 

780 v_new[: nn - 1] = 0.0 

781 gamma += v_new * beta[nn - 1] 

782 beta_scl = -beta_scl * (nn - 0.5) / (nn - 1.5) 

783 beta_scl[nn - 1] = 0.0 

784 gamma[nn - 1] += float(v_new[:n_beta].dot(beta_scl)) 

785 v_old = v.copy()