feanor_math/rings/
float_complex.rs

1use std::f64::EPSILON;
2use std::f64::consts::PI;
3
4use crate::divisibility::{DivisibilityRing, Domain};
5use crate::field::Field;
6use crate::homomorphism::CanHomFrom;
7use crate::impl_eq_based_self_iso;
8use crate::integer::{IntegerRing, IntegerRingStore};
9use crate::pid::{EuclideanRing, PrincipalIdealRing};
10use crate::ring::*;
11use crate::rings::approx_real::float::Real64;
12use crate::rings::rational::RationalFieldBase;
13
14/// An approximate implementation of the complex numbers `C`, using 64 bit floating
15/// point numbers.
16///
17/// # Warning
18///
19/// Since floating point numbers do not exactly represent the complex numbers, and this crate
20/// follows a mathematically precise approach, we cannot provide any function related to equality.
21/// In particular, `Complex64Base.eq_el(a, b)` is not supported, and will panic.
22/// Hence, this ring has only limited use within this crate, and is currently only used for
23/// floating-point FFTs.
24#[derive(Clone, Copy, PartialEq, Debug)]
25pub struct Complex64Base;
26
27/// An element of [`Complex64`].
28#[derive(Clone, Copy, Debug)]
29pub struct Complex64El(f64, f64);
30
31/// [`RingStore`] corresponding to [`Complex64Base`]
32pub type Complex64 = RingValue<Complex64Base>;
33
34impl Complex64 {
35    pub const RING: Self = RingValue::from(Complex64Base);
36    pub const I: Complex64El = Complex64El(0.0, 1.0);
37}
38
39impl Complex64Base {
40    pub fn abs(&self, Complex64El(re, im): Complex64El) -> f64 { (re * re + im * im).sqrt() }
41
42    pub fn conjugate(&self, Complex64El(re, im): Complex64El) -> Complex64El { Complex64El(re, -im) }
43
44    pub fn exp(&self, Complex64El(exp_re, exp_im): Complex64El) -> Complex64El {
45        let angle = exp_im;
46        let abs = exp_re.exp();
47        Complex64El(abs * angle.cos(), abs * angle.sin())
48    }
49
50    pub fn closest_gaussian_int(&self, Complex64El(re, im): Complex64El) -> (i64, i64) {
51        (re.round() as i64, im.round() as i64)
52    }
53
54    pub fn ln_main_branch(&self, Complex64El(re, im): Complex64El) -> Complex64El {
55        Complex64El(self.abs(Complex64El(re, im)).ln(), im.atan2(re))
56    }
57
58    pub fn is_absolute_approx_eq(&self, lhs: Complex64El, rhs: Complex64El, absolute_threshold: f64) -> bool {
59        self.abs(self.sub(lhs, rhs)) < absolute_threshold
60    }
61
62    pub fn is_relative_approx_eq(&self, lhs: Complex64El, rhs: Complex64El, relative_limit: f64) -> bool {
63        self.is_absolute_approx_eq(lhs, rhs, self.abs(lhs) * relative_limit)
64    }
65
66    pub fn is_approx_eq(&self, lhs: Complex64El, rhs: Complex64El, precision: u64) -> bool {
67        let scaled_precision = precision as f64 * EPSILON;
68        if self.is_absolute_approx_eq(lhs, self.zero(), scaled_precision) {
69            self.is_absolute_approx_eq(rhs, self.zero(), scaled_precision)
70        } else {
71            self.is_relative_approx_eq(lhs, rhs, scaled_precision)
72        }
73    }
74
75    pub fn from_f64(&self, x: f64) -> Complex64El { Complex64El(x, 0.0) }
76
77    pub fn root_of_unity(&self, i: i64, n: i64) -> Complex64El {
78        self.exp(self.mul(self.from_f64((i as f64 / n as f64) * (2.0 * PI)), Complex64::I))
79    }
80
81    pub fn re(&self, Complex64El(re, _im): Complex64El) -> f64 { re }
82
83    pub fn im(&self, Complex64El(_re, im): Complex64El) -> f64 { im }
84}
85
86impl Complex64 {
87    pub fn abs(&self, val: Complex64El) -> f64 { self.get_ring().abs(val) }
88
89    pub fn conjugate(&self, val: Complex64El) -> Complex64El { self.get_ring().conjugate(val) }
90
91    pub fn exp(&self, exp: Complex64El) -> Complex64El { self.get_ring().exp(exp) }
92
93    pub fn closest_gaussian_int(&self, val: Complex64El) -> (i64, i64) { self.get_ring().closest_gaussian_int(val) }
94
95    pub fn ln_main_branch(&self, val: Complex64El) -> Complex64El { self.get_ring().ln_main_branch(val) }
96
97    pub fn is_absolute_approx_eq(&self, lhs: Complex64El, rhs: Complex64El, absolute_threshold: f64) -> bool {
98        self.get_ring().is_absolute_approx_eq(lhs, rhs, absolute_threshold)
99    }
100
101    pub fn is_relative_approx_eq(&self, lhs: Complex64El, rhs: Complex64El, relative_limit: f64) -> bool {
102        self.get_ring().is_relative_approx_eq(lhs, rhs, relative_limit)
103    }
104
105    pub fn is_approx_eq(&self, lhs: Complex64El, rhs: Complex64El, precision: u64) -> bool {
106        self.get_ring().is_approx_eq(lhs, rhs, precision)
107    }
108
109    pub fn from_f64(&self, x: f64) -> Complex64El { self.get_ring().from_f64(x) }
110
111    pub fn root_of_unity(&self, i: i64, n: i64) -> Complex64El { self.get_ring().root_of_unity(i, n) }
112
113    pub fn re(&self, x: Complex64El) -> f64 { self.get_ring().re(x) }
114
115    pub fn im(&self, x: Complex64El) -> f64 { self.get_ring().im(x) }
116}
117
118impl RingBase for Complex64Base {
119    type Element = Complex64El;
120
121    fn clone_el(&self, val: &Self::Element) -> Self::Element { *val }
122
123    fn add_assign(&self, Complex64El(lhs_re, lhs_im): &mut Self::Element, Complex64El(rhs_re, rhs_im): Self::Element) {
124        *lhs_re += rhs_re;
125        *lhs_im += rhs_im;
126    }
127
128    fn negate_inplace(&self, Complex64El(re, im): &mut Self::Element) {
129        *re = -*re;
130        *im = -*im;
131    }
132
133    fn mul_assign(&self, Complex64El(lhs_re, lhs_im): &mut Self::Element, Complex64El(rhs_re, rhs_im): Self::Element) {
134        let new_im = *lhs_re * rhs_im + *lhs_im * rhs_re;
135        *lhs_re = *lhs_re * rhs_re - *lhs_im * rhs_im;
136        *lhs_im = new_im;
137    }
138
139    fn from_int(&self, value: i32) -> Self::Element { Complex64El(value as f64, 0.0) }
140
141    fn eq_el(&self, _: &Self::Element, _: &Self::Element) -> bool {
142        panic!("Cannot provide equality on approximate rings")
143    }
144
145    fn pow_gen<R: IntegerRingStore>(&self, x: Self::Element, power: &El<R>, integers: R) -> Self::Element
146    where
147        R::Type: IntegerRing,
148    {
149        self.exp(self.mul(
150            self.ln_main_branch(x),
151            Complex64El(integers.to_float_approx(power), 0.0),
152        ))
153    }
154
155    fn is_commutative(&self) -> bool { true }
156
157    fn is_noetherian(&self) -> bool { true }
158
159    fn is_approximate(&self) -> bool { true }
160
161    fn dbg_within<'a>(
162        &self,
163        Complex64El(re, im): &Self::Element,
164        out: &mut std::fmt::Formatter<'a>,
165        env: EnvBindingStrength,
166    ) -> std::fmt::Result {
167        if env >= EnvBindingStrength::Product {
168            write!(out, "({} + {}i)", re, im)
169        } else {
170            write!(out, "{} + {}i", re, im)
171        }
172    }
173
174    fn dbg<'a>(&self, value: &Self::Element, out: &mut std::fmt::Formatter<'a>) -> std::fmt::Result {
175        self.dbg_within(value, out, EnvBindingStrength::Weakest)
176    }
177
178    fn characteristic<I: IntegerRingStore + Copy>(&self, ZZ: I) -> Option<El<I>>
179    where
180        I::Type: IntegerRing,
181    {
182        Some(ZZ.zero())
183    }
184}
185
186impl_eq_based_self_iso! { Complex64Base }
187
188impl Domain for Complex64Base {}
189
190impl DivisibilityRing for Complex64Base {
191    fn checked_left_div(&self, lhs: &Self::Element, rhs: &Self::Element) -> Option<Self::Element> {
192        let abs_sqr = self.abs(*rhs) * self.abs(*rhs);
193        let Complex64El(res_re, res_im) = self.mul(*lhs, self.conjugate(*rhs));
194        return Some(Complex64El(res_re / abs_sqr, res_im / abs_sqr));
195    }
196}
197
198impl PrincipalIdealRing for Complex64Base {
199    fn checked_div_min(&self, lhs: &Self::Element, rhs: &Self::Element) -> Option<Self::Element> {
200        self.checked_left_div(lhs, rhs)
201    }
202
203    fn extended_ideal_gen(
204        &self,
205        _lhs: &Self::Element,
206        _rhs: &Self::Element,
207    ) -> (Self::Element, Self::Element, Self::Element) {
208        panic!("Since Complex64 is only approximate, this cannot be implemented properly")
209    }
210}
211
212impl EuclideanRing for Complex64Base {
213    fn euclidean_div_rem(&self, _lhs: Self::Element, _rhs: &Self::Element) -> (Self::Element, Self::Element) {
214        panic!("Since Complex64 is only approximate, this cannot be implemented properly")
215    }
216
217    fn euclidean_deg(&self, _: &Self::Element) -> Option<usize> {
218        panic!("Since Complex64 is only approximate, this cannot be implemented properly")
219    }
220}
221
222impl Field for Complex64Base {
223    fn div(&self, lhs: &Self::Element, rhs: &Self::Element) -> Self::Element {
224        self.checked_left_div(lhs, rhs).unwrap()
225    }
226}
227
228impl RingExtension for Complex64Base {
229    type BaseRing = Real64;
230
231    fn base_ring(&self) -> &Self::BaseRing { &Real64::RING }
232
233    fn from(&self, x: El<Self::BaseRing>) -> Self::Element { self.from_f64(x) }
234
235    fn mul_assign_base(&self, lhs: &mut Self::Element, rhs: &El<Self::BaseRing>) {
236        lhs.0 *= *rhs;
237        lhs.1 *= *rhs;
238    }
239}
240
241impl<I: ?Sized + IntegerRing> CanHomFrom<I> for Complex64Base {
242    type Homomorphism = ();
243
244    fn has_canonical_hom(&self, _from: &I) -> Option<Self::Homomorphism> { Some(()) }
245
246    fn map_in(&self, from: &I, el: <I as RingBase>::Element, hom: &Self::Homomorphism) -> Self::Element {
247        self.map_in_ref(from, &el, hom)
248    }
249
250    fn map_in_ref(&self, from: &I, el: &<I as RingBase>::Element, _hom: &Self::Homomorphism) -> Self::Element {
251        self.from_f64(from.to_float_approx(el))
252    }
253}
254
255impl<I> CanHomFrom<RationalFieldBase<I>> for Complex64Base
256where
257    I: IntegerRingStore,
258    I::Type: IntegerRing,
259{
260    type Homomorphism = <Self as CanHomFrom<I::Type>>::Homomorphism;
261
262    fn has_canonical_hom(&self, from: &RationalFieldBase<I>) -> Option<Self::Homomorphism> {
263        self.has_canonical_hom(from.base_ring().get_ring())
264    }
265
266    fn map_in(
267        &self,
268        from: &RationalFieldBase<I>,
269        el: <RationalFieldBase<I> as RingBase>::Element,
270        hom: &Self::Homomorphism,
271    ) -> Self::Element {
272        self.map_in_ref(from, &el, hom)
273    }
274
275    fn map_in_ref(
276        &self,
277        from: &RationalFieldBase<I>,
278        el: &<RationalFieldBase<I> as RingBase>::Element,
279        hom: &Self::Homomorphism,
280    ) -> Self::Element {
281        self.div(
282            &self.map_in_ref(from.base_ring().get_ring(), from.num(el), hom),
283            &self.map_in_ref(from.base_ring().get_ring(), from.den(el), hom),
284        )
285    }
286}
287
288#[test]
289fn test_pow() {
290    let CC = Complex64::RING;
291    let i = Complex64::I;
292    assert!(CC.is_approx_eq(CC.negate(i), CC.pow(i, 3), 1));
293    assert!(!CC.is_approx_eq(CC.negate(i), CC.pow(i, 1024 + 3), 1));
294    assert!(CC.is_approx_eq(CC.negate(i), CC.pow(i, 1024 + 3), 100));
295    assert!(CC.is_approx_eq(
296        CC.exp(CC.mul(CC.from_f64(PI / 4.0), i)),
297        CC.mul(CC.add(CC.one(), i), CC.from_f64(2.0f64.powf(-0.5))),
298        1
299    ));
300
301    let seventh_root_of_unity = CC.exp(CC.mul(i, CC.from_f64(2.0 * PI / 7.0)));
302    assert!(CC.is_approx_eq(CC.pow(seventh_root_of_unity, 7 * 100 + 1), seventh_root_of_unity, 1000));
303}
304
305#[test]
306fn test_mul() {
307    let CC = Complex64::RING;
308    let i = Complex64::I;
309    assert!(CC.is_approx_eq(CC.mul(i, i), CC.from_f64(-1.0), 1));
310    assert!(CC.is_approx_eq(CC.mul(i, CC.negate(i)), CC.from_f64(1.0), 1));
311    assert!(CC.is_approx_eq(CC.mul(CC.add(i, CC.one()), i), CC.sub(i, CC.one()), 1));
312}
feanor_math/rings/float_complex.rs

feanor_math/rings/
float_complex.rs
