/home/a220/proj/radnelac/src/common/math.rs

Source
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at https://mozilla.org/MPL/2.0/.

use crate::common::error::CalendarError;
use num_traits::AsPrimitive;
use num_traits::Bounded;
use num_traits::Euclid;
use num_traits::FromPrimitive;
use num_traits::NumAssign;
use num_traits::ToPrimitive;
use num_traits::Zero;
use std::cmp::PartialOrd;

// https://en.m.wikipedia.org/wiki/Double-precision_floating-point_format
// > Between 2^52=4,503,599,627,370,496 and 2^53=9,007,199,254,740,992 the
// > representable numbers are exactly the integers. For the next range,
// > from 2^53 to 2^54, everything is multiplied by 2, so the representable
// > numbers are the even ones, etc. Conversely, for the previous range from
// > 2^51 to 2^52, the spacing is 0.5, etc.
// >
// > The spacing as a fraction of the numbers in the range from 2^n to 2^n+1
// > is 2^n−52.

// We want to represent seconds as fractions of a day, and represent days
// since any calendar's epoch as whole numbers. We should avoid using floating
// point in the ranges where that would be inaccurate.
// 1/(24 * 60 * 60) = 0.000011574074074074073499346534
// 2 ** (36 - 52)   = 0.000015258789062500000000000000 (n = 36 is too large)
// 2 ** (35 - 52)   = 0.000007629394531250000000000000 (n = 35 is small, but risks off by 1 second)
// 2 ** (34 - 52)   = 0.000003814697265625000000000000 (n = 34 is probably small enough)
// 2 ** 34          = 17179869184
// 2 ** 35          = 34359738368
// 2 ** 36          = 68719476736
// Converted into years, it's still a lot of time:
// 2 ** 34 / 365.25 = 47035918.36824093

pub const EFFECTIVE_MAX: f64 = 17179869184.0;
pub const EFFECTIVE_MIN: f64 = -EFFECTIVE_MAX;
pub const EQ_SCALE: f64 = EFFECTIVE_MAX;
pub const EFFECTIVE_EPSILON: f64 = 0.000003814697265625;

pub trait TermNum:
    NumAssign
    + PartialOrd
    + ToPrimitive
    + FromPrimitive
    + AsPrimitive<f64>
    + AsPrimitive<i64>
    + AsPrimitive<Self>
    + Euclid
    + Bounded
    + Copy
{
    fn approx_eq(self, other: Self) -> bool {
        self == other
    }

    fn approx_floor(self) -> Self {
        self
    }

    fn floor_round(self) -> Self {
        self
    }

    fn is_a_number(self) -> bool {
        true
    }

    fn modulus(self, other: Self) -> Self {
        let x = self;
        let y = other;
        Euclid::rem_euclid(&x, &y)
    }

    fn approx_eq_iter<T: IntoIterator<Item = Self>>(x: T, y: T) -> bool {
        !x.into_iter()
            .zip(y.into_iter())
            .any768(|(zx, zy)| !zx.approx_eq(zy))
    }

    fn sign(self) -> Self {
        if self.is_zero() {
            Self::zero()
        } else {
            Self::one()
        }
    }

    fn gcd(self, other: Self) -> Self {
        //LISTING 1.22 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        let x = self;
        let y = other;
        if y.is_zero() {
            x
        } else {
            y.gcd(x.modulus(y))
        }
    }

    fn lcm(self, other: Self) -> Self {
        //LISTING 1.23 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        let x = self;
        let y = other;
        (x * y) / x.gcd(y)
    }

    fn interval_modulus(self, a: Self, b: Self) -> Self {
        //LISTING 1.24 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        let x = self;
        if a == b {
            x
        } else {
            a + (x - a).modulus(b - a)
        }
    }

    fn adjusted_remainder(self, b: Self) -> Self {
        //LISTING 1.28 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        let x = self;
        let m = x.modulus(b);
        if m == Self::zero() {
            b
        } else {
            m
        }
    }

    fn sum<U: TermNum>(f: impl Fn(U) -> Self, p: impl Fn(U) -> bool, k: U) -> Self {
        //LISTING 1.30 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        //Modified to use iteration instead of recursion
        let mut result = Self::zero();
        let mut i = k;
        while p(i) {
            result += f(i);
            i += U::one();
        }
        result
    }

    fn product<U: TermNum>(f: impl Fn(U) -> Self, p: impl Fn(U) -> bool, k: U) -> Self {
        //LISTING 1.31 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        //Modified to user iteration instead of recursion
        let mut result = Self::one();
        let mut i = k;
        while p(i) {
            result *= f(i);
            i += U::one();
        }
        result
    }

    fn validate_mixed_radix(a: &[Self], b: &[Self]) -> Result<(), CalendarError> {
        if a.len() != (b.len() + 1) {
            Err(CalendarError::MixedRadixWrongSize)
        } else if b.iter()55.6k.any55.6k(|&bx| bx.is_zero()) {
            Err(CalendarError::MixedRadixZeroBase)
        } else {
            Ok(())
        }
    }

    fn from_mixed_radix(a: &[Self], b: &[Self], k: usize) -> Result<f64, CalendarError> {
        //LISTING 1.41 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        let n = b.len();
        let as_f64 = <Self as AsPrimitive<f64>>::as_;
        match TermNum::validate_mixed_radix(a, b) {
            Ok(()) => (),
            Err(error) => return Err(error),
        };

        let sum_mul: f64 = TermNum::sum(
            |i| a[i]11.4k * TermNum::product11.4k(|j| b[j2.04k], |j| j < k, i11.4k),
            |i| i <= k,
            0,
        )
        .as_();

        let sum_div: f64 = TermNum::sum(
            |i| as_f64(a[i])26.5k / as_f6426.5k(TermNum::product26.5k(|j| b[j52.0k], |j| j < i, k26.5k)),
            |i| i <= n,
            k + 1,
        );

        Ok(sum_mul + sum_div)
    }

    fn to_mixed_radix(x: f64, b: &[Self], k: usize, a: &mut [Self]) -> Result<(), CalendarError> {
        //LISTING 1.42 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        let n = b.len();
        match TermNum::validate_mixed_radix(a, b) {
            Ok(()) => (),
            Err(error) => return Err(error),
        };

        for i in 0..(n + 1)45.7k {
            if i == 0 {
                let p045.7k: f6445.7k = TermNum::product45.7k(|j| b[j1.53k], |j| j < k, 0).as_45.7k();
                let q0 = Self::from_f64((x / p0).approx_floor() as f64);
                match q0 {
                    Some(q) => a[i] = q,
                    None => return Err(CalendarError::ImpossibleResult),
                }
            } else if i > 0 && i < k {
                let p1512: f64512 = TermNum::product512(|j| b[j512], |j| j < k, i512).as_512();
                let q1 = Self::from_f64((x / p1).approx_floor() as f64);
                match q1 {
                    Some(q) => a[i] = q.modulus(b[i - 1]),
                    None => return Err(CalendarError::ImpossibleResult),
                }
            } else if i >= k && i < n {
                let p289.5k: f6489.5k = TermNum::product89.5k(|j| b[j133k], |j| j < i, k89.5k).as_89.5k();
                let q2 = Self::from_f64((x * p2).approx_floor() as f64);
                match q2 {
                    Some(q) => a[i] = q.modulus(b[i - 1]),
                    None => return Err(CalendarError::ImpossibleResult),
                }
            } else {
                let p345.7k: f6445.7k = TermNum::product45.7k(|j| b[j134k], |j| j < n, k45.7k).as_45.7k();
                let q3 = x * p3;
                let m = q3.modulus(b[n - 1].as_());
                if m.approx_eq(b[n - 1].as_()) || m.approx_eq(0.0) {
                    a[i] = Self::zero();
                } else if m18.1k.fract18.1k().approx_eq18.1k(1.0) {
                    a[i] = match Self::from_f64(m.ceil()) {
                        Some(m3) => m3,
                        None => return Err(CalendarError::ImpossibleResult),
                    };
                } else {
                    a[i] = match Self::from_f64(m) {
                        Some(m3) => m3,
                        None => return Err(CalendarError::ImpossibleResult),
                    };
                }
            }
        }
        Ok(())
    }

    fn search_min(p: impl Fn(Self) -> bool, k: Self) -> Self {
        //LISTING 1.32 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        //Modified to use iteration instead of recursion
        let mut i = k;
        while !p(i) {
            i += Self::one()
        }
        i
    }

    fn search_max(p: impl Fn(Self) -> bool, k: Self) -> Self {
        //LISTING 1.33 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        //Modified to use iteration instead of recursion
        let mut i = k - Self::one();
        while p(i) {
            i += Self::one()
        }
        i
    }
}

