0
# Mathematical Utilities
1

2
Comprehensive mathematical functions including interpolation, algebra operations, distance calculations, array manipulation, geometric computations, statistics, optimization, and mathematical constants for colour science applications.
3

4
## Capabilities
5

6
### Interpolation Classes and Functions
7

8
Advanced interpolation capabilities with multiple algorithms and kernel-based methods.
9

10
#### Interpolator Classes
11

12
```python { .api }
13
class KernelInterpolator:
14
    """
15
    Kernel-based interpolation of a 1-D function using convolution operations.
16
    
17
    Parameters:
18
    - x: independent variable values
19
    - y: dependent variable values to interpolate
20
    - window: width of the window in samples on each side (default: 3)
21
    - kernel: kernel function to use (default: kernel_lanczos)
22
    - kernel_kwargs: arguments for the kernel function
23
    - padding_kwargs: arguments for numpy.pad when extending data
24
    - dtype: data type for internal conversions
25
    """
26
    def __init__(self, x: ArrayLike, y: ArrayLike, window: float = 3, kernel: Callable = None, kernel_kwargs: dict = None, padding_kwargs: dict = None, dtype: Type = None): ...
27
    
28
    def __call__(self, x: ArrayLike) -> NDArrayFloat: ...
29
    
30
    @property
31
    def x(self) -> NDArrayFloat: ...
32
    @property 
33
    def y(self) -> NDArrayFloat: ...
34
    @property
35
    def window(self) -> float: ...
36
    @property
37
    def kernel(self) -> Callable: ...
38
    @property
39
    def kernel_kwargs(self) -> dict: ...
40
    @property
41
    def padding_kwargs(self) -> dict: ...
42

43
class LinearInterpolator:
44
    """
45
    Linear interpolation of a 1-D function using numpy.interp.
46
    
47
    Parameters:
48
    - x: independent variable values
49
    - y: dependent variable values to interpolate
50
    - dtype: data type for internal conversions
51
    """
52
    def __init__(self, x: ArrayLike, y: ArrayLike, dtype: Type = None): ...
53
    
54
    def __call__(self, x: ArrayLike) -> NDArrayFloat: ...
55
    
56
    @property
57
    def x(self) -> NDArrayFloat: ...
58
    @property
59
    def y(self) -> NDArrayFloat: ...
60

61
class SpragueInterpolator:
62
    """
63
    Fifth-order polynomial interpolation using Sprague (1880) method.
64
    Recommended by CIE for uniformly spaced independent variables.
65
    
66
    Parameters:
67
    - x: independent variable values
68
    - y: dependent variable values to interpolate (minimum 6 points required)
69
    - dtype: data type for internal conversions
70
    """
71
    def __init__(self, x: ArrayLike, y: ArrayLike, dtype: Type = None): ...
72
    
73
    def __call__(self, x: ArrayLike) -> NDArrayFloat: ...
74
    
75
    @property
76
    def x(self) -> NDArrayFloat: ...
77
    @property
78
    def y(self) -> NDArrayFloat: ...
79

80
class CubicSplineInterpolator:
81
    """
82
    Cubic spline interpolation using scipy.interpolate.interp1d.
83
    
84
    Wrapper around scipy's cubic spline interpolation.
85
    """
86
    def __init__(self, *args, **kwargs): ...
87

88
class PchipInterpolator:
89
    """
90
    Piecewise Cubic Hermite Interpolating Polynomial interpolation.
91
    
92
    Wrapper around scipy.interpolate.PchipInterpolator.
93
    """
94
    def __init__(self, x: ArrayLike, y: ArrayLike, *args, **kwargs): ...
95
    
96
    @property
97
    def y(self) -> NDArrayFloat: ...
98

99
class NearestNeighbourInterpolator(KernelInterpolator):
100
    """
101
    Nearest-neighbour interpolation using kernel-based approach.
102
    
103
    Specialized KernelInterpolator using nearest-neighbour kernel.
104
    """
105
    def __init__(self, *args, **kwargs): ...
106

107
class NullInterpolator:
108
    """
109
    Null interpolation returning existing values within tolerances.
110
    
111
    Parameters:
112
    - x: independent variable values
113
    - y: dependent variable values
114
    - absolute_tolerance: absolute tolerance for matching
115
    - relative_tolerance: relative tolerance for matching
116
    - default: default value for points outside tolerances
117
    - dtype: data type for internal conversions
118
    """
119
    def __init__(self, x: ArrayLike, y: ArrayLike, absolute_tolerance: float = 1e-10, relative_tolerance: float = 1e-5, default: float = np.nan, dtype: Type = None): ...
120
    
121
    def __call__(self, x: ArrayLike) -> NDArrayFloat: ...
122
    
123
    @property
124
    def x(self) -> NDArrayFloat: ...
125
    @property
126
    def y(self) -> NDArrayFloat: ...
127
    @property
128
    def absolute_tolerance(self) -> float: ...
129
    @property
130
    def relative_tolerance(self) -> float: ...
131
    @property
132
    def default(self) -> float: ...
133
```
134

