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

scirs2_core/
numeric.rs