impl TermNum for usize {}
impl TermNum for u64 {}
impl TermNum for u32 {}
impl TermNum for u16 {}
impl TermNum for u8 {}

impl TermNum for i64 {
    fn sign(self) -> Self {
        if self.is_zero() {
            Self::zero()
        } else {
            self.signum()
        }
    }

    fn modulus(self, other: Self) -> Self {
        debug_assert!(other != Self::zero());
        if other > Self::zero() {
            let x = self;
            let y = other;
            Euclid::rem_euclid(&x, &y)
        } else {
            let x = -self;
            let y = -other;
            -Euclid::rem_euclid(&x, &y)
        }
    }
}

macro_rules! CastFn {
    ($n: ident, $t:ident) => {
        fn $n(self) -> Self {
            (self as $t).$n() as Self
        }
    };
    ($n: ident, $t:ident, $u: ident) => {
        fn $n(self, other: Self) -> $u {
            (self as $t).$n(other as $t) as $u
        }
    };
}

impl TermNum for i32 {
    CastFn!(sign, i64);
    CastFn!(modulus, i64, Self);
}

impl TermNum for i16 {
    CastFn!(sign, i64);
    CastFn!(modulus, i64, Self);
}

impl TermNum for i8 {
    CastFn!(sign, i64);
    CastFn!(modulus, i64, Self);
}

impl TermNum for f64 {
    fn sign(self) -> Self {
        //LISTING 1.16 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        //Modified to use `signum()`
        if self.is_zero() {
            Self::zero()
        } else {
            self.signum()
        }
    }

    fn approx_eq(self, other: Self) -> bool {
        let x = self;
        let y = other;
        if x == y {
            true
        } else if x.signum() != y.signum() && x != 0.00 && y != 0.00 {
            false
        } else {
            (x - y).abs() < (x.abs() / EQ_SCALE) || (x - y).abs() < EFFECTIVE_EPSILON134k
        }
    }

    fn approx_floor(self) -> Self {
        let x = self;
        let cx = x.ceil();
        if x.approx_eq(cx) {
            cx
        } else {
            x.floor()
        }
    }

    fn floor_round(self) -> Self {
        //LISTING 1.15 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        (self + 0.5).floor()
    }

    fn modulus(self, other: Self) -> Self {
        //LISTING 1.17 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
        let x = self;
        let y = other;
        debug_assert!(y != 0.0); //Cannot replace with NonZero - it's not supported for f64
        if x == 0.0 {
            0.0
        } else {
            x - (y * (x / y).floor())
        }
    }

    fn is_a_number(self) -> bool {
        !self.is_nan()
    }
}

impl TermNum for f32 {
    CastFn!(sign, f64);
    CastFn!(modulus, f64, Self);
    CastFn!(approx_eq, f64, bool);
    CastFn!(approx_floor, f64);
    CastFn!(floor_round, f64);

    fn is_a_number(self) -> bool {
        !self.is_nan()
    }
}

// TODO: binary search (listing 1.35)
// TODO: inverse f (listing 1.36)
// TODO: list_of_fixed_from_moments (listing 1.37)
// TODO: range (1.38)
// TODO: scan_range (1.39)
// TODO: positions_in_range (1.40)
// TODO: angles, minutes, degrees