135
#### Interpolation Kernel Functions
136

137
Mathematical kernels for interpolation algorithms.
138

139
```python { .api }
140
def kernel_nearest_neighbour(x: ArrayLike) -> NDArrayFloat:
141
    """
142
    Nearest-neighbour kernel returning 1 for |x| < 0.5, 0 otherwise.
143
    
144
    Parameters:
145
    - x: samples at which to evaluate the kernel
146
    
147
    Returns:
148
    Kernel values evaluated at given samples
149
    """
150

151
def kernel_linear(x: ArrayLike) -> NDArrayFloat:
152
    """
153
    Linear kernel returning (1 - |x|) for |x| < 1, 0 otherwise.
154
    
155
    Parameters:
156
    - x: samples at which to evaluate the kernel
157
    
158
    Returns:
159
    Kernel values evaluated at given samples
160
    """
161

162
def kernel_sinc(x: ArrayLike, a: float = 3) -> NDArrayFloat:
163
    """
164
    Sinc kernel for high-quality interpolation.
165
    
166
    Parameters:
167
    - x: samples at which to evaluate the kernel
168
    - a: kernel size parameter (must be >= 1)
169
    
170
    Returns:
171
    Kernel values evaluated at given samples
172
    """
173

174
def kernel_lanczos(x: ArrayLike, a: float = 3) -> NDArrayFloat:
175
    """
176
    Lanczos kernel for high-quality resampling.
177
    
178
    Parameters:
179
    - x: samples at which to evaluate the kernel
180
    - a: kernel size parameter (must be >= 1)
181
    
182
    Returns:
183
    Kernel values evaluated at given samples
184
    """
185

186
def kernel_cardinal_spline(x: ArrayLike, a: float = 0.5, b: float = 0.0) -> NDArrayFloat:
187
    """
188
    Cardinal spline kernel with configurable parameters.
189
    
190
    Notable parameterizations:
191
    - Catmull-Rom: (a=0.5, b=0)
192
    - Cubic B-Spline: (a=0, b=1) 
193
    - Mitchell-Netravalli: (a=1/3, b=1/3)
194
    
195
    Parameters:
196
    - x: samples at which to evaluate the kernel
197
    - a: control parameter a
198
    - b: control parameter b
199
    
200
    Returns:
201
    Kernel values evaluated at given samples
202
    """
203
```
204

205
#### Interpolation Support Functions
206

207
```python { .api }
208
def lagrange_coefficients(r: float, n: int = 4) -> NDArrayFloat:
209
    """
210
    Compute Lagrange coefficients for interpolation.
211
    
212
    Parameters:
213
    - r: point to get coefficients at
214
    - n: degree of the Lagrange coefficients
215
    
216
    Returns:
217
    Lagrange coefficients array
218
    """
219

220
def table_interpolation_trilinear(V_xyz: ArrayLike, table: ArrayLike) -> NDArrayFloat:
221
    """
222
    Trilinear interpolation using 4-dimensional lookup table.
223
    
224
    Parameters:
225
    - V_xyz: values to interpolate
226
    - table: 4-dimensional interpolation table (NxNxNx3)
227
    
228
    Returns:
229
    Interpolated values
230
    """
231

232
def table_interpolation_tetrahedral(V_xyz: ArrayLike, table: ArrayLike) -> NDArrayFloat:
233
    """
234
    Tetrahedral interpolation using 4-dimensional lookup table.
235
    
236
    Parameters:
237
    - V_xyz: values to interpolate  
238
    - table: 4-dimensional interpolation table (NxNxNx3)
239
    
240
    Returns:
241
    Interpolated values
242
    """
243

244
def table_interpolation(V_xyz: ArrayLike, table: ArrayLike, method: str = "Trilinear") -> NDArrayFloat:
245
    """
246
    Generic table interpolation with selectable method.
247
    
248
    Parameters:
249
    - V_xyz: values to interpolate
250
    - table: 4-dimensional interpolation table (NxNxNx3) 
251
    - method: interpolation method ("Trilinear" or "Tetrahedral")
252
    
253
    Returns:
254
    Interpolated values
255
    """
256

257
# Available table interpolation methods
258
TABLE_INTERPOLATION_METHODS = {
259
    "Trilinear": table_interpolation_trilinear,
260
    "Tetrahedral": table_interpolation_tetrahedral,
261
}
262
```
263

