0
# Polynomial Operations
1

2
CuPy provides comprehensive polynomial operations through the `cupy.polynomial` module and legacy polynomial functions, enabling mathematical computation with polynomials including arithmetic, evaluation, fitting, root finding, and advanced polynomial manipulations optimized for GPU acceleration.
3

4
## Capabilities
5

6
### Legacy Polynomial Functions
7

8
Legacy polynomial functions providing basic polynomial operations and manipulations.
9

10
```python { .api }
11
class poly1d:
12
    """
13
    One-dimensional polynomial class.
14
    
15
    A convenience class for dealing with polynomials, providing
16
    methods for polynomial arithmetic, evaluation, and analysis.
17
    """
18
    def __init__(self, c_or_r, r=False, variable=None):
19
        """
20
        Parameters:
21
            c_or_r: array_like - Polynomial coefficients or roots
22
            r: bool, optional - If True, c_or_r contains roots instead of coefficients
23
            variable: str, optional - Variable name for string representation
24
        """
25
    
26
    def __call__(self, val):
27
        """Evaluate polynomial at given values."""
28
    
29
    def __add__(self, other):
30
        """Add polynomials."""
31
    
32
    def __sub__(self, other):
33
        """Subtract polynomials."""
34
    
35
    def __mul__(self, other):
36
        """Multiply polynomials."""
37
    
38
    def __div__(self, other):
39
        """Divide polynomials."""
40
    
41
    def __pow__(self, val):
42
        """Raise polynomial to a power."""
43
    
44
    def deriv(self, m=1):
45
        """Return derivative of polynomial."""
46
    
47
    def integ(self, m=1, k=None):
48
        """Return integral of polynomial."""
49
    
50
    @property
51
    def coeffs(self):
52
        """Polynomial coefficients."""
53
    
54
    @property
55
    def order(self):
56
        """Order (degree) of polynomial."""
57
    
58
    @property
59
    def roots(self):
60
        """Roots of polynomial."""
61

62
def poly(seq_of_zeros):
63
    """
64
    Find the coefficients of a polynomial with the given sequence of roots.
65
    
66
    Parameters:
67
        seq_of_zeros: array_like - Sequence of polynomial roots
68
    
69
    Returns:
70
        ndarray: Polynomial coefficients
71
    """
72

73
def polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False):
74
    """
75
    Least squares polynomial fit.
76
    
77
    Parameters:
78
        x: array_like - x-coordinates of data points
79
        y: array_like - y-coordinates of data points  
80
        deg: int - Degree of fitting polynomial
81
        rcond: float, optional - Relative condition number for fit
82
        full: bool, optional - Switch determining nature of return value
83
        w: array_like, optional - Weights to apply to y-coordinates
84
        cov: bool or str, optional - Return covariance matrix
85
    
86
    Returns:
87
        ndarray or tuple: Polynomial coefficients (and additional info if requested)
88
    """
89

90
def polyval(p, x):
91
    """
92
    Evaluate a polynomial at specific values.
93
    
94
    Parameters:
95
        p: array_like - Polynomial coefficients in decreasing powers
96
        x: array_like - Values at which to evaluate polynomial
97
    
98
    Returns:
99
        ndarray: Polynomial evaluated at x
100
    """
101

102
def polyadd(a1, a2):
103
    """
104
    Add two polynomials.
105
    
106
    Parameters:
107
        a1, a2: array_like - Polynomial coefficients
108
    
109
    Returns:
110
        ndarray: Coefficients of sum polynomial
111
    """
112

113
def polysub(a1, a2):
114
    """
115
    Subtract second polynomial from first.
116
    
117
    Parameters:
118
        a1, a2: array_like - Polynomial coefficients
119
    
120
    Returns:
121
        ndarray: Coefficients of difference polynomial
122
    """
123

124
def polymul(a1, a2):
125
    """
126
    Multiply two polynomials.
127
    
128
    Parameters:
129
        a1, a2: array_like - Polynomial coefficients
130
    
131
    Returns:
132
        ndarray: Coefficients of product polynomial
133
    """
134

135
def polydiv(u, v):
136
    """
137
    Divide first polynomial by second.
138
    
139
    Parameters:
140
        u, v: array_like - Polynomial coefficients (dividend and divisor)
141
    
142
    Returns:
143
        tuple: (quotient, remainder) polynomial coefficients
144
    """
145

146
def roots(p):
147
    """
148
    Return the roots of a polynomial with coefficients given in p.
149
    
150
    Parameters:
151
        p: array_like - Polynomial coefficients
152
    
153
    Returns:
154
        ndarray: Roots of the polynomial
155
    """
156
```
157