#[cfg(test)]
mod tests {
    use super::*;
    use proptest::prop_assume;
    use proptest::proptest;

    #[test]
    fn modulus_basics() {
        assert_eq!((9.0).modulus(5.0), 4.0);
        assert_eq!((-9.0).modulus(5.0), 1.0);
        assert_eq!((9.0).modulus(-5.0), -1.0);
        assert_eq!((-9.0).modulus(-5.0), -4.0);
    }

    #[test]
    fn adjusted_remainder() {
        assert_eq!((13).adjusted_remainder(12), 1);
    }

    #[test]
    #[should_panic]
    fn modulus_zero() {
        (123.0).modulus(0.0);
    }

    #[test]
    fn gcd_wikipedia_examples() {
        //See https://en.wikipedia.org/wiki/Greatest_common_divisor
        assert_eq!(8.0.gcd(12.0), 4.0);
        assert_eq!(54.0.gcd(24.0), 6.0);
        assert_eq!(9.0.gcd(28.0), 1.0); //Coprime
        assert_eq!(24.0.gcd(60.0), 12.0);
        assert_eq!(42.0.gcd(56.0), 14.0);
    }

    #[test]
    fn lcm_wikipedia_examples() {
        //https://en.wikipedia.org/wiki/Least_common_multiple
        assert_eq!((5.0).lcm(4.0), 20.0);
        assert_eq!((6.0).lcm(4.0), 12.0);
    }

    #[test]
    fn sum_of_2x() {
        let y1 = TermNum::sum1(|x| x2 * 2.0, |i| i < 3.0, 1.0);
        assert_eq!(y, 3.0 * 2.0);
    }

    #[test]
    fn product_of_2x() {
        let y1 = TermNum::product1(|x| x2 * 2.0, |i| i < 3.0, 1.0);
        assert_eq!(y, 2.0 * 4.0);
    }

    #[test]
    fn search_min_sign() {
        let y1 = f64::search_min1(|i| i.sign() == 1.0, -10.0);
        assert_eq!(y, 1.0);
        let z = f64::search_min(|i| i.sign() == 1.0, 500.0);
        assert_eq!(z, 500.0);
    }

    #[test]
    fn search_max_sign() {
        let y1 = f64::search_max1(|i| i.sign() == -1.0, -10.0);
        assert_eq!(y, 0.0);
        let z = f64::search_max(|i| i.sign() == -1.0, 500.0);
        assert_eq!(z, 499.0);
    }

