logicaffeine_base/

hypernatural.rs

//! Infinite integers generated by the program's generator `H` (`dt = 1/H`) — the trip counts of
//! supertasks.
//!
//! A loop that steps `t` by `dt` from `0` while `t < T` runs `⌈T/dt⌉ = ⌈T·H⌉` times: an infinite
//! natural, exactly determined once the residue of `H` is. A [`Hypernatural`] is held as a
//! polynomial in `H` with rational coefficients on each residue class of `H` (a [`Depends`] of
//! polynomials), each one integer-valued on every admitted member of its class — so `⌈H/3⌉` is
//! `H/3`, `(H + 2)/3` or `(H + 1)/3` as `H ≡ 0, 1, 2 (mod 3)`, and `H/2` is an infinite integer
//! once `H` is declared even.
//!
//! Since `H` exceeds every standard natural, two values compare by the leading coefficient of their
//! difference, branch by branch. The residue of a value modulo `k` is periodic in `H`: an
//! integer-valued polynomial on an admitted class `H ≡ c (mod L)` is an integer combination
//! `Σ c_j·C(t, j)` of binomial coefficients in `t = (H − c)/L`. The term `c_j·C(t, j) mod k` reads
//! `C(t, j)` modulo `q_j = k / gcd(c_j, k)`, and by Kummer's theorem `C(t, j) mod p^e` has period
//! `p^(e + ⌊log_p j⌋)`; so `q_j` times the part of `lcm(1, …, j)` made of primes dividing `q_j` is a
//! period, and evaluating one period of representatives gives every residue exactly. `H^64 mod 2`
//! reads two classes and `C(H, 25) mod 2` the thirty-two Lucas's theorem names.

use crate::hyperreal::{check_power_degree, Hyperreal, PowerError};
use crate::limit_value::{Depends, ResidueError, ResidueType, MAX_MODULUS};
use crate::numeric::{BigInt, Rational};
use crate::poly::Poly;
use std::cmp::Ordering;

/// Why the integer part of a hyperreal is not an infinite integer.
#[derive(Clone, Debug, PartialEq, Eq)]
pub enum FloorError {
    /// A power of `1/dt` with an irrational coefficient: `⌊√2/dt⌋` is a Beatty sequence in `H`,
    /// no polynomial.
    Irrational,
    /// The residues the integer part depends on need more classes than one value holds.
    Residue(ResidueError),
}

impl FloorError {
    /// The refusal, with the infinitesimal written `inf`.
    pub fn describe(&self, inf: &str) -> String {
        match self {
            FloorError::Irrational => format!(
                "the integer part of a value with an irrational multiple of a power of 1/{inf} is \
                 no polynomial in 1/{inf}"
            ),
            FloorError::Residue(e) => e.describe(&format!("⌈1/{inf}⌉")),
        }
    }
}

/// The longest period of a polynomial's residues [`integer_on`] reads class by class; past it, the
/// polynomial is evaluated at consecutive members.
const DIRECT_PERIOD: u64 = 1 << 16;

/// Whether `p` is an integer at every `h = first + step·t` for natural `t`. A polynomial with integer
/// coefficients is, at once. Otherwise `p = q/D` with `q` integral and `D` the coefficients' common
/// denominator, and `q(h) mod D` depends only on `h mod D`, so it is read on the
/// `D/gcd(step, D)` residues the members take; when that period passes [`DIRECT_PERIOD`], `p` is
/// evaluated at `deg p + 1` consecutive members, which decides it by finite differences.
fn integer_on(p: &Poly, first: u64, step: u64) -> bool {
    let coeffs = p.coeffs();
    if coeffs.iter().all(Rational::is_integer) {
        return true;
    }
    let d = coeffs.iter().fold(BigInt::from_u64(1), |l, c| {
        let den = c.denominator();
        let g = l.gcd(den);
        l.div_rem(&g).expect("a gcd of positive denominators is positive").0.mul(den)
    });
    let modular = d.to_u64().and_then(|dd| {
        let period = dd / num_gcd(step % dd, dd).max(1);
        (period <= DIRECT_PERIOD).then_some((dd, period))
    });
    let Some((dd, period)) = modular else {
        let degree = p.degree().unwrap_or(0) as u64;
        return (0..=degree).all(|t| {
            let h = BigInt::from_u64(first).add(&BigInt::from_u64(step).mul(&BigInt::from_u64(t)));
            p.eval(&Rational::from_bigint(h)).is_integer()
        });
    };
    let residues: Vec<u128> = coeffs
        .iter()
        .map(|c| {
            let q = c.numerator().mul(&d.div_rem(c.denominator()).expect("D is a multiple of each denominator").0);
            let r = q.div_rem(&d).expect("D is positive").1;
            let r = if r.is_negative() { r.add(&d) } else { r };
            u128::from(r.to_u64().expect("a residue below D"))
        })
        .collect();
    let dd = u128::from(dd);
    (0..period).all(|t| {
        let h = (u128::from(first) + u128::from(step) * u128::from(t)) % dd;
        residues.iter().rev().fold(0u128, |acc, c| (acc * h + c) % dd) == 0
    })
}