264
### Function Extrapolation
265

266
Extrapolation capabilities for extending functions beyond their defined domain.
267

268
```python { .api }
269
class Extrapolator:
270
    """
271
    Extrapolate 1-D functions beyond their interpolation range.
272
    
273
    Supports linear and constant extrapolation methods with customizable boundary values.
274
    
275
    Parameters:
276
    - interpolator: interpolator object with x and y properties
277
    - method: extrapolation method ("Linear" or "Constant")
278
    - left: value to return for x < xi[0] (overrides method)
279
    - right: value to return for x > xi[-1] (overrides method)
280
    - dtype: data type for internal conversions
281
    """
282
    def __init__(self, interpolator: ProtocolInterpolator = None, method: str = "Linear", left: Real = None, right: Real = None, dtype: Type = None): ...
283
    
284
    def __call__(self, x: ArrayLike) -> NDArrayFloat: ...
285
    
286
    @property
287
    def interpolator(self) -> ProtocolInterpolator: ...
288
    @property
289
    def method(self) -> str: ...
290
    @property
291
    def left(self) -> Real: ...
292
    @property
293
    def right(self) -> Real: ...
294
```
295

296
### Safe Algebra Operations
297

298
Robust mathematical operations with configurable error handling.
299

300
#### Safe Division
301

302
```python { .api }
303
def sdiv(a: ArrayLike, b: ArrayLike) -> NDArrayFloat:
304
    """
305
    Safe division handling zero-division gracefully with configurable modes.
306
    
307
    Available modes (configured via sdiv_mode):
308
    - "Numpy": standard numpy zero-division handling
309
    - "Ignore": zero-division occurs silently
310
    - "Warning": zero-division with warnings
311
    - "Raise": raise exception on zero-division
312
    - "Ignore Zero Conversion": silent with NaN/inf converted to zero
313
    - "Warning Zero Conversion": warning with NaN/inf converted to zero
314
    - "Ignore Limit Conversion": silent with NaN/inf converted to limits
315
    - "Warning Limit Conversion": warning with NaN/inf converted to limits
316
    
317
    Parameters:
318
    - a: numerator array
319
    - b: denominator array
320
    
321
    Returns:
322
    Safely divided result
323
    """
324

325
def get_sdiv_mode() -> str:
326
    """Get current safe division mode."""
327

328
def set_sdiv_mode(mode: str) -> None:
329
    """Set safe division mode."""
330

331
class sdiv_mode:
332
    """
333
    Context manager for temporarily setting safe division mode.
334
    
335
    Parameters:
336
    - mode: temporary safe division mode
337
    """
338
    def __init__(self, mode: str = None): ...
339
    def __enter__(self) -> "sdiv_mode": ...
340
    def __exit__(self, *args): ...
341
```
342

343
#### Safe Power Functions
344

345
```python { .api }
346
def spow(a: ArrayLike, p: ArrayLike) -> NDArrayFloat:
347
    """
348
    Safe power function: sign(a) * |a|^p
349
    
350
    Avoids NaN generation when base is negative and exponent is fractional.
351
    
352
    Parameters:
353
    - a: base array
354
    - p: power/exponent array
355
    
356
    Returns:
357
    Safely computed power result
358
    """
359

360
def is_spow_enabled() -> bool:
361
    """Check if safe power function is enabled."""
362

363
def set_spow_enable(enable: bool) -> None:
364
    """Enable or disable safe power function."""
365

366
class spow_enable:
367
    """
368
    Context manager for temporarily setting safe power function state.
369
    
370
    Parameters:
371
    - enable: whether to enable safe power function
372
    """
373
    def __init__(self, enable: bool): ...
374
    def __enter__(self) -> "spow_enable": ...
375
    def __exit__(self, *args): ...
376
```
377