158
### Modern Polynomial Package
159

160
The `cupy.polynomial` package provides a comprehensive set of polynomial classes and functions.
161

162
```python { .api }
163
# Polynomial base classes and utilities
164
class Polynomial:
165
    """
166
    Base class for all polynomial classes.
167
    
168
    Provides common interface and functionality for polynomial
169
    operations across different polynomial bases.
170
    """
171
    def __init__(self, coef, domain=None, window=None): ...
172
    
173
    def __call__(self, x): ...
174
    def __add__(self, other): ...
175
    def __sub__(self, other): ...
176
    def __mul__(self, other): ...
177
    def __truediv__(self, other): ...
178
    def __pow__(self, other): ...
179
    
180
    def deriv(self, m=1): ...
181
    def integ(self, m=1, k=[]): ...
182
    def roots(self): ...
183
    
184
    @classmethod
185
    def fit(cls, x, y, deg, domain=None, rcond=None, full=False, w=None, window=None): ...
186
    
187
    @classmethod  
188
    def fromroots(cls, roots, domain=[], window=None): ...
189

190
# Power series polynomials (standard polynomials)
191
class Polynomial:
192
    """Power series polynomial in standard form."""
193

194
def polyval(x, c):
195
    """Evaluate polynomial at x."""
196

197
def polyval2d(x, y, c):
198
    """Evaluate 2-D polynomial at x, y."""
199

200
def polyval3d(x, y, z, c):
201
    """Evaluate 3-D polynomial at x, y, z."""
202

203
def polyvander(x, deg):
204
    """Generate Vandermonde matrix."""
205

206
def polyvander2d(x, y, deg):
207
    """Generate 2-D Vandermonde matrix."""
208

209
def polyvander3d(x, y, z, deg):
210
    """Generate 3-D Vandermonde matrix."""
211

212
def polyfit(x, y, deg, rcond=None, full=False, w=None):
213
    """Least squares fit of polynomial to data."""
214

215
def polycompanion(c):
216
    """Return companion matrix."""
217

218
def polyroots(c):
219
    """Return roots of polynomial."""
220

221
def polyfromroots(roots):
222
    """Generate polynomial from roots."""
223

224
# Chebyshev polynomials
225
class Chebyshev:
226
    """Chebyshev polynomial class."""
227

228
def chebval(x, c):
229
    """Evaluate Chebyshev series at x."""
230

231
def chebval2d(x, y, c):
232
    """Evaluate 2-D Chebyshev series."""
233

234
def chebval3d(x, y, z, c):
235
    """Evaluate 3-D Chebyshev series."""
236

237
def chebvander(x, deg):
238
    """Generate Chebyshev Vandermonde matrix."""
239

240
def chebfit(x, y, deg, rcond=None, full=False, w=None):
241
    """Fit Chebyshev series to data."""
242

243
def chebroots(c):
244
    """Return roots of Chebyshev polynomial."""
245

246
def chebfromroots(roots):
247
    """Generate Chebyshev polynomial from roots."""
248

249
# Legendre polynomials
250
class Legendre:
251
    """Legendre polynomial class."""
252

253
def legval(x, c):
254
    """Evaluate Legendre series at x."""
255

256
def legval2d(x, y, c):
257
    """Evaluate 2-D Legendre series."""
258

259
def legval3d(x, y, z, c):
260
    """Evaluate 3-D Legendre series."""
261

262
def legvander(x, deg):
263
    """Generate Legendre Vandermonde matrix."""
264

265
def legfit(x, y, deg, rcond=None, full=False, w=None):
266
    """Fit Legendre series to data."""
267

268
def legroots(c):
269
    """Return roots of Legendre polynomial."""
270

271
def legfromroots(roots):
272
    """Generate Legendre polynomial from roots."""
273

274
# Laguerre polynomials
275
class Laguerre:
276
    """Laguerre polynomial class."""
277

278
def lagval(x, c):
279
    """Evaluate Laguerre series at x."""
280

281
def lagval2d(x, y, c):
282
    """Evaluate 2-D Laguerre series."""
283

284
def lagval3d(x, y, z, c):
285
    """Evaluate 3-D Laguerre series."""
286

287
def lagvander(x, deg):
288
    """Generate Laguerre Vandermonde matrix."""
289

290
def lagfit(x, y, deg, rcond=None, full=False, w=None):
291
    """Fit Laguerre series to data."""
292

293
def lagroots(c):
294
    """Return roots of Laguerre polynomial."""
295

296
def lagfromroots(roots):
297
    """Generate Laguerre polynomial from roots."""
298

299
# Hermite polynomials (physicist's)
300
class HermiteE:
301
    """Hermite polynomial class (probabilist's)."""
302

303
def hermeval(x, c):
304
    """Evaluate Hermite series at x."""
305

306
def hermeval2d(x, y, c):
307
    """Evaluate 2-D Hermite series."""
308

309
def hermeval3d(x, y, z, c):
310
    """Evaluate 3-D Hermite series."""
311

312
def hermevander(x, deg):
313
    """Generate Hermite Vandermonde matrix."""
314

315
def hermefit(x, y, deg, rcond=None, full=False, w=None):
316
    """Fit Hermite series to data."""
317

318
def hermeroots(c):
319
    """Return roots of Hermite polynomial."""
320

321
def hermefromroots(roots):
322
    """Generate Hermite polynomial from roots."""
323

324
# Hermite polynomials (probabilist's)
325
class Hermite:
326
    """Hermite polynomial class (physicist's)."""
327

328
def hermval(x, c):
329
    """Evaluate Hermite series at x."""
330

331
def hermval2d(x, y, c):
332
    """Evaluate 2-D Hermite series."""
333

334
def hermval3d(x, y, z, c):
335
    """Evaluate 3-D Hermite series."""
336

337
def hermvander(x, deg):
338
    """Generate Hermite Vandermonde matrix."""
339

340
def hermfit(x, y, deg, rcond=None, full=False, w=None):
341
    """Fit Hermite series to data."""
342

343
def hermroots(c):
344
    """Return roots of Hermite polynomial."""
345

346
def hermfromroots(roots):
347
    """Generate Hermite polynomial from roots."""
348
```
349

