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

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: 

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 # Defensive: envelope[0] == 1 and a zero at index j-2 would already have 

374 # triggered the plateau at the previous j, so this never holds. 

375 if envelope[plateau_point - 1] == 0.0: # pragma: no cover - defensive, see above 

376 cutoff = plateau_point 

377 else: 

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

379 if j3 < j2: 

380 j2 = j3 + 1 

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

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

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

384 d = np.argmin(cc) 

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

386 return cutoff 

387 

388 

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

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

391 

392 This function implements an adaptive algorithm to determine the appropriate 

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

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

395 and checking if the resulting coefficients can be truncated. 

396 

397 Args: 

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

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

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

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

402 value from preferences. 

403 

404 Returns: 

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

406 

407 Warns: 

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

409 number of iterations. 

410 """ 

411 minpow2 = 4 # 17 points 

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

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

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

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

416 n = 2**k + 1 

417 points = cls._chebpts(n) 

418 values = fun(points) 

419 coeffs = cls._vals2coeffs(values) 

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

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

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

423 if vscale <= tol: 

424 coeffs = np.array([0.0]) 

425 break 

426 chplen = standard_chop(coeffs, tol=tol) 

427 if chplen < coeffs.size: 

428 coeffs = coeffs[:chplen] 

429 break 

430 if k == maxpow2: 

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

432 break 

433 return coeffs 

434 

435 

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

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

438 

439 This function performs multiplication of two Chebyshev series represented by 

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

441 

442 Args: 

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

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

445 

446 Returns: 

447 numpy.ndarray: Coefficients of the product series. 

448 

449 Note: 

450 The input series must have the same length. 

451 """ 

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

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

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

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

456 ak = ak[: fc.size] 

457 inputcfs = np.append(fc, gc) 

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

459 return out 

460 

461 

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

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

464 

465 This function calculates the barycentric weights used in the barycentric 

466 interpolation formula for Chebyshev points of the second kind. 

467 

468 Args: 

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

470 

471 Returns: 

472 numpy.ndarray: Array of barycentric weights. 

473 

474 Note: 

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

476 explicit formula with alternating signs. 

477 """ 

478 if n == 0: 

479 wts = np.array([]) 

480 elif n == 1: 

481 wts = np.array([1]) 

482 else: 

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

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

485 wts[0] = 0.5 * wts[0] 

486 return wts 

487 

488 

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

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

491 

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

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

494 together with the endpoints ±1. 

495 

496 Args: 

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

498 

499 Returns: 

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

501 

502 Note: 

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

504 """ 

505 if n == 1: 

506 pts = np.array([0.0]) 

507 else: 

508 nn = np.arange(n) 

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

510 return pts 

511 

512 

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

514 """Convert function values to Chebyshev coefficients. 

515 

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

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

518 It uses an FFT-based algorithm for efficiency. 

519 

520 Args: 

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

522 

523 Returns: 

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

525 

526 Note: 

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

528 which is implemented here using FFT for efficiency. 

529 """ 

530 n = vals.size 

531 if n <= 1: 

532 coeffs = vals 

533 return coeffs 

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

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

536 coeffs = ifft(tmp) 

537 coeffs = np.real(coeffs) 

538 elif np.isreal(1j * vals).all(): 

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

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

541 else: 

542 coeffs = ifft(tmp) 

543 coeffs = coeffs[:n] 

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

545 return coeffs 

546 

547 

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

549 """Convert Chebyshev coefficients to function values. 

550 

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

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

553 algorithm for efficiency. 

554 

555 Args: 

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

557 

558 Returns: 

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

560 

561 Note: 

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

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

564 """ 

565 n = coeffs.size 

566 if n <= 1: 

567 vals = coeffs 

568 return vals 

569 coeffs = coeffs.copy() 

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

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

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

573 vals = fft(tmp) 

574 vals = np.real(vals) 

575 elif np.isreal(1j * coeffs).all(): 

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

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

578 else: 

