malachite_float/basic/
extended.rs

1// Copyright © 2026 Mikhail Hogrefe
2//
3// This file is part of Malachite.
4//
5// Malachite is free software: you can redistribute it and/or modify it under the terms of the GNU
6// Lesser General Public License (LGPL) as published by the Free Software Foundation; either version
7// 3 of the License, or (at your option) any later version. See <https://www.gnu.org/licenses/>.
8
9use crate::Float;
10use crate::InnerFloat::{Finite, NaN};
11use crate::conversion::rational_from_float::RationalFromFloatError;
12use core::cmp::Ordering::{self, *};
13use core::cmp::max;
14use core::mem::swap;
15use core::ops::{Add, Mul, Shr, ShrAssign};
16use malachite_base::num::arithmetic::traits::{CeilingLogBase2, Parity, Sqrt, SqrtAssign};
17use malachite_base::num::basic::integers::PrimitiveInt;
18use malachite_base::num::basic::traits::Zero;
19use malachite_base::num::conversion::traits::{ExactFrom, SaturatingInto};
20use malachite_base::num::logic::traits::SignificantBits;
21use malachite_base::rounding_modes::RoundingMode::{self, *};
22use malachite_nz::natural::arithmetic::float_extras::float_can_round;
23use malachite_nz::natural::arithmetic::float_sub::exponent_shift_compare;
24use malachite_nz::platform::Limb;
25use malachite_q::Rational;
26
27pub_crate_test_struct! {
28#[derive(Clone)]
29ExtendedFloat {
30    pub(crate) x: Float,
31    exp: i64,
32}}
33
34impl ExtendedFloat {
35    fn is_valid(&self) -> bool {
36        if self.x == 0u32 && self.exp != 0 {
37            return false;
38        }
39        let exp = self.x.get_exponent();
40        exp.is_none() || exp == Some(0)
41    }
42
43    fn from_rational_prec_round(value: Rational, prec: u64, rm: RoundingMode) -> (Self, Ordering) {
44        if value == 0 {
45            return (
46                Self {
47                    x: Float::ZERO,
48                    exp: 0,
49                },
50                Equal,
51            );
52        }
53        let exp = value.floor_log_base_2_abs() + 1;
54        let (x, o) = Float::from_rational_prec_round(value >> exp, prec, rm);
55        let new_exp = x.get_exponent().unwrap();
56        (
57            Self {
58                x: x >> new_exp,
59                exp: i64::from(new_exp) + exp,
60            },
61            o,
62        )
63    }
64
65    pub(crate) fn from_rational_prec_round_ref(
66        value: &Rational,
67        prec: u64,
68        rm: RoundingMode,
69    ) -> (Self, Ordering) {
70        if *value == 0 {
71            return (
72                Self {
73                    x: Float::ZERO,
74                    exp: 0,
75                },
76                Equal,
77            );
78        }
79        let exp = value.floor_log_base_2_abs() + 1;
80        if exp >= i64::from(Float::MIN_EXPONENT) && exp <= i64::from(Float::MAX_EXPONENT) {
81            let (x, o) = Float::from_rational_prec_round_ref(value, prec, rm);
82            let exp = x.get_exponent().unwrap();
83            return (
84                Self {
85                    x: x >> exp,
86                    exp: i64::from(exp),
87                },
88                o,
89            );
90        }
91        let (x, o) = Float::from_rational_prec_round(value >> exp, prec, rm);
92        let new_exp = x.get_exponent().unwrap();
93        (
94            Self {
95                x: x >> new_exp,
96                exp: i64::from(new_exp) + exp,
97            },
98            o,
99        )
100    }
101
102    fn add_prec_ref_ref(&self, other: &Self, prec: u64) -> (Self, Ordering) {
103        assert!(self.is_valid());
104        assert!(other.is_valid());
105        assert!(self.x.is_normal());
106        assert!(other.x.is_normal());
107        Self::from_rational_prec_round(
108            Rational::exact_from(self) + Rational::exact_from(other),
109            prec,
110            Nearest,
111        )
112    }
113
114    fn sub_prec_ref_ref(&self, other: &Self, prec: u64) -> (Self, Ordering) {
115        assert!(self.is_valid());
116        assert!(other.is_valid());
117        assert!(self.x.is_normal());
118        assert!(other.x.is_normal());
119        Self::from_rational_prec_round(
120            Rational::exact_from(self) - Rational::exact_from(other),
121            prec,
122            Nearest,
123        )
124    }
125
126    fn sub_prec(self, other: Self, prec: u64) -> (Self, Ordering) {
127        assert!(self.is_valid());
128        assert!(other.is_valid());
129        assert!(self.x.is_normal());
130        assert!(other.x.is_normal());
131        Self::from_rational_prec_round(
132            Rational::exact_from(self) - Rational::exact_from(other),
133            prec,
134            Nearest,
135        )
136    }
137
138    fn mul_prec_ref_ref(&self, other: &Self, prec: u64) -> (Self, Ordering) {
139        assert!(self.is_valid());
140        assert!(other.is_valid());
141        assert!(self.x.is_normal());
142        assert!(other.x.is_normal());
143        let (mut product, o) = self.x.mul_prec_ref_ref(&other.x, prec);
144        let mut product_exp = self.exp + other.exp;
145        let extra_exp = product.get_exponent().unwrap();
146        product >>= extra_exp;
147        product_exp = product_exp.checked_add(i64::from(extra_exp)).unwrap();
148        (
149            Self {
150                x: product,
151                exp: product_exp,
152            },
153            o,
154        )
155    }
156
157    fn div_prec_val_ref(self, other: &Self, prec: u64) -> (Self, Ordering) {
158        assert!(self.is_valid());
159        assert!(other.is_valid());
160        assert!(self.x.is_normal());
161        assert!(other.x.is_normal());
162        let (mut quotient, o) = self.x.div_prec_ref_ref(&other.x, prec);
163        let mut quotient_exp = self.exp - other.exp;
164        let extra_exp = quotient.get_exponent().unwrap();
165        quotient >>= extra_exp;
166        quotient_exp = quotient_exp.checked_add(i64::from(extra_exp)).unwrap();
167        (
168            Self {
169                x: quotient,
170                exp: quotient_exp,
171            },
172            o,
173        )
174    }
175
176    fn div_prec_assign_ref(&mut self, other: &Self, prec: u64) -> Ordering {
177        let mut x = Self {
178            x: Float::ZERO,
179            exp: 0,
180        };
181        swap(self, &mut x);
182        let (q, o) = x.div_prec_val_ref(other, prec);
183        *self = q;
184        o
185    }
186
187    fn square_round_assign(&mut self, rm: RoundingMode) -> Ordering {
188        let mut x = Self {
189            x: Float::ZERO,
190            exp: 0,
191        };
192        swap(self, &mut x);
193        let (mut square, o) = x.x.square_round(rm);
194        let mut square_exp = x.exp << 1;
195        let extra_exp = square.get_exponent().unwrap();
196        square >>= extra_exp;
197        square_exp = square_exp.checked_add(i64::from(extra_exp)).unwrap();
198        *self = Self {
199            x: square,
200            exp: square_exp,
201        };
202        o
203    }
204
205    fn from_extended_float_prec_round(x: Self, prec: u64, rm: RoundingMode) -> (Self, Ordering) {
206        if let Ok(x) = Rational::try_from(&x) {
207            Self::from_rational_prec_round(x, prec, rm)
208        } else {
209            (x, Equal)
210        }
211    }
212
213    pub_crate_test! {from_extended_float_prec_round_ref(
214        x: &Self,
215        prec: u64,
216        rm: RoundingMode,
217    ) -> (Self, Ordering) {
218        if let Ok(x) = Rational::try_from(x) {
219            Self::from_rational_prec_round(x, prec, rm)
220        } else {
221            (x.clone(), Equal)
222        }
223    }}
224
225    fn shr_prec_round<T: PrimitiveInt>(
226        self,
227        bits: T,
228        prec: u64,
229        rm: RoundingMode,
230    ) -> (Self, Ordering) {
231        // assumes no overflow or underflow
232        let (out_x, o) = Self::from_extended_float_prec_round(self, prec, rm);
233        let out_exp =
234            i64::exact_from(i128::from(out_x.exp) - SaturatingInto::<i128>::saturating_into(bits));
235        (
236            Self {
237                x: out_x.x,
238                exp: out_exp,
239            },
240            o,
241        )
242    }
243
244    fn shr_round_assign<T: PrimitiveInt>(&mut self, bits: T, _rm: RoundingMode) -> Ordering {
245        // assumes no overflow or underflow
246        if self.x.is_normal() {
247            let out_exp = i64::exact_from(
248                i128::from(self.exp) - SaturatingInto::<i128>::saturating_into(bits),
249            );
250            self.exp = out_exp;
251        }
252        Equal
253    }
254
255    pub(crate) fn into_float_helper(
256        mut self,
257        prec: u64,
258        rm: RoundingMode,
259        previous_o: Ordering,
260    ) -> (Float, Ordering) {
261        let o = self
262            .x
263            .shl_prec_round_assign_helper(self.exp, prec, rm, previous_o);
264        (self.x, o)
265    }
266
267    pub(crate) fn increment(&mut self) {
268        self.x.increment();
269        if let Some(exp) = self.x.get_exponent()
270            && exp == 1
271        {
272            self.x >>= 1u32;
273            self.exp = 0;
274        }
275    }
276}
277
278impl From<Float> for ExtendedFloat {
279    fn from(value: Float) -> Self {
280        if let Some(exp) = value.get_exponent() {
281            Self {
282                x: value >> exp,
283                exp: i64::from(exp),
284            }
285        } else {
286            Self { x: value, exp: 0 }
287        }
288    }
289}
290
291impl TryFrom<ExtendedFloat> for Rational {
292    type Error = RationalFromFloatError;
293
294    fn try_from(value: ExtendedFloat) -> Result<Self, Self::Error> {
295        Self::try_from(value.x).map(|x| x << value.exp)
296    }
297}
298
299impl<'a> TryFrom<&'a ExtendedFloat> for Rational {
300    type Error = RationalFromFloatError;
301
302    fn try_from(value: &'a ExtendedFloat) -> Result<Self, Self::Error> {
303        Self::try_from(&value.x).map(|x| x << value.exp)
304    }
305}
306
307impl PartialEq for ExtendedFloat {
308    fn eq(&self, other: &Self) -> bool {
309        self.x == other.x && self.exp == other.exp
310    }
311}
312
313impl PartialOrd for ExtendedFloat {
314    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
315        assert!(self.is_valid());
316        assert!(other.is_valid());
317        let self_sign = self.x > 0u32;
318        let other_sign = other.x > 0u32;
319        match self_sign.cmp(&other_sign) {
320            Greater => Some(Greater),
321            Less => Some(Less),
322            Equal => match self.exp.cmp(&other.exp) {
323                Greater => Some(if self_sign { Greater } else { Less }),
324                Less => Some(if self_sign { Less } else { Greater }),
325                Equal => self.x.partial_cmp(&other.x),
326            },
327        }
328    }
329}
330
331impl Add<&ExtendedFloat> for &ExtendedFloat {
332    type Output = ExtendedFloat;
333
334    fn add(self, other: &ExtendedFloat) -> Self::Output {
335        let prec = max(self.x.significant_bits(), other.x.significant_bits());
336        self.add_prec_ref_ref(other, prec).0
337    }
338}
339
340impl Mul<&ExtendedFloat> for &ExtendedFloat {
341    type Output = ExtendedFloat;
342
343    fn mul(self, other: &ExtendedFloat) -> Self::Output {
344        let prec = max(self.x.significant_bits(), other.x.significant_bits());
345        self.mul_prec_ref_ref(other, prec).0
346    }
347}
348
349impl SqrtAssign for ExtendedFloat {
350    fn sqrt_assign(&mut self) {
351        if self.exp.odd() {
352            self.x <<= 1;
353            self.exp = self.exp.checked_sub(1).unwrap();
354        }
355        self.x.sqrt_assign();
356        self.exp >>= 1;
357        if let Some(new_exp) = self.x.get_exponent() {
358            self.exp = self.exp.checked_add(i64::from(new_exp)).unwrap();
359            self.x >>= new_exp;
360        }
361        assert!(self.is_valid());
362    }
363}
364
365impl Sqrt for ExtendedFloat {
366    type Output = Self;
367
368    fn sqrt(mut self) -> Self::Output {
369        self.sqrt_assign();
370        self
371    }
372}
373
374impl Shr<u32> for ExtendedFloat {
375    type Output = Self;
376
377    fn shr(mut self, bits: u32) -> Self::Output {
378        self.shr_round_assign(bits, Nearest);
379        self
380    }
381}
382
383impl ShrAssign<u32> for ExtendedFloat {
384    fn shr_assign(&mut self, bits: u32) {
385        self.shr_round_assign(bits, Nearest);
386    }
387}
388
389fn cmp2_helper_extended(b: &ExtendedFloat, c: &ExtendedFloat, cancel: &mut u64) -> Ordering {
390    match (&b.x, &c.x) {
391        (
392            Float(Finite {
393                precision: x_prec,
394                significand: x,
395                ..
396            }),
397            Float(Finite {
398                precision: y_prec,
399                significand: y,
400                ..
401            }),
402        ) => {
403            let (o, c) = exponent_shift_compare(
404                x.as_limbs_asc(),
405                b.exp,
406                *x_prec,
407                y.as_limbs_asc(),
408                c.exp,
409                *y_prec,
410            );
411            *cancel = c;
412            o
413        }
414        _ => panic!(),
415    }
416}
417
418pub(crate) fn agm_prec_round_normal_extended(
419    mut a: ExtendedFloat,
420    mut b: ExtendedFloat,
421    prec: u64,
422    rm: RoundingMode,
423) -> (ExtendedFloat, Ordering) {
424    if a.x < 0u32 || b.x < 0u32 {
425        return (
426            ExtendedFloat {
427                x: float_nan!(),
428                exp: 0,
429            },
430            Equal,
431        );
432    }
433    let q = prec;
434    let mut working_prec = q + q.ceiling_log_base_2() + 15;
435    // b (op2) and a (op1) are the 2 operands but we want b >= a
436    match a.partial_cmp(&b).unwrap() {
437        Equal => return ExtendedFloat::from_extended_float_prec_round(a, prec, rm),
438        Greater => swap(&mut a, &mut b),
439        _ => {}
440    }
441    let mut increment = Limb::WIDTH;
442    let mut v;
443    let mut scaleit;
444    loop {
445        let mut err: u64 = 0;
446        let mut u = a.mul_prec_ref_ref(&b, working_prec).0;
447        v = a.add_prec_ref_ref(&b, working_prec).0;
448        u.sqrt_assign();
449        v >>= 1u32;
450        scaleit = 0;
451        let mut n: u64 = 1;
452        let mut eq = 0;
453        while cmp2_helper_extended(&u, &v, &mut eq) != Equal && eq <= working_prec - 2 {
454            let mut vf;
455            vf = (&u + &v) >> 1;
456            // See proof in algorithms.tex
457            if eq > working_prec >> 2 {
458                // vf = V(k)
459                let low_p = (working_prec + 1) >> 1;
460                let mut w = v.sub_prec_ref_ref(&u, low_p).0; // e = V(k-1)-U(k-1)
461                w.square_round_assign(Nearest); // e = e^2
462                w.shr_round_assign(4, Nearest); // e*= (1/2)^2*1/4
463                w.div_prec_assign_ref(&vf, low_p); // 1/4*e^2/V(k)
464                let vf_exp = vf.exp;
465                v = vf.sub_prec(w, working_prec).0;
466                // 0 or 1
467                err = u64::exact_from(vf_exp - v.exp);
468                break;
469            }
470            let uf = &u * &v;
471            u = uf.sqrt();
472            swap(&mut v, &mut vf);
473            n += 1;
474        }
475        // the error on v is bounded by (18n+51) ulps, or twice if there was an exponent loss in the
476        // final subtraction
477        //
478        // 18n+51 should not overflow since n is about log(p)
479        err += (18 * n + 51).ceiling_log_base_2();
480        // we should have n+2 <= 2^(p/4) [see algorithms.tex]
481        if (n + 2).ceiling_log_base_2() <= working_prec >> 2
482            && float_can_round(v.x.significand_ref().unwrap(), working_prec - err, q, rm)
483        {
484            break;
485        }
486        working_prec += increment;
487        increment = working_prec >> 1;
488    }
489    v.shr_prec_round(scaleit, prec, rm)
490}
491
492pub_crate_test! {agm_prec_round_normal_ref_ref_extended<'a>(
493    mut a: &'a ExtendedFloat,
494    mut b: &'a ExtendedFloat,
495    prec: u64,
496    rm: RoundingMode,
497) -> (ExtendedFloat, Ordering) {
498    if a.x < 0u32 || b.x < 0u32 {
499        return (
500            ExtendedFloat {
501                x: float_nan!(),
502                exp: 0,
503            },
504            Equal,
505        );
506    }
507    let q = prec;
508    let mut working_prec = q + q.ceiling_log_base_2() + 15;
509    // b (op2) and a (op1) are the 2 operands but we want b >= a
510    match a.partial_cmp(b).unwrap() {
511        Equal => return ExtendedFloat::from_extended_float_prec_round_ref(a, prec, rm),
512        Greater => swap(&mut a, &mut b),
513        _ => {}
514    }
515    let mut increment = Limb::WIDTH;
516    let mut v;
517    let mut scaleit;
518    loop {
519        let mut err: u64 = 0;
520        let mut u = a.mul_prec_ref_ref(b, working_prec).0;
521        v = a.add_prec_ref_ref(b, working_prec).0;
522        u.sqrt_assign();
523        v >>= 1u32;
524        scaleit = 0;
525        let mut n: u64 = 1;
526        let mut eq = 0;
527        while cmp2_helper_extended(&u, &v, &mut eq) != Equal && eq <= working_prec - 2 {
528            let mut vf;
529            vf = (&u + &v) >> 1;
530            // See proof in algorithms.tex
531            if eq > working_prec >> 2 {
532                // vf = V(k)
533                let low_p = (working_prec + 1) >> 1;
534                let mut w = v.sub_prec_ref_ref(&u, low_p).0; // e = V(k-1)-U(k-1)
535                assert!(w.is_valid());
536                assert!(w.x.is_normal());
537                w.square_round_assign(Nearest); // e = e^2
538                assert!(w.is_valid());
539                assert!(w.x.is_normal());
540                w.shr_round_assign(4, Nearest); // e*= (1/2)^2*1/4
541                assert!(w.is_valid());
542                assert!(w.x.is_normal());
543                w.div_prec_assign_ref(&vf, low_p); // 1/4*e^2/V(k)
544                let vf_exp = vf.exp;
545                v = vf.sub_prec(w, working_prec).0;
546                // 0 or 1
547                err = u64::exact_from(vf_exp - v.exp);
548                break;
549            }
550            let uf = &u * &v;
551            u = uf.sqrt();
552            swap(&mut v, &mut vf);
553            n += 1;
554        }
555        // the error on v is bounded by (18n+51) ulps, or twice if there was an exponent loss in the
556        // final subtraction
557        //
558        // 18n+51 should not overflow since n is about log(p)
559        err += (18 * n + 51).ceiling_log_base_2();
560        // we should have n+2 <= 2^(p/4) [see algorithms.tex]
561        if (n + 2).ceiling_log_base_2() <= working_prec >> 2
562            && float_can_round(v.x.significand_ref().unwrap(), working_prec - err, q, rm)
563        {
564            break;
565        }
566        working_prec += increment;
567        increment = working_prec >> 1;
568    }
569    v.shr_prec_round(scaleit, prec, rm)
570}}
malachite_float/basic/extended.rs

malachite_float/basic/
extended.rs