350
## Usage Examples
351

352
### Basic Polynomial Operations
353

354
```python
355
import cupy as cp
356
import cupy.polynomial as poly
357

358
# Create and manipulate polynomials using poly1d
359
p1 = cp.poly1d([1, 2, 3])  # x^2 + 2x + 3
360
p2 = cp.poly1d([2, 1])     # 2x + 1
361

362
print("p1:", p1)
363
print("p2:", p2)
364

365
# Polynomial arithmetic
366
sum_poly = p1 + p2
367
diff_poly = p1 - p2  
368
prod_poly = p1 * p2
369
print("Sum:", sum_poly)
370
print("Product:", prod_poly)
371

372
# Evaluate polynomials
373
x_values = cp.linspace(-5, 5, 100)
374
y1 = p1(x_values)
375
y2 = p2(x_values)
376
print("p1 evaluated at x=2:", p1(2))
377
print("p2 evaluated at x=2:", p2(2))
378

379
# Find roots
380
roots_p1 = p1.roots
381
print("Roots of p1:", roots_p1)
382

383
# Derivatives and integrals
384
p1_deriv = p1.deriv()
385
p1_integ = p1.integ()
386
print("Derivative of p1:", p1_deriv)
387
print("Integral of p1:", p1_integ)
388
```
389

390
### Polynomial Fitting
391

