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

malachite_float/basic/
extended.rs
