scirs2_core/
numeric.rs

1//! Numeric traits and utilities for SciRS2
2//!
3//! This module provides unified access to numeric traits and complex numbers.
4//!
5//! # SciRS2 POLICY Compliance
6//!
7//! This is the **only** module in the SciRS2 ecosystem that should directly
8//! import external numeric libraries. All other crates must use this module.
9//!
10//! # Re-exports
11//!
12//! - `num_traits`: All numeric traits (Float, Zero, One, etc.)
13//! - `num_complex`: Complex number types
14//! - `num_integer`: Integer-specific functions (gcd, lcm, etc.)
15
16use crate::error::{CoreError, CoreResult, ErrorContext};
17use std::fmt::Debug;
18use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
19
20// ============================================================================
21// Complete Re-exports for SciRS2 POLICY Compliance
22// ============================================================================
23
24/// Complete re-export of num-traits for SciRS2 ecosystem
25pub use num_traits::{
26    self,
27    bounds::Bounded,
28    cast::{AsPrimitive, FromPrimitive, NumCast, ToPrimitive},
29    float::{Float, FloatConst},
30    identities::{One, Zero},
31    int::PrimInt,
32    ops::{
33        checked::*,
34        inv::Inv,
35        mul_add::{MulAdd, MulAddAssign},
36        saturating::*,
37        wrapping::*,
38    },
39    pow::{checked_pow, Pow},
40    real::Real,
41    sign::{Signed, Unsigned},
42    Num, NumAssign, NumAssignOps, NumAssignRef, NumOps, NumRef,
43};
44
45/// Complete re-export of num-complex for SciRS2 ecosystem
46pub use num_complex::{self, Complex, Complex32, Complex64, ComplexFloat};
47
48/// Complete re-export of num-integer for SciRS2 ecosystem
49pub use num_integer::{
50    self, div_ceil, div_floor, div_mod_floor, gcd, gcd_lcm, lcm, Average, ExtendedGcd, Integer,
51    Roots,
52};
53
54// ============================================================================
55// SciRS2 Custom Traits and Utilities
56// ============================================================================
57
58/// A trait for numeric types that can be used in scientific calculations
59///
60/// This trait combines common numeric traits required for scientific computing
61/// operations, providing a unified trait bound for generic code.
62pub trait ScientificNumber:
63    Num
64    + Clone
65    + Copy
66    + PartialOrd
67    + Debug
68    + Zero
69    + One
70    + Add<Output = Self>
71    + Sub<Output = Self>
72    + Mul<Output = Self>
73    + Div<Output = Self>
74    + AddAssign
75    + SubAssign
76    + MulAssign
77    + DivAssign
78    + NumCast
79{
80    /// Absolute value
81    #[must_use]
82    fn abs(self) -> Self;
83
84    /// Square root
85    #[must_use]
86    fn sqrt(self) -> Self;
87
88    /// Square
89    #[must_use]
90    fn square(self) -> Self {
91        self * self
92    }
93
94    /// Maximum of two values
95    #[must_use]
96    fn max(self, other: Self) -> Self;
97
98    /// Minimum of two values
99    #[must_use]
100    fn min(self, other: Self) -> Self;
101
102    /// Check if the value is finite
103    #[must_use]
104    fn is_finite(self) -> bool;
105
106    /// Convert to f64
107    #[must_use]
108    fn to_f64(self) -> Option<f64>;
109
110    /// Convert from f64
111    #[must_use]
112    fn from_f64(value: f64) -> Option<Self>;
113
114    /// Convert from little-endian bytes
115    #[must_use]
116    fn from_le_bytes(bytes: &[u8]) -> Self;
117
118    /// Convert from big-endian bytes
119    #[must_use]
120    fn from_be_bytes(bytes: &[u8]) -> Self;
121
122    /// Convert to little-endian bytes
123    #[must_use]
124    fn to_le_bytes(self) -> Vec<u8>;
125
126    /// Convert to big-endian bytes
127    #[must_use]
128    fn to_be_bytes(self) -> Vec<u8>;
129}
130
131/// A trait for real-valued floating point types
132pub trait RealNumber: ScientificNumber + Float {
133    /// Returns the machine epsilon (the difference between 1.0 and the least value greater than 1.0)
134    #[must_use]
135    fn epsilon() -> Self;
136
137    /// Exponential function (e^x)
138    #[must_use]
139    fn exp(self) -> Self;
140
141    /// Natural logarithm (ln(x))
142    #[must_use]
143    fn ln(self) -> Self;
144
145    /// Base-10 logarithm
146    #[must_use]
147    fn log10(self) -> Self;
148
149    /// Base-2 logarithm
150    #[must_use]
151    fn log2(self) -> Self;
152
153    /// Sine function
154    #[must_use]
155    fn sin(self) -> Self;
156
157    /// Cosine function
158    #[must_use]
159    fn cos(self) -> Self;
160
161    /// Tangent function
162    #[must_use]
163    fn tan(self) -> Self;
164
165    /// Hyperbolic sine
166    #[must_use]
167    fn sinh(self) -> Self;
168
169    /// Hyperbolic cosine
170    #[must_use]
171    fn cosh(self) -> Self;
172
173    /// Hyperbolic tangent
174    #[must_use]
175    fn tanh(self) -> Self;
176
177    /// Inverse sine
178    #[must_use]
179    fn asin(self) -> Self;
180
181    /// Inverse cosine
182    #[must_use]
183    fn acos(self) -> Self;
184
185    /// Inverse tangent
186    #[must_use]
187    fn atan(self) -> Self;
188
189    /// Inverse tangent of y/x with correct quadrant
190    #[must_use]
191    fn atan2(self, other: Self) -> Self;
192
193    /// Power function
194    #[must_use]
195    fn powf(self, n: Self) -> Self;
196
197    /// Integer power function
198    #[must_use]
199    fn powi(self, n: i32) -> Self;
200
201    /// Factorial function (approximation for non-integers)
202    #[must_use]
203    fn factorial(self) -> Self;
204}
205
206/// A trait for complex number types
207pub trait ComplexNumber: ScientificNumber {
208    /// The real part of the complex number
209    type RealPart: RealNumber;
210
211    /// Returns the real part of the complex number
212    #[must_use]
213    fn re(&self) -> Self::RealPart;
214
215    /// Returns the imaginary part of the complex number
216    #[must_use]
217    fn im(&self) -> Self::RealPart;
218
219    /// Create a new complex number from real and imaginary parts
220    #[must_use]
221    fn from_parts(re: Self::RealPart, im: Self::RealPart) -> Self;
222
223    /// Returns the complex conjugate
224    #[must_use]
225    fn conj(self) -> Self;
226
227    /// Returns the magnitude (absolute value)
228    #[must_use]
229    fn abs(self) -> Self::RealPart;
230
231    /// Returns the argument (phase)
232    #[must_use]
233    fn arg(self) -> Self::RealPart;
234
235    /// Returns the complex number in exponential form (r, theta)
236    #[must_use]
237    fn to_polar(self) -> (Self::RealPart, Self::RealPart);
238
239    /// Creates a complex number from polar coordinates
240    #[must_use]
241    fn from_polar(r: Self::RealPart, theta: Self::RealPart) -> Self;
242
243    /// Exponential function
244    #[must_use]
245    fn exp(self) -> Self;
246
247    /// Natural logarithm
248    #[must_use]
249    fn ln(self) -> Self;
250
251    /// Power function with complex exponent
252    #[must_use]
253    fn powc(self, exp: Self) -> Self;
254
255    /// Power function with real exponent
256    #[must_use]
257    fn powf(self, exp: Self::RealPart) -> Self;
258
259    /// Square root
260    #[must_use]
261    fn sqrt(self) -> Self;
262}
263
264/// A trait for integers that can be used in scientific calculations
265pub trait ScientificInteger: ScientificNumber + Eq {
266    /// Greatest common divisor
267    #[must_use]
268    fn gcd(self, other: Self) -> Self;
269
270    /// Least common multiple
271    #[must_use]
272    fn lcm(self, other: Self) -> Self;
273
274    /// Check if the number is prime
275    #[must_use]
276    fn is_prime(self) -> bool;
277
278    /// Check if the number is even
279    #[must_use]
280    fn is_even(self) -> bool;
281
282    /// Check if the number is odd
283    #[must_use]
284    fn is_odd(self) -> bool;
285
286    /// Modular exponentiation (self^exp mod modulus)
287    #[must_use]
288    fn mod_pow(self, exp: Self, modulus: Self) -> Self;
289
290    /// Factorial
291    fn factorial(self) -> CoreResult<Self>;
292
293    /// Binomial coefficient (n choose k)
294    #[must_use]
295    fn binomial(self, k: Self) -> Self;
296}
297
298// Implement ScientificNumber for f32
299impl ScientificNumber for f32 {
300    fn abs(self) -> Self {
301        self.abs()
302    }
303
304    fn sqrt(self) -> Self {
305        self.sqrt()
306    }
307
308    fn max(self, other: Self) -> Self {
309        self.max(other)
310    }
311
312    fn min(self, other: Self) -> Self {
313        self.min(other)
314    }
315
316    fn is_finite(self) -> bool {
317        self.is_finite()
318    }
319
320    fn to_f64(self) -> Option<f64> {
321        Some(self as f64)
322    }
323
324    fn from_f64(value: f64) -> Option<Self> {
325        if value.is_finite() && value <= Self::MAX as f64 && value >= Self::MIN as f64 {
326            Some(value as f32)
327        } else {
328            None
329        }
330    }
331
332    fn from_le_bytes(bytes: &[u8]) -> Self {
333        let mut array = [0u8; 4];
334        array.copy_from_slice(&bytes[..4]);
335        f32::from_le_bytes(array)
336    }
337
338    fn from_be_bytes(bytes: &[u8]) -> Self {
339        let mut array = [0u8; 4];
340        array.copy_from_slice(&bytes[..4]);
341        f32::from_be_bytes(array)
342    }
343
344    fn to_le_bytes(self) -> Vec<u8> {
345        self.to_le_bytes().to_vec()
346    }
347
348    fn to_be_bytes(self) -> Vec<u8> {
349        self.to_be_bytes().to_vec()
350    }
351}
352
353// Implement RealNumber for f32
354impl RealNumber for f32 {
355    fn epsilon() -> Self {
356        Self::EPSILON
357    }
358
359    fn exp(self) -> Self {
360        self.exp()
361    }
362
363    fn ln(self) -> Self {
364        self.ln()
365    }
366
367    fn log10(self) -> Self {
368        self.log10()
369    }
370
371    fn log2(self) -> Self {
372        self.log2()
373    }
374
375    fn sin(self) -> Self {
376        self.sin()
377    }
378
379    fn cos(self) -> Self {
380        self.cos()
381    }
382
383    fn tan(self) -> Self {
384        self.tan()
385    }
386
387    fn sinh(self) -> Self {
388        self.sinh()
389    }
390
391    fn cosh(self) -> Self {
392        self.cosh()
393    }
394
395    fn tanh(self) -> Self {
396        self.tanh()
397    }
398
399    fn asin(self) -> Self {
400        self.asin()
401    }
402
403    fn acos(self) -> Self {
404        self.acos()
405    }
406
407    fn atan(self) -> Self {
408        self.atan()
409    }
410
411    fn atan2(self, other: Self) -> Self {
412        self.atan2(other)
413    }
414
415    fn powf(self, n: Self) -> Self {
416        self.powf(n)
417    }
418
419    fn powi(self, n: i32) -> Self {
420        self.powi(n)
421    }
422
423    fn factorial(self) -> Self {
424        if self < 0.0 {
425            return Self::NAN;
426        }
427
428        // Use Stirling's approximation for non-integers or large values
429        if self != self.trunc() || self > 100.0 {
430            const SQRT_TWO_PI: f32 = 2.506_628_3;
431            return SQRT_TWO_PI * self.powf(self + 0.5) * (-self).exp();
432        }
433
434        let mut result = 1.0;
435        let n = self as u32;
436        for i in 2..=n {
437            result *= i as f32;
438        }
439
440        result
441    }
442}
443
444// Implement ScientificNumber for f64
445impl ScientificNumber for f64 {
446    fn abs(self) -> Self {
447        self.abs()
448    }
449
450    fn sqrt(self) -> Self {
451        self.sqrt()
452    }
453
454    fn max(self, other: Self) -> Self {
455        self.max(other)
456    }
457
458    fn min(self, other: Self) -> Self {
459        self.min(other)
460    }
461
462    fn is_finite(self) -> bool {
463        self.is_finite()
464    }
465
466    fn to_f64(self) -> Option<f64> {
467        Some(self)
468    }
469
470    fn from_f64(value: f64) -> Option<Self> {
471        Some(value)
472    }
473
474    fn from_le_bytes(bytes: &[u8]) -> Self {
475        let mut array = [0u8; 8];
476        array.copy_from_slice(&bytes[..8]);
477        f64::from_le_bytes(array)
478    }
479
480    fn from_be_bytes(bytes: &[u8]) -> Self {
481        let mut array = [0u8; 8];
482        array.copy_from_slice(&bytes[..8]);
483        f64::from_be_bytes(array)
484    }
485
486    fn to_le_bytes(self) -> Vec<u8> {
487        self.to_le_bytes().to_vec()
488    }
489
490    fn to_be_bytes(self) -> Vec<u8> {
491        self.to_be_bytes().to_vec()
492    }
493}
494
495// Implement RealNumber for f64
496impl RealNumber for f64 {
497    fn epsilon() -> Self {
498        Self::EPSILON
499    }
500
501    fn exp(self) -> Self {
502        self.exp()
503    }
504
505    fn ln(self) -> Self {
506        self.ln()
507    }
508
509    fn log10(self) -> Self {
510        self.log10()
511    }
512
513    fn log2(self) -> Self {
514        self.log2()
515    }
516
517    fn sin(self) -> Self {
518        self.sin()
519    }
520
521    fn cos(self) -> Self {
522        self.cos()
523    }
524
525    fn tan(self) -> Self {
526        self.tan()
527    }
528
529    fn sinh(self) -> Self {
530        self.sinh()
531    }
532
533    fn cosh(self) -> Self {
534        self.cosh()
535    }
536
537    fn tanh(self) -> Self {
538        self.tanh()
539    }
540
541    fn asin(self) -> Self {
542        self.asin()
543    }
544
545    fn acos(self) -> Self {
546        self.acos()
547    }
548
549    fn atan(self) -> Self {
550        self.atan()
551    }
552
553    fn atan2(self, other: Self) -> Self {
554        self.atan2(other)
555    }
556
557    fn powf(self, n: Self) -> Self {
558        self.powf(n)
559    }
560
561    fn powi(self, n: i32) -> Self {
562        self.powi(n)
563    }
564
565    fn factorial(self) -> Self {
566        if self < 0.0 {
567            return Self::NAN;
568        }
569
570        // Use Stirling's approximation for non-integers or large values
571        if self != self.trunc() || self > 170.0 {
572            const SQRT_TWO_PI: f64 = 2.506_628_274_631_000_2;
573            return SQRT_TWO_PI * self.powf(self + 0.5) * (-self).exp();
574        }
575
576        let mut result = 1.0;
577        let n = self as u32;
578        for i in 2..=n {
579            result *= i as f64;
580        }
581
582        result
583    }
584}
585
586// Implement ScientificNumber for i32
587impl ScientificNumber for i32 {
588    fn abs(self) -> Self {
589        self.abs()
590    }
591
592    fn sqrt(self) -> Self {
593        (self as f64).sqrt() as i32
594    }
595
596    fn max(self, other: Self) -> Self {
597        std::cmp::max(self, other)
598    }
599
600    fn min(self, other: Self) -> Self {
601        std::cmp::min(self, other)
602    }
603
604    fn is_finite(self) -> bool {
605        true // integers are always finite
606    }
607
608    fn to_f64(self) -> Option<f64> {
609        Some(self as f64)
610    }
611
612    fn from_f64(value: f64) -> Option<Self> {
613        if value.is_finite() && value <= i32::MAX as f64 && value >= i32::MIN as f64 {
614            Some(value as i32)
615        } else {
616            None
617        }
618    }
619
620    fn from_le_bytes(bytes: &[u8]) -> Self {
621        let mut array = [0u8; 4];
622        array.copy_from_slice(&bytes[..4]);
623        i32::from_le_bytes(array)
624    }
625
626    fn from_be_bytes(bytes: &[u8]) -> Self {
627        let mut array = [0u8; 4];
628        array.copy_from_slice(&bytes[..4]);
629        i32::from_be_bytes(array)
630    }
631
632    fn to_le_bytes(self) -> Vec<u8> {
633        self.to_le_bytes().to_vec()
634    }
635
636    fn to_be_bytes(self) -> Vec<u8> {
637        self.to_be_bytes().to_vec()
638    }
639}
640
641// Implement ScientificInteger for i32
642impl ScientificInteger for i32 {
643    fn gcd(self, other: Self) -> Self {
644        let mut a = self.abs();
645        let mut b = other.abs();
646
647        if a == 0 {
648            return b;
649        }
650        if b == 0 {
651            return a;
652        }
653
654        while b != 0 {
655            let temp = b;
656            b = a % b;
657            a = temp;
658        }
659
660        a
661    }
662
663    fn lcm(self, other: Self) -> Self {
664        if self == 0 || other == 0 {
665            return 0;
666        }
667
668        let gcd = self.gcd(other);
669
670        (self / gcd) * other
671    }
672
673    fn is_prime(self) -> bool {
674        if self <= 1 {
675            return false;
676        }
677        if self <= 3 {
678            return true;
679        }
680        if self % 2 == 0 || self % 3 == 0 {
681            return false;
682        }
683
684        let mut i = 5;
685        while i * i <= self {
686            if self % i == 0 || self % (i + 2) == 0 {
687                return false;
688            }
689            i += 6;
690        }
691
692        true
693    }
694
695    fn is_even(self) -> bool {
696        self % 2 == 0
697    }
698
699    fn is_odd(self) -> bool {
700        self % 2 != 0
701    }
702
703    fn mod_pow(self, exp: Self, modulus: Self) -> Self {
704        if modulus == 1 {
705            return 0;
706        }
707
708        let mut base = self % modulus;
709        let mut result = 1;
710        let mut exponent = exp;
711
712        while exponent > 0 {
713            if exponent % 2 == 1 {
714                result = (result * base) % modulus;
715            }
716            exponent >>= 1;
717            base = (base * base) % modulus;
718        }
719
720        result
721    }
722
723    fn factorial(self) -> CoreResult<Self> {
724        if self < 0 {
725            return Err(CoreError::ValueError(ErrorContext::new(
726                "Factorial not defined for negative numbers".to_string(),
727            )));
728        }
729
730        let mut result = 1;
731        for i in 2..=self {
732            result *= i;
733        }
734
735        Ok(result)
736    }
737
738    fn binomial(self, k: Self) -> Self {
739        if k < 0 || k > self {
740            return 0;
741        }
742
743        let k = std::cmp::min(k, self - k);
744        let mut result = 1;
745
746        for i in 0..k {
747            result = result * (self - i) / (i + 1);
748        }
749
750        result
751    }
752}
753
754/// Conversion between different numeric types
755pub trait NumericConversion<T> {
756    /// Try to convert to the target type
757    fn try_convert(self) -> Option<T>;
758
759    /// Convert to the target type, or a default value if conversion fails
760    fn convert_or(self, default: T) -> T;
761
762    /// Convert to the target type, or panic if conversion fails
763    fn convert(self) -> T;
764}
765
766/// Implement NumericConversion for all types that implement NumCast
767impl<F, T> NumericConversion<T> for F
768where
769    F: NumCast,
770    T: NumCast,
771{
772    fn try_convert(self) -> Option<T> {
773        num_traits::cast(self)
774    }
775
776    fn convert_or(self, default: T) -> T {
777        self.try_convert().unwrap_or(default)
778    }
779
780    fn convert(self) -> T {
781        self.try_convert().expect("Numeric conversion failed")
782    }
783}
784
785/// Trait for converting between degrees and radians
786pub trait AngleConversion: Sized {
787    /// Convert from degrees to radians
788    fn to_radians(&self) -> CoreResult<Self>
789    where
790        Self: std::marker::Sized;
791
792    /// Convert from radians to degrees
793    fn to_degrees(&self) -> CoreResult<Self>
794    where
795        Self: std::marker::Sized;
796}
797
798/// Implement AngleConversion for all RealNumber types
799impl<T: RealNumber> AngleConversion for T {
800    fn to_radians(&self) -> CoreResult<Self> {
801        let pi = T::from_f64(std::f64::consts::PI).ok_or_else(|| {
802            CoreError::ValueError(ErrorContext::new(
803                "Failed to convert PI constant to target type".to_string(),
804            ))
805        })?;
806        let one_eighty = T::from_f64(180.0).ok_or_else(|| {
807            CoreError::ValueError(ErrorContext::new(
808                "Failed to convert 180.0 to target type".to_string(),
809            ))
810        })?;
811        Ok(*self * pi / one_eighty)
812    }
813
814    fn to_degrees(&self) -> CoreResult<Self> {
815        let pi = T::from_f64(std::f64::consts::PI).ok_or_else(|| {
816            CoreError::ValueError(ErrorContext::new(
817                "Failed to convert PI constant to target type".to_string(),
818            ))
819        })?;
820        let one_eighty = T::from_f64(180.0).ok_or_else(|| {
821            CoreError::ValueError(ErrorContext::new(
822                "Failed to convert 180.0 to target type".to_string(),
823            ))
824        })?;
825        Ok(*self * one_eighty / pi)
826    }
827}
828
829/// Type-safe representation of a unitless quantity
830#[derive(Clone, Copy, PartialEq, Debug)]
831pub struct Scalar<T: ScientificNumber>(pub T);
832
833impl<T: ScientificNumber> Scalar<T> {
834    /// Create a new scalar
835    pub fn new(value: T) -> Self {
836        Scalar(value)
837    }
838
839    /// Get the underlying value
840    pub fn value(&self) -> T {
841        self.0
842    }
843}
844
845impl<T: ScientificNumber> From<T> for Scalar<T> {
846    fn from(value: T) -> Self {
847        Scalar(value)
848    }
849}
850
851impl<T: ScientificNumber> Add for Scalar<T> {
852    type Output = Self;
853
854    fn add(self, other: Self) -> Self::Output {
855        Scalar(self.0 + other.0)
856    }
857}
858
859impl<T: ScientificNumber> Sub for Scalar<T> {
860    type Output = Self;
861
862    fn sub(self, other: Self) -> Self::Output {
863        Scalar(self.0 - other.0)
864    }
865}
866
867impl<T: ScientificNumber> Mul for Scalar<T> {
868    type Output = Self;
869
870    fn mul(self, other: Self) -> Self::Output {
871        Scalar(self.0 * other.0)
872    }
873}
874
875impl<T: ScientificNumber> Div for Scalar<T> {
876    type Output = Self;
877
878    fn div(self, other: Self) -> Self::Output {
879        Scalar(self.0 / other.0)
880    }
881}
882
883impl<T: ScientificNumber + Neg<Output = T>> Neg for Scalar<T> {
884    type Output = Self;
885
886    fn neg(self) -> Self::Output {
887        Scalar(self.0.neg())
888    }
889}
890
891/// Automated precision tracking for numerical computations
892pub mod precision_tracking;
893
894/// Specialized numeric types for scientific domains
895pub mod scientific_types;
896
897/// Arbitrary precision numerical computation support
898#[cfg(feature = "arbitrary-precision")]
899pub mod arbitrary_precision;
900
901/// Numerical stability improvements
902pub mod stability;
903
904/// Stable numerical algorithms
905pub mod stable_algorithms;
906
907/// Advanced-optimized SIMD operations for numerical computations
908///
909/// This module provides vectorized implementations of common mathematical operations
910/// using the highest available SIMD instruction sets for maximum performance.
911pub mod advanced_simd {
912    #[allow(unused_imports)]
913    use super::*;
914
915    #[cfg(target_arch = "aarch64")]
916    use std::arch::aarch64::*;
917    #[cfg(target_arch = "x86_64")]
918    use std::arch::x86_64::*;
919
920    /// Optimized vectorized addition for f32 arrays
921    #[inline]
922    pub fn add_f32_advanced(a: &[f32], b: &[f32], result: &mut [f32]) {
923        debug_assert_eq!(a.len(), b.len());
924        debug_assert_eq!(a.len(), result.len());
925
926        let len = a.len();
927        let mut i = 0;
928
929        // AVX2 path - process 8 elements at a time
930        #[cfg(target_arch = "x86_64")]
931        {
932            if is_x86_feature_detected!("avx2") {
933                unsafe {
934                    while i + 8 <= len {
935                        let va = _mm256_loadu_ps(a.as_ptr().add(i));
936                        let vb = _mm256_loadu_ps(b.as_ptr().add(i));
937                        let vr = _mm256_add_ps(va, vb);
938                        _mm256_storeu_ps(result.as_mut_ptr().add(i), vr);
939                        i += 8;
940                    }
941                }
942            }
943            // SSE path - process 4 elements at a time
944            else if is_x86_feature_detected!("sse") {
945                unsafe {
946                    while i + 4 <= len {
947                        let va = _mm_loadu_ps(a.as_ptr().add(i));
948                        let vb = _mm_loadu_ps(b.as_ptr().add(i));
949                        let vr = _mm_add_ps(va, vb);
950                        _mm_storeu_ps(result.as_mut_ptr().add(i), vr);
951                        i += 4;
952                    }
953                }
954            }
955        }
956
957        // ARM NEON path - process 4 elements at a time
958        #[cfg(target_arch = "aarch64")]
959        {
960            if std::arch::is_aarch64_feature_detected!("neon") {
961                unsafe {
962                    while i + 4 <= len {
963                        let va = vld1q_f32(a.as_ptr().add(i));
964                        let vb = vld1q_f32(b.as_ptr().add(i));
965                        let vr = vaddq_f32(va, vb);
966                        vst1q_f32(result.as_mut_ptr().add(i), vr);
967                        i += 4;
968                    }
969                }
970            }
971        }
972
973        // Scalar fallback for remaining elements
974        while i < len {
975            result[i] = a[i] + b[i];
976            i += 1;
977        }
978    }
979
980    /// Optimized vectorized multiplication for f32 arrays
981    #[inline]
982    pub fn mul_f32_advanced(a: &[f32], b: &[f32], result: &mut [f32]) {
983        debug_assert_eq!(a.len(), b.len());
984        debug_assert_eq!(a.len(), result.len());
985
986        let len = a.len();
987        let mut i = 0;
988
989        // AVX2 + FMA path for maximum performance
990        #[cfg(target_arch = "x86_64")]
991        {
992            if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") {
993                unsafe {
994                    while i + 8 <= len {
995                        let va = _mm256_loadu_ps(a.as_ptr().add(i));
996                        let vb = _mm256_loadu_ps(b.as_ptr().add(i));
997                        let vr = _mm256_mul_ps(va, vb);
998                        _mm256_storeu_ps(result.as_mut_ptr().add(i), vr);
999                        i += 8;
1000                    }
1001                }
1002            } else if is_x86_feature_detected!("sse") {
1003                unsafe {
1004                    while i + 4 <= len {
1005                        let va = _mm_loadu_ps(a.as_ptr().add(i));
1006                        let vb = _mm_loadu_ps(b.as_ptr().add(i));
1007                        let vr = _mm_mul_ps(va, vb);
1008                        _mm_storeu_ps(result.as_mut_ptr().add(i), vr);
1009                        i += 4;
1010                    }
1011                }
1012            }
1013        }
1014
1015        // ARM NEON path
1016        #[cfg(target_arch = "aarch64")]
1017        {
1018            if std::arch::is_aarch64_feature_detected!("neon") {
1019                unsafe {
1020                    while i + 4 <= len {
1021                        let va = vld1q_f32(a.as_ptr().add(i));
1022                        let vb = vld1q_f32(b.as_ptr().add(i));
1023                        let vr = vmulq_f32(va, vb);
1024                        vst1q_f32(result.as_mut_ptr().add(i), vr);
1025                        i += 4;
1026                    }
1027                }
1028            }
1029        }
1030
1031        // Scalar fallback
1032        while i < len {
1033            result[i] = a[i] * b[i];
1034            i += 1;
1035        }
1036    }
1037
1038    /// Optimized fused multiply-add (a * b + c) for f32 arrays
1039    #[inline]
1040    pub fn fma_f32_advanced(a: &[f32], b: &[f32], c: &[f32], result: &mut [f32]) {
1041        debug_assert_eq!(a.len(), b.len());
1042        debug_assert_eq!(a.len(), c.len());
1043        debug_assert_eq!(a.len(), result.len());
1044
1045        let len = a.len();
1046        let mut i = 0;
1047
1048        // AVX2 + FMA path for optimal performance
1049        #[cfg(target_arch = "x86_64")]
1050        {
1051            if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") {
1052                unsafe {
1053                    while i + 8 <= len {
1054                        let va = _mm256_loadu_ps(a.as_ptr().add(i));
1055                        let vb = _mm256_loadu_ps(b.as_ptr().add(i));
1056                        let vc = _mm256_loadu_ps(c.as_ptr().add(i));
1057                        let vr = _mm256_fmadd_ps(va, vb, vc);
1058                        _mm256_storeu_ps(result.as_mut_ptr().add(i), vr);
1059                        i += 8;
1060                    }
1061                }
1062            } else if is_x86_feature_detected!("sse") {
1063                unsafe {
1064                    while i + 4 <= len {
1065                        let va = _mm_loadu_ps(a.as_ptr().add(i));
1066                        let vb = _mm_loadu_ps(b.as_ptr().add(i));
1067                        let vc = _mm_loadu_ps(c.as_ptr().add(i));
1068                        let vr = _mm_add_ps(_mm_mul_ps(va, vb), vc);
1069                        _mm_storeu_ps(result.as_mut_ptr().add(i), vr);
1070                        i += 4;
1071                    }
1072                }
1073            }
1074        }
1075
1076        // ARM NEON path with FMA
1077        #[cfg(target_arch = "aarch64")]
1078        {
1079            if std::arch::is_aarch64_feature_detected!("neon") {
1080                unsafe {
1081                    while i + 4 <= len {
1082                        let va = vld1q_f32(a.as_ptr().add(i));
1083                        let vb = vld1q_f32(b.as_ptr().add(i));
1084                        let vc = vld1q_f32(c.as_ptr().add(i));
1085                        let vr = vfmaq_f32(vc, va, vb);
1086                        vst1q_f32(result.as_mut_ptr().add(i), vr);
1087                        i += 4;
1088                    }
1089                }
1090            }
1091        }
1092
1093        // Scalar fallback
1094        while i < len {
1095            result[i] = a[i] * b[i] + c[0];
1096            i += 1;
1097        }
1098    }
1099
1100    /// Optimized vectorized dot product for f32 arrays
1101    #[inline]
1102    pub fn dot_product_f32_advanced(a: &[f32], b: &[f32]) -> f32 {
1103        debug_assert_eq!(a.len(), b.len());
1104
1105        let len = a.len();
1106        let mut i = 0;
1107        let mut sum = 0.0f32;
1108
1109        // AVX2 + FMA path for maximum throughput
1110        #[cfg(target_arch = "x86_64")]
1111        {
1112            if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") {
1113                unsafe {
1114                    let mut acc = _mm256_setzero_ps();
1115                    while i + 8 <= len {
1116                        let va = _mm256_loadu_ps(a.as_ptr().add(i));
1117                        let vb = _mm256_loadu_ps(b.as_ptr().add(i));
1118                        acc = _mm256_fmadd_ps(va, vb, acc);
1119                        i += 8;
1120                    }
1121                    // Horizontal sum of 8 floats
1122                    let hi = _mm256_extractf128_ps(acc, 1);
1123                    let lo = _mm256_castps256_ps128(acc);
1124                    let sum4 = _mm_add_ps(hi, lo);
1125                    let sum2 = _mm_add_ps(sum4, _mm_movehl_ps(sum4, sum4));
1126                    let sum1 = _mm_add_ss(sum2, _mm_shuffle_ps(sum2, sum2, 1));
1127                    sum = _mm_cvtss_f32(sum1);
1128                }
1129            } else if is_x86_feature_detected!("sse") {
1130                unsafe {
1131                    let mut acc = _mm_setzero_ps();
1132                    while i + 4 <= len {
1133                        let va = _mm_loadu_ps(a.as_ptr().add(i));
1134                        let vb = _mm_loadu_ps(b.as_ptr().add(i));
1135                        acc = _mm_add_ps(acc, _mm_mul_ps(va, vb));
1136                        i += 4;
1137                    }
1138                    // Horizontal sum of 4 floats
1139                    let sum2 = _mm_add_ps(acc, _mm_movehl_ps(acc, acc));
1140                    let sum1 = _mm_add_ss(sum2, _mm_shuffle_ps(sum2, sum2, 1));
1141                    sum = _mm_cvtss_f32(sum1);
1142                }
1143            }
1144        }
1145
1146        // ARM NEON path with accumulation
1147        #[cfg(target_arch = "aarch64")]
1148        {
1149            if std::arch::is_aarch64_feature_detected!("neon") {
1150                unsafe {
1151                    let mut acc = vdupq_n_f32(0.0);
1152                    while i + 4 <= len {
1153                        let va = vld1q_f32(a.as_ptr().add(i));
1154                        let vb = vld1q_f32(b.as_ptr().add(i));
1155                        acc = vfmaq_f32(acc, va, vb);
1156                        i += 4;
1157                    }
1158                    // Horizontal sum
1159                    sum = vaddvq_f32(acc);
1160                }
1161            }
1162        }
1163
1164        // Scalar accumulation for remaining elements
1165        while i < len {
1166            sum += a[i] * b[i];
1167            i += 1;
1168        }
1169
1170        sum
1171    }
1172
1173    /// Optimized vectorized sum reduction for f32 arrays
1174    #[inline]
1175    pub fn sum_f32_advanced(data: &[f32]) -> f32 {
1176        let len = data.len();
1177        let mut i = 0;
1178        let mut sum = 0.0f32;
1179
1180        // AVX2 path
1181        #[cfg(target_arch = "x86_64")]
1182        {
1183            if is_x86_feature_detected!("avx2") {
1184                unsafe {
1185                    let mut acc = _mm256_setzero_ps();
1186                    while i + 8 <= len {
1187                        let v = _mm256_loadu_ps(data.as_ptr().add(i));
1188                        acc = _mm256_add_ps(acc, v);
1189                        i += 8;
1190                    }
1191                    // Horizontal sum
1192                    let hi = _mm256_extractf128_ps(acc, 1);
1193                    let lo = _mm256_castps256_ps128(acc);
1194                    let sum4 = _mm_add_ps(hi, lo);
1195                    let sum2 = _mm_add_ps(sum4, _mm_movehl_ps(sum4, sum4));
1196                    let sum1 = _mm_add_ss(sum2, _mm_shuffle_ps(sum2, sum2, 1));
1197                    sum = _mm_cvtss_f32(sum1);
1198                }
1199            } else if is_x86_feature_detected!("sse") {
1200                unsafe {
1201                    let mut acc = _mm_setzero_ps();
1202                    while i + 4 <= len {
1203                        let v = _mm_loadu_ps(data.as_ptr().add(i));
1204                        acc = _mm_add_ps(acc, v);
1205                        i += 4;
1206                    }
1207                    let sum2 = _mm_add_ps(acc, _mm_movehl_ps(acc, acc));
1208                    let sum1 = _mm_add_ss(sum2, _mm_shuffle_ps(sum2, sum2, 1));
1209                    sum = _mm_cvtss_f32(sum1);
1210                }
1211            }
1212        }
1213
1214        // ARM NEON path
1215        #[cfg(target_arch = "aarch64")]
1216        {
1217            if std::arch::is_aarch64_feature_detected!("neon") {
1218                unsafe {
1219                    let mut acc = vdupq_n_f32(0.0);
1220                    while i + 4 <= len {
1221                        let v = vld1q_f32(data.as_ptr().add(i));
1222                        acc = vaddq_f32(acc, v);
1223                        i += 4;
1224                    }
1225                    sum = vaddvq_f32(acc);
1226                }
1227            }
1228        }
1229
1230        // Scalar accumulation
1231        while i < len {
1232            sum += data[i];
1233            i += 1;
1234        }
1235
1236        sum
1237    }
1238}
1239
1240#[cfg(test)]
1241mod tests {
1242    use super::*;
1243
1244    #[test]
1245    fn test_scientific_number_f32() {
1246        let a: f32 = 3.0;
1247        let b: f32 = 4.0;
1248
1249        assert_eq!(a.abs(), 3.0);
1250        assert_eq!(a.sqrt(), (3.0_f32).sqrt());
1251        assert_eq!(a.square(), 9.0);
1252        assert_eq!(a.max(b), 4.0);
1253        assert_eq!(a.min(b), 3.0);
1254        assert!(a.is_finite());
1255        assert_eq!(a.to_f64(), Some(3.0));
1256        assert_eq!(<f32 as ScientificNumber>::from_f64(3.5), Some(3.5));
1257    }
1258
1259    #[test]
1260    fn test_scientific_number_f64() {
1261        let a: f64 = 3.0;
1262        let b: f64 = 4.0;
1263
1264        assert_eq!(a.abs(), 3.0);
1265        assert_eq!(a.sqrt(), (3.0_f64).sqrt());
1266        assert_eq!(a.square(), 9.0);
1267        assert_eq!(a.max(b), 4.0);
1268        assert_eq!(a.min(b), 3.0);
1269        assert!(a.is_finite());
1270        assert_eq!(a.to_f64(), Some(3.0));
1271        assert_eq!(<f64 as ScientificNumber>::from_f64(3.5), Some(3.5));
1272    }
1273
1274    #[test]
1275    fn test_real_number_f32() {
1276        let a: f32 = 3.0;
1277
1278        assert_eq!(<f32 as RealNumber>::epsilon(), f32::EPSILON);
1279        assert_eq!(a.exp(), (3.0_f32).exp());
1280        assert_eq!(a.ln(), (3.0_f32).ln());
1281        assert_eq!(a.sin(), (3.0_f32).sin());
1282        assert_eq!(a.cos(), (3.0_f32).cos());
1283        assert_eq!(a.powf(2.0), 9.0);
1284    }
1285
1286    #[test]
1287    fn test_real_number_f64() {
1288        let a: f64 = 3.0;
1289
1290        assert_eq!(<f64 as RealNumber>::epsilon(), f64::EPSILON);
1291        assert_eq!(a.exp(), (3.0_f64).exp());
1292        assert_eq!(a.ln(), (3.0_f64).ln());
1293        assert_eq!(a.sin(), (3.0_f64).sin());
1294        assert_eq!(a.cos(), (3.0_f64).cos());
1295        assert_eq!(a.powf(2.0), 9.0);
1296    }
1297
1298    #[test]
1299    fn test_scientific_integer_i32() {
1300        let a: i32 = 12;
1301        let b: i32 = 8;
1302
1303        assert_eq!(a.gcd(b), 4);
1304        assert_eq!(a.lcm(b), 24);
1305        assert!(!a.is_prime());
1306        assert!(11_i32.is_prime());
1307        assert!(a.is_even());
1308        assert!(!a.is_odd());
1309        assert_eq!(a.mod_pow(2, 10), 4); // 12^2 mod 10 = 4
1310        assert_eq!(5_i32.factorial().unwrap(), 120);
1311        assert_eq!(5_i32.binomial(2), 10); // 5 choose 2 = 10
1312    }
1313
1314    #[test]
1315    fn test_numeric_conversion() {
1316        let a: f64 = 3.5;
1317        let b: i32 = a.try_convert().expect("3.5 should convert to i32 as 3");
1318        assert_eq!(b, 3);
1319
1320        let c: f32 = 100.5;
1321        let d: f64 = c.convert();
1322        assert_eq!(d, 100.5);
1323
1324        let e: i32 = 100;
1325        let f: f32 = e.convert();
1326        assert_eq!(f, 100.0);
1327    }
1328
1329    #[test]
1330    fn test_angle_conversion() {
1331        let degrees: f64 = 180.0;
1332        let radians = <f64 as AngleConversion>::to_radians(&degrees).unwrap();
1333        assert!((radians - std::f64::consts::PI).abs() < 1e-10);
1334
1335        let radians: f64 = std::f64::consts::PI / 2.0;
1336        let degrees = <f64 as AngleConversion>::to_degrees(&radians).unwrap();
1337        assert!((degrees - 90.0).abs() < 1e-10);
1338    }
1339
1340    #[test]
1341    fn test_scalar() {
1342        let a = Scalar::new(3.0_f64);
1343        let b = Scalar::new(4.0_f64);
1344
1345        assert_eq!((a + b).value(), 7.0);
1346        assert_eq!((a - b).value(), -1.0);
1347        assert_eq!((a * b).value(), 12.0);
1348        assert_eq!((a / b).value(), 0.75);
1349        assert_eq!((-a).value(), -3.0);
1350
1351        let c: Scalar<f64> = 5.0.into();
1352        assert_eq!(c.value(), 5.0);
1353    }
1354}
scirs2_core/numeric.rs

scirs2_core/
numeric.rs
