rustymath_core/linalg/
scalar.rs

1//! The item types of the matrices that `factor` decomposes: `f32`, `f64`, and complex numbers with parts of either type. Each decomposition is written once for all of them. Where a real algorithm transposes, a complex one also conjugates. Real quantities, such as norms, the sizes of pivots and singular values, have the parts' type.
2use crate::{Native, float::{FUSED, sealed::Sealed}};
3use num_complex::Complex;
4use num_traits::{One, Zero};
5use std::ops::{Add, Mul, Neg, Sub};
6
7/// A real or complex number that the decompositions take: `f32`, `f64`, or a complex number with parts of either type. rustymath implements it for those four types, and no other crate can. For a real type each method is the plain operation, and the generic algorithms give the same real results as algorithms written for reals alone. The trait has no `/`. Quotients go through `quotient` and `reciprocal`, which keep a complex quotient from overflowing where its result fits.
8pub trait Scalar:
9    Sealed + Copy + PartialEq + Zero + One + Add<Output = Self> + Sub<Output = Self> + Mul<Output = Self> + Neg<Output = Self> + Send + Sync
10{
11    /// The type of the real and imaginary parts.
12    type Real: Native;
13    /// Whether the type is complex, with a real and an imaginary part.
14    const COMPLEX: bool;
15    fn conj(self) -> Self;
16    fn re(self) -> Self::Real;
17    fn im(self) -> Self::Real;
18    fn from_real(x: Self::Real) -> Self;
19    /// `|re| + |im|`, within a factor of √2 of the modulus and cheaper. Pivot searches and scaling compare it.
20    fn abs1(self) -> Self::Real;
21    /// `|self|`, computed without overflow where the result fits, as `hypot` computes it.
22    fn modulus(self) -> Self::Real;
23    fn mul_real(self, x: Self::Real) -> Self;
24    fn div_real(self, x: Self::Real) -> Self;
25    /// `self × b + c`, each multiply-add rounding once where `FUSED` is true, as `Core::mla` does.
26    fn mul_acc(self, b: Self, c: Self) -> Self;
27    fn reciprocal(self) -> Self;
28    /// `self / d`. A complex quotient halves operands near overflow and scales up operands near underflow, by powers of two, as Baudin and Smith's form of Smith's algorithm does. It then applies Smith's formula.
29    fn quotient(self, d: Self) -> Self;
30}
31
32impl<T: Native> Scalar for T {
33    type Real = T;
34    const COMPLEX: bool = false;
35    #[inline(always)]
36    fn conj(self) -> T { self }
37    #[inline(always)]
38    fn re(self) -> T { self }
39    #[inline(always)]
40    fn im(self) -> T { T::zero() }
41    #[inline(always)]
42    fn from_real(x: T) -> T { x }
43    #[inline(always)]
44    fn abs1(self) -> T { self.abs() }
45    #[inline(always)]
46    fn modulus(self) -> T { self.abs() }
47    #[inline(always)]
48    fn mul_real(self, x: T) -> T { self * x }
49    #[inline(always)]
50    fn div_real(self, x: T) -> T { self / x }
51    #[inline(always)]
52    fn mul_acc(self, b: T, c: T) -> T { self.mla::<FUSED>(b, c) }
53    #[inline(always)]
54    fn reciprocal(self) -> T { self.recip() }
55    #[inline(always)]
56    fn quotient(self, d: T) -> T { self / d }
57}
58
59impl<T: Native> Scalar for Complex<T> {
60    type Real = T;
61    const COMPLEX: bool = true;
62    #[inline(always)]
63    fn conj(self) -> Self { Complex::new(self.re, -self.im) }
64    #[inline(always)]
65    fn re(self) -> T { self.re }
66    #[inline(always)]
67    fn im(self) -> T { self.im }
68    #[inline(always)]
69    fn from_real(x: T) -> Self { Complex::new(x, T::zero()) }
70    #[inline(always)]
71    fn abs1(self) -> T { self.re.abs() + self.im.abs() }
72    #[inline(always)]
73    fn modulus(self) -> T { crate::complex::modulus::<T, FUSED>(self.re, self.im) }
74    #[inline(always)]
75    fn mul_real(self, x: T) -> Self { Complex::new(self.re * x, self.im * x) }
76    #[inline(always)]
77    fn div_real(self, x: T) -> Self { Complex::new(self.re / x, self.im / x) }
78    #[inline(always)]
79    fn mul_acc(self, b: Self, c: Self) -> Self {
80        let re = self.re.mla::<FUSED>(b.re, (-self.im).mla::<FUSED>(b.im, c.re));
81        Complex::new(re, self.re.mla::<FUSED>(b.im, self.im.mla::<FUSED>(b.re, c.im)))
82    }
83    #[inline(always)]
84    fn reciprocal(self) -> Self { Self::one().quotient(self) }
85    #[inline(always)]
86    fn quotient(self, d: Self) -> Self {
87        let (two, eps) = (T::one() + T::one(), T::epsilon());
88        let (big, small, up) = (T::max_value() / two, T::min_positive_value() * two / eps, two / (eps * eps));
89        let size = |z: Self| z.re.abs().max(z.im.abs());
90        let (xs, ys) = (size(self), size(d));
91        let (mut x, mut y, mut scale) = (self, d, T::one());
92        if xs >= big { (x, scale) = (x.div_real(two), scale * two) }
93        if ys >= big { (y, scale) = (y.div_real(two), scale / two) }
94        if xs <= small { (x, scale) = (x.mul_real(up), scale / up) }
95        if ys <= small { (y, scale) = (y.mul_real(up), scale * up) }
96        // (a + bi) / (c + di) for |d| ≤ |c|, as its real and imaginary parts. Where `d / c` underflows to 0, the parts use `b / c` and `a / c` instead.
97        let smith = |a: T, b: T, c: T, d: T| {
98            let r = d / c;
99            let s = c + d * r;
100            if r != T::zero() { ((a + b * r) / s, (b - a * r) / s) } else { ((a + d * (b / c)) / s, (b - d * (a / c)) / s) }
101        };
102        let (re, im) = if y.re.abs() >= y.im.abs() { smith(x.re, x.im, y.re, y.im) } else { let (re, im) = smith(x.im, x.re, y.im, y.re); (re, -im) };
103        Complex::new(re, im).mul_real(scale)
104    }
105}
rustymath_core/linalg/scalar.rs

rustymath_core/linalg/
scalar.rs