579 vals = fft(tmp) 

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

581 return vals 

582 

583 

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

585 """Convert Chebyshev coefficients to Legendre coefficients. 

586 

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

588 coefficients such that:: 

589 

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

591 

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

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

594 

595 Args: 

596 c (array-like): Chebyshev coefficients. 

597 

598 Returns: 

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

600 """ 

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

602 n = c.size 

603 if n <= 1: 

604 return c.copy() 

605 

606 # Build Legendre coefficients via the recurrence: 

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

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

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

610 # - M[j,col-2] 

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

612 leg_coeffs = np.zeros(n) 

613 

614 prev_prev = np.zeros(n) 

615 prev_prev[0] = 1.0 # T_0 = P_0 

616 leg_coeffs += c[0] * prev_prev 

617 

618 prev = np.zeros(n) 

619 prev[1] = 1.0 # T_1 = P_1 

620 leg_coeffs += c[1] * prev 

621 

622 j = np.arange(n) 

623 for col in range(2, n): 

624 curr = np.zeros(n) 

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

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

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

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

629 curr -= prev_prev 

630 leg_coeffs += c[col] * curr 

631 prev_prev = prev 

632 prev = curr 

633 

634 return leg_coeffs 

635 

636 

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

638 """Convert Legendre coefficients to Chebyshev coefficients. 

639 

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

641 coefficients such that:: 

642 

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

644 

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

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

647 

648 Args: 

649 c (array-like): Legendre coefficients. 

650 

651 Returns: 

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

653 """ 

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

655 n = c.size 

656 if n == 0: 

657 return np.zeros(0) 

658 if n == 1: 

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

660 

661 # Build Chebyshev coefficients via the Legendre recurrence. 

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

663 # Multiplication by x in Chebyshev basis: 

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

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

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

667 result = np.zeros(n) 

668 

669 prev_prev = np.zeros(n) 

670 prev_prev[0] = 1.0 # P_0 = T_0 

671 result += c[0] * prev_prev 

672 

673 prev = np.zeros(n) 

674 prev[1] = 1.0 # P_1 = T_1 

675 result += c[1] * prev 

676 

677 for j in range(2, n): 

678 # x * prev in Chebyshev basis 

679 xprev = np.zeros(n) 

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

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

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

683 

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

685 result += c[j] * curr 

686 prev_prev = prev 

687 prev = curr 

688 

689 return result 

690 

691 

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

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

694 

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

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

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

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

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

700 

701 The algorithm is based on: 

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

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

704 

705 Args: 

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

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

708 

709 Returns: 

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

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

712 respectively. 

713 """ 

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

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

716 

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

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

719 a, b = b, a 

720 

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

722 mn = na + nb 

723 

724 # Pad a to length mn 

725 alpha = np.zeros(mn) 

726 alpha[:na] = a 

727 

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

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

730 # S[0, 0] = 1 

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

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

733 k = np.arange(mn) 

734 main = np.zeros(mn) 

735 main[0] = 1.0 

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

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

738 

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

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

741 res = main * v 

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

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

744 return cast(np.ndarray, res) 

745 

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

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

748 

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

750 """ 

751 n_beta = len(beta) 

752 # Save / restore main[0] for S 

753 save_main0 = main[0] 

754 main[0] = s00 

755 

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

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

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

759 

760 # First column 

761 v_new = _s_apply(alpha_arg) 

762 v = v_new.copy() 

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

764 beta_scl = scl * beta 

765 beta_scl[0] = 0.0 

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

767 

768 if n_beta > 1: 

769 # Second column 

770 v_new = _s_apply(v) + sgn * v 

771 v_old = v.copy() 

772 v = v_new.copy() 

773 v_new[0] = 0.0 

774 gamma += beta[1] * v_new 

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

776 beta_scl[1] = 0.0 

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

778 

779 # Remaining columns 

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

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

782 v_new[: nn - 1] = 0.0 

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

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

785 beta_scl[nn - 1] = 0.0 

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