378
### Vector and Matrix Operations
379

380
Linear algebra operations for colour science computations.
381

382
```python { .api }
383
def normalise_vector(a: ArrayLike) -> NDArrayFloat:
384
    """
385
    Normalize vector to unit length using L2 norm.
386
    
387
    Parameters:
388
    - a: vector to normalize
389
    
390
    Returns:
391
    Normalized vector
392
    """
393

394
def normalise_maximum(a: ArrayLike, axis: int = None, factor: float = 1, clip: bool = True) -> NDArrayFloat:
395
    """
396
    Normalize array values by maximum value with optional clipping.
397
    
398
    Parameters:
399
    - a: array to normalize
400
    - axis: normalization axis
401
    - factor: normalization factor
402
    - clip: whether to clip values to domain [0, factor]
403
    
404
    Returns:
405
    Maximum-normalized array
406
    """
407

408
def vecmul(m: ArrayLike, v: ArrayLike) -> NDArrayFloat:
409
    """
410
    Batched matrix-vector multiplication with broadcasting.
411
    
412
    Equivalent to np.einsum('...ij,...j->...i', m, v).
413
    
414
    Parameters:
415
    - m: matrix array 
416
    - v: vector array
417
    
418
    Returns:
419
    Matrix-vector multiplication result
420
    """
421

422
def is_identity(a: ArrayLike) -> bool:
423
    """
424
    Test whether array is an identity matrix.
425
    
426
    Parameters:
427
    - a: array to test
428
    
429
    Returns:
430
    True if array is identity matrix
431
    """
432

433
def eigen_decomposition(a: ArrayLike, eigen_w_v_count: int = None, descending_order: bool = True, covariance_matrix: bool = False) -> tuple:
434
    """
435
    Compute eigenvalues and eigenvectors with optional ordering.
436
    
437
    Parameters:
438
    - a: array to decompose
439
    - eigen_w_v_count: number of eigenvalues/eigenvectors to return
440
    - descending_order: whether to return in descending eigenvalue order
441
    - covariance_matrix: whether to use covariance matrix A^T * A
442
    
443
    Returns:
444
    Tuple of (eigenvalues, eigenvectors)
445
    """
446
```
447

448
### Distance Metrics
449

450
Distance calculations for colour difference and similarity analysis.
451

452
```python { .api }
453
def euclidean_distance(a: ArrayLike, b: ArrayLike) -> NDArrayFloat:
454
    """
455
    Euclidean (L2) distance between point arrays.
456
    
457
    For 2D: sqrt((xa - xb)² + (ya - yb)²)
458
    
459
    Parameters:
460
    - a: first point array
461
    - b: second point array
462
    
463
    Returns:
464
    Euclidean distance
465
    """
466

467
def manhattan_distance(a: ArrayLike, b: ArrayLike) -> NDArrayFloat:
468
    """
469
    Manhattan (L1, City-Block) distance between point arrays.
470
    
471
    For 2D: |xa - xb| + |ya - yb|
472
    
473
    Parameters:
474
    - a: first point array
475
    - b: second point array
476
    
477
    Returns:
478
    Manhattan distance
479
    """
480
```
481

482
### Linear Transformations and Interpolation
483

484
Mathematical transformations and interpolation functions.
485

