cinrs_core/complex.rs
1//! Complex arithmetic, for folding a constant expression.
2//!
3//! `static double _Complex z = (1.0 + 2.0i) * (3.0 - 4.0i);` has to become a
4//! literal: a Rust `static` initialiser is a constant expression, and
5//! [`crate::sema`] reduces every arithmetic one to a value before code
6//! generation ever sees it. So the front end needs the same arithmetic the
7//! *runtime* has, and this is it.
8//!
9//! # Why it is written twice
10//!
11//! The generated code calls `cinrs_rt::complex`, a `no_std` crate built around
12//! `num_complex::Complex`; this crate is the procedural macro's own front end
13//! and has no business depending on the runtime — a build that said
14//! `default-features = false` would then still compile `num-complex` for a
15//! feature it had switched off. The two copies are kept honest by a test
16//! rather than by a shared crate: `tests/complex.rs` folds a table of
17//! constants through *this* code and computes the same products and quotients
18//! through the runtime at run time, and fails if any bit differs.
19//!
20//! See `cinrs_rt::complex` for what the algorithms are and why: the naive
21//! product with Annex G.5.1's infinity recovery, Smith's algorithm for the
22//! `double` quotient, and the closed form in the wider format for the `float`
23//! one.
24
25/// The two parts of a complex value, always carried as `f64`.
26///
27/// A `float _Complex` constant is carried the same way and rounded to `f32`
28/// where it is stored, exactly as [`crate::ir::ConstValue::Float`] carries a
29/// `float`.
30pub type Parts = (f64, f64);
31
32/// Generates the product and the `f64` quotient for one component width.
33macro_rules! narrow_ops {
34 ($t:ty, $mul:ident) => {
35 /// The naive product with Annex G.5.1's recovery, in this width.
36 fn $mul(mut a: $t, mut b: $t, mut c: $t, mut d: $t) -> ($t, $t) {
37 let (ac, bd, ad, bc) = (a * c, b * d, a * d, b * c);
38 let mut x = ac - bd;
39 let mut y = ad + bc;
40 if x.is_nan() && y.is_nan() {
41 let mut recalc = false;
42 let unit = |v: $t| if v.is_infinite() { 1.0 } else { 0.0 as $t }.copysign(v);
43 let tame = |v: $t| {
44 if v.is_nan() {
45 (0.0 as $t).copysign(v)
46 } else {
47 v
48 }
49 };
50 if a.is_infinite() || b.is_infinite() {
51 a = unit(a);
52 b = unit(b);
53 c = tame(c);
54 d = tame(d);
55 recalc = true;
56 }
57 if c.is_infinite() || d.is_infinite() {
58 c = unit(c);
59 d = unit(d);
60 a = tame(a);
61 b = tame(b);
62 recalc = true;
63 }
64 if !recalc
65 && (ac.is_infinite()
66 || bd.is_infinite()
67 || ad.is_infinite()
68 || bc.is_infinite())
69 {
70 a = tame(a);
71 b = tame(b);
72 c = tame(c);
73 d = tame(d);
74 recalc = true;
75 }
76 if recalc {
77 let inf = <$t>::INFINITY;
78 x = inf * (a * c - b * d);
79 y = inf * (a * d + b * c);
80 }
81 }
82 (x, y)
83 }
84 };
85}
86
87narrow_ops!(f32, mul_parts_f32);
88narrow_ops!(f64, mul_parts_f64);
89
90/// The product of two `double _Complex` constants.
91pub fn mul((a, b): Parts, (c, d): Parts) -> Parts {
92 mul_parts_f64(a, b, c, d)
93}
94
95/// The product of two `float _Complex` constants, computed in single
96/// precision.
97///
98/// The width matters: an intermediate that overflows to infinity in `float`
99/// makes `(10²⁰ + 10²⁰i)²` a NaN there and a finite zero in `double`, and the
100/// runtime computes it in `float`.
101pub fn mul_f32((a, b): Parts, (c, d): Parts) -> Parts {
102 let (x, y) = mul_parts_f32(a as f32, b as f32, c as f32, d as f32);
103 (f64::from(x), f64::from(y))
104}
105
106/// The quotient of two `double _Complex` constants: Smith's algorithm with
107/// Annex G.5.1's recovery.
108pub fn div((a, b): Parts, (c, d): Parts) -> Parts {
109 let (x, y) = if abs(c) < abs(d) {
110 let ratio = c / d;
111 let denom = c * ratio + d;
112 if abs(ratio) > f64::MIN_POSITIVE {
113 ((a * ratio + b) / denom, (b * ratio - a) / denom)
114 } else {
115 (((a / d) * c + b) / denom, ((b / d) * c - a) / denom)
116 }
117 } else {
118 let ratio = d / c;
119 let denom = d * ratio + c;
120 if abs(ratio) > f64::MIN_POSITIVE {
121 ((b * ratio + a) / denom, (b - a * ratio) / denom)
122 } else {
123 (((b / c) * d + a) / denom, (b - (a / c) * d) / denom)
124 }
125 };
126 recover((a, b), (c, d), (x, y))
127}
128
129/// The quotient of two `float _Complex` constants: the closed form in
130/// `double`, which cannot overflow for `float` operands, with the same
131/// recovery and a single rounding at the end.
132pub fn div_f32((a, b): Parts, (c, d): Parts) -> Parts {
133 let (a, b) = (f64::from(a as f32), f64::from(b as f32));
134 let (c, d) = (f64::from(c as f32), f64::from(d as f32));
135 let denom = c * c + d * d;
136 let quotient = ((a * c + b * d) / denom, (b * c - a * d) / denom);
137 let (x, y) = recover((a, b), (c, d), quotient);
138 (f64::from(x as f32), f64::from(y as f32))
139}
140
141/// Annex G.5.1's recovery for a quotient that came out NaN + iNaN.
142fn recover((a, b): Parts, (c, d): Parts, (x, y): Parts) -> Parts {
143 if !(x.is_nan() && y.is_nan()) {
144 return (x, y);
145 }
146 let inf = f64::INFINITY;
147 let unit = |v: f64| if v.is_infinite() { 1.0 } else { 0.0f64 }.copysign(v);
148 if c == 0.0 && d == 0.0 && (!a.is_nan() || !b.is_nan()) {
149 let scale = inf.copysign(c);
150 return (scale * a, scale * b);
151 }
152 if (a.is_infinite() || b.is_infinite()) && c.is_finite() && d.is_finite() {
153 let (a, b) = (unit(a), unit(b));
154 return (inf * (a * c + b * d), inf * (b * c - a * d));
155 }
156 if (c.is_infinite() || d.is_infinite()) && a.is_finite() && b.is_finite() {
157 let (c, d) = (unit(c), unit(d));
158 return (0.0 * (a * c + b * d), 0.0 * (b * c - a * d));
159 }
160 (x, y)
161}
162
163fn abs(x: f64) -> f64 {
164 f64::from_bits(x.to_bits() & !(1u64 << 63))
165}