fn num_gcd(a: u64, b: u64) -> u64 {
    if b == 0 {
        a
    } else {
        num_gcd(b, a % b)
    }
}

/// `⌊y⌋` for a limited hyperreal: the floor of its standard part, one less when the standard part
/// is an integer approached from below.
fn floor_limited(y: &Hyperreal) -> BigInt {
    let st = y.standard_part().expect("a limited value has a standard part");
    let fl = st.floor();
    let on_integer = st == crate::algebraic::Algebraic::from_rational(Rational::from_bigint(fl.clone()));
    if on_integer && y.sub(&Hyperreal::from_algebraic(st)).signum() < 0 {
        fl.sub(&BigInt::from_u64(1))
    } else {
        fl
    }
}

/// An integer polynomial in the generator, by residue class. See the module documentation.
#[derive(Clone, Debug, PartialEq, Eq, Hash)]
pub struct Hypernatural {
    value: Depends<Poly>,
}

/// `lcm(a, b)` of positive integers.
fn lcm(a: &BigInt, b: &BigInt) -> BigInt {
    a.mul(b).div_rem(&a.gcd(b)).expect("a positive gcd").0
}

/// The part of `n` made of the primes that divide `q`.
fn smooth_part(n: &BigInt, q: &BigInt) -> BigInt {
    let one = BigInt::from_u64(1);
    let mut rest = n.clone();
    loop {
        let g = rest.gcd(q);
        if g == one {
            break;
        }
        rest = rest.div_rem(&g).expect("a positive gcd").0;
    }
    n.div_rem(&rest).expect("a positive remainder").0
}

/// The integer value of `f` at the integer `x` (`f` is integer-valued there).
fn eval_integer(f: &Poly, x: u64) -> BigInt {
    f.eval(&Rational::from_bigint(BigInt::from_u64(x)))
        .to_bigint()
        .expect("a hypernatural branch is integer-valued on its class")
}

/// `p` in descending powers of `generator`: `2·H + 1`, `H^2 + 3·H`, `(1/3)·H + 2/3`.
fn poly_text(p: &Poly, generator: &str) -> String {
    let terms: Vec<(bool, String)> = (0..p.coeffs().len())
        .rev()
        .filter(|&i| !p.coeff(i).is_zero())
        .map(|i| {
            let c = p.coeff(i);
            let negative = c.is_negative();
            let m = c.abs();
            let power = match i {
                0 => String::new(),
                1 => generator.to_string(),
                k => format!("{generator}^{k}"),
            };
            let body = if i == 0 {
                m.to_string()
            } else if m == Rational::one() {
                power
            } else if m.is_integer() {
                format!("{m}·{power}")
            } else {
                format!("({m})·{power}")
            };
            (negative, body)
        })
        .collect();
    join_signed(terms)
}

