dashu_cmplx/
div.rs

1//! Complex division and reciprocal (near-correctly rounded via Smith's method + guard re-round).
2
3use crate::cbig::CBig;
4use crate::repr::{combine_parts, exact, riemann, CfpResult, Context};
5use core::ops::{Div, DivAssign};
6use dashu_base::{AbsOrd, Inverse};
7use dashu_float::round::Round;
8use dashu_float::{FBig, FpError};
9use dashu_int::Word;
10
11/// Guard digits (base-B) for `div`/`inv`. The naive complex-division error is `~(3+√5)·u`; a fixed
12/// guard comfortably absorbs it for well-conditioned denominators.
13const DIV_GUARD: usize = 14;
14
15impl<R: Round> Context<R> {
16    /// Reciprocal `1/z = conj(z)/|z|²` under this context (context layer).
17    pub fn inv<const B: Word>(&self, z: &CBig<R, B>) -> CfpResult<R, B> {
18        if z.is_infinite() {
19            return Ok(exact(FBig::ZERO, FBig::ZERO)); // 1/∞ = 0
20        }
21        if z.is_zero() {
22            return Ok(riemann(*self)); // 1/0 = ∞
23        }
24        let gctx = self.guard(DIV_GUARD);
25        let p = self.precision();
26        let (x, y) = (z.re(), z.im());
27        // n = x² + y²
28        let x2 = gctx.sqr(x)?.value();
29        let y2 = gctx.sqr(y)?.value();
30        let n = gctx.add(x2.repr(), y2.repr())?.value();
31        // 1/z = (x/n) + i(-y/n)
32        let re = gctx.div(x, n.repr())?.value().with_precision(p);
33        let neg_y = -y.clone();
34        let im = gctx.div(&neg_y, n.repr())?.value().with_precision(p);
35        Ok(combine_parts(re, im))
36    }
37
38    /// Divide two complex numbers under this context (context layer), using Smith's overflow-safe
39    /// method: the branch `|u| >= |v|` avoids forming `|denominator|²`.
40    pub fn div<const B: Word>(&self, z: &CBig<R, B>, w: &CBig<R, B>) -> CfpResult<R, B> {
41        if let Some(special) = div_special(z, w) {
42            return special;
43        }
44        let gctx = self.guard(DIV_GUARD);
45        let p = self.precision();
46        let (x, y) = (z.re(), z.im());
47        let (u, v) = (w.re(), w.im());
48        // Determine |u| >= |v| via abs_cmp on temporary FBig views (Smith's method).
49        let u_ge_v = {
50            let fu = FBig::from_repr(u.clone(), gctx);
51            let fv = FBig::from_repr(v.clone(), gctx);
52            fu.abs_cmp(&fv).is_ge()
53        };
54
55        // r, d depend on which of |u|, |v| is larger
56        let (r, d) = if u_ge_v {
57            // r = v/u, d = u + r·v
58            let r = gctx.div(v, u)?.value();
59            let rv = gctx.mul(r.repr(), v)?.value();
60            let d = gctx.add(u, rv.repr())?.value();
61            (r, d)
62        } else {
63            // r = u/v, d = v + r·u
64            let r = gctx.div(u, v)?.value();
65            let ru = gctx.mul(r.repr(), u)?.value();
66            let d = gctx.add(v, ru.repr())?.value();
67            (r, d)
68        };
69
70        let (re, im) = if u_ge_v {
71            // re = (x + r·y)/d, im = (y - r·x)/d
72            let ry = gctx.mul(r.repr(), y)?.value();
73            let rx = gctx.mul(r.repr(), x)?.value();
74            let num_re = gctx.add(x, ry.repr())?.value();
75            let num_im = gctx.sub(y, rx.repr())?.value();
76            (
77                gctx.div(num_re.repr(), d.repr())?.value().with_precision(p),
78                gctx.div(num_im.repr(), d.repr())?.value().with_precision(p),
79            )
80        } else {
81            // re = (r·x + y)/d, im = (r·y - x)/d
82            let rx = gctx.mul(r.repr(), x)?.value();
83            let ry = gctx.mul(r.repr(), y)?.value();
84            let num_re = gctx.add(rx.repr(), y)?.value();
85            let num_im = gctx.sub(ry.repr(), x)?.value();
86            (
87                gctx.div(num_re.repr(), d.repr())?.value().with_precision(p),
88                gctx.div(num_im.repr(), d.repr())?.value().with_precision(p),
89            )
90        };
91        Ok(combine_parts(re, im))
92    }
93
94    /// Divide a complex number by a real scalar (context layer): `(x+iy)/s = (x/s) + i(y/s)`.
95    pub fn div_real<const B: Word>(&self, z: &CBig<R, B>, s: &FBig<R, B>) -> CfpResult<R, B> {
96        if z.is_infinite() || s.repr().is_infinite() {
97            if z.is_infinite() && s.repr().is_infinite() {
98                return Err(FpError::Indeterminate); // ∞/∞
99            }
100            if s.repr().is_infinite() {
101                return Ok(exact(FBig::ZERO, FBig::ZERO)); // finite/∞ = 0
102            }
103            // z infinite, s finite nonzero → ∞
104            return Ok(riemann(*self));
105        }
106        if s.repr().is_pos_zero() || s.repr().is_neg_zero() {
107            if z.is_zero() {
108                return Err(FpError::Indeterminate); // 0/0 (incl. ±0)
109            }
110            return Ok(riemann(*self)); // z/±0 (z≠0) = ∞
111        }
112        let gctx = self.guard(DIV_GUARD);
113        let p = self.precision();
114        let re = gctx.div(z.re(), s.repr())?.value().with_precision(p);
115        let im = gctx.div(z.im(), s.repr())?.value().with_precision(p);
116        Ok(combine_parts(re, im))
117    }
118}
119
120/// Annex-G short-circuit for `z / w`.
121fn div_special<R: Round, const B: Word>(z: &CBig<R, B>, w: &CBig<R, B>) -> Option<CfpResult<R, B>> {
122    let (zi, wi) = (z.is_infinite(), w.is_infinite());
123    let (zz, wz) = (z.is_zero(), w.is_zero());
124    let ctx = Context::max(z.context(), w.context());
125    if (zi && wi) || (zz && wz) {
126        Some(Err(FpError::Indeterminate)) // ∞/∞ or 0/0
127    } else if wi {
128        Some(Ok(exact(FBig::ZERO, FBig::ZERO))) // (finite or 0) / ∞ = 0
129    } else if wz || zi {
130        Some(Ok(riemann(ctx))) // (nonzero or ∞) / 0, or ∞ / finite = ∞
131    } else if zz {
132        Some(Ok(exact(FBig::ZERO, FBig::ZERO))) // 0 / finite = 0
133    } else {
134        None
135    }
136}
137
138impl<R: Round, const B: Word> Inverse for CBig<R, B> {
139    type Output = CBig<R, B>;
140
141    #[inline]
142    fn inv(self) -> Self::Output {
143        self.context().unwrap_cfp(self.context().inv(&self))
144    }
145}
146
147impl<R: Round, const B: Word> Inverse for &CBig<R, B> {
148    type Output = CBig<R, B>;
149
150    #[inline]
151    fn inv(self) -> Self::Output {
152        self.context().unwrap_cfp(self.context().inv(self))
153    }
154}
155
156// CBig / CBig operators
157crate::helper_macros::impl_cbig_binop!(Div, div, DivAssign, div_assign);
158
159// --- scalar division by a real FBig (mixed-type operators) ---
160
161// CBig / FBig (componentwise, via the shared scalar macro).
162crate::helper_macros::impl_cbig_scalar_binop!(Div, div, div_real);
163
164// FBig / CBig = (s + 0i) / z, reusing complex division.
165impl<R: Round, const B: Word> Div<&CBig<R, B>> for &FBig<R, B> {
166    type Output = CBig<R, B>;
167    #[inline]
168    fn div(self, rhs: &CBig<R, B>) -> CBig<R, B> {
169        let s = CBig::from(self.clone());
170        let ctx = Context::max(s.context(), rhs.context());
171        ctx.unwrap_cfp(ctx.div(&s, rhs))
172    }
173}
174impl<R: Round, const B: Word> Div<CBig<R, B>> for &FBig<R, B> {
175    type Output = CBig<R, B>;
176    #[inline]
177    fn div(self, rhs: CBig<R, B>) -> CBig<R, B> {
178        let this = self.clone();
179        let s = CBig::from(this);
180        let ctx = Context::max(s.context(), rhs.context());
181        ctx.unwrap_cfp(ctx.div(&s, &rhs))
182    }
183}
184impl<R: Round, const B: Word> Div<&CBig<R, B>> for FBig<R, B> {
185    type Output = CBig<R, B>;
186    #[inline]
187    fn div(self, rhs: &CBig<R, B>) -> CBig<R, B> {
188        let s = CBig::from(self);
189        let ctx = Context::max(s.context(), rhs.context());
190        ctx.unwrap_cfp(ctx.div(&s, rhs))
191    }
192}
193impl<R: Round, const B: Word> Div<CBig<R, B>> for FBig<R, B> {
194    type Output = CBig<R, B>;
195    #[inline]
196    fn div(self, rhs: CBig<R, B>) -> CBig<R, B> {
197        let s = CBig::from(self);
198        let ctx = Context::max(s.context(), rhs.context());
199        ctx.unwrap_cfp(ctx.div(&s, &rhs))
200    }
201}
202
203#[cfg(test)]
204mod tests {
205    use super::*;
206    use dashu_float::round::mode;
207
208    type C = CBig<mode::HalfAway, 10>;
209    type F = FBig<mode::HalfAway, 10>;
210
211    fn c(re: i32, im: i32) -> C {
212        let mk = |v: i32| -> F { F::from(v).with_precision(53).value() };
213        C::from_parts(mk(re), mk(im))
214    }
215
216    #[test]
217    fn div_inverse() {
218        // z / z = 1 for z != 0
219        let z = c(3, 4);
220        let q = &z / &z;
221        assert_eq!(q.re().significand(), &1.into());
222        assert!(q.im().significand().is_zero());
223    }
224
225    #[test]
226    fn div_basic() {
227        // (6+8i)/(3+4i) = 2  (since 6+8i = 2·(3+4i))
228        let z = c(6, 8);
229        let w = c(3, 4);
230        let q = &z / &w;
231        assert_eq!(q.re().significand(), &2.into());
232        assert!(q.im().significand().is_zero());
233    }
234
235    #[test]
236    fn inv_basic() {
237        // 1/(3+4i) = (3-4i)/25 = 0.12 - 0.16i — use precision 53 for exact-ish
238        type F = FBig<mode::HalfAway, 10>;
239        let mk = |v: i32| -> F { F::from(v).with_precision(53).value() };
240        let z = C::from_parts(mk(3), mk(4));
241        let r = (&z).inv();
242        // (3-4i)/25: re = 3/25, im = -4/25
243        assert_eq!(r.context().precision(), 53);
244        // re ≈ 0.12, im ≈ -0.16; check via multiplying back: z·inv(z) = 1
245        let one = &z * &r;
246        assert_eq!(one.re().significand(), &1.into());
247        assert!(one.im().significand().is_zero());
248    }
249
250    #[test]
251    fn scalar_div_by_real() {
252        let z = c(6, 8);
253        let s = FBig::<mode::HalfAway, 10>::from(2);
254        let q = &z / &s;
255        assert_eq!(q.re().significand(), &3.into());
256        assert_eq!(q.im().significand(), &4.into());
257    }
258
259    #[test]
260    fn div_real_by_negative_zero_matches_positive_zero() {
261        use dashu_float::Repr;
262        // Dividing a nonzero complex by the scalar -0 must give the same Riemann infinity as +0,
263        // not fall through to the general path (which would divide by -0 componentwise).
264        let z = c(3, 4);
265        let neg_zero = F::from_repr(Repr::neg_zero(), z.context().float());
266        let q_neg = &z / &neg_zero;
267        let q_pos = &z / &F::ZERO;
268        assert!(q_neg.is_infinite());
269        assert!(q_neg == q_pos); // both are +∞ + i·0
270    }
271}
dashu_cmplx/div.rs

dashu_cmplx/
div.rs