    proptest! {
        #[test]
        fn mixed_radix_time(ahr in 0..24,amn in 0..59,asc in 0..59) {
            let ahr = ahr as f64;
            let amn = amn as f64;
            let asc = asc as f64;
            let a = [ahr, amn, asc];
            let b = [60.0, 60.0];
            let seconds = TermNum::from_mixed_radix(&a, &b, 2).unwrap();
            let minutes = TermNum::from_mixed_radix(&a, &b, 1).unwrap();
            let hours = TermNum::from_mixed_radix(&a, &b, 0).unwrap();
            let expected_seconds = (ahr * 3600.0) + (amn* 60.0) + asc;
            let expected_minutes = (ahr * 60.0) + amn + (asc / 60.0);
            let expected_hours = ahr + (amn / 60.0) + (asc / 3600.0);
            assert!(seconds.approx_eq(expected_seconds));
            assert!(minutes.approx_eq(expected_minutes));
            assert!(hours.approx_eq(expected_hours));

            let mut a_seconds = [0.0, 0.0, 0.0];
            let mut a_minutes = [0.0, 0.0, 0.0];
            let mut a_hours = [0.0, 0.0, 0.0];
            TermNum::to_mixed_radix(seconds, &b, 2, &mut a_seconds).unwrap();
            TermNum::to_mixed_radix(minutes, &b, 1, &mut a_minutes).unwrap();
            TermNum::to_mixed_radix(hours, &b, 0, &mut a_hours).unwrap();

            println!("a: {a:?}, a_hours: {a_hours:?}, hours: {hours}");

            assert!(TermNum::approx_eq_iter(a_seconds, a));
            assert!(TermNum::approx_eq_iter(a_minutes, a));
            assert!(TermNum::approx_eq_iter(a_hours, a));
        }

        #[test]
        fn mixed_radix_time_i(ahr in 0..24,amn in 0..59,asc in 0..59) {
            let ahr = ahr as i32;
            let amn = amn as i32;
            let asc = asc as i32;
            let a = [ahr, amn, asc];
            let b = [60, 60];
            let seconds = TermNum::from_mixed_radix(&a, &b, 2).unwrap();
            let minutes = TermNum::from_mixed_radix(&a, &b, 1).unwrap();
            let hours = TermNum::from_mixed_radix(&a, &b, 0).unwrap();

            let ahr = ahr as f64;
            let amn = amn as f64;
            let asc = asc as f64;
            let expected_seconds = (ahr * 3600.0) + (amn* 60.0) + asc;
            let expected_minutes = (ahr * 60.0) + amn + (asc / 60.0);
            let expected_hours = ahr + (amn / 60.0) + (asc / 3600.0);
            assert!(seconds.approx_eq(expected_seconds));
            assert!(minutes.approx_eq(expected_minutes));
            assert!(hours.approx_eq(expected_hours));

            let mut a_seconds = [0, 0, 0];
            let mut a_minutes = [0, 0, 0];
            let mut a_hours = [0, 0, 0];
            TermNum::to_mixed_radix(seconds, &b, 2, &mut a_seconds).unwrap();
            TermNum::to_mixed_radix(minutes, &b, 1, &mut a_minutes).unwrap();
            TermNum::to_mixed_radix(hours, &b, 0, &mut a_hours).unwrap();

            println!("a: {a:?}, a_hours: {a_hours:?}, hours: {hours}");

            assert_eq!(&a_seconds, &a);
            assert_eq!(&a_minutes, &a);
            assert_eq!(&a_hours, &a);
        }


        #[test]
        fn modulus_positivity(x in EFFECTIVE_MIN..EFFECTIVE_MAX, y in 0.0..EFFECTIVE_MAX) {
            assert!((x as f64).modulus(y as f64) >= 0.0);
        }

        #[test]
        fn modulus_i_positivity(x: i32, y in 1..i32::MAX) {
            assert!(x.modulus(y) >= 0);
        }


        #[test]
        fn modulus_negative_x(x in 0.0..EFFECTIVE_MAX, y in 0.0..EFFECTIVE_MAX) {
            prop_assume!(y != 0.0);
            prop_assume!(x.modulus(y) != 0.0);
            let a0 = (-x).modulus(y);
            let a1 = y - x.modulus(y);
            assert!(a0.approx_eq(a1));
        }

        #[test]
        fn modulus_i_negative_x(x in 0..i32::MAX, y in 1..i32::MAX) {
            prop_assume!(x.modulus(y) != 0);
            let a0 = (-x).modulus(y);
            let a1 = y - x.modulus(y);
            assert_eq!(a0, a1);
        }

        #[test]
        fn modulus_mult(
            x in -131072.0..131072.0,
            y in -131072.0..131072.0,
            z in -131072.0..131072.0) {
            //LISTING 1.19 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
            //Using sqrt(EFFECTIVE_MAX) as limit
            let x = x as f64;
            let y = y as f64;
            let z = z as f64;
            prop_assume!(y != 0.0);
            prop_assume!(z != 0.0);
            let a: f64 = x.modulus(y);
            let az: f64 = (x*z).modulus(y*z);
            println!("x={}; y={}; z={}; a={}; a*z= {}; az= {};", x, y, z, a, a*z, az);
            assert!((a * z).approx_eq(az));
        }

        #[test]
        fn modulus_i_mult(
            x in i16::MIN..i16::MAX,
            y in i16::MIN..i16::MAX,
            z in i16::MIN..i16::MAX) {
            //LISTING 1.19 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
            let x = x as i32;
            let y = y as i32;
            let z = z as i32;
            prop_assume!(y != 0);
            prop_assume!(z != 0);
            let a = x.modulus(y);
            let az = (x*z).modulus(y*z);
            println!("x={}; y={}; z={}; a={}; a*z= {}; az= {};", x, y, z, a, a*z, az);
            assert_eq!(a * z, az);
        }

        #[test]
        fn modulus_mult_minus_1(x in 0.0..EFFECTIVE_MAX, y in 0.0..EFFECTIVE_MAX) {
            //LISTING 1.19 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
            prop_assume!(y != 0.0);
            let a0 = (-x).modulus(-y);
            let a1 = -(x.modulus(y));
            assert_eq!(a0, a1);
        }

        #[test]
        fn modulus_i_mult_minus_1(x in 0..i32::MAX, y in 1..i32::MAX) {
            //LISTING 1.19 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
            let a0 = (-x).modulus(-y);
            let a1 = -(x.modulus(y));
            assert_eq!(a0, a1);
        }

        #[test]
        fn modulus_i_multiple_of_y(x: i32, y: i32) {
            //LISTING 1.20 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
            prop_assume!(y != 0);
            let a = (x as i64) - (x.modulus(y) as i64);
            assert_eq!(a % (y as i64), 0);
        }

        #[test]
        fn modulus_bounds(x in EFFECTIVE_MIN..EFFECTIVE_MAX, y in EFFECTIVE_MIN..EFFECTIVE_MAX) {
            //LISTING 1.21 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
            prop_assume!(y != 0.0);
            let a = x.modulus(y) * y.sign();
            assert!(0.0 <= a && a < y.abs());
        }
        #[test]
        fn modulus_i_bounds(x: i32, y: i32) {
            //LISTING 1.21 (*Calendrical Calculations: The Ultimate Edition* by Reingold & Dershowitz.)
            prop_assume!(y != 0);
            let a = x.modulus(y) * (y.sign());
            assert!(0 <= a && a < y.abs());
        }
    }
}

Coverage Report

Created: 2025-08-13 21:02