486
```python { .api }
487
def linear_conversion(a: ArrayLike, old_range: ArrayLike, new_range: ArrayLike) -> NDArrayFloat:
488
    """
489
    Linear conversion between value ranges.
490
    
491
    Parameters:
492
    - a: array to convert
493
    - old_range: [old_min, old_max] current range
494
    - new_range: [new_min, new_max] target range
495
    
496
    Returns:
497
    Linearly converted values
498
    """
499

500
def linstep_function(x: ArrayLike, a: ArrayLike = 0, b: ArrayLike = 1, clip: bool = False) -> NDArrayFloat:
501
    """
502
    Linear interpolation between two values: (1-x)*a + x*b
503
    
504
    Parameters:
505
    - x: interpolation parameter
506
    - a: start value  
507
    - b: end value
508
    - clip: whether to clip result to [a, b]
509
    
510
    Returns:
511
    Linearly interpolated result
512
    """
513

514
# Alias for linstep_function
515
lerp = linstep_function
516

517
def smoothstep_function(x: ArrayLike, a: ArrayLike = 0, b: ArrayLike = 1, clip: bool = False) -> NDArrayFloat:
518
    """
519
    Smooth sigmoid-like interpolation function.
520
    
521
    Uses polynomial: t² * (3 - 2t) where t = (x-a)/(b-a)
522
    
523
    Parameters:
524
    - x: input array
525
    - a: low input domain limit (left edge)
526
    - b: high input domain limit (right edge)
527
    - clip: whether to scale and clip input to [a, b]
528
    
529
    Returns:
530
    Smoothly interpolated result
531
    """
532

533
# Alias for smoothstep_function  
534
smooth = smoothstep_function
535
```
536

537
### Coordinate System Transformations
538

539
Geometric transformations between coordinate systems.
540

541
```python { .api }
542
def cartesian_to_spherical(a: ArrayLike) -> NDArrayFloat:
543
    """
544
    Transform Cartesian coordinates (x,y,z) to spherical (ρ,θ,φ).
545
    
546
    Parameters:
547
    - a: Cartesian coordinates [x, y, z]
548
    
549
    Returns:
550
    Spherical coordinates [ρ, θ, φ] where:
551
    - ρ: radial distance [0, +∞)
552
    - θ: inclination [0, π] radians
553
    - φ: azimuth [-π, π] radians
554
    """
555

556
def spherical_to_cartesian(a: ArrayLike) -> NDArrayFloat:
557
    """
558
    Transform spherical coordinates (ρ,θ,φ) to Cartesian (x,y,z).
559
    
560
    Parameters:
561
    - a: Spherical coordinates [ρ, θ, φ]
562
    
563
    Returns:
564
    Cartesian coordinates [x, y, z]
565
    """
566

567
def cartesian_to_polar(a: ArrayLike) -> NDArrayFloat:
568
    """
569
    Transform Cartesian coordinates (x,y) to polar (ρ,φ).
570
    
571
    Parameters:
572
    - a: Cartesian coordinates [x, y]
573
    
574
    Returns:
575
    Polar coordinates [ρ, φ] where:
576
    - ρ: radial coordinate [0, +∞)
577
    - φ: angular coordinate [-π, π] radians
578
    """
579

580
def polar_to_cartesian(a: ArrayLike) -> NDArrayFloat:
581
    """
582
    Transform polar coordinates (ρ,φ) to Cartesian (x,y).
583
    
584
    Parameters:
585
    - a: Polar coordinates [ρ, φ]
586
    
587
    Returns:
588
    Cartesian coordinates [x, y]
589
    """
590

591
def cartesian_to_cylindrical(a: ArrayLike) -> NDArrayFloat:
592
    """
593
    Transform Cartesian coordinates (x,y,z) to cylindrical (ρ,φ,z).
594
    
595
    Parameters:
596
    - a: Cartesian coordinates [x, y, z]
597
    
598
    Returns:
599
    Cylindrical coordinates [ρ, φ, z] where:
600
    - ρ: radial distance [0, +∞) 
601
    - φ: azimuth [-π, π] radians
602
    - z: height [0, +∞)
603
    """
604

605
def cylindrical_to_cartesian(a: ArrayLike) -> NDArrayFloat:
606
    """
607
    Transform cylindrical coordinates (ρ,φ,z) to Cartesian (x,y,z).
608
    
609
    Parameters:
610
    - a: Cylindrical coordinates [ρ, φ, z]
611
    
612
    Returns:
613
    Cartesian coordinates [x, y, z]
614
    """
615
```
616

617
### Statistics and Regression
618

619
Statistical analysis and optimization functions.
620

621
```python { .api }
622
def least_square_mapping_MoorePenrose(y: ArrayLike, x: ArrayLike) -> NDArrayFloat:
623
    """
624
    Least-squares mapping using Moore-Penrose inverse.
625
    
626
    Computes optimal linear mapping from y to x variables.
627
    
628
    Parameters:
629
    - y: dependent variable (known values)
630
    - x: independent variable (target values)
631
    
632
    Returns:
633
    Least-squares mapping matrix
634
    """
635
```
636