392
```python
393
# Generate noisy data
394
x_data = cp.linspace(0, 10, 100)
395
true_coeffs = [0.5, -2, 1, 3]  # 0.5x^3 - 2x^2 + x + 3
396
y_true = cp.polyval(true_coeffs, x_data)
397
noise = cp.random.normal(0, 0.5, x_data.shape)
398
y_noisy = y_true + noise
399

400
# Fit polynomial of different degrees
401
degrees = [2, 3, 4, 5]
402
fitted_polys = {}
403

404
for degree in degrees:
405
    coeffs = cp.polyfit(x_data, y_noisy, degree)
406
    fitted_polys[degree] = coeffs
407
    
408
    # Evaluate fitted polynomial
409
    y_fitted = cp.polyval(coeffs, x_data)
410
    mse = cp.mean((y_noisy - y_fitted) ** 2)
411
    print(f"Degree {degree}: MSE = {mse:.4f}")
412
    print(f"  Coefficients: {coeffs}")
413

414
# Best fit analysis
415
best_degree = min(fitted_polys.keys(), 
416
                 key=lambda d: cp.mean((y_noisy - cp.polyval(fitted_polys[d], x_data)) ** 2))
417
print(f"Best degree: {best_degree}")
418

419
# Weighted fitting for heteroscedastic noise
420
weights = 1.0 / (1.0 + cp.abs(x_data))  # Less weight for larger x
421
weighted_coeffs = cp.polyfit(x_data, y_noisy, 3, w=weights)
422
print("Weighted fit coefficients:", weighted_coeffs)
423

424
# Robust fitting with full output
425
coeffs, residuals, rank, singular_vals, rcond = cp.polyfit(
426
    x_data, y_noisy, 3, full=True
427
)
428
print("Fit diagnostics:")
429
print(f"  Residuals: {residuals}")
430
print(f"  Rank: {rank}")
431
print(f"  Condition number: {1/rcond:.2e}")
432
```
433

434
### Advanced Polynomial Types
435

436
```python
437
# Chebyshev polynomials for better numerical properties
438
x_cheb = cp.linspace(-1, 1, 100)
439
y_data = cp.exp(x_cheb) + 0.1 * cp.random.normal(size=x_cheb.shape)
440

441
# Fit using Chebyshev polynomials
442
cheb_coeffs = poly.chebyshev.chebfit(x_cheb, y_data, 5)
443
y_cheb_fitted = poly.chebyshev.chebval(x_cheb, cheb_coeffs)
444

445
print("Chebyshev coefficients:", cheb_coeffs)
446
print("Chebyshev fit MSE:", cp.mean((y_data - y_cheb_fitted) ** 2))
447

448
# Compare with standard polynomial fit
449
std_coeffs = cp.polyfit(x_cheb, y_data, 5)
450
y_std_fitted = cp.polyval(std_coeffs, x_cheb)
451
print("Standard polynomial MSE:", cp.mean((y_data - y_std_fitted) ** 2))
452

453
# Legendre polynomial approximation
454
leg_coeffs = poly.legendre.legfit(x_cheb, y_data, 5)
455
y_leg_fitted = poly.legendre.legval(x_cheb, leg_coeffs)
456
print("Legendre fit MSE:", cp.mean((y_data - y_leg_fitted) ** 2))
457

458
# Hermite polynomials for Gaussian-weighted data
459
x_herm = cp.linspace(-3, 3, 100)
460
gaussian_weight = cp.exp(-x_herm**2 / 2)
461
y_weighted = cp.sin(x_herm) * gaussian_weight + 0.05 * cp.random.normal(size=x_herm.shape)
462

463
herm_coeffs = poly.hermite_e.hermefit(x_herm, y_weighted, 6)
464
y_herm_fitted = poly.hermite_e.hermeval(x_herm, herm_coeffs)
465
print("Hermite fit MSE:", cp.mean((y_weighted - y_herm_fitted) ** 2))
466
```
467

468
### Multi-dimensional Polynomial Operations
469

