rustymath-core 0.1.56

Numerical algorithms, formulas and compact array storage.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105

//! 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.
use crate::{Native, float::{FUSED, sealed::Sealed}};
use num_complex::Complex;
use num_traits::{One, Zero};
use std::ops::{Add, Mul, Neg, Sub};

/// 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.
pub trait Scalar:
    Sealed + Copy + PartialEq + Zero + One + Add<Output = Self> + Sub<Output = Self> + Mul<Output = Self> + Neg<Output = Self> + Send + Sync
{
    /// The type of the real and imaginary parts.
    type Real: Native;
    /// Whether the type is complex, with a real and an imaginary part.
    const COMPLEX: bool;
    fn conj(self) -> Self;
    fn re(self) -> Self::Real;
    fn im(self) -> Self::Real;
    fn from_real(x: Self::Real) -> Self;
    /// `|re| + |im|`, within a factor of √2 of the modulus and cheaper. Pivot searches and scaling compare it.
    fn abs1(self) -> Self::Real;
    /// `|self|`, computed without overflow where the result fits, as `hypot` computes it.
    fn modulus(self) -> Self::Real;
    fn mul_real(self, x: Self::Real) -> Self;
    fn div_real(self, x: Self::Real) -> Self;
    /// `self × b + c`, each multiply-add rounding once where `FUSED` is true, as `Core::mla` does.
    fn mul_acc(self, b: Self, c: Self) -> Self;
    fn reciprocal(self) -> Self;
    /// `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.
    fn quotient(self, d: Self) -> Self;
}

impl<T: Native> Scalar for T {
    type Real = T;
    const COMPLEX: bool = false;
    #[inline(always)]
    fn conj(self) -> T { self }
    #[inline(always)]
    fn re(self) -> T { self }
    #[inline(always)]
    fn im(self) -> T { T::zero() }
    #[inline(always)]
    fn from_real(x: T) -> T { x }
    #[inline(always)]
    fn abs1(self) -> T { self.abs() }
    #[inline(always)]
    fn modulus(self) -> T { self.abs() }
    #[inline(always)]
    fn mul_real(self, x: T) -> T { self * x }
    #[inline(always)]
    fn div_real(self, x: T) -> T { self / x }
    #[inline(always)]
    fn mul_acc(self, b: T, c: T) -> T { self.mla::<FUSED>(b, c) }
    #[inline(always)]
    fn reciprocal(self) -> T { self.recip() }
    #[inline(always)]
    fn quotient(self, d: T) -> T { self / d }
}

impl<T: Native> Scalar for Complex<T> {
    type Real = T;
    const COMPLEX: bool = true;
    #[inline(always)]
    fn conj(self) -> Self { Complex::new(self.re, -self.im) }
    #[inline(always)]
    fn re(self) -> T { self.re }
    #[inline(always)]
    fn im(self) -> T { self.im }
    #[inline(always)]
    fn from_real(x: T) -> Self { Complex::new(x, T::zero()) }
    #[inline(always)]
    fn abs1(self) -> T { self.re.abs() + self.im.abs() }
    #[inline(always)]
    fn modulus(self) -> T { crate::complex::modulus::<T, FUSED>(self.re, self.im) }
    #[inline(always)]
    fn mul_real(self, x: T) -> Self { Complex::new(self.re * x, self.im * x) }
    #[inline(always)]
    fn div_real(self, x: T) -> Self { Complex::new(self.re / x, self.im / x) }
    #[inline(always)]
    fn mul_acc(self, b: Self, c: Self) -> Self {
        let re = self.re.mla::<FUSED>(b.re, (-self.im).mla::<FUSED>(b.im, c.re));
        Complex::new(re, self.re.mla::<FUSED>(b.im, self.im.mla::<FUSED>(b.re, c.im)))
    }
    #[inline(always)]
    fn reciprocal(self) -> Self { Self::one().quotient(self) }
    #[inline(always)]
    fn quotient(self, d: Self) -> Self {
        let (two, eps) = (T::one() + T::one(), T::epsilon());
        let (big, small, up) = (T::max_value() / two, T::min_positive_value() * two / eps, two / (eps * eps));
        let size = |z: Self| z.re.abs().max(z.im.abs());
        let (xs, ys) = (size(self), size(d));
        let (mut x, mut y, mut scale) = (self, d, T::one());
        if xs >= big { (x, scale) = (x.div_real(two), scale * two) }
        if ys >= big { (y, scale) = (y.div_real(two), scale / two) }
        if xs <= small { (x, scale) = (x.mul_real(up), scale / up) }
        if ys <= small { (y, scale) = (y.mul_real(up), scale * up) }
        // (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.
        let smith = |a: T, b: T, c: T, d: T| {
            let r = d / c;
            let s = c + d * r;
            if r != T::zero() { ((a + b * r) / s, (b - a * r) / s) } else { ((a + d * (b / c)) / s, (b - d * (a / c)) / s) }
        };
        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) };
        Complex::new(re, im).mul_real(scale)
    }
}