637
### Random Number Generation
638

639
Utilities for generating random data for testing and Monte Carlo methods.
640

641
```python { .api }
642
def random_triplet_generator(size: int, limits: ArrayLike = [[0,1], [0,1], [0,1]], random_state: np.random.RandomState = None) -> NDArrayFloat:
643
    """
644
    Generate random triplets within specified limits.
645
    
646
    Parameters:
647
    - size: number of triplets to generate
648
    - limits: [[x_min, x_max], [y_min, y_max], [z_min, z_max]] limits for each axis
649
    - random_state: random number generator state for reproducibility
650
    
651
    Returns:
652
    Array of random triplets shape (size, 3)
653
    """
654

655
# Global random state for reproducible results
656
RANDOM_STATE = np.random.RandomState()
657
```
658

659
### Mathematical Constants
660

661
Physical and mathematical constants for colour science calculations.
662

663
```python { .api }
664
# CIE Constants
665
CONSTANT_K_M: float  # CIE luminous efficacy constant (683 lm/W)
666
CONSTANT_KP_M: float  # CIE scotopic luminous efficacy constant (1700 lm/W)
667

668
# CODATA Physical Constants  
669
CONSTANT_AVOGADRO: float      # Avogadro constant (6.02214076e23 mol⁻¹)
670
CONSTANT_BOLTZMANN: float     # Boltzmann constant (1.380649e-23 J/K)
671
CONSTANT_LIGHT_SPEED: float   # Speed of light in vacuum (299792458 m/s)
672
CONSTANT_PLANCK: float        # Planck constant (6.62607015e-34 J⋅s)
673

674
# Numerical Constants
675
EPSILON: float                    # Machine epsilon for floating point comparisons
676
DTYPE_INT_DEFAULT: type          # Default integer data type (np.int64)  
677
DTYPE_FLOAT_DEFAULT: type        # Default float data type (np.float64)
678
TOLERANCE_ABSOLUTE_DEFAULT: float # Default absolute tolerance (1e-10)
679
TOLERANCE_RELATIVE_DEFAULT: float # Default relative tolerance (1e-5)
680
TOLERANCE_ABSOLUTE_TESTS: float  # Absolute tolerance for tests (1e-7)
681
TOLERANCE_RELATIVE_TESTS: float  # Relative tolerance for tests (1e-7)
682

683
# Utility Constants
684
THRESHOLD_INTEGER: float         # Threshold for integer detection
685
PATTERN_FLOATING_POINT_NUMBER: str # Regex pattern for floating point numbers
686
```
687

688
## Usage Examples
689

690
### Basic Interpolation
691

692
```python
693
import numpy as np
694
from colour.algebra import LinearInterpolator, SpragueInterpolator
695

696
# Linear interpolation
697
x = np.array([1, 2, 3, 4, 5])
698
y = np.array([1, 4, 9, 16, 25])  # x²
699
interp = LinearInterpolator(x, y)
700

701
# Interpolate at intermediate points
702
result = interp([1.5, 2.5, 3.5])
703
print(result)  # [2.5, 6.5, 12.5]
704

705
# High-quality Sprague interpolation
706
sprague = SpragueInterpolator(x, y)
707
result = sprague([1.5, 2.5, 3.5])
708
print(result)  # More accurate cubic interpolation
709
```
710

711
### Safe Mathematical Operations
712

713
```python
714
from colour.algebra import sdiv, spow, sdiv_mode
715

716
# Safe division handling zero denominators
717
a = np.array([1, 2, 3])
718
b = np.array([2, 0, 4])
719

720
# Default mode converts inf/nan to zero
721
result = sdiv(a, b)
722
print(result)  # [0.5, 0.0, 0.75]
723

724
# Change division behavior
725
with sdiv_mode("Warning"):
726
    result = sdiv(a, b)  # Will warn about zero division
727

728
# Safe power for negative bases
729
result = spow(-2, 0.5)  # Returns -√2 instead of NaN
730
print(result)  # -1.4142135...
731
```
732

733
### Vector Operations and Distance Metrics
734