470
```python
471
# 2D polynomial fitting
472
x2d = cp.random.rand(200) * 4 - 2
473
y2d = cp.random.rand(200) * 4 - 2
474
z_true = 2*x2d**2 + 3*y2d**2 - x2d*y2d + x2d + 2*y2d + 1
475
z_noisy = z_true + 0.1 * cp.random.normal(size=z_true.shape)
476

477
# Fit 2D polynomial surface
478
degree_2d = [2, 2]  # Degree in x and y directions
479
coeffs_2d = poly.polynomial.polyfit2d(x2d, y2d, z_noisy, degree_2d)
480

481
# Evaluate on regular grid
482
x_grid = cp.linspace(-2, 2, 50)
483
y_grid = cp.linspace(-2, 2, 50)
484
X_grid, Y_grid = cp.meshgrid(x_grid, y_grid)
485
Z_fitted = poly.polynomial.polyval2d(X_grid, Y_grid, coeffs_2d)
486

487
print("2D polynomial coefficients shape:", coeffs_2d.shape)
488
print("Fitted surface shape:", Z_fitted.shape)
489

490
# 3D polynomial evaluation
491
x3d = cp.linspace(-1, 1, 20)
492
y3d = cp.linspace(-1, 1, 20) 
493
z3d = cp.linspace(-1, 1, 20)
494
X3d, Y3d, Z3d = cp.meshgrid(x3d, y3d, z3d, indexing='ij')
495

496
# Simple 3D polynomial coefficients
497
coeffs_3d = cp.zeros((3, 3, 3))
498
coeffs_3d[1, 0, 0] = 1  # x term
499
coeffs_3d[0, 1, 0] = 1  # y term  
500
coeffs_3d[0, 0, 1] = 1  # z term
501
coeffs_3d[2, 0, 0] = 0.5  # x^2 term
502

503
# Evaluate 3D polynomial
504
result_3d = poly.polynomial.polyval3d(X3d, Y3d, Z3d, coeffs_3d)
505
print("3D polynomial result shape:", result_3d.shape)
506
```
507

508
### Root Finding and Analysis
509

510
```python
511
# Complex polynomial root finding
512
# Polynomial: x^4 - 2x^3 + 3x^2 - 2x + 1
513
coeffs = [1, -2, 3, -2, 1]
514
roots = cp.roots(coeffs)
515

516
print("Polynomial coefficients:", coeffs)
517
print("Roots:", roots)
518
print("Root types:", [type(r).__name__ for r in roots])
519

520
# Verify roots by evaluation
521
for i, root in enumerate(roots):
522
    value = cp.polyval(coeffs, root)
523
    print(f"Root {i}: {root}, P(root) = {value}")
524

525
# Create polynomial from roots
526
reconstructed_coeffs = cp.poly(roots)
527
print("Reconstructed coefficients:", reconstructed_coeffs)
528
print("Original coefficients (normalized):", cp.array(coeffs) / coeffs[0])
529

530
# Companion matrix approach for root finding
531
companion = poly.polynomial.polycompanion([1, -2, 3, -2, 1])
532
eigenvals = cp.linalg.eigvals(companion)
533
print("Roots from companion matrix:", eigenvals)
534

535
# Analysis of polynomial behavior
536
def analyze_polynomial(coeffs, x_range=(-5, 5), n_points=1000):
537
    """Analyze polynomial behavior over a range."""
538
    x = cp.linspace(x_range[0], x_range[1], n_points)
539
    y = cp.polyval(coeffs, x)
540
    
541
    # Find extrema (approximate)
542
    dy = cp.gradient(y, x[1] - x[0])
543
    extrema_indices = cp.where(cp.abs(dy) < cp.std(dy) * 0.1)[0]
544
    
545
    analysis = {
546
        'min_value': cp.min(y),
547
        'max_value': cp.max(y),
548
        'range': cp.max(y) - cp.min(y),
549
        'approximate_extrema': len(extrema_indices),
550
        'zeros_in_range': cp.sum(cp.abs(y) < cp.std(y) * 0.01)
551
    }
552
    
553
    return analysis
554

555
# Analyze different polynomials
556
test_polynomials = [
557
    [1, 0, -1],           # x^2 - 1
558
    [1, -3, 3, -1],       # (x-1)^3  
559
    [1, 0, 0, 0, -1],     # x^4 - 1
560
    [1, -2, 1]            # (x-1)^2
561
]
562

563
for i, coeffs in enumerate(test_polynomials):
564
    analysis = analyze_polynomial(coeffs)
565
    print(f"\nPolynomial {i+1}: {coeffs}")
566
    for key, value in analysis.items():
567
        print(f"  {key}: {value}")
568
```
569

570
### Specialized Polynomial Applications
571