/// `p(1/dt)` written with the infinitesimal's name: `2/dt + 1`, `1/dt^2`, `1/(3·dt)`.
fn reciprocal_text(p: &Poly, inf: &str) -> String {
    let terms: Vec<(bool, String)> = (0..p.coeffs().len())
        .rev()
        .filter(|&i| !p.coeff(i).is_zero())
        .map(|i| {
            let c = p.coeff(i);
            let negative = c.is_negative();
            let m = c.abs();
            if i == 0 {
                return (negative, m.to_string());
            }
            let power = if i == 1 {
                inf.to_string()
            } else {
                format!("{inf}^{i}")
            };
            let (num, den) = (m.numerator().clone(), m.denominator().clone());
            let body = if den == BigInt::from_i64(1) {
                format!("{num}/{power}")
            } else {
                format!("{num}/({den}·{power})")
            };
            (negative, body)
        })
        .collect();
    join_signed(terms)
}

fn join_signed(terms: Vec<(bool, String)>) -> String {
    if terms.is_empty() {
        return "0".to_string();
    }
    let mut out = String::new();
    for (i, (negative, body)) in terms.into_iter().enumerate() {
        match (i, negative) {
            (0, true) => out.push_str(&format!("-{body}")),
            (0, false) => out.push_str(&body),
            (_, true) => out.push_str(&format!(" - {body}")),
            (_, false) => out.push_str(&format!(" + {body}")),
        }
    }
    out
}

impl Hypernatural {
    /// `H` itself.
    pub fn generator(ty: &ResidueType) -> Hypernatural {
        Hypernatural {
            value: Depends::constant(Poly::x(), ty),
        }
    }

    /// A standard integer.
    pub fn from_integer(n: BigInt, ty: &ResidueType) -> Hypernatural {
        Hypernatural {
            value: Depends::constant(Poly::constant(Rational::from_bigint(n)), ty),
        }
    }

    /// The polynomial `p` in `H` as an infinite integer, when it is an integer on every class the
    /// type admits: a polynomial of degree `d` that is an integer at `d + 1` consecutive integers is
    /// one at every integer, and likewise along one residue class.
    pub fn from_polynomial(p: Poly, ty: &ResidueType) -> Option<Hypernatural> {
        let m = ty.modulus();
        if !integer_on(&p, 0, 1) && !ty.residues().into_iter().all(|r| integer_on(&p, r, m)) {
            return None;
        }
        Some(Hypernatural {
            value: Depends::constant(p, ty),
        })
    }

    /// The value `p` on each class, when each `p` is an integer on the admitted members of its
    /// class.
    pub fn from_branches(value: Depends<Poly>) -> Option<Hypernatural> {
        let candidate = Hypernatural { value };
        let (l, classes) = candidate.admitted_classes();
        for (c, p) in classes {
            if !integer_on(p, c, l) {
                return None;
            }
        }
        Some(candidate)
    }

    /// `⌊x⌋` for a hyperreal `x` under the residue type `ty`. With `x = P(H) + r` — `P` its
    /// principal part with rational coefficients, `r` limited — and `D` the common denominator of
    /// `P`, the fractional part of `P(H)` is `(D·P(H) mod D)/D` on each class, so
    /// `⌊x⌋ = P(H) − frac + ⌊frac + r⌋` there.
    pub fn floor_of(x: &Hyperreal, ty: &ResidueType) -> Result<Hypernatural, FloorError> {
        let (principal, rest) = x.principal_part();
        let mut coefficients = vec![Rational::zero()];
        for a in &principal {
            coefficients.push(a.to_rational().ok_or(FloorError::Irrational)?);
        }
        let p = Poly::from_coeffs(coefficients);
        let denominator = p
            .coeffs()
            .iter()
            .fold(BigInt::from_u64(1), |acc, c| lcm(&acc, c.denominator()));
        let scaled = Hypernatural::from_polynomial(p.scale(&Rational::from_bigint(denominator.clone())), ty)
            .expect("an integer polynomial is an integer everywhere");
        let fractions = scaled.rem(&denominator).map_err(FloorError::Residue)?;
        let value = fractions.map(|f| {
            let frac = Rational::new(f.clone(), denominator.clone()).expect("a positive denominator");
            let carry = floor_limited(&rest.add(&Hyperreal::from_rational(frac.clone())));
            p.sub(&Poly::constant(frac))
                .add(&Poly::constant(Rational::from_bigint(carry)))
        });
        Ok(Hypernatural::from_branches(value).expect("each branch is an integer on its class"))
    }