Line	Count	Source
1		// This Source Code Form is subject to the terms of the Mozilla Public
2		// License, v. 2.0. If a copy of the MPL was not distributed with this
3		// file, You can obtain one at https://mozilla.org/MPL/2.0/.
4
5		use crate::common::error::CalendarError;
6		use num_traits::AsPrimitive;
7		use num_traits::Bounded;
8		use num_traits::Euclid;
9		use num_traits::FromPrimitive;
10		use num_traits::NumAssign;
11		use num_traits::ToPrimitive;
12		use num_traits::Zero;
13		use std::cmp::PartialOrd;
14
15		// https://en.m.wikipedia.org/wiki/Double-precision_floating-point_format
16		// > Between 2^52=4,503,599,627,370,496 and 2^53=9,007,199,254,740,992 the
17		// > representable numbers are exactly the integers. For the next range,
18		// > from 2^53 to 2^54, everything is multiplied by 2, so the representable
19		// > numbers are the even ones, etc. Conversely, for the previous range from
20		// > 2^51 to 2^52, the spacing is 0.5, etc.
21		// >
22		// > The spacing as a fraction of the numbers in the range from 2^n to 2^n+1
23		// > is 2^n−52.
24
25		// We want to represent seconds as fractions of a day, and represent days
26		// since any calendar's epoch as whole numbers. We should avoid using floating
27		// point in the ranges where that would be inaccurate.
28		// 1/(24 * 60 * 60) = 0.000011574074074074073499346534
29		// 2 ** (36 - 52) = 0.000015258789062500000000000000 (n = 36 is too large)
30		// 2 ** (35 - 52) = 0.000007629394531250000000000000 (n = 35 is small, but risks off by 1 second)
31		// 2 ** (34 - 52) = 0.000003814697265625000000000000 (n = 34 is probably small enough)
32		// 2 ** 34 = 17179869184
33		// 2 ** 35 = 34359738368
34		// 2 ** 36 = 68719476736
35		// Converted into years, it's still a lot of time:
36		// 2 ** 34 / 365.25 = 47035918.36824093
37
38		pub const EFFECTIVE_MAX: f64 = 17179869184.0;
39		pub const EFFECTIVE_MIN: f64 = -EFFECTIVE_MAX;
40		pub const EQ_SCALE: f64 = EFFECTIVE_MAX;
41		pub const EFFECTIVE_EPSILON: f64 = 0.000003814697265625;
42
43		pub trait TermNum:
44		NumAssign
45		+ PartialOrd
46		+ ToPrimitive
47		+ FromPrimitive
48		+ AsPrimitive<f64>
49		+ AsPrimitive<i64>
50		+ AsPrimitive<Self>
51		+ Euclid
52		+ Bounded
53		+ Copy
54		{
55	0	fn approx_eq(self, other: Self) -> bool {
56	0	self == other
57	0	}
58
59	0	fn approx_floor(self) -> Self {
60	0	self
61	0	}
62
63	0	fn floor_round(self) -> Self {
64	0	self
65	0	}
66
67	256	fn is_a_number(self) -> bool {
68	256	true
69	256	}
70
71	65.0k	fn modulus(self, other: Self) -> Self {
72	65.0k	let x = self;
73	65.0k	let y = other;
74	65.0k	Euclid::rem_euclid(&x, &y)
75	65.0k	}
76
77	768	fn approx_eq_iter<T: IntoIterator<Item = Self>>(x: T, y: T) -> bool {
78	768	!x.into_iter()
79	768	.zip(y.into_iter())
80	2.30k	. any768 (\|(zx, zy)\| !zx.approx_eq(zy))
81	768	}
82
83	0	fn sign(self) -> Self {
84	0	if self.is_zero() {
85	0	Self::zero()
86		} else {
87	0	Self::one()
88		}
89	0	}
90
91	25	fn gcd(self, other: Self) -> Self {
92		//LISTING 1.22 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
93	25	let x = self;
94	25	let y = other;
95	25	if y.is_zero() {
96	7	x
97		} else {
98	18	y.gcd(x.modulus(y))
99		}
100	25	}
101
102	2	fn lcm(self, other: Self) -> Self {
103		//LISTING 1.23 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
104	2	let x = self;
105	2	let y = other;
106	2	(x * y) / x.gcd(y)
107	2	}
108
109	12.0k	fn interval_modulus(self, a: Self, b: Self) -> Self {
110		//LISTING 1.24 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
111	12.0k	let x = self;
112	12.0k	if a == b {
113	0	x
114		} else {
115	12.0k	a + (x - a).modulus(b - a)
116		}
117	12.0k	}
118
119	106k	fn adjusted_remainder(self, b: Self) -> Self {
120		//LISTING 1.28 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
121	106k	let x = self;
122	106k	let m = x.modulus(b);
123	106k	if m == Self::zero() {
124	3.60k	b
125		} else {
126	102k	m
127		}
128	106k	}
129
130	19.7k	fn sum<U: TermNum>(f: impl Fn(U) -> Self, p: impl Fn(U) -> bool, k: U) -> Self {
131		//LISTING 1.30 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
132		//Modified to use iteration instead of recursion
133	19.7k	let mut result = Self::zero();
134	19.7k	let mut i = k;
135	57.6k	while p(i) {
136	37.9k	result += f(i);
137	37.9k	i += U::one();
138	37.9k	}
139	19.7k	result
140	19.7k	}
141
142	219k	fn product<U: TermNum>(f: impl Fn(U) -> Self, p: impl Fn(U) -> bool, k: U) -> Self {
143		//LISTING 1.31 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
144		//Modified to user iteration instead of recursion
145	219k	let mut result = Self::one();
146	219k	let mut i = k;
147	543k	while p(i) {
148	323k	result *= f(i);
149	323k	i += U::one();
150	323k	}
151	219k	result
152	219k	}
153
154	55.6k	fn validate_mixed_radix(a: &[Self], b: &[Self]) -> Result<(), CalendarError> {
155	55.6k	if a.len() != (b.len() + 1) {
156	0	Err(CalendarError::MixedRadixWrongSize)
157	163k	} else if b.iter()55.6k . any55.6k (\|&bx\| bx.is_zero()) {
158	0	Err(CalendarError::MixedRadixZeroBase)
159		} else {
160	55.6k	Ok(())
161		}
162	55.6k	}
163
164	9.87k	fn from_mixed_radix(a: &[Self], b: &[Self], k: usize) -> Result<f64, CalendarError> {
165		//LISTING 1.41 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
166	9.87k	let n = b.len();
167	9.87k	let as_f64 = <Self as AsPrimitive<f64>>::as_;
168	9.87k	match TermNum::validate_mixed_radix(a, b) {
169	9.87k	Ok(()) => (),
170	0	Err(error) => return Err(error),
171		};
172
173	9.87k	let sum_mul: f64 = TermNum::sum(
174	13.4k	\|i\| a[i]11.4k * TermNum::product11.4k (\|j\| b[ j2.04k ], \|j\| j < k, i11.4k ),
175	21.2k	\|i\| i <= k,
176		0,
177		)
178	9.87k	.as_();
179
180	9.87k	let sum_div: f64 = TermNum::sum(
181	78.5k	\|i\| as_f64(a[i])26.5k / as_f6426.5k ( TermNum::product26.5k (\|j\| b[ j52.0k ], \|j\| j < i, k26.5k )),
182	36.4k	\|i\| i <= n,
183	9.87k	k + 1,
184		);
185
186	9.87k	Ok(sum_mul + sum_div)
187	9.87k	}
188
189	45.7k	fn to_mixed_radix(x: f64, b: &[Self], k: usize, a: &mut [Self]) -> Result<(), CalendarError> {
190		//LISTING 1.42 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
191	45.7k	let n = b.len();
192	45.7k	match TermNum::validate_mixed_radix(a, b) {
193	45.7k	Ok(()) => (),
194	0	Err(error) => return Err(error),
195		};
196
197	181k	for i in 0.. (n + 1)45.7k {
198	181k	if i == 0 {
199	47.3k	let p045.7k : f6445.7k = TermNum::product45.7k (\|j\| b[ j1.53k ], \|j\| j < k, 0). as_45.7k ();
200	45.7k	let q0 = Self::from_f64((x / p0).approx_floor() as f64);
201	45.7k	match q0 {
202	45.7k	Some(q) => a[i] = q,
203	0	None => return Err(CalendarError::ImpossibleResult),
204		}
205	135k	} else if i > 0 && i < k {
206	1.02k	let p1512 : f64512 = TermNum::product512 (\|j\| b[ j512 ], \|j\| j < k, i512 ). as_512 ();
207	512	let q1 = Self::from_f64((x / p1).approx_floor() as f64);
208	512	match q1 {
209	512	Some(q) => a[i] = q.modulus(b[i - 1]),
210	0	None => return Err(CalendarError::ImpossibleResult),
211		}
212	135k	} else if i >= k && i < n {
213	222k	let p289.5k : f6489.5k = TermNum::product89.5k (\|j\| b[ j133k ], \|j\| j < i, k89.5k ). as_89.5k ();
214	89.5k	let q2 = Self::from_f64((x * p2).approx_floor() as f64);
215	89.5k	match q2 {
216	89.5k	Some(q) => a[i] = q.modulus(b[i - 1]),
217	0	None => return Err(CalendarError::ImpossibleResult),
218		}
219		} else {
220	180k	let p345.7k : f6445.7k = TermNum::product45.7k (\|j\| b[ j134k ], \|j\| j < n, k45.7k ). as_45.7k ();
221	45.7k	let q3 = x * p3;
222	45.7k	let m = q3.modulus(b[n - 1].as_());
223	45.7k	if m.approx_eq(b[n - 1].as_()) \|\| m.approx_eq(0.0) {
224	27.6k	a[i] = Self::zero();
225	27.6k	} else if m18.1k . fract18.1k (). approx_eq18.1k (1.0) {
226	74	a[i] = match Self::from_f64(m.ceil()) {
227	74	Some(m3) => m3,
228	0	None => return Err(CalendarError::ImpossibleResult),
229		};
230		} else {
231	18.0k	a[i] = match Self::from_f64(m) {
232	18.0k	Some(m3) => m3,
233	0	None => return Err(CalendarError::ImpossibleResult),
234		};
235		}
236		}
237		}
238	45.7k	Ok(())
239	45.7k	}
240
241	2	fn search_min(p: impl Fn(Self) -> bool, k: Self) -> Self {
242		//LISTING 1.32 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
243		//Modified to use iteration instead of recursion
244	2	let mut i = k;
245	13	while !p(i) {
246	11	i += Self::one()
247		}
248	2	i
249	2	}
250
251	2	fn search_max(p: impl Fn(Self) -> bool, k: Self) -> Self {
252		//LISTING 1.33 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
253		//Modified to use iteration instead of recursion
254	2	let mut i = k - Self::one();
255	13	while p(i) {
256	11	i += Self::one()
257		}
258	2	i
259	2	}
260		}
261
262		impl TermNum for usize {}
263		impl TermNum for u64 {}
264		impl TermNum for u32 {}
265		impl TermNum for u16 {}
266		impl TermNum for u8 {}
267
268		impl TermNum for i64 {
269	256	fn sign(self) -> Self {
270	256	if self.is_zero() {
271	0	Self::zero()
272		} else {
273	256	self.signum()
274		}
275	256	}
276
277	2.20M	fn modulus(self, other: Self) -> Self {
278	2.20M	debug_assert!(other != Self::zero());
279	2.20M	if other > Self::zero() {
280	2.18M	let x = self;
281	2.18M	let y = other;
282	2.18M	Euclid::rem_euclid(&x, &y)
283		} else {
284	12.7k	let x = -self;
285	12.7k	let y = -other;
286	12.7k	-Euclid::rem_euclid(&x, &y)
287		}
288	2.20M	}
289		}
290
291		macro_rules! CastFn {
292		($n: ident, $t:ident) => {
293	256	fn $n(self) -> Self {
294	256	(self as $t).$n() as Self
295	256	}
296		};
297		($n: ident, $t:ident, $u: ident) => {
298	1.66M	fn $n(self, other: Self) -> $u {
299	0	(self as $t).$n(other as $t) as $u
300	1.66M	}
301		};
302		}
303
304		impl TermNum for i32 {
305		CastFn!(sign, i64);
306		CastFn!(modulus, i64, Self);
307		}
308
309		impl TermNum for i16 {
310		CastFn!(sign, i64);
311		CastFn!(modulus, i64, Self);
312		}
313
314		impl TermNum for i8 {
315		CastFn!(sign, i64);
316		CastFn!(modulus, i64, Self);
317		}
318
319		impl TermNum for f64 {
320	282	fn sign(self) -> Self {
321		//LISTING 1.16 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
322		//Modified to use `signum()`
323	282	if self.is_zero() {
324	2	Self::zero()
325		} else {
326	280	self.signum()
327		}
328	282	}
329
330	255k	fn approx_eq(self, other: Self) -> bool {
331	255k	let x = self;
332	255k	let y = other;
333	255k	if x == y {
334	119k	true
335	135k	} else if x.signum() != y.signum() && x != 0.00 && y != 0.00 {
336	0	false
337		} else {
338	135k	(x - y).abs() < (x.abs() / EQ_SCALE) \|\| (x - y).abs() < EFFECTIVE_EPSILON134k
339		}
340	255k	}
341
342	135k	fn approx_floor(self) -> Self {
343	135k	let x = self;
344	135k	let cx = x.ceil();
345	135k	if x.approx_eq(cx) {
346	82.8k	cx
347		} else {
348	52.9k	x.floor()
349		}
350	135k	}
351
352	0	fn floor_round(self) -> Self {
353		//LISTING 1.15 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
354	0	(self + 0.5).floor()
355	0	}
356
357	157k	fn modulus(self, other: Self) -> Self {
358		//LISTING 1.17 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
359	157k	let x = self;
360	157k	let y = other;
361	157k	debug_assert!(y != 0.0); //Cannot replace with NonZero - it's not supported for f64
362	157k	if x == 0.0 {
363	83.5k	0.0
364		} else {
365	74.0k	x - (y * (x / y).floor())
366		}
367	157k	}
368
369	1.63M	fn is_a_number(self) -> bool {
370	1.63M	!self.is_nan()
371	1.63M	}
372		}
373
374		impl TermNum for f32 {
375		CastFn!(sign, f64);
376		CastFn!(modulus, f64, Self);
377		CastFn!(approx_eq, f64, bool);
378		CastFn!(approx_floor, f64);
379		CastFn!(floor_round, f64);
380
381	0	fn is_a_number(self) -> bool {
382	0	!self.is_nan()
383	0	}
384		}
385
386		// TODO: binary search (listing 1.35)
387		// TODO: inverse f (listing 1.36)
388		// TODO: list_of_fixed_from_moments (listing 1.37)
389		// TODO: range (1.38)
390		// TODO: scan_range (1.39)
391		// TODO: positions_in_range (1.40)
392		// TODO: angles, minutes, degrees
393
394		#[cfg(test)]
395		mod tests {
396		use super::*;
397		use proptest::prop_assume;
398		use proptest::proptest;
399
400		#[test]
401	1	fn modulus_basics() {
402	1	assert_eq!((9.0).modulus(5.0), 4.0);
403	1	assert_eq!((-9.0).modulus(5.0), 1.0);
404	1	assert_eq!((9.0).modulus(-5.0), -1.0);
405	1	assert_eq!((-9.0).modulus(-5.0), -4.0);
406	1	}
407
408		#[test]
409	1	fn adjusted_remainder() {
410	1	assert_eq!((13).adjusted_remainder(12), 1);
411	1	}
412
413		#[test]
414		#[should_panic]
415	1	fn modulus_zero() {
416	1	(123.0).modulus(0.0);
417	1	}
418
419		#[test]
420	1	fn gcd_wikipedia_examples() {
421		//See https://en.wikipedia.org/wiki/Greatest_common_divisor
422	1	assert_eq!(8.0.gcd(12.0), 4.0);
423	1	assert_eq!(54.0.gcd(24.0), 6.0);
424	1	assert_eq!(9.0.gcd(28.0), 1.0); //Coprime
425	1	assert_eq!(24.0.gcd(60.0), 12.0);
426	1	assert_eq!(42.0.gcd(56.0), 14.0);
427	1	}
428
429		#[test]
430	1	fn lcm_wikipedia_examples() {
431		//https://en.wikipedia.org/wiki/Least_common_multiple
432	1	assert_eq!((5.0).lcm(4.0), 20.0);
433	1	assert_eq!((6.0).lcm(4.0), 12.0);
434	1	}
435
436		#[test]
437	1	fn sum_of_2x() {
438	3	let y1 = TermNum::sum1 (\|x\| x2 * 2.0, \|i\| i < 3.0, 1.0);
439	1	assert_eq!(y, 3.0 * 2.0);
440	1	}
441
442		#[test]
443	1	fn product_of_2x() {
444	3	let y1 = TermNum::product1 (\|x\| x2 * 2.0, \|i\| i < 3.0, 1.0);
445	1	assert_eq!(y, 2.0 * 4.0);
446	1	}
447
448		#[test]
449	1	fn search_min_sign() {
450	12	let y1 = f64::search_min1 (\|i\| i.sign() == 1.0, -10.0);
451	1	assert_eq!(y, 1.0);
452	1	let z = f64::search_min(\|i\| i.sign() == 1.0, 500.0);
453	1	assert_eq!(z, 500.0);
454	1	}
455
456		#[test]
457	1	fn search_max_sign() {
458	12	let y1 = f64::search_max1 (\|i\| i.sign() == -1.0, -10.0);
459	1	assert_eq!(y, 0.0);
460	1	let z = f64::search_max(\|i\| i.sign() == -1.0, 500.0);
461	1	assert_eq!(z, 499.0);
462	1	}
463
464		proptest! {
465		#[test]
466		fn mixed_radix_time(ahr in 0..24,amn in 0..59,asc in 0..59) {
467		let ahr = ahr as f64;
468		let amn = amn as f64;
469		let asc = asc as f64;
470		let a = [ahr, amn, asc];
471		let b = [60.0, 60.0];
472		let seconds = TermNum::from_mixed_radix(&a, &b, 2).unwrap();
473		let minutes = TermNum::from_mixed_radix(&a, &b, 1).unwrap();
474		let hours = TermNum::from_mixed_radix(&a, &b, 0).unwrap();
475		let expected_seconds = (ahr * 3600.0) + (amn* 60.0) + asc;
476		let expected_minutes = (ahr * 60.0) + amn + (asc / 60.0);
477		let expected_hours = ahr + (amn / 60.0) + (asc / 3600.0);
478		assert!(seconds.approx_eq(expected_seconds));
479		assert!(minutes.approx_eq(expected_minutes));
480		assert!(hours.approx_eq(expected_hours));
481
482		let mut a_seconds = [0.0, 0.0, 0.0];
483		let mut a_minutes = [0.0, 0.0, 0.0];
484		let mut a_hours = [0.0, 0.0, 0.0];
485		TermNum::to_mixed_radix(seconds, &b, 2, &mut a_seconds).unwrap();
486		TermNum::to_mixed_radix(minutes, &b, 1, &mut a_minutes).unwrap();
487		TermNum::to_mixed_radix(hours, &b, 0, &mut a_hours).unwrap();
488
489		println!("a: {a:?}, a_hours: {a_hours:?}, hours: {hours}");
490
491		assert!(TermNum::approx_eq_iter(a_seconds, a));
492		assert!(TermNum::approx_eq_iter(a_minutes, a));
493		assert!(TermNum::approx_eq_iter(a_hours, a));
494		}
495
496		#[test]
497		fn mixed_radix_time_i(ahr in 0..24,amn in 0..59,asc in 0..59) {
498		let ahr = ahr as i32;
499		let amn = amn as i32;
500		let asc = asc as i32;
501		let a = [ahr, amn, asc];
502		let b = [60, 60];
503		let seconds = TermNum::from_mixed_radix(&a, &b, 2).unwrap();
504		let minutes = TermNum::from_mixed_radix(&a, &b, 1).unwrap();
505		let hours = TermNum::from_mixed_radix(&a, &b, 0).unwrap();
506
507		let ahr = ahr as f64;
508		let amn = amn as f64;
509		let asc = asc as f64;
510		let expected_seconds = (ahr * 3600.0) + (amn* 60.0) + asc;
511		let expected_minutes = (ahr * 60.0) + amn + (asc / 60.0);
512		let expected_hours = ahr + (amn / 60.0) + (asc / 3600.0);
513		assert!(seconds.approx_eq(expected_seconds));
514		assert!(minutes.approx_eq(expected_minutes));
515		assert!(hours.approx_eq(expected_hours));
516
517		let mut a_seconds = [0, 0, 0];
518		let mut a_minutes = [0, 0, 0];
519		let mut a_hours = [0, 0, 0];
520		TermNum::to_mixed_radix(seconds, &b, 2, &mut a_seconds).unwrap();
521		TermNum::to_mixed_radix(minutes, &b, 1, &mut a_minutes).unwrap();
522		TermNum::to_mixed_radix(hours, &b, 0, &mut a_hours).unwrap();
523
524		println!("a: {a:?}, a_hours: {a_hours:?}, hours: {hours}");
525
526		assert_eq!(&a_seconds, &a);
527		assert_eq!(&a_minutes, &a);
528		assert_eq!(&a_hours, &a);
529		}
530
531
532		#[test]
533		fn modulus_positivity(x in EFFECTIVE_MIN..EFFECTIVE_MAX, y in 0.0..EFFECTIVE_MAX) {
534		assert!((x as f64).modulus(y as f64) >= 0.0);
535		}
536
537		#[test]
538		fn modulus_i_positivity(x: i32, y in 1..i32::MAX) {
539		assert!(x.modulus(y) >= 0);
540		}
541
542
543		#[test]
544		fn modulus_negative_x(x in 0.0..EFFECTIVE_MAX, y in 0.0..EFFECTIVE_MAX) {
545		prop_assume!(y != 0.0);
546		prop_assume!(x.modulus(y) != 0.0);
547		let a0 = (-x).modulus(y);
548		let a1 = y - x.modulus(y);
549		assert!(a0.approx_eq(a1));
550		}
551
552		#[test]
553		fn modulus_i_negative_x(x in 0..i32::MAX, y in 1..i32::MAX) {
554		prop_assume!(x.modulus(y) != 0);
555		let a0 = (-x).modulus(y);
556		let a1 = y - x.modulus(y);
557		assert_eq!(a0, a1);
558		}
559
560		#[test]
561		fn modulus_mult(
562		x in -131072.0..131072.0,
563		y in -131072.0..131072.0,
564		z in -131072.0..131072.0) {
565		//LISTING 1.19 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
566		//Using sqrt(EFFECTIVE_MAX) as limit
567		let x = x as f64;
568		let y = y as f64;
569		let z = z as f64;
570		prop_assume!(y != 0.0);
571		prop_assume!(z != 0.0);
572		let a: f64 = x.modulus(y);
573		let az: f64 = (xz).modulus(yz);
574		println!("x={}; y={}; z={}; a={}; az= {}; az= {};", x, y, z, a, az, az);
575		assert!((a * z).approx_eq(az));
576		}
577
578		#[test]
579		fn modulus_i_mult(
580		x in i16::MIN..i16::MAX,
581		y in i16::MIN..i16::MAX,
582		z in i16::MIN..i16::MAX) {
583		//LISTING 1.19 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
584		let x = x as i32;
585		let y = y as i32;
586		let z = z as i32;
587		prop_assume!(y != 0);
588		prop_assume!(z != 0);
589		let a = x.modulus(y);
590		let az = (xz).modulus(yz);
591		println!("x={}; y={}; z={}; a={}; az= {}; az= {};", x, y, z, a, az, az);
592		assert_eq!(a * z, az);
593		}
594
595		#[test]
596		fn modulus_mult_minus_1(x in 0.0..EFFECTIVE_MAX, y in 0.0..EFFECTIVE_MAX) {
597		//LISTING 1.19 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
598		prop_assume!(y != 0.0);
599		let a0 = (-x).modulus(-y);
600		let a1 = -(x.modulus(y));
601		assert_eq!(a0, a1);
602		}
603
604		#[test]
605		fn modulus_i_mult_minus_1(x in 0..i32::MAX, y in 1..i32::MAX) {
606		//LISTING 1.19 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
607		let a0 = (-x).modulus(-y);
608		let a1 = -(x.modulus(y));
609		assert_eq!(a0, a1);
610		}
611
612		#[test]
613		fn modulus_i_multiple_of_y(x: i32, y: i32) {
614		//LISTING 1.20 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
615		prop_assume!(y != 0);
616		let a = (x as i64) - (x.modulus(y) as i64);
617		assert_eq!(a % (y as i64), 0);
618		}
619
620		#[test]
621		fn modulus_bounds(x in EFFECTIVE_MIN..EFFECTIVE_MAX, y in EFFECTIVE_MIN..EFFECTIVE_MAX) {
622		//LISTING 1.21 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
623		prop_assume!(y != 0.0);
624		let a = x.modulus(y) * y.sign();
625		assert!(0.0 <= a && a < y.abs());
626		}
627		#[test]
628		fn modulus_i_bounds(x: i32, y: i32) {
629		//LISTING 1.21 (Calendrical Calculations: The Ultimate Edition by Reingold & Dershowitz.)
630		prop_assume!(y != 0);
631		let a = x.modulus(y) * (y.sign());
632		assert!(0 <= a && a < y.abs());
633		}
634		}
635		}