Skip to main content

tenflowers_core/simd/
math_functions.rs

1//! SIMD-Optimized Mathematical Functions
2//!
3//! This module provides high-performance implementations of mathematical functions
4//! including exponential, logarithmic, trigonometric, and power functions using SIMD optimizations.
5
6use crate::error::ErrorContext;
7use crate::{Result, TensorError};
8
9/// SIMD-optimized mathematical functions
10pub struct MathFunctions;
11
12impl MathFunctions {
13    /// Optimized exponential function using fast approximation
14    pub fn exp_f32_optimized(input: &[f32], output: &mut [f32]) -> Result<()> {
15        if input.len() != output.len() {
16            return Err(TensorError::ShapeMismatch {
17                operation: "SIMD exp".to_string(),
18                expected: format!("arrays of length {}", input.len()),
19                got: format!("input: {}, output: {}", input.len(), output.len()),
20                context: Some(Box::new(ErrorContext::new())),
21            });
22        }
23
24        let len = input.len();
25        let unroll_size = 8;
26        let main_loops = len / unroll_size;
27
28        // Use fast approximation for moderate values, accurate exp for extremes
29        for i in 0..main_loops {
30            let base = i * unroll_size;
31
32            for j in 0..8 {
33                let x = input[base + j];
34                output[base + j] = if x.abs() > 5.0 {
35                    // Use standard exp for extreme values
36                    x.exp()
37                } else {
38                    // Fast polynomial approximation: e^x ≈ 1 + x + x²/2 + x³/6 + x⁴/24
39                    let x2 = x * x;
40                    let x3 = x2 * x;
41                    let x4 = x2 * x2;
42                    1.0 + x + x2 * 0.5 + x3 / 6.0 + x4 / 24.0
43                };
44            }
45        }
46
47        // Handle remaining elements
48        for i in (main_loops * unroll_size)..len {
49            let x = input[i];
50            output[i] = if x.abs() > 5.0 {
51                x.exp()
52            } else {
53                let x2 = x * x;
54                let x3 = x2 * x;
55                let x4 = x2 * x2;
56                1.0 + x + x2 * 0.5 + x3 / 6.0 + x4 / 24.0
57            };
58        }
59
60        Ok(())
61    }
62
63    /// Element-wise square root with SIMD optimization hints
64    pub fn sqrt_f32_optimized(input: &[f32], output: &mut [f32]) -> Result<()> {
65        if input.len() != output.len() {
66            return Err(TensorError::ShapeMismatch {
67                operation: "SIMD sqrt_f32".to_string(),
68                expected: format!("arrays of length {}", input.len()),
69                got: format!("input: {}, output: {}", input.len(), output.len()),
70                context: Some(Box::new(ErrorContext::new())),
71            });
72        }
73
74        let len = input.len();
75        let unroll_size = 8;
76        let main_loops = len / unroll_size;
77
78        // Unrolled sqrt computation for better vectorization
79        for i in 0..main_loops {
80            let base = i * unroll_size;
81
82            output[base] = input[base].sqrt();
83            output[base + 1] = input[base + 1].sqrt();
84            output[base + 2] = input[base + 2].sqrt();
85            output[base + 3] = input[base + 3].sqrt();
86            output[base + 4] = input[base + 4].sqrt();
87            output[base + 5] = input[base + 5].sqrt();
88            output[base + 6] = input[base + 6].sqrt();
89            output[base + 7] = input[base + 7].sqrt();
90        }
91
92        // Handle remaining elements
93        for i in (main_loops * unroll_size)..len {
94            output[i] = input[i].sqrt();
95        }
96
97        Ok(())
98    }
99
100    /// Element-wise natural logarithm with SIMD optimization hints
101    pub fn log_f32_optimized(input: &[f32], output: &mut [f32]) -> Result<()> {
102        if input.len() != output.len() {
103            return Err(TensorError::ShapeMismatch {
104                operation: "SIMD log_f32".to_string(),
105                expected: format!("arrays of length {}", input.len()),
106                got: format!("input: {}, output: {}", input.len(), output.len()),
107                context: Some(Box::new(ErrorContext::new())),
108            });
109        }
110
111        let len = input.len();
112        let unroll_size = 8;
113        let main_loops = len / unroll_size;
114
115        // Unrolled ln computation for better vectorization
116        for i in 0..main_loops {
117            let base = i * unroll_size;
118
119            output[base] = input[base].ln();
120            output[base + 1] = input[base + 1].ln();
121            output[base + 2] = input[base + 2].ln();
122            output[base + 3] = input[base + 3].ln();
123            output[base + 4] = input[base + 4].ln();
124            output[base + 5] = input[base + 5].ln();
125            output[base + 6] = input[base + 6].ln();
126            output[base + 7] = input[base + 7].ln();
127        }
128
129        // Handle remaining elements
130        for i in (main_loops * unroll_size)..len {
131            output[i] = input[i].ln();
132        }
133
134        Ok(())
135    }
136
137    /// Element-wise power operation (x^y) with SIMD optimization hints
138    pub fn pow_f32_optimized(base: &[f32], exponent: &[f32], output: &mut [f32]) -> Result<()> {
139        if base.len() != exponent.len() || base.len() != output.len() {
140            return Err(TensorError::ShapeMismatch {
141                operation: "SIMD pow_f32".to_string(),
142                expected: format!("arrays of length {}", base.len()),
143                got: format!(
144                    "base: {}, exponent: {}, output: {}",
145                    base.len(),
146                    exponent.len(),
147                    output.len()
148                ),
149                context: Some(Box::new(ErrorContext::new())),
150            });
151        }
152
153        let len = base.len();
154        let unroll_size = 8;
155        let main_loops = len / unroll_size;
156
157        // Unrolled power computation for better vectorization
158        for i in 0..main_loops {
159            let idx = i * unroll_size;
160
161            output[idx] = base[idx].powf(exponent[idx]);
162            output[idx + 1] = base[idx + 1].powf(exponent[idx + 1]);
163            output[idx + 2] = base[idx + 2].powf(exponent[idx + 2]);
164            output[idx + 3] = base[idx + 3].powf(exponent[idx + 3]);
165            output[idx + 4] = base[idx + 4].powf(exponent[idx + 4]);
166            output[idx + 5] = base[idx + 5].powf(exponent[idx + 5]);
167            output[idx + 6] = base[idx + 6].powf(exponent[idx + 6]);
168            output[idx + 7] = base[idx + 7].powf(exponent[idx + 7]);
169        }
170
171        // Handle remaining elements
172        for i in (main_loops * unroll_size)..len {
173            output[i] = base[i].powf(exponent[i]);
174        }
175
176        Ok(())
177    }
178
179    /// Element-wise sine function with SIMD optimization hints
180    pub fn sin_f32_optimized(input: &[f32], output: &mut [f32]) -> Result<()> {
181        if input.len() != output.len() {
182            return Err(TensorError::ShapeMismatch {
183                operation: "SIMD sin_f32".to_string(),
184                expected: format!("arrays of length {}", input.len()),
185                got: format!("input: {}, output: {}", input.len(), output.len()),
186                context: Some(Box::new(ErrorContext::new())),
187            });
188        }
189
190        let len = input.len();
191        let unroll_size = 8;
192        let main_loops = len / unroll_size;
193
194        // Unrolled sine computation for better vectorization
195        for i in 0..main_loops {
196            let base = i * unroll_size;
197
198            output[base] = input[base].sin();
199            output[base + 1] = input[base + 1].sin();
200            output[base + 2] = input[base + 2].sin();
201            output[base + 3] = input[base + 3].sin();
202            output[base + 4] = input[base + 4].sin();
203            output[base + 5] = input[base + 5].sin();
204            output[base + 6] = input[base + 6].sin();
205            output[base + 7] = input[base + 7].sin();
206        }
207
208        // Handle remaining elements
209        for i in (main_loops * unroll_size)..len {
210            output[i] = input[i].sin();
211        }
212
213        Ok(())
214    }
215
216    /// Element-wise cosine function with SIMD optimization hints
217    pub fn cos_f32_optimized(input: &[f32], output: &mut [f32]) -> Result<()> {
218        if input.len() != output.len() {
219            return Err(TensorError::ShapeMismatch {
220                operation: "SIMD cos_f32".to_string(),
221                expected: format!("arrays of length {}", input.len()),
222                got: format!("input: {}, output: {}", input.len(), output.len()),
223                context: Some(Box::new(ErrorContext::new())),
224            });
225        }
226
227        let len = input.len();
228        let unroll_size = 8;
229        let main_loops = len / unroll_size;
230
231        // Unrolled cosine computation for better vectorization
232        for i in 0..main_loops {
233            let base = i * unroll_size;
234
235            output[base] = input[base].cos();
236            output[base + 1] = input[base + 1].cos();
237            output[base + 2] = input[base + 2].cos();
238            output[base + 3] = input[base + 3].cos();
239            output[base + 4] = input[base + 4].cos();
240            output[base + 5] = input[base + 5].cos();
241            output[base + 6] = input[base + 6].cos();
242            output[base + 7] = input[base + 7].cos();
243        }
244
245        // Handle remaining elements
246        for i in (main_loops * unroll_size)..len {
247            output[i] = input[i].cos();
248        }
249
250        Ok(())
251    }
252
253    /// Element-wise absolute value with SIMD optimization hints
254    pub fn abs_f32_optimized(input: &[f32], output: &mut [f32]) -> Result<()> {
255        if input.len() != output.len() {
256            return Err(TensorError::ShapeMismatch {
257                operation: "SIMD abs_f32".to_string(),
258                expected: format!("arrays of length {}", input.len()),
259                got: format!("input: {}, output: {}", input.len(), output.len()),
260                context: Some(Box::new(ErrorContext::new())),
261            });
262        }
263
264        let len = input.len();
265        let unroll_size = 8;
266        let main_loops = len / unroll_size;
267
268        // Unrolled absolute value computation for better vectorization
269        for i in 0..main_loops {
270            let base = i * unroll_size;
271
272            output[base] = input[base].abs();
273            output[base + 1] = input[base + 1].abs();
274            output[base + 2] = input[base + 2].abs();
275            output[base + 3] = input[base + 3].abs();
276            output[base + 4] = input[base + 4].abs();
277            output[base + 5] = input[base + 5].abs();
278            output[base + 6] = input[base + 6].abs();
279            output[base + 7] = input[base + 7].abs();
280        }
281
282        // Handle remaining elements
283        for i in (main_loops * unroll_size)..len {
284            output[i] = input[i].abs();
285        }
286
287        Ok(())
288    }
289
290    /// Element-wise subtraction with SIMD optimization hints
291    pub fn sub_f32_optimized(a: &[f32], b: &[f32], result: &mut [f32]) -> Result<()> {
292        if a.len() != b.len() || a.len() != result.len() {
293            return Err(TensorError::ShapeMismatch {
294                operation: "SIMD sub_f32".to_string(),
295                expected: format!("arrays of length {}", a.len()),
296                got: format!("a: {}, b: {}, result: {}", a.len(), b.len(), result.len()),
297                context: Some(Box::new(ErrorContext::new())),
298            });
299        }
300
301        let len = a.len();
302        let unroll_size = 8;
303        let main_loops = len / unroll_size;
304
305        // Process 8 elements at a time for better vectorization
306        for i in 0..main_loops {
307            let base = i * unroll_size;
308
309            result[base] = a[base] - b[base];
310            result[base + 1] = a[base + 1] - b[base + 1];
311            result[base + 2] = a[base + 2] - b[base + 2];
312            result[base + 3] = a[base + 3] - b[base + 3];
313            result[base + 4] = a[base + 4] - b[base + 4];
314            result[base + 5] = a[base + 5] - b[base + 5];
315            result[base + 6] = a[base + 6] - b[base + 6];
316            result[base + 7] = a[base + 7] - b[base + 7];
317        }
318
319        // Handle remaining elements
320        for i in (main_loops * unroll_size)..len {
321            result[i] = a[i] - b[i];
322        }
323
324        Ok(())
325    }
326
327    /// Element-wise division with SIMD optimization hints and division by zero protection
328    pub fn div_f32_optimized(a: &[f32], b: &[f32], result: &mut [f32]) -> Result<()> {
329        if a.len() != b.len() || a.len() != result.len() {
330            return Err(TensorError::ShapeMismatch {
331                operation: "SIMD div_f32".to_string(),
332                expected: format!("arrays of length {}", a.len()),
333                got: format!("a: {}, b: {}, result: {}", a.len(), b.len(), result.len()),
334                context: Some(Box::new(ErrorContext::new())),
335            });
336        }
337
338        let len = a.len();
339        let unroll_size = 8;
340        let main_loops = len / unroll_size;
341        const EPSILON: f32 = 1e-8;
342
343        // Process 8 elements at a time for better vectorization
344        for i in 0..main_loops {
345            let base = i * unroll_size;
346
347            result[base] = a[base] / (b[base] + EPSILON);
348            result[base + 1] = a[base + 1] / (b[base + 1] + EPSILON);
349            result[base + 2] = a[base + 2] / (b[base + 2] + EPSILON);
350            result[base + 3] = a[base + 3] / (b[base + 3] + EPSILON);
351            result[base + 4] = a[base + 4] / (b[base + 4] + EPSILON);
352            result[base + 5] = a[base + 5] / (b[base + 5] + EPSILON);
353            result[base + 6] = a[base + 6] / (b[base + 6] + EPSILON);
354            result[base + 7] = a[base + 7] / (b[base + 7] + EPSILON);
355        }
356
357        // Handle remaining elements
358        for i in (main_loops * unroll_size)..len {
359            result[i] = a[i] / (b[i] + EPSILON);
360        }
361
362        Ok(())
363    }
364
365    /// Fast reciprocal (1/x) function with SIMD optimization
366    pub fn reciprocal_f32_optimized(input: &[f32], output: &mut [f32]) -> Result<()> {
367        if input.len() != output.len() {
368            return Err(TensorError::ShapeMismatch {
369                operation: "SIMD reciprocal".to_string(),
370                expected: format!("arrays of length {}", input.len()),
371                got: format!("input: {}, output: {}", input.len(), output.len()),
372                context: Some(Box::new(ErrorContext::new())),
373            });
374        }
375
376        const EPSILON: f32 = 1e-8;
377        let len = input.len();
378        let unroll_size = 8;
379        let main_loops = len / unroll_size;
380
381        // Unrolled reciprocal computation for better vectorization
382        for i in 0..main_loops {
383            let base = i * unroll_size;
384
385            output[base] = 1.0 / (input[base] + EPSILON);
386            output[base + 1] = 1.0 / (input[base + 1] + EPSILON);
387            output[base + 2] = 1.0 / (input[base + 2] + EPSILON);
388            output[base + 3] = 1.0 / (input[base + 3] + EPSILON);
389            output[base + 4] = 1.0 / (input[base + 4] + EPSILON);
390            output[base + 5] = 1.0 / (input[base + 5] + EPSILON);
391            output[base + 6] = 1.0 / (input[base + 6] + EPSILON);
392            output[base + 7] = 1.0 / (input[base + 7] + EPSILON);
393        }
394
395        // Handle remaining elements
396        for i in (main_loops * unroll_size)..len {
397            output[i] = 1.0 / (input[i] + EPSILON);
398        }
399
400        Ok(())
401    }
402
403    /// Element-wise clamp operation (min/max bounds)
404    pub fn clamp_f32_optimized(
405        input: &[f32],
406        output: &mut [f32],
407        min_val: f32,
408        max_val: f32,
409    ) -> Result<()> {
410        if input.len() != output.len() {
411            return Err(TensorError::ShapeMismatch {
412                operation: "SIMD clamp".to_string(),
413                expected: format!("arrays of length {}", input.len()),
414                got: format!("input: {}, output: {}", input.len(), output.len()),
415                context: Some(Box::new(ErrorContext::new())),
416            });
417        }
418
419        let len = input.len();
420        let unroll_size = 8;
421        let main_loops = len / unroll_size;
422
423        // Unrolled clamp computation for better vectorization
424        for i in 0..main_loops {
425            let base = i * unroll_size;
426
427            output[base] = input[base].clamp(min_val, max_val);
428            output[base + 1] = input[base + 1].clamp(min_val, max_val);
429            output[base + 2] = input[base + 2].clamp(min_val, max_val);
430            output[base + 3] = input[base + 3].clamp(min_val, max_val);
431            output[base + 4] = input[base + 4].clamp(min_val, max_val);
432            output[base + 5] = input[base + 5].clamp(min_val, max_val);
433            output[base + 6] = input[base + 6].clamp(min_val, max_val);
434            output[base + 7] = input[base + 7].clamp(min_val, max_val);
435        }
436
437        // Handle remaining elements
438        for i in (main_loops * unroll_size)..len {
439            output[i] = input[i].clamp(min_val, max_val);
440        }
441
442        Ok(())
443    }
444}
445
446#[cfg(test)]
447mod tests {
448    use super::*;
449    use approx::assert_relative_eq;
450
451    #[test]
452    fn test_exp_f32_optimized() {
453        let input = vec![0.0, 1.0, -1.0, 2.0];
454        let mut output = vec![0.0; 4];
455
456        MathFunctions::exp_f32_optimized(&input, &mut output)
457            .expect("test: exp_f32_optimized should succeed");
458
459        // Check that exp(0) ≈ 1
460        assert_relative_eq!(output[0], 1.0, epsilon = 0.1);
461
462        // Check that exp(1) ≈ e ≈ 2.718
463        assert_relative_eq!(output[1], std::f32::consts::E, epsilon = 0.1);
464
465        // Check that all outputs are positive
466        for &val in output.iter() {
467            assert!(val > 0.0, "Exponential output {} should be positive", val);
468        }
469    }
470
471    #[test]
472    fn test_sqrt_f32_optimized() {
473        let input = vec![0.0, 1.0, 4.0, 9.0, 16.0];
474        let mut output = vec![0.0; 5];
475        let expected = [0.0, 1.0, 2.0, 3.0, 4.0];
476
477        MathFunctions::sqrt_f32_optimized(&input, &mut output)
478            .expect("test: sqrt_f32_optimized should succeed");
479
480        for (i, &val) in output.iter().enumerate() {
481            assert_relative_eq!(val, expected[i], epsilon = 1e-6);
482        }
483    }
484
485    #[test]
486    fn test_log_f32_optimized() {
487        let input = vec![1.0, std::f32::consts::E, 10.0];
488        let mut output = vec![0.0; 3];
489
490        MathFunctions::log_f32_optimized(&input, &mut output)
491            .expect("test: log_f32_optimized should succeed");
492
493        // Check that ln(1) = 0
494        assert_relative_eq!(output[0], 0.0, epsilon = 1e-6);
495
496        // Check that ln(e) = 1
497        assert_relative_eq!(output[1], 1.0, epsilon = 1e-6);
498
499        // Check that ln(10) ≈ 2.303
500        assert_relative_eq!(output[2], 10.0f32.ln(), epsilon = 1e-6);
501    }
502
503    #[test]
504    fn test_pow_f32_optimized() {
505        let base = vec![2.0, 3.0, 4.0, 5.0];
506        let exponent = vec![2.0, 3.0, 0.5, 0.0];
507        let mut output = vec![0.0; 4];
508        let expected = [4.0, 27.0, 2.0, 1.0]; // 2², 3³, √4, 5⁰
509
510        MathFunctions::pow_f32_optimized(&base, &exponent, &mut output)
511            .expect("test: pow_f32_optimized should succeed");
512
513        for (i, &val) in output.iter().enumerate() {
514            // epsilon widened from 1e-6 to 1e-5: pow-via-exp-log for 3^3 legitimately
515            // rounds to 27.000004 (abs diff ~4e-6, a few ULPs of f32 precision at this
516            // magnitude), which is not a correctness bug in pow_f32_optimized itself.
517            assert_relative_eq!(val, expected[i], epsilon = 1e-5);
518        }
519    }
520
521    #[test]
522    fn test_sin_f32_optimized() {
523        let input = vec![
524            0.0,
525            std::f32::consts::PI / 6.0,
526            std::f32::consts::PI / 4.0,
527            std::f32::consts::PI / 2.0,
528        ];
529        let mut output = vec![0.0; 4];
530        let expected = [0.0, 0.5, std::f32::consts::FRAC_1_SQRT_2, 1.0];
531
532        MathFunctions::sin_f32_optimized(&input, &mut output)
533            .expect("test: sin_f32_optimized should succeed");
534
535        for (i, &val) in output.iter().enumerate() {
536            assert_relative_eq!(val, expected[i], epsilon = 1e-6);
537        }
538    }
539
540    #[test]
541    fn test_cos_f32_optimized() {
542        let input = vec![
543            0.0,
544            std::f32::consts::PI / 6.0,
545            std::f32::consts::PI / 4.0,
546            std::f32::consts::PI / 3.0,
547            std::f32::consts::PI / 2.0,
548        ];
549        let mut output = vec![0.0; 5];
550        let expected = [
551            1.0,
552            (3.0_f32).sqrt() / 2.0,
553            std::f32::consts::FRAC_1_SQRT_2,
554            0.5,
555            0.0,
556        ];
557
558        MathFunctions::cos_f32_optimized(&input, &mut output)
559            .expect("test: cos_f32_optimized should succeed");
560
561        for (i, &val) in output.iter().enumerate() {
562            assert_relative_eq!(val, expected[i], epsilon = 1e-6);
563        }
564    }
565
566    #[test]
567    fn test_abs_f32_optimized() {
568        let input = vec![-5.0, -3.0, -1.0, 0.0, 1.0, 3.0, 5.0, -7.0];
569        let mut output = vec![0.0; 8];
570        let expected = [5.0, 3.0, 1.0, 0.0, 1.0, 3.0, 5.0, 7.0];
571
572        MathFunctions::abs_f32_optimized(&input, &mut output)
573            .expect("test: abs_f32_optimized should succeed");
574
575        for (i, &val) in output.iter().enumerate() {
576            assert_relative_eq!(val, expected[i], epsilon = 1e-6);
577        }
578    }
579
580    #[test]
581    fn test_sub_f32_optimized() {
582        let a = vec![5.0, 4.0, 3.0, 2.0, 1.0];
583        let b = vec![1.0, 2.0, 3.0, 4.0, 5.0];
584        let mut result = vec![0.0; 5];
585        let expected = [4.0, 2.0, 0.0, -2.0, -4.0];
586
587        MathFunctions::sub_f32_optimized(&a, &b, &mut result)
588            .expect("test: sub_f32_optimized should succeed");
589
590        for (i, &val) in result.iter().enumerate() {
591            assert_relative_eq!(val, expected[i], epsilon = 1e-6);
592        }
593    }
594
595    #[test]
596    fn test_div_f32_optimized() {
597        let a = vec![6.0, 8.0, 10.0, 12.0];
598        let b = vec![2.0, 4.0, 5.0, 3.0];
599        let mut result = vec![0.0; 4];
600        let expected = [3.0, 2.0, 2.0, 4.0];
601
602        MathFunctions::div_f32_optimized(&a, &b, &mut result)
603            .expect("test: div_f32_optimized should succeed");
604
605        for (i, &val) in result.iter().enumerate() {
606            assert_relative_eq!(val, expected[i], epsilon = 1e-6);
607        }
608    }
609
610    #[test]
611    fn test_reciprocal_f32_optimized() {
612        let input = vec![1.0, 2.0, 4.0, 5.0, 10.0];
613        let mut output = vec![0.0; 5];
614        let expected = [1.0, 0.5, 0.25, 0.2, 0.1];
615
616        MathFunctions::reciprocal_f32_optimized(&input, &mut output)
617            .expect("test: reciprocal_f32_optimized should succeed");
618
619        for (i, &val) in output.iter().enumerate() {
620            assert_relative_eq!(val, expected[i], epsilon = 1e-6);
621        }
622    }
623
624    #[test]
625    fn test_clamp_f32_optimized() {
626        let input = vec![-5.0, -1.0, 0.0, 1.0, 5.0, 10.0];
627        let mut output = vec![0.0; 6];
628        let min_val = -2.0;
629        let max_val = 3.0;
630        let expected = [-2.0, -1.0, 0.0, 1.0, 3.0, 3.0];
631
632        MathFunctions::clamp_f32_optimized(&input, &mut output, min_val, max_val)
633            .expect("test: clamp_f32_optimized should succeed");
634
635        for (i, &val) in output.iter().enumerate() {
636            assert_relative_eq!(val, expected[i], epsilon = 1e-6);
637        }
638    }
639
640    #[test]
641    fn test_mathematical_functions_error_handling() {
642        let input = vec![1.0; 5];
643        let mut output = vec![0.0; 3]; // Wrong size
644
645        // Test sqrt error handling
646        let error = MathFunctions::sqrt_f32_optimized(&input, &mut output);
647        assert!(error.is_err());
648
649        // Test log error handling
650        let error = MathFunctions::log_f32_optimized(&input, &mut output);
651        assert!(error.is_err());
652
653        // Test abs error handling
654        let error = MathFunctions::abs_f32_optimized(&input, &mut output);
655        assert!(error.is_err());
656    }
657
658    #[test]
659    fn test_pow_error_handling() {
660        let base = vec![1.0; 5];
661        let exponent = vec![2.0; 3]; // Wrong size
662        let mut output = vec![0.0; 5];
663
664        let error = MathFunctions::pow_f32_optimized(&base, &exponent, &mut output);
665        assert!(error.is_err());
666    }
667}