sva-formula 0.7.13

Laws in t and f: the normal form, the rule table, and the observations a law answers without sampling
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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

// Concern: the complex amplitude every atom factor is measured in | Non-concern: what an amplitude multiplies (spectral_sum/atom.rs) | IO: (C64, C64) -> C64

use std::ops::{Add, Div, Mul, Neg, Sub};

#[derive(Clone, Copy, Debug)]
pub struct C64 {
    pub re: f64,
    pub im: f64,
}

/// Bitwise after canonicalizing `-0.0`, so equality is transitive.
impl PartialEq for C64 {
    fn eq(&self, other: &C64) -> bool {
        self.bits() == other.bits()
    }
}
impl Eq for C64 {}

impl C64 {
    pub const ZERO: C64 = C64 { re: 0.0, im: 0.0 };
    pub const ONE: C64 = C64 { re: 1.0, im: 0.0 };
    pub const I: C64 = C64 { re: 0.0, im: 1.0 };

    pub fn new(re: f64, im: f64) -> C64 {
        C64 { re, im }
    }

    pub fn real(re: f64) -> C64 {
        C64 { re, im: 0.0 }
    }

    pub fn bits(self) -> (u64, u64) {
        (canonical(self.re), canonical(self.im))
    }

    pub fn is_finite(self) -> bool {
        self.re.is_finite() && self.im.is_finite()
    }

    pub fn is_zero(self) -> bool {
        self.bits() == C64::ZERO.bits()
    }

    pub fn is_real(self) -> bool {
        canonical(self.im) == canonical(0.0)
    }

    pub fn conj(self) -> C64 {
        C64::new(self.re, -self.im)
    }

    pub fn norm_sqr(self) -> f64 {
        self.re * self.re + self.im * self.im
    }

    pub fn abs(self) -> f64 {
        self.re.hypot(self.im)
    }

    /// A component alone on its axis inverts by one division, correctly rounded; only a
    /// value off both axes takes the conjugate over the squared norm.
    pub fn inv(self) -> C64 {
        let d = self.norm_sqr();
        match (self.re == 0.0, self.im == 0.0) {
            (false, true) => C64::new(1.0 / self.re, -self.im / d),
            (true, false) => C64::new(self.re / d, -1.0 / self.im),
            _ => C64::new(self.re / d, -self.im / d),
        }
    }

    pub fn exp(self) -> C64 {
        let r = self.re.exp();
        C64::new(r * self.im.cos(), r * self.im.sin())
    }

    pub fn powi(self, n: u32) -> C64 {
        let mut out = C64::ONE;
        for _ in 0..n {
            out = out * self;
        }
        out
    }

    pub fn scale(self, k: f64) -> C64 {
        C64::new(self.re * k, self.im * k)
    }

    pub fn over(self, k: f64) -> C64 {
        C64::new(self.re / k, self.im / k)
    }
}

pub fn canonical(x: f64) -> u64 {
    if x == 0.0 {
        0f64.to_bits()
    } else {
        x.to_bits()
    }
}

impl Add for C64 {
    type Output = C64;
    fn add(self, o: C64) -> C64 {
        C64::new(self.re + o.re, self.im + o.im)
    }
}
impl Sub for C64 {
    type Output = C64;
    fn sub(self, o: C64) -> C64 {
        C64::new(self.re - o.re, self.im - o.im)
    }
}
impl Mul for C64 {
    type Output = C64;
    fn mul(self, o: C64) -> C64 {
        C64::new(
            self.re * o.re - self.im * o.im,
            self.re * o.im + self.im * o.re,
        )
    }
}
impl Div for C64 {
    type Output = C64;
    fn div(self, o: C64) -> C64 {
        if o.im == 0.0 && o.re != 0.0 {
            return C64::new(self.re / o.re, self.im / o.re);
        }
        if o.re == 0.0 && o.im != 0.0 {
            return C64::new(self.im / o.im, -self.re / o.im);
        }
        let d = o.norm_sqr();
        C64::new(
            (self.re * o.re + self.im * o.im) / d,
            (self.im * o.re - self.re * o.im) / d,
        )
    }
}
impl Neg for C64 {
    type Output = C64;
    fn neg(self) -> C64 {
        C64::new(-self.re, -self.im)
    }
}