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/// Trait for types that can be used as sparse matrix elements
299///
300/// This trait captures the minimal requirements for sparse matrix element types,
301/// supporting both integer and floating-point types. It provides the basic
302/// arithmetic operations and identity elements needed for sparse matrix operations.
303///
304/// # Design Philosophy
305///
306/// Unlike `Float` or `ScientificNumber`, this trait is intentionally minimal to support
307/// the widest range of numeric types in sparse matrices, including:
308/// - Integer types (u8, i32, u64, etc.) for graph adjacency matrices, binary operators, etc.
309/// - Floating-point types (f32, f64) for numerical computation
310///
311/// # Examples
312///
313/// ```
314/// use scirs2_core::numeric::SparseElement;
315///
316/// // Integer sparse matrix element
317/// let a: u8 = 1;
318/// let b: u8 = 2;
319/// assert_eq!(a + b, 3);
320/// assert_eq!(u8::sparse_zero(), 0);
321/// assert_eq!(u8::sparse_one(), 1);
322/// assert!(u8::sparse_zero().is_zero());
323///
324/// // Float sparse matrix element
325/// let x: f64 = 1.0;
326/// let y: f64 = 2.0;
327/// assert_eq!(x + y, 3.0);
328/// assert_eq!(f64::sparse_zero(), 0.0);
329/// assert_eq!(f64::sparse_one(), 1.0);
330/// assert!(f64::sparse_zero().is_zero());
331/// ```
332pub trait SparseElement:
333    Copy + PartialEq + Add<Output = Self> + Sub<Output = Self> + Mul<Output = Self> + Default + Debug
334{
335    /// The zero element (additive identity) for sparse matrix operations
336    ///
337    /// This method is prefixed with `sparse_` to avoid ambiguity with `num_traits::Zero::zero()`.
338    #[must_use]
339    fn sparse_zero() -> Self;
340
341    /// The one element (multiplicative identity) for sparse matrix operations
342    ///
343    /// This method is prefixed with `sparse_` to avoid ambiguity with `num_traits::One::one()`.
344    #[must_use]
345    fn sparse_one() -> Self;
346
347    /// The zero element (additive identity)
348    ///
349    /// # Deprecation Notice
350    ///
351    /// This method is deprecated to avoid ambiguity with `num_traits::Zero::zero()`.
352    /// Use `sparse_zero()` instead.
353    #[deprecated(
354        since = "0.1.0",
355        note = "Use `sparse_zero()` instead to avoid ambiguity with num_traits::Zero"
356    )]
357    #[must_use]
358    fn zero() -> Self {
359        Self::sparse_zero()
360    }
361
362    /// The one element (multiplicative identity)
363    ///
364    /// # Deprecation Notice
365    ///
366    /// This method is deprecated to avoid ambiguity with `num_traits::One::one()`.
367    /// Use `sparse_one()` instead.
368    #[deprecated(
369        since = "0.1.0",
370        note = "Use `sparse_one()` instead to avoid ambiguity with num_traits::One"
371    )]
372    #[must_use]
373    fn one() -> Self {
374        Self::sparse_one()
375    }
376
377    /// Check if value is zero
378    #[must_use]
379    fn is_zero(&self) -> bool {
380        *self == Self::sparse_zero()
381    }
382}
383
384// Implement SparseElement for unsigned integer types
385impl SparseElement for u8 {
386    fn sparse_zero() -> Self {
387        0
388    }
389    fn sparse_one() -> Self {
390        1
391    }
392}
393
394impl SparseElement for u16 {
395    fn sparse_zero() -> Self {
396        0
397    }
398    fn sparse_one() -> Self {
399        1
400    }
401}
402
403impl SparseElement for u32 {
404    fn sparse_zero() -> Self {
405        0
406    }
407    fn sparse_one() -> Self {
408        1
409    }
410}
411
412impl SparseElement for u64 {
413    fn sparse_zero() -> Self {
414        0
415    }
416    fn sparse_one() -> Self {
417        1
418    }
419}
420
421impl SparseElement for u128 {
422    fn sparse_zero() -> Self {
423        0
424    }
425    fn sparse_one() -> Self {
426        1
427    }
428}
429
430impl SparseElement for usize {
431    fn sparse_zero() -> Self {
432        0
433    }
434    fn sparse_one() -> Self {
435        1
436    }
437}
438
439// Implement SparseElement for signed integer types
440impl SparseElement for i8 {
441    fn sparse_zero() -> Self {
442        0
443    }
444    fn sparse_one() -> Self {
445        1
446    }
447}
448
449impl SparseElement for i16 {
450    fn sparse_zero() -> Self {
451        0
452    }
453    fn sparse_one() -> Self {
454        1
455    }
456}
457
458impl SparseElement for i32 {
459    fn sparse_zero() -> Self {
460        0
461    }
462    fn sparse_one() -> Self {
463        1
464    }
465}
466
467impl SparseElement for i64 {
468    fn sparse_zero() -> Self {
469        0
470    }
471    fn sparse_one() -> Self {
472        1
473    }
474}
475
476impl SparseElement for i128 {
477    fn sparse_zero() -> Self {
478        0
479    }
480    fn sparse_one() -> Self {
481        1
482    }
483}
484
485impl SparseElement for isize {
486    fn sparse_zero() -> Self {
487        0
488    }
489    fn sparse_one() -> Self {
490        1
491    }
492}
493
494// Implement SparseElement for floating-point types
495impl SparseElement for f32 {
496    fn sparse_zero() -> Self {
497        0.0
498    }
499    fn sparse_one() -> Self {
500        1.0
501    }
502}
503
504impl SparseElement for f64 {
505    fn sparse_zero() -> Self {
506        0.0
507    }
508    fn sparse_one() -> Self {
509        1.0
510    }
511}
512
513// Implement ScientificNumber for f32
514impl ScientificNumber for f32 {
515    fn abs(self) -> Self {
516        self.abs()
517    }
518
519    fn sqrt(self) -> Self {
520        self.sqrt()
521    }
522
523    fn max(self, other: Self) -> Self {
524        self.max(other)
525    }
526
527    fn min(self, other: Self) -> Self {
528        self.min(other)
529    }
530
531    fn is_finite(self) -> bool {
532        self.is_finite()
533    }
534
535    fn to_f64(self) -> Option<f64> {
536        Some(self as f64)
537    }
538
539    fn from_f64(value: f64) -> Option<Self> {
540        if value.is_finite() && value <= Self::MAX as f64 && value >= Self::MIN as f64 {
541            Some(value as f32)
542        } else {
543            None
544        }
545    }
546
547    fn from_le_bytes(bytes: &[u8]) -> Self {
548        let mut array = [0u8; 4];
549        array.copy_from_slice(&bytes[..4]);
550        f32::from_le_bytes(array)
551    }
552
553    fn from_be_bytes(bytes: &[u8]) -> Self {
554        let mut array = [0u8; 4];
555        array.copy_from_slice(&bytes[..4]);
556        f32::from_be_bytes(array)
557    }
558
559    fn to_le_bytes(self) -> Vec<u8> {
560        self.to_le_bytes().to_vec()
561    }
562
563    fn to_be_bytes(self) -> Vec<u8> {
564        self.to_be_bytes().to_vec()
565    }
566}
567
568// Implement RealNumber for f32
569impl RealNumber for f32 {
570    fn epsilon() -> Self {
571        Self::EPSILON
572    }
573
574    fn exp(self) -> Self {
575        self.exp()
576    }
577
578    fn ln(self) -> Self {
579        self.ln()
580    }
581
582    fn log10(self) -> Self {
583        self.log10()
584    }
585
586    fn log2(self) -> Self {
587        self.log2()
588    }
589
590    fn sin(self) -> Self {
591        self.sin()
592    }
593
594    fn cos(self) -> Self {
595        self.cos()
596    }
597
598    fn tan(self) -> Self {
599        self.tan()
600    }
601
602    fn sinh(self) -> Self {
603        self.sinh()
604    }
605
606    fn cosh(self) -> Self {
607        self.cosh()
608    }
609
610    fn tanh(self) -> Self {
611        self.tanh()
612    }
613
614    fn asin(self) -> Self {
615        self.asin()
616    }
617
618    fn acos(self) -> Self {
619        self.acos()
620    }
621
622    fn atan(self) -> Self {
623        self.atan()
624    }
625
626    fn atan2(self, other: Self) -> Self {
627        self.atan2(other)
628    }
629
630    fn powf(self, n: Self) -> Self {
631        self.powf(n)
632    }
633
634    fn powi(self, n: i32) -> Self {
635        self.powi(n)
636    }
637
638    fn factorial(self) -> Self {
639        if self < 0.0 {
640            return Self::NAN;
641        }
642
643        // Use Stirling's approximation for non-integers or large values
644        if self != self.trunc() || self > 100.0 {
645            const SQRT_TWO_PI: f32 = 2.506_628_3;
646            return SQRT_TWO_PI * self.powf(self + 0.5) * (-self).exp();
647        }
648
649        let mut result = 1.0;
650        let n = self as u32;
651        for i in 2..=n {
652            result *= i as f32;
653        }
654
655        result
656    }
657}
658
659// Implement ScientificNumber for f64
660impl ScientificNumber for f64 {
661    fn abs(self) -> Self {
662        self.abs()
663    }
664
665    fn sqrt(self) -> Self {
666        self.sqrt()
667    }
668
669    fn max(self, other: Self) -> Self {
670        self.max(other)
671    }
672
673    fn min(self, other: Self) -> Self {
674        self.min(other)
675    }
676
677    fn is_finite(self) -> bool {
678        self.is_finite()
679    }
680
681    fn to_f64(self) -> Option<f64> {
682        Some(self)
683    }
684
685    fn from_f64(value: f64) -> Option<Self> {
686        Some(value)
687    }
688
689    fn from_le_bytes(bytes: &[u8]) -> Self {
690        let mut array = [0u8; 8];
691        array.copy_from_slice(&bytes[..8]);
692        f64::from_le_bytes(array)
693    }
694
695    fn from_be_bytes(bytes: &[u8]) -> Self {
696        let mut array = [0u8; 8];
697        array.copy_from_slice(&bytes[..8]);
698        f64::from_be_bytes(array)
699    }
700
701    fn to_le_bytes(self) -> Vec<u8> {
702        self.to_le_bytes().to_vec()
703    }
704
705    fn to_be_bytes(self) -> Vec<u8> {
706        self.to_be_bytes().to_vec()
707    }
708}
709
710// Implement RealNumber for f64
711impl RealNumber for f64 {
712    fn epsilon() -> Self {
713        Self::EPSILON
714    }
715
716    fn exp(self) -> Self {
717        self.exp()
718    }
719
720    fn ln(self) -> Self {
721        self.ln()
722    }
723
724    fn log10(self) -> Self {
725        self.log10()
726    }
727
728    fn log2(self) -> Self {
729        self.log2()
730    }
731
732    fn sin(self) -> Self {
733        self.sin()
734    }
735
736    fn cos(self) -> Self {
737        self.cos()
738    }
739
740    fn tan(self) -> Self {
741        self.tan()
742    }
743
744    fn sinh(self) -> Self {
745        self.sinh()
746    }
747
748    fn cosh(self) -> Self {
749        self.cosh()
750    }
751
752    fn tanh(self) -> Self {
753        self.tanh()
754    }
755
756    fn asin(self) -> Self {
757        self.asin()
758    }
759
760    fn acos(self) -> Self {
761        self.acos()
762    }
763
764    fn atan(self) -> Self {
765        self.atan()
766    }
767
768    fn atan2(self, other: Self) -> Self {
769        self.atan2(other)
770    }
771
772    fn powf(self, n: Self) -> Self {
773        self.powf(n)
774    }
775
776    fn powi(self, n: i32) -> Self {
777        self.powi(n)
778    }
779
780    fn factorial(self) -> Self {
781        if self < 0.0 {
782            return Self::NAN;
783        }
784
785        // Use Stirling's approximation for non-integers or large values
786        if self != self.trunc() || self > 170.0 {
787            const SQRT_TWO_PI: f64 = 2.506_628_274_631_000_2;
788            return SQRT_TWO_PI * self.powf(self + 0.5) * (-self).exp();
789        }
790
791        let mut result = 1.0;
792        let n = self as u32;
793        for i in 2..=n {
794            result *= i as f64;
795        }
796
797        result
798    }
799}
800
801// Implement ScientificNumber for i32
802impl ScientificNumber for i32 {
803    fn abs(self) -> Self {
804        self.abs()
805    }
806
807    fn sqrt(self) -> Self {
808        (self as f64).sqrt() as i32
809    }
810
811    fn max(self, other: Self) -> Self {
812        std::cmp::max(self, other)
813    }
814
815    fn min(self, other: Self) -> Self {
816        std::cmp::min(self, other)
817    }
818
819    fn is_finite(self) -> bool {
820        true // integers are always finite
821    }
822
823    fn to_f64(self) -> Option<f64> {
824        Some(self as f64)
825    }
826
827    fn from_f64(value: f64) -> Option<Self> {
828        if value.is_finite() && value <= i32::MAX as f64 && value >= i32::MIN as f64 {
829            Some(value as i32)
830        } else {
831            None
832        }
833    }
834
835    fn from_le_bytes(bytes: &[u8]) -> Self {
836        let mut array = [0u8; 4];
837        array.copy_from_slice(&bytes[..4]);
838        i32::from_le_bytes(array)
839    }
840
841    fn from_be_bytes(bytes: &[u8]) -> Self {
842        let mut array = [0u8; 4];
843        array.copy_from_slice(&bytes[..4]);
844        i32::from_be_bytes(array)
845    }
846
847    fn to_le_bytes(self) -> Vec<u8> {
848        self.to_le_bytes().to_vec()
849    }
850
851    fn to_be_bytes(self) -> Vec<u8> {
852        self.to_be_bytes().to_vec()
853    }
854}
855
856// Implement ScientificInteger for i32
857impl ScientificInteger for i32 {
858    fn gcd(self, other: Self) -> Self {
859        let mut a = self.abs();
860        let mut b = other.abs();
861
862        if a == 0 {
863            return b;
864        }
865        if b == 0 {
866            return a;
867        }
868
869        while b != 0 {
870            let temp = b;
871            b = a % b;
872            a = temp;
873        }
874
875        a
876    }
877
878    fn lcm(self, other: Self) -> Self {
879        if self == 0 || other == 0 {
880            return 0;
881        }
882
883        let gcd = self.gcd(other);
884
885        (self / gcd) * other
886    }
887
888    fn is_prime(self) -> bool {
889        if self <= 1 {
890            return false;
891        }
892        if self <= 3 {
893            return true;
894        }
895        if self % 2 == 0 || self % 3 == 0 {
896            return false;
897        }
898
899        let mut i = 5;
900        while i * i <= self {
901            if self % i == 0 || self % (i + 2) == 0 {
902                return false;
903            }
904            i += 6;
905        }
906
907        true
908    }
909
910    fn is_even(self) -> bool {
911        self % 2 == 0
912    }
913
914    fn is_odd(self) -> bool {
915        self % 2 != 0
916    }
917
918    fn mod_pow(self, exp: Self, modulus: Self) -> Self {
919        if modulus == 1 {
920            return 0;
921        }
922
923        let mut base = self % modulus;
924        let mut result = 1;
925        let mut exponent = exp;
926
927        while exponent > 0 {
928            if exponent % 2 == 1 {
929                result = (result * base) % modulus;
930            }
931            exponent >>= 1;
932            base = (base * base) % modulus;
933        }
934
935        result
936    }
937
938    fn factorial(self) -> CoreResult<Self> {
939        if self < 0 {
940            return Err(CoreError::ValueError(ErrorContext::new(
941                "Factorial not defined for negative numbers".to_string(),
942            )));
943        }
944
945        let mut result = 1;
946        for i in 2..=self {
947            result *= i;
948        }
949
950        Ok(result)
951    }
952
953    fn binomial(self, k: Self) -> Self {
954        if k < 0 || k > self {
955            return 0;
956        }
957
958        let k = std::cmp::min(k, self - k);
959        let mut result = 1;
960
961        for i in 0..k {
962            result = result * (self - i) / (i + 1);
963        }
964
965        result
966    }
967}
968
969/// Conversion between different numeric types
970pub trait NumericConversion<T> {
971    /// Try to convert to the target type
972    fn try_convert(self) -> Option<T>;
973
974    /// Convert to the target type, or a default value if conversion fails
975    fn convert_or(self, default: T) -> T;
976
977    /// Convert to the target type, or panic if conversion fails
978    fn convert(self) -> T;
979}
980
981/// Implement NumericConversion for all types that implement NumCast
982impl<F, T> NumericConversion<T> for F
983where
984    F: NumCast,
985    T: NumCast,
986{
987    fn try_convert(self) -> Option<T> {
988        num_traits::cast(self)
989    }
990
991    fn convert_or(self, default: T) -> T {
992        self.try_convert().unwrap_or(default)
993    }
994
995    fn convert(self) -> T {
996        self.try_convert().expect("Numeric conversion failed")
997    }
998}
999
1000/// Trait for converting between degrees and radians
1001pub trait AngleConversion: Sized {
1002    /// Convert from degrees to radians
1003    fn to_radians(&self) -> CoreResult<Self>
1004    where
1005        Self: std::marker::Sized;
1006
1007    /// Convert from radians to degrees
1008    fn to_degrees(&self) -> CoreResult<Self>
1009    where
1010        Self: std::marker::Sized;
1011}
1012
1013/// Implement AngleConversion for all RealNumber types
1014impl<T: RealNumber> AngleConversion for T {
1015    fn to_radians(&self) -> CoreResult<Self> {
1016        let pi = T::from_f64(std::f64::consts::PI).ok_or_else(|| {
1017            CoreError::ValueError(ErrorContext::new(
1018                "Failed to convert PI constant to target type".to_string(),
1019            ))
1020        })?;
1021        let one_eighty = T::from_f64(180.0).ok_or_else(|| {
1022            CoreError::ValueError(ErrorContext::new(
1023                "Failed to convert 180.0 to target type".to_string(),
1024            ))
1025        })?;
1026        Ok(*self * pi / one_eighty)
1027    }
1028
1029    fn to_degrees(&self) -> CoreResult<Self> {
1030        let pi = T::from_f64(std::f64::consts::PI).ok_or_else(|| {
1031            CoreError::ValueError(ErrorContext::new(
1032                "Failed to convert PI constant to target type".to_string(),
1033            ))
1034        })?;
1035        let one_eighty = T::from_f64(180.0).ok_or_else(|| {
1036            CoreError::ValueError(ErrorContext::new(
1037                "Failed to convert 180.0 to target type".to_string(),
1038            ))
1039        })?;
1040        Ok(*self * one_eighty / pi)
1041    }
1042}
1043
1044/// Type-safe representation of a unitless quantity
1045#[derive(Clone, Copy, PartialEq, Debug)]
1046pub struct Scalar<T: ScientificNumber>(pub T);
1047
1048impl<T: ScientificNumber> Scalar<T> {
1049    /// Create a new scalar
1050    pub fn new(value: T) -> Self {
1051        Scalar(value)
1052    }
1053
1054    /// Get the underlying value
1055    pub fn value(&self) -> T {
1056        self.0
1057    }
1058}
1059
1060impl<T: ScientificNumber> From<T> for Scalar<T> {
1061    fn from(value: T) -> Self {
1062        Scalar(value)
1063    }
1064}
1065
1066impl<T: ScientificNumber> Add for Scalar<T> {
1067    type Output = Self;
1068
1069    fn add(self, other: Self) -> Self::Output {
1070        Scalar(self.0 + other.0)
1071    }
1072}
1073
1074impl<T: ScientificNumber> Sub for Scalar<T> {
1075    type Output = Self;
1076
1077    fn sub(self, other: Self) -> Self::Output {
1078        Scalar(self.0 - other.0)
1079    }
1080}
1081
1082impl<T: ScientificNumber> Mul for Scalar<T> {
1083    type Output = Self;
1084
1085    fn mul(self, other: Self) -> Self::Output {
1086        Scalar(self.0 * other.0)
1087    }
1088}
1089
1090impl<T: ScientificNumber> Div for Scalar<T> {
1091    type Output = Self;
1092
1093    fn div(self, other: Self) -> Self::Output {
1094        Scalar(self.0 / other.0)
1095    }
1096}
1097
1098impl<T: ScientificNumber + Neg<Output = T>> Neg for Scalar<T> {
1099    type Output = Self;
1100
1101    fn neg(self) -> Self::Output {
1102        Scalar(self.0.neg())
1103    }
1104}
1105
1106/// Numeric utilities: approximate equality, ULP distance, relative error, safe clamping
1107pub mod utilities;
1108
1109/// Automated precision tracking for numerical computations
1110pub mod precision_tracking;
1111
1112/// Specialized numeric types for scientific domains
1113pub mod scientific_types;
1114
1115/// Arbitrary precision numerical computation support
1116#[cfg(feature = "arbitrary-precision")]
1117pub mod arbitrary_precision;
1118
1119/// Numerical stability improvements
1120pub mod stability;
1121
1122/// Stable numerical algorithms
1123pub mod stable_algorithms;
1124
1125/// Extended numeric stability toolkit (gradients, error helpers, array-based activations)
1126pub mod stability_toolkit;
1127
1128/// Advanced-optimized SIMD operations for numerical computations
1129///
1130/// This module provides vectorized implementations of common mathematical operations
1131/// using the highest available SIMD instruction sets for maximum performance.
1132pub mod advanced_simd {
1133    #[allow(unused_imports)]
1134    use super::*;
1135
1136    #[cfg(target_arch = "aarch64")]
1137    use std::arch::aarch64::*;
1138    #[cfg(target_arch = "x86_64")]
1139    use std::arch::x86_64::*;
1140
1141    /// Optimized vectorized addition for f32 arrays
1142    #[inline]
1143    pub fn add_f32_advanced(a: &[f32], b: &[f32], result: &mut [f32]) {
1144        debug_assert_eq!(a.len(), b.len());
1145        debug_assert_eq!(a.len(), result.len());
1146
1147        let len = a.len();
1148        let mut i = 0;
1149
1150        // AVX2 path - process 8 elements at a time
1151        #[cfg(target_arch = "x86_64")]
1152        {
1153            if is_x86_feature_detected!("avx2") {
1154                unsafe {
1155                    while i + 8 <= len {
1156                        let va = _mm256_loadu_ps(a.as_ptr().add(i));
1157                        let vb = _mm256_loadu_ps(b.as_ptr().add(i));
1158                        let vr = _mm256_add_ps(va, vb);
1159                        _mm256_storeu_ps(result.as_mut_ptr().add(i), vr);
1160                        i += 8;
1161                    }
1162                }
1163            }
1164            // SSE path - process 4 elements at a time
1165            else if is_x86_feature_detected!("sse") {
1166                unsafe {
1167                    while i + 4 <= len {
1168                        let va = _mm_loadu_ps(a.as_ptr().add(i));
1169                        let vb = _mm_loadu_ps(b.as_ptr().add(i));
1170                        let vr = _mm_add_ps(va, vb);
1171                        _mm_storeu_ps(result.as_mut_ptr().add(i), vr);
1172                        i += 4;
1173                    }
1174                }
1175            }
1176        }
1177
1178        // ARM NEON path - process 4 elements at a time
1179        #[cfg(target_arch = "aarch64")]
1180        {
1181            if std::arch::is_aarch64_feature_detected!("neon") {
1182                unsafe {
1183                    while i + 4 <= len {
1184                        let va = vld1q_f32(a.as_ptr().add(i));
1185                        let vb = vld1q_f32(b.as_ptr().add(i));
1186                        let vr = vaddq_f32(va, vb);
1187                        vst1q_f32(result.as_mut_ptr().add(i), vr);
1188                        i += 4;
1189                    }
1190                }
1191            }
1192        }
1193
1194        // Scalar fallback for remaining elements
1195        while i < len {
1196            result[i] = a[i] + b[i];
1197            i += 1;
1198        }
1199    }
1200
1201    /// Optimized vectorized multiplication for f32 arrays
1202    #[inline]
1203    pub fn mul_f32_advanced(a: &[f32], b: &[f32], result: &mut [f32]) {
1204        debug_assert_eq!(a.len(), b.len());
1205        debug_assert_eq!(a.len(), result.len());
1206
1207        let len = a.len();
1208        let mut i = 0;
1209
1210        // AVX2 + FMA path for maximum performance
1211        #[cfg(target_arch = "x86_64")]
1212        {
1213            if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") {
1214                unsafe {
1215                    while i + 8 <= len {
1216                        let va = _mm256_loadu_ps(a.as_ptr().add(i));
1217                        let vb = _mm256_loadu_ps(b.as_ptr().add(i));
1218                        let vr = _mm256_mul_ps(va, vb);
1219                        _mm256_storeu_ps(result.as_mut_ptr().add(i), vr);
1220                        i += 8;
1221                    }
1222                }
1223            } else if is_x86_feature_detected!("sse") {
1224                unsafe {
1225                    while i + 4 <= len {
1226                        let va = _mm_loadu_ps(a.as_ptr().add(i));
1227                        let vb = _mm_loadu_ps(b.as_ptr().add(i));
1228                        let vr = _mm_mul_ps(va, vb);
1229                        _mm_storeu_ps(result.as_mut_ptr().add(i), vr);
1230                        i += 4;
1231                    }
1232                }
1233            }
1234        }
1235
1236        // ARM NEON path
1237        #[cfg(target_arch = "aarch64")]
1238        {
1239            if std::arch::is_aarch64_feature_detected!("neon") {
1240                unsafe {
1241                    while i + 4 <= len {
1242                        let va = vld1q_f32(a.as_ptr().add(i));
1243                        let vb = vld1q_f32(b.as_ptr().add(i));
1244                        let vr = vmulq_f32(va, vb);
1245                        vst1q_f32(result.as_mut_ptr().add(i), vr);
1246                        i += 4;
1247                    }
1248                }
1249            }
1250        }
1251
1252        // Scalar fallback
1253        while i < len {
1254            result[i] = a[i] * b[i];
1255            i += 1;
1256        }
1257    }
1258
1259    /// Optimized fused multiply-add (a * b + c) for f32 arrays
1260    #[inline]
1261    pub fn fma_f32_advanced(a: &[f32], b: &[f32], c: &[f32], result: &mut [f32]) {
1262        debug_assert_eq!(a.len(), b.len());
1263        debug_assert_eq!(a.len(), c.len());
1264        debug_assert_eq!(a.len(), result.len());
1265
1266        let len = a.len();
1267        let mut i = 0;
1268
1269        // AVX2 + FMA path for optimal performance
1270        #[cfg(target_arch = "x86_64")]
1271        {
1272            if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") {
1273                unsafe {
1274                    while i + 8 <= len {
1275                        let va = _mm256_loadu_ps(a.as_ptr().add(i));
1276                        let vb = _mm256_loadu_ps(b.as_ptr().add(i));
1277                        let vc = _mm256_loadu_ps(c.as_ptr().add(i));
1278                        let vr = _mm256_fmadd_ps(va, vb, vc);
1279                        _mm256_storeu_ps(result.as_mut_ptr().add(i), vr);
1280                        i += 8;
1281                    }
1282                }
1283            } else if is_x86_feature_detected!("sse") {
1284                unsafe {
1285                    while i + 4 <= len {
1286                        let va = _mm_loadu_ps(a.as_ptr().add(i));
1287                        let vb = _mm_loadu_ps(b.as_ptr().add(i));
1288                        let vc = _mm_loadu_ps(c.as_ptr().add(i));
1289                        let vr = _mm_add_ps(_mm_mul_ps(va, vb), vc);
1290                        _mm_storeu_ps(result.as_mut_ptr().add(i), vr);
1291                        i += 4;
1292                    }
1293                }
1294            }
1295        }
1296
1297        // ARM NEON path with FMA
1298        #[cfg(target_arch = "aarch64")]
1299        {
1300            if std::arch::is_aarch64_feature_detected!("neon") {
1301                unsafe {
1302                    while i + 4 <= len {
1303                        let va = vld1q_f32(a.as_ptr().add(i));
1304                        let vb = vld1q_f32(b.as_ptr().add(i));
1305                        let vc = vld1q_f32(c.as_ptr().add(i));
1306                        let vr = vfmaq_f32(vc, va, vb);
1307                        vst1q_f32(result.as_mut_ptr().add(i), vr);
1308                        i += 4;
1309                    }
1310                }
1311            }
1312        }
1313
1314        // Scalar fallback
1315        while i < len {
1316            result[i] = a[i] * b[i] + c[0];
1317            i += 1;
1318        }
1319    }
1320
1321    /// Optimized vectorized dot product for f32 arrays
1322    #[inline]
1323    pub fn dot_product_f32_advanced(a: &[f32], b: &[f32]) -> f32 {
1324        debug_assert_eq!(a.len(), b.len());
1325
1326        let len = a.len();
1327        let mut i = 0;
1328        let mut sum = 0.0f32;
1329
1330        // AVX2 + FMA path for maximum throughput
1331        #[cfg(target_arch = "x86_64")]
1332        {
1333            if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") {
1334                unsafe {
1335                    let mut acc = _mm256_setzero_ps();
1336                    while i + 8 <= len {
1337                        let va = _mm256_loadu_ps(a.as_ptr().add(i));
1338                        let vb = _mm256_loadu_ps(b.as_ptr().add(i));
1339                        acc = _mm256_fmadd_ps(va, vb, acc);
1340                        i += 8;
1341                    }
1342                    // Horizontal sum of 8 floats
1343                    let hi = _mm256_extractf128_ps(acc, 1);
1344                    let lo = _mm256_castps256_ps128(acc);
1345                    let sum4 = _mm_add_ps(hi, lo);
1346                    let sum2 = _mm_add_ps(sum4, _mm_movehl_ps(sum4, sum4));
1347                    let sum1 = _mm_add_ss(sum2, _mm_shuffle_ps(sum2, sum2, 1));
1348                    sum = _mm_cvtss_f32(sum1);
1349                }
1350            } else if is_x86_feature_detected!("sse") {
1351                unsafe {
1352                    let mut acc = _mm_setzero_ps();
1353                    while i + 4 <= len {
1354                        let va = _mm_loadu_ps(a.as_ptr().add(i));
1355                        let vb = _mm_loadu_ps(b.as_ptr().add(i));
1356                        acc = _mm_add_ps(acc, _mm_mul_ps(va, vb));
1357                        i += 4;
1358                    }
1359                    // Horizontal sum of 4 floats
1360                    let sum2 = _mm_add_ps(acc, _mm_movehl_ps(acc, acc));
1361                    let sum1 = _mm_add_ss(sum2, _mm_shuffle_ps(sum2, sum2, 1));
1362                    sum = _mm_cvtss_f32(sum1);
1363                }
1364            }
1365        }
1366
1367        // ARM NEON path with accumulation
1368        #[cfg(target_arch = "aarch64")]
1369        {
1370            if std::arch::is_aarch64_feature_detected!("neon") {
1371                unsafe {
1372                    let mut acc = vdupq_n_f32(0.0);
1373                    while i + 4 <= len {
1374                        let va = vld1q_f32(a.as_ptr().add(i));
1375                        let vb = vld1q_f32(b.as_ptr().add(i));
1376                        acc = vfmaq_f32(acc, va, vb);
1377                        i += 4;
1378                    }
1379                    // Horizontal sum
1380                    sum = vaddvq_f32(acc);
1381                }
1382            }
1383        }
1384
1385        // Scalar accumulation for remaining elements
1386        while i < len {
1387            sum += a[i] * b[i];
1388            i += 1;
1389        }
1390
1391        sum
1392    }
1393
1394    /// Optimized vectorized sum reduction for f32 arrays
1395    #[inline]
1396    pub fn sum_f32_advanced(data: &[f32]) -> f32 {
1397        let len = data.len();
1398        let mut i = 0;
1399        let mut sum = 0.0f32;
1400
1401        // AVX2 path
1402        #[cfg(target_arch = "x86_64")]
1403        {
1404            if is_x86_feature_detected!("avx2") {
1405                unsafe {
1406                    let mut acc = _mm256_setzero_ps();
1407                    while i + 8 <= len {
1408                        let v = _mm256_loadu_ps(data.as_ptr().add(i));
1409                        acc = _mm256_add_ps(acc, v);
1410                        i += 8;
1411                    }
1412                    // Horizontal sum
1413                    let hi = _mm256_extractf128_ps(acc, 1);
1414                    let lo = _mm256_castps256_ps128(acc);
1415                    let sum4 = _mm_add_ps(hi, lo);
1416                    let sum2 = _mm_add_ps(sum4, _mm_movehl_ps(sum4, sum4));
1417                    let sum1 = _mm_add_ss(sum2, _mm_shuffle_ps(sum2, sum2, 1));
1418                    sum = _mm_cvtss_f32(sum1);
1419                }
1420            } else if is_x86_feature_detected!("sse") {
1421                unsafe {
1422                    let mut acc = _mm_setzero_ps();
1423                    while i + 4 <= len {
1424                        let v = _mm_loadu_ps(data.as_ptr().add(i));
1425                        acc = _mm_add_ps(acc, v);
1426                        i += 4;
1427                    }
1428                    let sum2 = _mm_add_ps(acc, _mm_movehl_ps(acc, acc));
1429                    let sum1 = _mm_add_ss(sum2, _mm_shuffle_ps(sum2, sum2, 1));
1430                    sum = _mm_cvtss_f32(sum1);
1431                }
1432            }
1433        }
1434
1435        // ARM NEON path
1436        #[cfg(target_arch = "aarch64")]
1437        {
1438            if std::arch::is_aarch64_feature_detected!("neon") {
1439                unsafe {
1440                    let mut acc = vdupq_n_f32(0.0);
1441                    while i + 4 <= len {
1442                        let v = vld1q_f32(data.as_ptr().add(i));
1443                        acc = vaddq_f32(acc, v);
1444                        i += 4;
1445                    }
1446                    sum = vaddvq_f32(acc);
1447                }
1448            }
1449        }
1450
1451        // Scalar accumulation
1452        while i < len {
1453            sum += data[i];
1454            i += 1;
1455        }
1456
1457        sum
1458    }
1459}
1460
1461#[cfg(test)]
1462mod tests {
1463    use super::*;
1464
1465    #[test]
1466    fn test_scientific_number_f32() {
1467        let a: f32 = 3.0;
1468        let b: f32 = 4.0;
1469
1470        assert_eq!(a.abs(), 3.0);
1471        assert_eq!(a.sqrt(), (3.0_f32).sqrt());
1472        assert_eq!(a.square(), 9.0);
1473        assert_eq!(a.max(b), 4.0);
1474        assert_eq!(a.min(b), 3.0);
1475        assert!(a.is_finite());
1476        assert_eq!(a.to_f64(), Some(3.0));
1477        assert_eq!(<f32 as ScientificNumber>::from_f64(3.5), Some(3.5));
1478    }
1479
1480    #[test]
1481    fn test_scientific_number_f64() {
1482        let a: f64 = 3.0;
1483        let b: f64 = 4.0;
1484
1485        assert_eq!(a.abs(), 3.0);
1486        assert_eq!(a.sqrt(), (3.0_f64).sqrt());
1487        assert_eq!(a.square(), 9.0);
1488        assert_eq!(a.max(b), 4.0);
1489        assert_eq!(a.min(b), 3.0);
1490        assert!(a.is_finite());
1491        assert_eq!(a.to_f64(), Some(3.0));
1492        assert_eq!(<f64 as ScientificNumber>::from_f64(3.5), Some(3.5));
1493    }
1494
1495    #[test]
1496    fn test_real_number_f32() {
1497        let a: f32 = 3.0;
1498
1499        assert_eq!(<f32 as RealNumber>::epsilon(), f32::EPSILON);
1500        assert_eq!(a.exp(), (3.0_f32).exp());
1501        assert_eq!(a.ln(), (3.0_f32).ln());
1502        assert_eq!(a.sin(), (3.0_f32).sin());
1503        assert_eq!(a.cos(), (3.0_f32).cos());
1504        assert_eq!(a.powf(2.0), 9.0);
1505    }
1506
1507    #[test]
1508    fn test_real_number_f64() {
1509        let a: f64 = 3.0;
1510
1511        assert_eq!(<f64 as RealNumber>::epsilon(), f64::EPSILON);
1512        assert_eq!(a.exp(), (3.0_f64).exp());
1513        assert_eq!(a.ln(), (3.0_f64).ln());
1514        assert_eq!(a.sin(), (3.0_f64).sin());
1515        assert_eq!(a.cos(), (3.0_f64).cos());
1516        assert_eq!(a.powf(2.0), 9.0);
1517    }
1518
1519    #[test]
1520    fn test_scientific_integer_i32() {
1521        let a: i32 = 12;
1522        let b: i32 = 8;
1523
1524        assert_eq!(a.gcd(b), 4);
1525        assert_eq!(a.lcm(b), 24);
1526        assert!(!a.is_prime());
1527        assert!(11_i32.is_prime());
1528        assert!(a.is_even());
1529        assert!(!a.is_odd());
1530        assert_eq!(a.mod_pow(2, 10), 4); // 12^2 mod 10 = 4
1531        assert_eq!(5_i32.factorial().expect("Operation failed"), 120);
1532        assert_eq!(5_i32.binomial(2), 10); // 5 choose 2 = 10
1533    }
1534
1535    #[test]
1536    fn test_numeric_conversion() {
1537        let a: f64 = 3.5;
1538        let b: i32 = a.try_convert().expect("3.5 should convert to i32 as 3");
1539        assert_eq!(b, 3);
1540
1541        let c: f32 = 100.5;
1542        let d: f64 = c.convert();
1543        assert_eq!(d, 100.5);
1544
1545        let e: i32 = 100;
1546        let f: f32 = e.convert();
1547        assert_eq!(f, 100.0);
1548    }
1549
1550    #[test]
1551    fn test_angle_conversion() {
1552        let degrees: f64 = 180.0;
1553        let radians = <f64 as AngleConversion>::to_radians(&degrees).expect("Operation failed");
1554        assert!((radians - std::f64::consts::PI).abs() < 1e-10);
1555
1556        let radians: f64 = std::f64::consts::PI / 2.0;
1557        let degrees = <f64 as AngleConversion>::to_degrees(&radians).expect("Operation failed");
1558        assert!((degrees - 90.0).abs() < 1e-10);
1559    }
1560
1561    #[test]
1562    fn test_scalar() {
1563        let a = Scalar::new(3.0_f64);
1564        let b = Scalar::new(4.0_f64);
1565
1566        assert_eq!((a + b).value(), 7.0);
1567        assert_eq!((a - b).value(), -1.0);
1568        assert_eq!((a * b).value(), 12.0);
1569        assert_eq!((a / b).value(), 0.75);
1570        assert_eq!((-a).value(), -3.0);
1571
1572        let c: Scalar<f64> = 5.0.into();
1573        assert_eq!(c.value(), 5.0);
1574    }
1575}
scirs2_core/numeric.rs

scirs2_core/
numeric.rs