572
```python
573
# Interpolation using different polynomial bases
574
def compare_interpolation_methods(x_data, y_data, eval_points):
575
    """Compare different polynomial interpolation methods."""
576
    n_points = len(x_data)
577
    degree = n_points - 1
578
    
579
    methods = {}
580
    
581
    # Standard polynomial interpolation
582
    std_coeffs = cp.polyfit(x_data, y_data, degree)
583
    methods['Standard'] = cp.polyval(std_coeffs, eval_points)
584
    
585
    # Chebyshev interpolation
586
    cheb_coeffs = poly.chebyshev.chebfit(x_data, y_data, degree)
587
    methods['Chebyshev'] = poly.chebyshev.chebval(eval_points, cheb_coeffs)
588
    
589
    # Legendre interpolation
590
    leg_coeffs = poly.legendre.legfit(x_data, y_data, degree)
591
    methods['Legendre'] = poly.legendre.legval(eval_points, leg_coeffs)
592
    
593
    return methods
594

595
# Test interpolation with Runge function (challenging for high-degree polynomials)
596
def runge_function(x):
597
    return 1 / (1 + 25 * x**2)
598

599
# Interpolation points
600
n_interp = 10
601
x_interp = cp.linspace(-1, 1, n_interp)
602
y_interp = runge_function(x_interp)
603

604
# Evaluation points
605
x_eval = cp.linspace(-1, 1, 200)
606
y_true = runge_function(x_eval)
607

608
# Compare methods
609
methods = compare_interpolation_methods(x_interp, y_interp, x_eval)
610

611
print("Interpolation comparison (MSE vs true function):")
612
for method_name, y_approx in methods.items():
613
    mse = cp.mean((y_true - y_approx) ** 2)
614
    print(f"  {method_name}: {mse:.6f}")
615

616
# Orthogonal polynomial series expansion
617
def orthogonal_series_approximation(func, x_range, n_terms, poly_type='chebyshev'):
618
    """Approximate a function using orthogonal polynomial series."""
619
    x = cp.linspace(x_range[0], x_range[1], 1000)
620
    y = func(x)
621
    
622
    if poly_type == 'chebyshev':
623
        coeffs = poly.chebyshev.chebfit(x, y, n_terms-1)
624
        y_approx = poly.chebyshev.chebval(x, coeffs)
625
    elif poly_type == 'legendre':
626
        coeffs = poly.legendre.legfit(x, y, n_terms-1)
627
        y_approx = poly.legendre.legval(x, coeffs)
628
    elif poly_type == 'hermite':
629
        coeffs = poly.hermite.hermfit(x, y, n_terms-1)
630
        y_approx = poly.hermite.hermval(x, coeffs)
631
    
632
    mse = cp.mean((y - y_approx) ** 2)
633
    return coeffs, y_approx, mse
634

635
# Test with exponential function
636
exponential = lambda x: cp.exp(x)
637
x_range = (-2, 2)
638

639
for n_terms in [5, 10, 15, 20]:
640
    for poly_type in ['chebyshev', 'legendre']:
641
        coeffs, y_approx, mse = orthogonal_series_approximation(
642
            exponential, x_range, n_terms, poly_type
643
        )
644
        print(f"{poly_type.capitalize()} ({n_terms} terms): MSE = {mse:.8f}")
645

646
# Polynomial differentiation and integration
647
p = cp.poly1d([1, -2, 1, 3])  # x^3 - 2x^2 + x + 3
648
print(f"Original polynomial: {p}")
649

650
# Multiple derivatives
651
for k in range(1, 4):
652
    p_deriv = p.deriv(k)
653
    print(f"Derivative {k}: {p_deriv}")
654

655
# Multiple integrals
656
for k in range(1, 3):
657
    p_integ = p.integ(k)
658
    print(f"Integral {k}: {p_integ}")
659

660
# Definite integration
661
x_vals = cp.linspace(0, 2, 100)
662
y_vals = p(x_vals)
663
numerical_integral = cp.trapz(y_vals, x_vals)
664
analytical_integral = p.integ()(2) - p.integ()(0)
665
print(f"Numerical integral (0 to 2): {numerical_integral}")
666
print(f"Analytical integral (0 to 2): {analytical_integral}")
667
```
668

