feanor_math/rings/
float_complex.rs

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

feanor_math/rings/
float_complex.rs