735
```python
736
from colour.algebra import normalise_vector, euclidean_distance, vecmul
737

738
# Normalize vectors
739
vector = np.array([3, 4, 5])
740
normalized = normalise_vector(vector)
741
print(normalized)  # [0.424, 0.566, 0.707] (unit length)
742

743
# Calculate distances
744
point_a = np.array([1, 2, 3])
745
point_b = np.array([4, 6, 8])
746
distance = euclidean_distance(point_a, point_b)
747
print(distance)  # 7.071...
748

749
# Matrix-vector multiplication with broadcasting
750
matrix = np.random.random((100, 3, 3))  # 100 3x3 matrices
751
vectors = np.random.random((100, 3))     # 100 3D vectors
752
result = vecmul(matrix, vectors)         # 100 transformed vectors
753
print(result.shape)  # (100, 3)
754
```
755

756
### Coordinate Transformations
757

758
```python
759
from colour.algebra import (
760
    cartesian_to_spherical, spherical_to_cartesian,
761
    cartesian_to_polar, polar_to_cartesian
762
)
763

764
# 3D coordinate transformations
765
cartesian = np.array([3, 1, 6])
766
spherical = cartesian_to_spherical(cartesian)
767
print(spherical)  # [6.782, 0.485, 0.322] (ρ, θ, φ)
768

769
# Convert back
770
back_to_cartesian = spherical_to_cartesian(spherical)
771
print(back_to_cartesian)  # [3, 1, 6] (approximately)
772

773
# 2D coordinate transformations  
774
cartesian_2d = np.array([3, 4])
775
polar = cartesian_to_polar(cartesian_2d)
776
print(polar)  # [5.0, 0.927] (ρ, φ)
777
```
778

779
### Function Extrapolation
780

781
```python
782
from colour.algebra import LinearInterpolator, Extrapolator
783

784
# Create base interpolator
785
x = np.array([2, 3, 4])
786
y = np.array([4, 9, 16])
787
interp = LinearInterpolator(x, y)
788

789
# Linear extrapolation
790
extrap = Extrapolator(interp, method="Linear")
791
result = extrap([0, 1, 5, 6])  # Extrapolate beyond [2, 4] range
792
print(result)  # [-11, -4, 25, 32] (linear extrapolation)
793

794
# Constant extrapolation with custom boundaries
795
extrap_const = Extrapolator(interp, method="Constant", left=0, right=100)
796
result = extrap_const([0, 1, 5, 6])
797
print(result)  # [0, 0, 100, 100] (constant boundaries)
798
```
799

800
### High-Quality Table Interpolation
801

802
```python
803
from colour.algebra import table_interpolation
804
import colour
805

806
# Load a 3D LUT for color transformation
807
lut_path = "path/to/color_lut.cube"
808
lut = colour.read_LUT(lut_path)
809

810
# Interpolate RGB values using the LUT
811
rgb_values = np.array([
812
    [0.5, 0.3, 0.8],
813
    [0.2, 0.7, 0.1],
814
    [0.9, 0.1, 0.4]
815
])
816

817
# Trilinear interpolation (default)
818
result_trilinear = table_interpolation(rgb_values, lut.table, method="Trilinear")
819

820
# Tetrahedral interpolation (more accurate)
821
result_tetrahedral = table_interpolation(rgb_values, lut.table, method="Tetrahedral")
822

823
print("Trilinear:", result_trilinear)
824
print("Tetrahedral:", result_tetrahedral)
825
```
826

827
## Key Features
828

829
### Performance Optimizations
830

831
- **Vectorized Operations**: All functions support NumPy array broadcasting for efficient batch processing
832
- **Memory Efficiency**: In-place operations where possible to minimize memory allocation
833
- **Numerical Stability**: Safe mathematical operations prevent overflow, underflow, and invalid values
834

835
### Error Handling
836

837
- **Configurable Modes**: Safe division and power functions with multiple error handling strategies
838
- **Graceful Degradation**: Functions continue operation with sensible defaults when encountering edge cases
839
- **Comprehensive Validation**: Input validation with informative error messages
840

841
### Flexibility
842

843
- **Multiple Algorithms**: Choice of interpolation methods for different quality/performance trade-offs
844
- **Customizable Parameters**: Extensive parameterization for fine-tuning mathematical operations
845
- **Modular Design**: Mix and match components to create custom mathematical pipelines
846

847
The mathematical utilities provide the computational foundation for all colour science operations in the package, offering both high-level convenience functions and low-level building blocks for advanced applications.