669
### Performance Optimization
670

671
```python
672
# Efficient polynomial evaluation using Horner's method
673
def horner_evaluation(coeffs, x):
674
    """Evaluate polynomial using Horner's method for numerical stability."""
675
    result = cp.zeros_like(x, dtype=cp.result_type(coeffs.dtype, x.dtype))
676
    for coeff in coeffs:
677
        result = result * x + coeff
678
    return result
679

680
# Compare with standard evaluation
681
coeffs_high_degree = cp.random.rand(50)  # High-degree polynomial
682
x_test = cp.linspace(-10, 10, 10000)
683

684
# Time both methods
685
import time
686

687
start_time = time.time()
688
result_standard = cp.polyval(coeffs_high_degree, x_test)
689
standard_time = time.time() - start_time
690

691
start_time = time.time()
692
result_horner = horner_evaluation(coeffs_high_degree, x_test)
693
horner_time = time.time() - start_time
694

695
print(f"Standard evaluation time: {standard_time:.4f}s")
696
print(f"Horner evaluation time: {horner_time:.4f}s")
697
print(f"Results match: {cp.allclose(result_standard, result_horner)}")
698

699
# Vectorized polynomial operations for multiple polynomials
700
def batch_polynomial_evaluation(poly_coeffs_list, x_values):
701
    """Evaluate multiple polynomials efficiently."""
702
    results = []
703
    for coeffs in poly_coeffs_list:
704
        results.append(cp.polyval(coeffs, x_values))
705
    return cp.stack(results)
706

707
# Create batch of polynomials
708
n_polynomials = 100
709
poly_list = [cp.random.rand(cp.random.randint(3, 10)) for _ in range(n_polynomials)]
710
x_batch = cp.linspace(-5, 5, 1000)
711

712
# Batch evaluation
713
batch_results = batch_polynomial_evaluation(poly_list, x_batch)
714
print(f"Batch evaluation shape: {batch_results.shape}")
715
print(f"Memory usage: {batch_results.nbytes / 1024**2:.2f} MB")
716

717
# Memory-efficient polynomial fitting for large datasets
718
def streaming_polynomial_fit(data_generator, degree, max_batch_size=10000):
719
    """Fit polynomial to streaming data in batches."""
720
    x_accumulator = []
721
    y_accumulator = []
722
    
723
    for x_batch, y_batch in data_generator:
724
        x_accumulator.extend(x_batch)
725
        y_accumulator.extend(y_batch)
726
        
727
        # Process when batch is full
728
        if len(x_accumulator) >= max_batch_size:
729
            x_array = cp.array(x_accumulator)
730
            y_array = cp.array(y_accumulator)
731
            
732
            # Fit polynomial to current batch
733
            coeffs = cp.polyfit(x_array, y_array, degree)
734
            yield coeffs
735
            
736
            # Clear accumulators
737
            x_accumulator = []
738
            y_accumulator = []
739
    
740
    # Process remaining data
741
    if x_accumulator:
742
        x_array = cp.array(x_accumulator)
743
        y_array = cp.array(y_accumulator)
744
        coeffs = cp.polyfit(x_array, y_array, degree)
745
        yield coeffs
746

747
# Example streaming data generator
748
def synthetic_data_generator(n_batches=10, batch_size=5000):
749
    for i in range(n_batches):
750
        x = cp.random.rand(batch_size) * 10 - 5
751
        y = 2*x**3 - 3*x**2 + x + 1 + 0.1*cp.random.randn(batch_size)
752
        yield cp.asnumpy(x), cp.asnumpy(y)
753

754
# Process streaming data
755
coeffs_list = list(streaming_polynomial_fit(synthetic_data_generator(), degree=3))
756
print(f"Processed {len(coeffs_list)} batches")
757
print("Average coefficients:", cp.mean(cp.stack(coeffs_list), axis=0))
758
```
759

760
Polynomial operations in CuPy provide comprehensive mathematical tools for polynomial manipulation, fitting, evaluation, and analysis, supporting both legacy interfaces and modern orthogonal polynomial bases with GPU acceleration for high-performance computational mathematics and numerical analysis applications.