    /// `⌈x⌉ = −⌊−x⌋`.
    pub fn ceil_of(x: &Hyperreal, ty: &ResidueType) -> Result<Hypernatural, FloorError> {
        Ok(Self::floor_of(&x.neg(), ty)?.neg())
    }

    pub fn depends(&self) -> &Depends<Poly> {
        &self.value
    }

    pub fn add(&self, o: &Hypernatural) -> Result<Hypernatural, ResidueError> {
        Ok(Hypernatural {
            value: self.value.zip_with(&o.value, |a, b| a.add(b))?,
        })
    }

    pub fn sub(&self, o: &Hypernatural) -> Result<Hypernatural, ResidueError> {
        Ok(Hypernatural {
            value: self.value.zip_with(&o.value, |a, b| a.sub(b))?,
        })
    }

    pub fn mul(&self, o: &Hypernatural) -> Result<Hypernatural, ResidueError> {
        Ok(Hypernatural {
            value: self.value.zip_with(&o.value, |a, b| a.mul(b))?,
        })
    }

    /// The highest degree in `H` over every class.
    pub fn degree(&self) -> usize {
        self.value
            .branches()
            .map(|(_, p)| p.degree().unwrap_or(0))
            .max()
            .unwrap_or(0)
    }

    /// `self^k` by repeated squaring, refused when its degree would pass
    /// [`MAX_POWER_DEGREE`](crate::hyperreal::MAX_POWER_DEGREE).
    pub fn pow(&self, k: u64) -> Result<Hypernatural, PowerError> {
        check_power_degree(self.degree(), k)?;
        let square = |a: &Hypernatural, b: &Hypernatural| {
            a.mul(b).expect("a value meets its own residue classes")
        };
        let mut acc = Hypernatural::from_integer(BigInt::from_u64(1), self.value.residue_type());
        let mut base = self.clone();
        let mut e = k;
        while e > 0 {
            if e & 1 == 1 {
                acc = square(&acc, &base);
            }
            e >>= 1;
            if e > 0 {
                base = square(&base, &base);
            }
        }
        Ok(acc)
    }

    pub fn neg(&self) -> Hypernatural {
        Hypernatural {
            value: self.value.map(|a| a.neg()),
        }
    }

    /// Each class modulo `L = lcm(modulus, type modulus)` the type admits, with its branch: the
    /// arithmetic progressions `c + L·t` along which a branch is integer-valued. A branch held at a
    /// coarser modulus than the type's is integral only on the admitted members of its class.
    fn admitted_classes(&self) -> (u64, Vec<(u64, &Poly)>) {
        let l = lcm(
            &BigInt::from_u64(self.value.modulus()),
            &BigInt::from_u64(self.value.residue_type().modulus()),
        )
        .to_u64()
        .expect("a value is held within the bound of its own type");
        let classes = (0..l).filter_map(|c| self.value.at(c).map(|f| (c, f))).collect();
        (l, classes)
    }

    /// A period in `H` of `self mod k`. See the module documentation.
    fn residue_period(&self, k: &BigInt) -> BigInt {
        let one = BigInt::from_u64(1);
        let (l, classes) = self.admitted_classes();
        let mut period = one.clone();
        for (c, f) in classes {
            let degree = f.degree().unwrap_or(0) as u64;
            let mut differences: Vec<BigInt> =
                (0..=degree).map(|t| eval_integer(f, c + l * t)).collect();
            let mut lcm_to_j = one.clone();
            for j in 1..=degree {
                differences = differences.windows(2).map(|w| w[1].sub(&w[0])).collect();
                lcm_to_j = lcm(&lcm_to_j, &BigInt::from_u64(j));
                let q = k.div_rem(&differences[0].gcd(k)).expect("k > 0").0;
                if q != one {
                    period = lcm(&period, &q.mul(&smooth_part(&lcm_to_j, &q)));
                }
            }
        }
        period.mul(&BigInt::from_u64(l))
    }

