1use crate::error::ErrorContext;
7use crate::{Result, TensorError};
8
9pub struct MathFunctions;
11
12impl MathFunctions {
13 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 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 x.exp()
37 } else {
38 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 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 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 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 for i in (main_loops * unroll_size)..len {
94 output[i] = input[i].sqrt();
95 }
96
97 Ok(())
98 }
99
100 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 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 for i in (main_loops * unroll_size)..len {
131 output[i] = input[i].ln();
132 }
133
134 Ok(())
135 }
136
137 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 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 for i in (main_loops * unroll_size)..len {
173 output[i] = base[i].powf(exponent[i]);
174 }
175
176 Ok(())
177 }
178
179 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 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 for i in (main_loops * unroll_size)..len {
210 output[i] = input[i].sin();
211 }
212
213 Ok(())
214 }
215
216 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 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 for i in (main_loops * unroll_size)..len {
247 output[i] = input[i].cos();
248 }
249
250 Ok(())
251 }
252
253 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 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 for i in (main_loops * unroll_size)..len {
284 output[i] = input[i].abs();
285 }
286
287 Ok(())
288 }
289
290 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 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 for i in (main_loops * unroll_size)..len {
321 result[i] = a[i] - b[i];
322 }
323
324 Ok(())
325 }
326
327 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 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 for i in (main_loops * unroll_size)..len {
359 result[i] = a[i] / (b[i] + EPSILON);
360 }
361
362 Ok(())
363 }
364
365 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 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 for i in (main_loops * unroll_size)..len {
397 output[i] = 1.0 / (input[i] + EPSILON);
398 }
399
400 Ok(())
401 }
402
403 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 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 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 assert_relative_eq!(output[0], 1.0, epsilon = 0.1);
461
462 assert_relative_eq!(output[1], std::f32::consts::E, epsilon = 0.1);
464
465 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 assert_relative_eq!(output[0], 0.0, epsilon = 1e-6);
495
496 assert_relative_eq!(output[1], 1.0, epsilon = 1e-6);
498
499 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]; 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 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]; let error = MathFunctions::sqrt_f32_optimized(&input, &mut output);
647 assert!(error.is_err());
648
649 let error = MathFunctions::log_f32_optimized(&input, &mut output);
651 assert!(error.is_err());
652
653 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]; 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}