feanor_math/rings/approx_real/
float.rs

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

feanor_math/rings/approx_real/
float.rs