    /// `self mod k` (in `0..k`) on every residue class of `H`; refused when reading it needs more
    /// than [`MAX_MODULUS`] classes.
    pub fn rem(&self, k: &BigInt) -> Result<Depends<BigInt>, ResidueError> {
        assert!(k.is_positive(), "a modulus is positive");
        let period = self.residue_period(k);
        let ty = self.value.residue_type();
        let classes = lcm(&period, &BigInt::from_u64(ty.modulus()));
        let period = match period.to_u64() {
            Some(p) if classes <= BigInt::from_u64(MAX_MODULUS) => p,
            _ => return Err(ResidueError::TooFine { classes }),
        };
        let m = self.value.modulus();
        let raw = self.value.raw_branches();
        Depends::from_fn(period, ty, move |x| {
            let f = raw[(x % m) as usize]
                .as_ref()
                .expect("an admitted residue has a branch");
            let r = eval_integer(f, x).div_rem(k).expect("k > 0").1;
            if r.is_negative() {
                r.add(k)
            } else {
                r
            }
        })
    }

    /// `⌊t·self⌋` for a rational `t`.
    pub fn scaled_floor(&self, t: &Rational) -> Result<Hypernatural, ResidueError> {
        let p = Rational::from_bigint(t.numerator().clone());
        let q = t.denominator().clone();
        let scaled = Hypernatural {
            value: self.value.map(|f| f.scale(&p)),
        };
        let rem = scaled.rem(&q)?;
        let inv = Rational::from_bigint(q).recip().expect("q > 0");
        Ok(Hypernatural {
            value: scaled.value.zip_with(&rem, |f, r| {
                f.sub(&Poly::constant(Rational::from_bigint(r.clone()))).scale(&inv)
            })?,
        })
    }

    /// `⌈t·self⌉ = −⌊−t·self⌋`.
    pub fn scaled_ceil(&self, t: &Rational) -> Result<Hypernatural, ResidueError> {
        Ok(self.scaled_floor(&t.negated())?.neg())
    }

    /// The order against another value on every class (`H` exceeds every standard natural, so the
    /// leading coefficient of the difference decides).
    pub fn compare(&self, o: &Hypernatural) -> Result<Depends<Ordering>, ResidueError> {
        self.value.zip_with(&o.value, |a, b| {
            let d = a.sub(b);
            match d.leading() {
                None => Ordering::Equal,
                Some(c) if c.is_negative() => Ordering::Less,
                Some(_) => Ordering::Greater,
            }
        })
    }

    pub fn is_infinite(&self) -> Depends<bool> {
        self.value.map(|p| p.degree().unwrap_or(0) > 0)
    }

    pub fn is_negative(&self) -> Depends<bool> {
        self.value.map(|p| p.leading().is_some_and(Rational::is_negative))
    }

    /// The value at the standard point `H = h`, when the type admits it.
    pub fn eval_at(&self, h: u64) -> Option<BigInt> {
        self.value.at(h).map(|f| eval_integer(f, h))
    }

    /// `H ↦ 1/dt`.
    pub fn to_hyperreal(&self) -> Depends<Hyperreal> {
        self.value.map(|p| Hyperreal::from_reciprocal_powers(p.coeffs()))
    }

    /// The value as a polynomial in `generator`, with its dependence on the residue when there is
    /// one.
    pub fn render(&self, generator: &str) -> String {
        self.value.render(generator, |p| poly_text(p, generator))
    }

    /// The value written in the infinitesimal `dt = 1/H`: `2/dt + 1`.
    pub fn render_reciprocal(&self, inf: &str) -> String {
        let generator = format!("⌈1/{inf}⌉");
        self.value.render(&generator, |p| reciprocal_text(p, inf))
    }
}

impl crate::exact_codec::ExactCodec for Hypernatural {
    /// The residue-indexed polynomials in the generator.
    fn write(&self, out: &mut Vec<u8>) {
        self.value.write(out);
    }

    /// Refused unless every branch is integer-valued on the admitted members of its class: a
    /// polynomial of degree `d` that is an integer at `d + 1` consecutive members of an arithmetic
    /// progression is one at every member.
    fn read(buf: &[u8], pos: &mut usize) -> Option<Self> {
        Hypernatural::from_branches(crate::exact_codec::ExactCodec::read(buf, pos)?)
    }
}

#[cfg(test)]
#[path = "../../logicaffeine_tests/unit/logicaffeine_base/hypernatural/tests.rs"]
mod tests;