feanor_math/rings/approx_real/
float.rs

1use core::f64;
2use std::f64::EPSILON;
3
4use crate::algorithms::convolution::KaratsubaHint;
5use crate::algorithms::matmul::StrassenHint;
6use crate::ordered::OrderedRing;
7use crate::pid::{EuclideanRing, PrincipalIdealRing};
8use crate::field::Field;
9use crate::integer::{int_cast, IntegerRing, IntegerRingStore};
10use crate::primitive_int::StaticRing;
11use crate::rings::approx_real::{ApproxRealField, SqrtRing};
12use crate::{impl_eq_based_self_iso, ring::*};
13use crate::homomorphism::*;
14use crate::divisibility::{DivisibilityRing, Domain};
15use crate::rings::rational::{RationalField, RationalFieldBase};
16
17///
18/// An approximate implementation of the real numbers `R`, using 64 bit floating
19/// point numbers.
20/// 
21/// # Warning
22/// 
23/// Since floating point numbers do not exactly represent the real numbers, and this crate follows
24/// a mathematically precise approach, we cannot provide any function related to equality.
25/// In particular, `Real64Base.eq_el(a, b)` is not supported, and will panic. 
26/// Hence, this ring has only limited use within this crate, and is currently only used for
27/// floating-point FFTs and some approximate computations in the LLL algorithm. 
28/// 
29#[derive(Clone, Copy, PartialEq, Debug)]
30pub struct Real64Base;
31
32///
33/// [`RingStore`] corresponding to [`Real64Base`]
34/// 
35pub type Real64 = RingValue<Real64Base>;
36
37impl Real64 {
38
39    ///
40    /// The singleton ring instance of [`Real64`].
41    /// 
42    pub const RING: RingValue<Real64Base> = RingValue::from(Real64Base);
43}
44
45impl Real64Base {
46
47    pub fn is_absolute_approx_eq(&self, lhs: <Self as RingBase>::Element, rhs: <Self as RingBase>::Element, absolute_threshold: f64) -> bool {
48        (lhs - rhs).abs() < absolute_threshold
49    }
50
51    pub fn is_relative_approx_eq(&self, lhs: <Self as RingBase>::Element, rhs: <Self as RingBase>::Element, relative_threshold: f64) -> bool {
52        self.is_absolute_approx_eq(lhs, rhs, (lhs.abs() + rhs.abs()) * relative_threshold)
53    }
54
55    pub fn is_approx_eq(&self, lhs: <Self as RingBase>::Element, rhs: <Self as RingBase>::Element, precision: u64) -> bool {
56        let scaled_precision = precision as f64 * EPSILON;
57        if self.is_absolute_approx_eq(lhs, self.zero(), scaled_precision) {
58            self.is_absolute_approx_eq(rhs, self.zero(), scaled_precision)
59        } else {
60            self.is_relative_approx_eq(lhs, rhs, scaled_precision)
61        }
62    }
63}
64
65impl RingBase for Real64Base {
66 
67    type Element = f64;
68    
69    fn clone_el(&self, val: &Self::Element) -> Self::Element {
70        *val
71    }
72
73    fn add_assign(&self, lhs: &mut Self::Element, rhs: Self::Element) {
74        *lhs += rhs;
75    }
76
77    fn negate_inplace(&self, x: &mut Self::Element) {
78        *x = -*x;
79    }
80
81    fn mul_assign(&self, lhs: &mut Self::Element, rhs: Self::Element) {
82        *lhs *= rhs;
83    }
84
85    fn from_int(&self, value: i32) -> Self::Element {
86        value as f64
87    }
88    
89    fn eq_el(&self, _: &Self::Element, _: &Self::Element) -> bool {
90        panic!("Cannot provide equality on approximate rings")
91    }
92
93    fn pow_gen<R: IntegerRingStore>(&self, x: Self::Element, power: &El<R>, integers: R) -> Self::Element 
94        where R::Type: IntegerRing
95    {
96        if integers.get_ring().representable_bits().is_some() && integers.get_ring().representable_bits().unwrap() < i32::BITS as usize {
97            x.powi(int_cast(integers.clone_el(power), &StaticRing::<i32>::RING, integers))
98        } else {
99            x.powf(integers.to_float_approx(power))
100        }
101    }
102
103    fn is_commutative(&self) -> bool { true }
104
105    fn is_noetherian(&self) -> bool { true }
106
107    fn is_approximate(&self) -> bool { true }
108
109    fn dbg_within<'a>(&self, x: &Self::Element, out: &mut std::fmt::Formatter<'a>, _: EnvBindingStrength) -> std::fmt::Result {
110        write!(out, "{}", x)
111    }
112    
113    fn characteristic<I: IntegerRingStore + Copy>(&self, ZZ: I) -> Option<El<I>>
114        where I::Type: IntegerRing
115    {
116        Some(ZZ.zero())
117    }
118}
119
120impl_eq_based_self_iso!{ Real64Base }
121
122impl Domain for Real64Base {}
123
124impl DivisibilityRing for Real64Base {
125
126    fn checked_left_div(&self, lhs: &Self::Element, rhs: &Self::Element) -> Option<Self::Element> {
127        assert!(*rhs != 0.);
128        return Some(*lhs / *rhs);
129    }
130}
131
132impl PrincipalIdealRing for Real64Base {
133
134    fn checked_div_min(&self, lhs: &Self::Element, rhs: &Self::Element) -> Option<Self::Element> {
135        self.checked_left_div(lhs, rhs)
136    }
137    
138    fn extended_ideal_gen(&self, _lhs: &Self::Element, _rhs: &Self::Element) -> (Self::Element, Self::Element, Self::Element) {
139        panic!("Since Complex64 is only approximate, this cannot be implemented properly")
140    }
141}
142
143impl EuclideanRing for Real64Base {
144
145    fn euclidean_div_rem(&self, _lhs: Self::Element, _rhs: &Self::Element) -> (Self::Element, Self::Element) {
146        panic!("Since Complex64 is only approximate, this cannot be implemented properly")
147    }
148
149    fn euclidean_deg(&self, _: &Self::Element) -> Option<usize> {
150        panic!("Since Complex64 is only approximate, this cannot be implemented properly")
151    }
152}
153
154impl StrassenHint for Real64Base {
155    fn strassen_threshold(&self) -> usize {
156        // disable Strassen's algorithm, as it is very numerically unstable
157        usize::MAX
158    }
159}
160
161impl KaratsubaHint for Real64Base {
162    fn karatsuba_threshold(&self) -> usize {
163        // disable Karatsuba's algorithm, as it is very numerically unstable
164        usize::MAX
165    }
166}
167
168impl Field for Real64Base {
169
170    fn div(&self, lhs: &Self::Element, rhs: &Self::Element) -> Self::Element {
171        self.checked_left_div(lhs, rhs).unwrap()
172    }
173}
174
175impl OrderedRing for Real64Base {
176
177    fn cmp(&self, lhs: &Self::Element, rhs: &Self::Element) -> std::cmp::Ordering {
178        f64::partial_cmp(lhs, rhs).unwrap()
179    }
180}
181
182impl<I> CanHomFrom<I> for Real64Base 
183    where I: ?Sized + IntegerRing
184{
185    type Homomorphism = ();
186
187    fn has_canonical_hom(&self, _from: &I) -> Option<Self::Homomorphism> {
188        Some(())
189    }
190
191    fn map_in(&self, from: &I, el: <I as RingBase>::Element, _hom: &Self::Homomorphism) -> Self::Element {
192        from.to_float_approx(&el)
193    }
194
195    fn map_in_ref(&self, from: &I, el: &<I as RingBase>::Element, _hom: &Self::Homomorphism) -> Self::Element {
196        from.to_float_approx(el)
197    }
198}
199
200impl<I> CanHomFrom<RationalFieldBase<I>> for Real64Base 
201    where I: IntegerRingStore,
202        I::Type: IntegerRing
203{
204    type Homomorphism = ();
205
206    fn has_canonical_hom(&self, _from: &RationalFieldBase<I>) -> Option<Self::Homomorphism> {
207        Some(())
208    }
209
210    fn map_in(&self, from: &RationalFieldBase<I>, el: El<RationalField<I>>, hom: &Self::Homomorphism) -> Self::Element {
211        self.map_in_ref(from, &el, hom)
212    }
213
214    fn map_in_ref(&self, from: &RationalFieldBase<I>, el: &El<RationalField<I>>, _hom: &Self::Homomorphism) -> Self::Element {
215        from.base_ring().to_float_approx(from.num(el)) / from.base_ring().to_float_approx(from.den(el))
216    }
217}
218
219impl ApproxRealField for Real64Base {
220
221    fn epsilon(&self) -> &Self::Element {
222        &f64::EPSILON
223    }
224
225    fn infinity(&self) -> Self::Element {
226        f64::INFINITY
227    }
228
229    fn round_to_integer<I>(&self, ZZ: I, x: Self::Element) -> Option<El<I>>
230        where I: RingStore, I::Type: IntegerRing
231    {
232        ZZ.from_float_approx(x.round())
233    }
234}
235
236impl SqrtRing for Real64Base {
237
238    fn sqrt(&self, x: Self::Element) -> Self::Element {
239        x.sqrt()
240    }
241}
feanor_math/rings/approx_real/float.rs

feanor_math/rings/approx_real/
float.rs
