tokitai-operator 0.1.0

//! p-adic valuation domain and the `Padic` value type.
//!
//! `PadicDomain` is the p-adic precision-bounded domain used by the
//! p-adic valuation-skip planner. The `Padic` value type is the
//! `p`-adic integer with a fixed digit-count precision.
//!
//! Used by the p-adic valuation-skip release-gate
//! (`scripts/check_padic_valuation_skip_theorem.sh`) and by the
//! paper-claim evidence path.
//!
//! See `docs/padic_valuation_skip_theorem.md` and
//! `docs/padic_planner_resource.md` for the policy and the
//! `PadicPlanningResource` type.
//!
use super::{
    AlgebraicStructure, Claim, Contract, ContractId, ContractSet, Evidence, Scope, ValuationBound,
};
use super::{Domain, DomainId, NonArchimedeanDomain, PrecisionModel, ValuedDomain};
use crate::{Error, Result};
use std::collections::BTreeMap;

#[derive(Debug, Clone, PartialEq, Eq, Hash)]
pub enum PadicRepresentation {
    CanonicalResidue,
    Expansion,
    Ball,
    Lazy,
}

#[derive(Debug, Clone, PartialEq, Eq, Hash)]
pub enum PadicErrorModel {
    ExactModulo,
    Ball,
}

#[derive(Debug, Clone, PartialEq, Eq, Hash)]
pub struct PadicMeta {
    pub prime: u64,
    pub precision: u32,
    pub valuation_bound: Option<ValuationBound>,
    pub representation: PadicRepresentation,
    pub error: PadicErrorModel,
}

#[derive(Debug, Clone, PartialEq, Eq)]
pub struct PadicDomain {
    pub meta: PadicMeta,
}

#[derive(Debug, Clone, PartialEq, Eq, Hash)]
pub struct DynamicPadicMeta {
    pub prime: u64,
    pub precision: u32,
}

#[derive(Debug, Clone, PartialEq, Eq)]
pub struct DynamicPadicDomain {
    pub meta: DynamicPadicMeta,
}

#[derive(Debug, Clone, PartialEq, Eq, Hash)]
pub struct DynamicPadic {
    pub meta: DynamicPadicMeta,
    digits: Vec<u64>,
}

#[derive(Debug, Clone, PartialEq, Eq)]
pub struct DynamicPadicMatrix {
    pub domain: DynamicPadicDomain,
    pub rows: usize,
    pub cols: usize,
    pub data: Vec<DynamicPadic>,
}

#[derive(Debug, Clone, PartialEq, Eq)]
pub struct ValuationSkipReport {
    pub result: Padic,
    pub skipped_terms: usize,
    pub evaluated_terms: usize,
}

#[derive(Debug, Clone, PartialEq, Eq)]
pub struct PadicMatrix {
    pub domain: PadicDomain,
    pub rows: usize,
    pub cols: usize,
    pub data: Vec<Padic>,
}

#[derive(Debug, Clone, PartialEq, Eq)]
pub struct PadicMatrixMetadata {
    pub prime: u64,
    pub precision: u32,
    pub rows: usize,
    pub cols: usize,
    pub residues: Vec<u128>,
    pub valuations: Vec<u32>,
    pub unit_residues: Vec<u128>,
    pub valuation_histogram: BTreeMap<u32, usize>,
    pub bucket_fingerprint: String,
}

#[derive(Debug, Clone, PartialEq, Eq)]
pub struct PadicOutputCertificate {
    pub row: usize,
    pub col: usize,
    pub evaluated_product_count: usize,
    pub skipped_product_count: usize,
    pub min_skipped_valuation: Option<u32>,
    pub precision_cutoff: u32,
    pub precision_safety_margin: Option<u32>,
}

#[derive(Debug, Clone, PartialEq, Eq)]
pub struct PadicCertifiedMatmulReport {
    pub output: PadicMatrix,
    pub dense_oracle: PadicMatrix,
    pub lhs_metadata: PadicMatrixMetadata,
    pub rhs_metadata: PadicMatrixMetadata,
    pub output_certificates: Vec<PadicOutputCertificate>,
    pub skipped_products: usize,
    pub evaluated_products: usize,
    pub dense_oracle_matches: bool,
}

impl PadicDomain {
    pub fn new(prime: u64, precision: u32) -> Result<Self> {
        if prime < 2 {
            return Err(Error::domain("p-adic prime must be at least 2"));
        }
        if !is_prime(prime) {
            return Err(Error::domain(format!("{prime} is not prime")));
        }
        let _ = checked_modulus(prime, precision)?;
        Ok(Self {
            meta: PadicMeta {
                prime,
                precision,
                valuation_bound: None,
                representation: PadicRepresentation::CanonicalResidue,
                error: PadicErrorModel::ExactModulo,
            },
        })
    }

    pub fn modulus(&self) -> u128 {
        checked_modulus(self.meta.prime, self.meta.precision)
            .expect("PadicDomain::new validates modulus")
    }

    pub fn element(&self, value: u128) -> Padic {
        Padic {
            meta: self.meta.clone(),
            residue: value % self.modulus(),
        }
    }

    pub fn zero(&self) -> Padic {
        self.element(0)
    }

    pub fn one(&self) -> Padic {
        self.element(1)
    }

    pub fn add(&self, lhs: &Padic, rhs: &Padic) -> Result<Padic> {
        self.ensure_same_domain(lhs)?;
        self.ensure_same_domain(rhs)?;
        Ok(self.element((lhs.residue + rhs.residue) % self.modulus()))
    }

    pub fn sub(&self, lhs: &Padic, rhs: &Padic) -> Result<Padic> {
        self.ensure_same_domain(lhs)?;
        self.ensure_same_domain(rhs)?;
        let modulus = self.modulus();
        Ok(self.element((lhs.residue + modulus - rhs.residue) % modulus))
    }

    pub fn mul(&self, lhs: &Padic, rhs: &Padic) -> Result<Padic> {
        self.ensure_same_domain(lhs)?;
        self.ensure_same_domain(rhs)?;
        Ok(self.element((lhs.residue * rhs.residue) % self.modulus()))
    }

    pub fn neg(&self, value: &Padic) -> Result<Padic> {
        self.ensure_same_domain(value)?;
        if value.residue == 0 {
            Ok(value.clone())
        } else {
            Ok(self.element(self.modulus() - value.residue))
        }
    }

    pub fn valuation_of(&self, value: &Padic) -> Result<Option<u32>> {
        self.ensure_same_domain(value)?;
        Ok(valuation_residue(
            value.residue,
            self.meta.prime,
            self.meta.precision,
        ))
    }

    pub fn is_unit(&self, value: &Padic) -> Result<bool> {
        self.ensure_same_domain(value)?;
        Ok(value.residue % self.meta.prime as u128 != 0)
    }

    pub fn inverse(&self, value: &Padic) -> Result<Padic> {
        self.ensure_same_domain(value)?;
        if !self.is_unit(value)? {
            return Err(Error::domain(
                "only p-adic units are invertible in fixed residue representation",
            ));
        }
        let modulus = self.modulus() as i128;
        let residue = value.residue as i128;
        let inverse = mod_inverse(residue, modulus)
            .ok_or_else(|| Error::domain("failed to compute modular inverse"))?;
        Ok(self.element(inverse as u128))
    }

    pub fn div(&self, lhs: &Padic, rhs: &Padic) -> Result<Padic> {
        let inv = self.inverse(rhs)?;
        self.mul(lhs, &inv)
    }

    pub fn truncate(&self, value: &Padic, precision: u32) -> Result<Padic> {
        self.ensure_same_domain(value)?;
        if precision > self.meta.precision {
            return Err(Error::domain(format!(
                "cannot truncate from precision {} to higher precision {precision}",
                self.meta.precision
            )));
        }
        let truncated_domain = Self::new(self.meta.prime, precision)?;
        Ok(truncated_domain.element(value.residue))
    }

    /// Reduce a p-adic value to a smaller target precision by taking its
    /// residue modulo p^target_precision. The result is the canonical
    /// element(target_residue) on the lower-precision domain.
    pub fn reduce_to_precision(&self, value: &Padic, target_precision: u32) -> Result<Padic> {
        self.ensure_same_domain(value)?;
        if target_precision > self.meta.precision {
            return Err(Error::domain(format!(
                "cannot reduce from precision {} to higher precision {target_precision}",
                self.meta.precision
            )));
        }
        let target_modulus = checked_modulus(self.meta.prime, target_precision)
            .expect("target precision is already validated");
        let reduced_residue = value.residue % target_modulus;
        let target_domain = Self::new(self.meta.prime, target_precision)?;
        Ok(target_domain.element(reduced_residue))
    }

    pub fn sum_products(&self, lhs: &[Padic], rhs: &[Padic]) -> Result<Padic> {
        if lhs.len() != rhs.len() {
            return Err(Error::domain(format!(
                "sum_products length mismatch: left={}, right={}",
                lhs.len(),
                rhs.len()
            )));
        }
        let mut acc = self.zero();
        for (a, b) in lhs.iter().zip(rhs.iter()) {
            let product = self.mul(a, b)?;
            acc = self.add(&acc, &product)?;
        }
        Ok(acc)
    }

    pub fn vector_add(&self, lhs: &[Padic], rhs: &[Padic]) -> Result<Vec<Padic>> {
        if lhs.len() != rhs.len() {
            return Err(Error::domain(format!(
                "p-adic vector_add length mismatch: left={}, right={}",
                lhs.len(),
                rhs.len()
            )));
        }
        lhs.iter()
            .zip(rhs.iter())
            .map(|(left, right)| self.add(left, right))
            .collect()
    }

    pub fn vector_hadamard(&self, lhs: &[Padic], rhs: &[Padic]) -> Result<Vec<Padic>> {
        if lhs.len() != rhs.len() {
            return Err(Error::domain(format!(
                "p-adic vector_hadamard length mismatch: left={}, right={}",
                lhs.len(),
                rhs.len()
            )));
        }
        lhs.iter()
            .zip(rhs.iter())
            .map(|(left, right)| self.mul(left, right))
            .collect()
    }

    pub fn dot_product(&self, lhs: &[Padic], rhs: &[Padic]) -> Result<Padic> {
        self.sum_products(lhs, rhs)
    }

    pub fn matrix_vector_mul(&self, matrix: &PadicMatrix, vector: &[Padic]) -> Result<Vec<Padic>> {
        matrix.ensure_domain_and_shape(self, "matrix")?;
        if matrix.cols != vector.len() {
            return Err(Error::domain(format!(
                "p-adic matrix_vector_mul dimension mismatch: matrix cols={}, vector len={}",
                matrix.cols,
                vector.len()
            )));
        }
        for value in vector {
            self.ensure_same_domain(value)?;
        }
        let mut output = Vec::with_capacity(matrix.rows);
        for row in 0..matrix.rows {
            let start = row * matrix.cols;
            let end = start + matrix.cols;
            output.push(self.dot_product(&matrix.data[start..end], vector)?);
        }
        Ok(output)
    }

    pub fn sum_products_with_valuation_skip(
        &self,
        lhs: &[Padic],
        rhs: &[Padic],
    ) -> Result<ValuationSkipReport> {
        if lhs.len() != rhs.len() {
            return Err(Error::domain(format!(
                "sum_products length mismatch: left={}, right={}",
                lhs.len(),
                rhs.len()
            )));
        }
        let mut acc = self.zero();
        let mut skipped_terms = 0;
        let mut evaluated_terms = 0;
        for (a, b) in lhs.iter().zip(rhs.iter()) {
            let va = self.valuation_of(a)?.unwrap_or(self.meta.precision);
            let vb = self.valuation_of(b)?.unwrap_or(self.meta.precision);
            if va.saturating_add(vb) >= self.meta.precision {
                skipped_terms += 1;
                continue;
            }
            let product = self.mul(a, b)?;
            acc = self.add(&acc, &product)?;
            evaluated_terms += 1;
        }
        Ok(ValuationSkipReport {
            result: acc,
            skipped_terms,
            evaluated_terms,
        })
    }

    pub fn matrix(&self, rows: usize, cols: usize, data: Vec<Padic>) -> Result<PadicMatrix> {
        PadicMatrix::new(self.clone(), rows, cols, data)
    }

    pub fn matrix_metadata(
        &self,
        rows: usize,
        cols: usize,
        data: &[Padic],
    ) -> Result<PadicMatrixMetadata> {
        if data.len() != rows.saturating_mul(cols) {
            return Err(Error::domain(format!(
                "p-adic matrix data length mismatch: expected {}, got {}",
                rows.saturating_mul(cols),
                data.len()
            )));
        }

        let mut residues = Vec::with_capacity(data.len());
        let mut valuations = Vec::with_capacity(data.len());
        let mut unit_residues = Vec::with_capacity(data.len());
        let mut valuation_histogram = BTreeMap::new();
        for value in data {
            self.ensure_same_domain(value)?;
            let residue = value.residue % self.modulus();
            let valuation = self.valuation_of(value)?.unwrap_or(self.meta.precision);
            let unit_residue = unit_residue_at_valuation(
                residue,
                valuation,
                self.meta.prime,
                self.meta.precision,
            )?;
            residues.push(residue);
            valuations.push(valuation);
            unit_residues.push(unit_residue);
            *valuation_histogram.entry(valuation).or_insert(0) += 1;
        }
        let bucket_fingerprint = valuation_bucket_fingerprint(
            self.meta.prime,
            self.meta.precision,
            rows,
            cols,
            &valuations,
            &valuation_histogram,
        );
        Ok(PadicMatrixMetadata {
            prime: self.meta.prime,
            precision: self.meta.precision,
            rows,
            cols,
            residues,
            valuations,
            unit_residues,
            valuation_histogram,
            bucket_fingerprint,
        })
    }

    pub fn dense_matrix_mul(&self, lhs: &PadicMatrix, rhs: &PadicMatrix) -> Result<PadicMatrix> {
        lhs.ensure_domain_and_shape(self, "left")?;
        rhs.ensure_domain_and_shape(self, "right")?;
        if lhs.cols != rhs.rows {
            return Err(Error::domain(format!(
                "p-adic matrix inner dimension mismatch: left={}, right={}",
                lhs.cols, rhs.rows
            )));
        }

        let mut data = Vec::with_capacity(lhs.rows * rhs.cols);
        for row in 0..lhs.rows {
            for col in 0..rhs.cols {
                let mut acc = self.zero();
                for inner in 0..lhs.cols {
                    let product = self.mul(lhs.get(row, inner)?, rhs.get(inner, col)?)?;
                    acc = self.add(&acc, &product)?;
                }
                data.push(acc);
            }
        }
        self.matrix(lhs.rows, rhs.cols, data)
    }

    pub fn certified_valuation_sparse_matrix_mul(
        &self,
        lhs: &PadicMatrix,
        rhs: &PadicMatrix,
    ) -> Result<PadicCertifiedMatmulReport> {
        lhs.ensure_domain_and_shape(self, "left")?;
        rhs.ensure_domain_and_shape(self, "right")?;
        if lhs.cols != rhs.rows {
            return Err(Error::domain(format!(
                "p-adic matrix inner dimension mismatch: left={}, right={}",
                lhs.cols, rhs.rows
            )));
        }
        let lhs_metadata = lhs.metadata()?;
        let rhs_metadata = rhs.metadata()?;
        let dense_oracle = self.dense_matrix_mul(lhs, rhs)?;

        let mut data = Vec::with_capacity(lhs.rows * rhs.cols);
        let mut output_certificates = Vec::with_capacity(lhs.rows * rhs.cols);
        let mut skipped_products = 0;
        let mut evaluated_products = 0;
        for row in 0..lhs.rows {
            for col in 0..rhs.cols {
                let mut acc = self.zero();
                let mut skipped_product_count = 0;
                let mut evaluated_product_count = 0;
                let mut min_skipped_valuation: Option<u32> = None;
                for inner in 0..lhs.cols {
                    let lhs_index = row * lhs.cols + inner;
                    let rhs_index = inner * rhs.cols + col;
                    let lhs_valuation = lhs_metadata.valuations[lhs_index];
                    let rhs_valuation = rhs_metadata.valuations[rhs_index];
                    let product_valuation_bound = lhs_valuation.saturating_add(rhs_valuation);
                    if product_valuation_bound >= self.meta.precision {
                        skipped_product_count += 1;
                        skipped_products += 1;
                        min_skipped_valuation = Some(
                            min_skipped_valuation
                                .map(|current| current.min(product_valuation_bound))
                                .unwrap_or(product_valuation_bound),
                        );
                        continue;
                    }

                    let product = self.mul(lhs.get(row, inner)?, rhs.get(inner, col)?)?;
                    acc = self.add(&acc, &product)?;
                    evaluated_product_count += 1;
                    evaluated_products += 1;
                }
                let precision_safety_margin = min_skipped_valuation
                    .map(|valuation| valuation.saturating_sub(self.meta.precision));
                output_certificates.push(PadicOutputCertificate {
                    row,
                    col,
                    evaluated_product_count,
                    skipped_product_count,
                    min_skipped_valuation,
                    precision_cutoff: self.meta.precision,
                    precision_safety_margin,
                });
                data.push(acc);
            }
        }

        let output = self.matrix(lhs.rows, rhs.cols, data)?;
        let dense_oracle_matches = output
            .data
            .iter()
            .zip(dense_oracle.data.iter())
            .all(|(candidate, reference)| candidate.residue == reference.residue);
        if !dense_oracle_matches {
            return Err(Error::verification(
                "certified valuation-sparse p-adic matmul disagrees with dense CPU oracle",
            ));
        }
        Ok(PadicCertifiedMatmulReport {
            output,
            dense_oracle,
            lhs_metadata,
            rhs_metadata,
            output_certificates,
            skipped_products,
            evaluated_products,
            dense_oracle_matches,
        })
    }

    fn ensure_same_domain(&self, value: &Padic) -> Result<()> {
        if value.meta.prime == self.meta.prime && value.meta.precision == self.meta.precision {
            Ok(())
        } else {
            Err(Error::domain(format!(
                "p-adic domain mismatch: expected p={} precision={}, got p={} precision={}",
                self.meta.prime, self.meta.precision, value.meta.prime, value.meta.precision
            )))
        }
    }
}

impl DynamicPadicDomain {
    pub fn new(prime: u64, precision: u32) -> Result<Self> {
        if prime < 2 {
            return Err(Error::domain("p-adic prime must be at least 2"));
        }
        if !is_prime(prime) {
            return Err(Error::domain(format!("{prime} is not prime")));
        }
        Ok(Self {
            meta: DynamicPadicMeta { prime, precision },
        })
    }

    pub fn from_u128(&self, mut value: u128) -> DynamicPadic {
        let mut digits = Vec::with_capacity(self.meta.precision as usize);
        let prime = self.meta.prime as u128;
        for _ in 0..self.meta.precision {
            digits.push((value % prime) as u64);
            value /= prime;
        }
        DynamicPadic {
            meta: self.meta.clone(),
            digits,
        }
    }

    pub fn from_digits(&self, digits: &[u64]) -> Result<DynamicPadic> {
        let precision = self.meta.precision as usize;
        let mut canonical = Vec::with_capacity(precision);
        for digit in digits.iter().take(precision) {
            if *digit >= self.meta.prime {
                return Err(Error::domain(format!(
                    "dynamic p-adic digit {digit} is not canonical for p={}",
                    self.meta.prime
                )));
            }
            canonical.push(*digit);
        }
        canonical.resize(precision, 0);
        Ok(DynamicPadic {
            meta: self.meta.clone(),
            digits: canonical,
        })
    }

    pub fn zero(&self) -> DynamicPadic {
        self.from_u128(0)
    }

    pub fn one(&self) -> DynamicPadic {
        self.from_u128(1)
    }

    pub fn add(&self, lhs: &DynamicPadic, rhs: &DynamicPadic) -> Result<DynamicPadic> {
        self.ensure_same_dynamic_domain(lhs)?;
        self.ensure_same_dynamic_domain(rhs)?;
        let prime = self.meta.prime as u128;
        let mut carry = 0u128;
        let mut digits = Vec::with_capacity(self.meta.precision as usize);
        for (left, right) in lhs.digits.iter().zip(rhs.digits.iter()) {
            let total = *left as u128 + *right as u128 + carry;
            digits.push((total % prime) as u64);
            carry = total / prime;
        }
        Ok(DynamicPadic {
            meta: self.meta.clone(),
            digits,
        })
    }

    pub fn sub(&self, lhs: &DynamicPadic, rhs: &DynamicPadic) -> Result<DynamicPadic> {
        self.ensure_same_dynamic_domain(lhs)?;
        self.ensure_same_dynamic_domain(rhs)?;
        let prime = self.meta.prime as i128;
        let mut borrow = 0i128;
        let mut digits = Vec::with_capacity(self.meta.precision as usize);
        for (left, right) in lhs.digits.iter().zip(rhs.digits.iter()) {
            let mut diff = *left as i128 - *right as i128 - borrow;
            if diff < 0 {
                diff += prime;
                borrow = 1;
            } else {
                borrow = 0;
            }
            digits.push(diff as u64);
        }
        Ok(DynamicPadic {
            meta: self.meta.clone(),
            digits,
        })
    }

    pub fn mul(&self, lhs: &DynamicPadic, rhs: &DynamicPadic) -> Result<DynamicPadic> {
        self.ensure_same_dynamic_domain(lhs)?;
        self.ensure_same_dynamic_domain(rhs)?;
        let precision = self.meta.precision as usize;
        let prime = self.meta.prime as u128;
        let mut accum = vec![0u128; precision];
        for (left_index, left_digit) in lhs.digits.iter().enumerate() {
            for (right_index, right_digit) in rhs.digits.iter().enumerate() {
                let output_index = left_index + right_index;
                if output_index >= precision {
                    break;
                }
                let product = (*left_digit as u128)
                    .checked_mul(*right_digit as u128)
                    .ok_or_else(|| Error::domain("dynamic p-adic multiplication overflow"))?;
                accum[output_index] = accum[output_index]
                    .checked_add(product)
                    .ok_or_else(|| Error::domain("dynamic p-adic multiplication overflow"))?;
            }
        }

        let mut carry = 0u128;
        let mut digits = Vec::with_capacity(precision);
        for value in accum {
            let total = value
                .checked_add(carry)
                .ok_or_else(|| Error::domain("dynamic p-adic carry overflow"))?;
            digits.push((total % prime) as u64);
            carry = total / prime;
        }
        Ok(DynamicPadic {
            meta: self.meta.clone(),
            digits,
        })
    }

    pub fn neg(&self, value: &DynamicPadic) -> Result<DynamicPadic> {
        self.sub(&self.zero(), value)
    }

    pub fn valuation_of(&self, value: &DynamicPadic) -> Result<Option<u32>> {
        self.ensure_same_dynamic_domain(value)?;
        Ok(value
            .digits
            .iter()
            .position(|digit| *digit != 0)
            .map(|index| index as u32))
    }

    pub fn is_unit(&self, value: &DynamicPadic) -> Result<bool> {
        self.ensure_same_dynamic_domain(value)?;
        Ok(value.digits.first().copied().unwrap_or(0) != 0)
    }

    pub fn inverse(&self, value: &DynamicPadic) -> Result<DynamicPadic> {
        self.ensure_same_dynamic_domain(value)?;
        if !self.is_unit(value)? {
            return Err(Error::domain(
                "only dynamic p-adic units are invertible in bounded digit representation",
            ));
        }

        let first_digit = value.digits[0] as i128;
        let prime = self.meta.prime as i128;
        let first_inverse = mod_inverse(first_digit, prime)
            .ok_or_else(|| Error::domain("failed to compute dynamic p-adic unit inverse"))?;
        let mut inverse = self.from_digits(&[first_inverse as u64])?;
        let two = self.from_u128(2);
        let mut correct_digits = 1u32;
        while correct_digits < self.meta.precision {
            let value_times_inverse = self.mul(value, &inverse)?;
            let correction = self.sub(&two, &value_times_inverse)?;
            inverse = self.mul(&inverse, &correction)?;
            correct_digits = correct_digits.saturating_mul(2);
        }
        Ok(inverse)
    }

    pub fn div(&self, lhs: &DynamicPadic, rhs: &DynamicPadic) -> Result<DynamicPadic> {
        let inverse = self.inverse(rhs)?;
        self.mul(lhs, &inverse)
    }

    pub fn truncate(&self, value: &DynamicPadic, precision: u32) -> Result<DynamicPadic> {
        self.ensure_same_dynamic_domain(value)?;
        if precision > self.meta.precision {
            return Err(Error::domain(format!(
                "cannot truncate dynamic p-adic from precision {} to higher precision {precision}",
                self.meta.precision
            )));
        }
        let truncated_domain = Self::new(self.meta.prime, precision)?;
        truncated_domain.from_digits(&value.digits[..precision as usize])
    }

    pub fn equal_at_precision(
        &self,
        lhs: &DynamicPadic,
        rhs: &DynamicPadic,
        precision: u32,
    ) -> Result<bool> {
        self.ensure_same_dynamic_domain(lhs)?;
        self.ensure_same_dynamic_domain(rhs)?;
        if precision > self.meta.precision {
            return Err(Error::domain(format!(
                "cannot compare dynamic p-adics at precision {precision}; domain precision is {}",
                self.meta.precision
            )));
        }
        Ok(lhs.digits[..precision as usize] == rhs.digits[..precision as usize])
    }

    pub fn sum_products(&self, lhs: &[DynamicPadic], rhs: &[DynamicPadic]) -> Result<DynamicPadic> {
        if lhs.len() != rhs.len() {
            return Err(Error::domain(format!(
                "dynamic p-adic sum_products length mismatch: left={}, right={}",
                lhs.len(),
                rhs.len()
            )));
        }
        let mut acc = self.zero();
        for (left, right) in lhs.iter().zip(rhs.iter()) {
            let product = self.mul(left, right)?;
            acc = self.add(&acc, &product)?;
        }
        Ok(acc)
    }

    pub fn dot_product(&self, lhs: &[DynamicPadic], rhs: &[DynamicPadic]) -> Result<DynamicPadic> {
        self.sum_products(lhs, rhs)
    }

    pub fn vector_add(
        &self,
        lhs: &[DynamicPadic],
        rhs: &[DynamicPadic],
    ) -> Result<Vec<DynamicPadic>> {
        if lhs.len() != rhs.len() {
            return Err(Error::domain(format!(
                "dynamic p-adic vector_add length mismatch: left={}, right={}",
                lhs.len(),
                rhs.len()
            )));
        }
        lhs.iter()
            .zip(rhs.iter())
            .map(|(left, right)| self.add(left, right))
            .collect()
    }

    pub fn vector_hadamard(
        &self,
        lhs: &[DynamicPadic],
        rhs: &[DynamicPadic],
    ) -> Result<Vec<DynamicPadic>> {
        if lhs.len() != rhs.len() {
            return Err(Error::domain(format!(
                "dynamic p-adic vector_hadamard length mismatch: left={}, right={}",
                lhs.len(),
                rhs.len()
            )));
        }
        lhs.iter()
            .zip(rhs.iter())
            .map(|(left, right)| self.mul(left, right))
            .collect()
    }

    pub fn matrix(
        &self,
        rows: usize,
        cols: usize,
        data: Vec<DynamicPadic>,
    ) -> Result<DynamicPadicMatrix> {
        DynamicPadicMatrix::new(self.clone(), rows, cols, data)
    }

    pub fn identity_matrix(&self, size: usize) -> Result<DynamicPadicMatrix> {
        let mut data = Vec::with_capacity(size.saturating_mul(size));
        for row in 0..size {
            for col in 0..size {
                if row == col {
                    data.push(self.one());
                } else {
                    data.push(self.zero());
                }
            }
        }
        self.matrix(size, size, data)
    }

    pub fn matrix_add(
        &self,
        lhs: &DynamicPadicMatrix,
        rhs: &DynamicPadicMatrix,
    ) -> Result<DynamicPadicMatrix> {
        lhs.ensure_domain_and_shape(self, "left")?;
        rhs.ensure_domain_and_shape(self, "right")?;
        ensure_dynamic_matrix_same_shape(lhs, rhs, "matrix_add")?;
        let data = lhs
            .data
            .iter()
            .zip(rhs.data.iter())
            .map(|(left, right)| self.add(left, right))
            .collect::<Result<Vec<_>>>()?;
        self.matrix(lhs.rows, lhs.cols, data)
    }

    pub fn matrix_hadamard(
        &self,
        lhs: &DynamicPadicMatrix,
        rhs: &DynamicPadicMatrix,
    ) -> Result<DynamicPadicMatrix> {
        lhs.ensure_domain_and_shape(self, "left")?;
        rhs.ensure_domain_and_shape(self, "right")?;
        ensure_dynamic_matrix_same_shape(lhs, rhs, "matrix_hadamard")?;
        let data = lhs
            .data
            .iter()
            .zip(rhs.data.iter())
            .map(|(left, right)| self.mul(left, right))
            .collect::<Result<Vec<_>>>()?;
        self.matrix(lhs.rows, lhs.cols, data)
    }

    pub fn transpose(&self, matrix: &DynamicPadicMatrix) -> Result<DynamicPadicMatrix> {
        matrix.ensure_domain_and_shape(self, "matrix")?;
        let mut data = Vec::with_capacity(matrix.data.len());
        for col in 0..matrix.cols {
            for row in 0..matrix.rows {
                data.push(matrix.get(row, col)?.clone());
            }
        }
        self.matrix(matrix.cols, matrix.rows, data)
    }

    pub fn matrix_vector_mul(
        &self,
        matrix: &DynamicPadicMatrix,
        vector: &[DynamicPadic],
    ) -> Result<Vec<DynamicPadic>> {
        matrix.ensure_domain_and_shape(self, "matrix")?;
        if matrix.cols != vector.len() {
            return Err(Error::domain(format!(
                "dynamic p-adic matrix_vector_mul dimension mismatch: matrix cols={}, vector len={}",
                matrix.cols,
                vector.len()
            )));
        }
        for value in vector {
            self.ensure_same_dynamic_domain(value)?;
        }
        let mut output = Vec::with_capacity(matrix.rows);
        for row in 0..matrix.rows {
            let start = row * matrix.cols;
            let end = start + matrix.cols;
            output.push(self.dot_product(&matrix.data[start..end], vector)?);
        }
        Ok(output)
    }

    pub fn dense_matrix_mul(
        &self,
        lhs: &DynamicPadicMatrix,
        rhs: &DynamicPadicMatrix,
    ) -> Result<DynamicPadicMatrix> {
        lhs.ensure_domain_and_shape(self, "left")?;
        rhs.ensure_domain_and_shape(self, "right")?;
        if lhs.cols != rhs.rows {
            return Err(Error::domain(format!(
                "dynamic p-adic matrix inner dimension mismatch: left={}, right={}",
                lhs.cols, rhs.rows
            )));
        }

        let mut data = Vec::with_capacity(lhs.rows * rhs.cols);
        for row in 0..lhs.rows {
            for col in 0..rhs.cols {
                let mut acc = self.zero();
                for inner in 0..lhs.cols {
                    let product = self.mul(lhs.get(row, inner)?, rhs.get(inner, col)?)?;
                    acc = self.add(&acc, &product)?;
                }
                data.push(acc);
            }
        }
        self.matrix(lhs.rows, rhs.cols, data)
    }

    pub fn to_fixed_domain(&self) -> Result<PadicDomain> {
        PadicDomain::new(self.meta.prime, self.meta.precision)
    }

    fn ensure_same_dynamic_domain(&self, value: &DynamicPadic) -> Result<()> {
        if value.meta == self.meta && value.digits.len() == self.meta.precision as usize {
            Ok(())
        } else {
            Err(Error::domain(format!(
                "dynamic p-adic domain mismatch: expected p={} precision={}, got p={} precision={}",
                self.meta.prime, self.meta.precision, value.meta.prime, value.meta.precision
            )))
        }
    }
}

impl DynamicPadic {
    pub fn digits(&self) -> &[u64] {
        &self.digits
    }

    pub fn to_u128_checked(&self) -> Result<u128> {
        let mut value = 0u128;
        let mut place = 1u128;
        for digit in &self.digits {
            let term = place
                .checked_mul(*digit as u128)
                .ok_or_else(|| Error::domain("dynamic p-adic value does not fit in u128"))?;
            value = value
                .checked_add(term)
                .ok_or_else(|| Error::domain("dynamic p-adic value does not fit in u128"))?;
            place = place
                .checked_mul(self.meta.prime as u128)
                .ok_or_else(|| Error::domain("dynamic p-adic modulus does not fit in u128"))?;
        }
        Ok(value)
    }
}

impl DynamicPadicMatrix {
    pub fn new(
        domain: DynamicPadicDomain,
        rows: usize,
        cols: usize,
        data: Vec<DynamicPadic>,
    ) -> Result<Self> {
        if data.len() != rows.saturating_mul(cols) {
            return Err(Error::domain(format!(
                "dynamic p-adic matrix data length mismatch: expected {}, got {}",
                rows.saturating_mul(cols),
                data.len()
            )));
        }
        for value in &data {
            domain.ensure_same_dynamic_domain(value)?;
        }
        Ok(Self {
            domain,
            rows,
            cols,
            data,
        })
    }

    pub fn get(&self, row: usize, col: usize) -> Result<&DynamicPadic> {
        if row >= self.rows || col >= self.cols {
            return Err(Error::domain(format!(
                "dynamic p-adic matrix index out of bounds: row={row}, col={col}, shape={}x{}",
                self.rows, self.cols
            )));
        }
        Ok(&self.data[row * self.cols + col])
    }

    fn ensure_domain_and_shape(&self, domain: &DynamicPadicDomain, side: &str) -> Result<()> {
        if self.domain.meta != domain.meta {
            return Err(Error::domain(format!(
                "{side} dynamic p-adic matrix domain mismatch: expected p={} precision={}, got p={} precision={}",
                domain.meta.prime,
                domain.meta.precision,
                self.domain.meta.prime,
                self.domain.meta.precision
            )));
        }
        if self.data.len() != self.rows.saturating_mul(self.cols) {
            return Err(Error::domain(format!(
                "{side} dynamic p-adic matrix shape/data mismatch: shape={}x{}, data={}",
                self.rows,
                self.cols,
                self.data.len()
            )));
        }
        Ok(())
    }
}

fn ensure_dynamic_matrix_same_shape(
    lhs: &DynamicPadicMatrix,
    rhs: &DynamicPadicMatrix,
    op: &str,
) -> Result<()> {
    if lhs.rows != rhs.rows || lhs.cols != rhs.cols {
        return Err(Error::domain(format!(
            "dynamic p-adic {op} shape mismatch: left={}x{}, right={}x{}",
            lhs.rows, lhs.cols, rhs.rows, rhs.cols
        )));
    }
    Ok(())
}

impl PadicMatrix {
    pub fn new(domain: PadicDomain, rows: usize, cols: usize, data: Vec<Padic>) -> Result<Self> {
        if data.len() != rows.saturating_mul(cols) {
            return Err(Error::domain(format!(
                "p-adic matrix data length mismatch: expected {}, got {}",
                rows.saturating_mul(cols),
                data.len()
            )));
        }
        for value in &data {
            domain.ensure_same_domain(value)?;
        }
        Ok(Self {
            domain,
            rows,
            cols,
            data,
        })
    }

    pub fn metadata(&self) -> Result<PadicMatrixMetadata> {
        self.domain
            .matrix_metadata(self.rows, self.cols, &self.data)
    }

    pub fn get(&self, row: usize, col: usize) -> Result<&Padic> {
        if row >= self.rows || col >= self.cols {
            return Err(Error::domain(format!(
                "p-adic matrix index out of bounds: row={row}, col={col}, shape={}x{}",
                self.rows, self.cols
            )));
        }
        Ok(&self.data[row * self.cols + col])
    }

    fn ensure_domain_and_shape(&self, domain: &PadicDomain, side: &str) -> Result<()> {
        if self.domain.meta.prime != domain.meta.prime
            || self.domain.meta.precision != domain.meta.precision
        {
            return Err(Error::domain(format!(
                "{side} p-adic matrix domain mismatch: expected p={} precision={}, got p={} precision={}",
                domain.meta.prime,
                domain.meta.precision,
                self.domain.meta.prime,
                self.domain.meta.precision
            )));
        }
        if self.data.len() != self.rows.saturating_mul(self.cols) {
            return Err(Error::domain(format!(
                "{side} p-adic matrix shape/data mismatch: shape={}x{}, data={}",
                self.rows,
                self.cols,
                self.data.len()
            )));
        }
        Ok(())
    }
}

#[derive(Debug, Clone, PartialEq, Eq, Hash)]
pub struct Padic {
    pub meta: PadicMeta,
    pub residue: u128,
}

impl std::fmt::Display for Padic {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(
            f,
            "{{prime:{}, precision:{}, residue:{}}}",
            self.meta.prime, self.meta.precision, self.residue
        )
    }
}

impl Domain for DynamicPadicDomain {
    type Element = Vec<u64>;

    fn id(&self) -> DomainId {
        DomainId::new(format!(
            "Q_{}[dynamic:{}]",
            self.meta.prime, self.meta.precision
        ))
    }

    fn structure(&self) -> AlgebraicStructure {
        AlgebraicStructure::Field
    }

    fn precision_model(&self) -> PrecisionModel {
        PrecisionModel::Adaptive {
            min_digits: 0,
            max_digits: Some(self.meta.precision),
        }
    }

    fn contracts(&self) -> ContractSet {
        let scope = Scope::Domain(self.id().0);
        ContractSet::from_iter([
            Contract::new(
                ContractId(120),
                Claim::NonArchimedeanNorm,
                scope.clone(),
                Evidence::Axiom,
            ),
            Contract::new(
                ContractId(121),
                Claim::ValuationIdentity("v_p(x * y) = v_p(x) + v_p(y)".to_string()),
                scope.clone(),
                Evidence::Axiom,
            ),
            Contract::new(
                ContractId(122),
                Claim::PrecisionPreserved,
                scope,
                Evidence::Axiom,
            ),
        ])
    }
}

impl Domain for PadicDomain {
    type Element = u128;

    fn id(&self) -> DomainId {
        DomainId::new(format!("Q_{}[{}]", self.meta.prime, self.meta.precision))
    }

    fn structure(&self) -> AlgebraicStructure {
        AlgebraicStructure::Field
    }

    fn precision_model(&self) -> PrecisionModel {
        PrecisionModel::Fixed {
            digits: self.meta.precision,
        }
    }

    fn contracts(&self) -> ContractSet {
        let scope = Scope::Domain(self.id().0);
        ContractSet::from_iter([
            Contract::new(
                ContractId(100),
                Claim::NonArchimedeanNorm,
                scope.clone(),
                Evidence::Axiom,
            ),
            Contract::new(
                ContractId(101),
                Claim::ValuationIdentity("v_p(x * y) = v_p(x) + v_p(y)".to_string()),
                scope.clone(),
                Evidence::Axiom,
            ),
            Contract::new(
                ContractId(102),
                Claim::ValuationIdentity("v_p(x + y) >= min(v_p(x), v_p(y))".to_string()),
                scope.clone(),
                Evidence::Axiom,
            ),
            Contract::new(
                ContractId(103),
                Claim::PrecisionPreserved,
                scope,
                Evidence::Axiom,
            ),
        ])
    }
}

impl ValuedDomain for PadicDomain {
    type ValueGroup = Result<Option<u32>>;

    fn valuation(&self, value: &Self::Element) -> Self::ValueGroup {
        Ok(valuation_residue(
            *value % self.modulus(),
            self.meta.prime,
            self.meta.precision,
        ))
    }
}

impl NonArchimedeanDomain for PadicDomain {}

fn valuation_residue(value: u128, prime: u64, precision: u32) -> Option<u32> {
    if value == 0 {
        return None;
    }

    let mut n = value;
    let mut valuation = 0;
    let p = prime as u128;
    while valuation < precision && n % p == 0 {
        valuation += 1;
        n /= p;
    }
    Some(valuation)
}

fn checked_modulus(prime: u64, precision: u32) -> Result<u128> {
    let mut value = 1u128;
    for _ in 0..precision {
        value = value
            .checked_mul(prime as u128)
            .ok_or_else(|| Error::domain("p^N does not fit in u128"))?;
    }
    Ok(value)
}

fn unit_residue_at_valuation(
    residue: u128,
    valuation: u32,
    prime: u64,
    precision: u32,
) -> Result<u128> {
    if residue == 0 || valuation >= precision {
        return Ok(0);
    }
    let divisor = checked_modulus(prime, valuation)?;
    let unit_modulus = checked_modulus(prime, precision - valuation)?;
    Ok((residue / divisor) % unit_modulus)
}

fn valuation_bucket_fingerprint(
    prime: u64,
    precision: u32,
    rows: usize,
    cols: usize,
    valuations: &[u32],
    histogram: &BTreeMap<u32, usize>,
) -> String {
    let mut material =
        format!("padic-matrix-buckets-v1;p={prime};k={precision};shape={rows}x{cols};");
    material.push_str("valuations=");
    for (index, valuation) in valuations.iter().enumerate() {
        if index > 0 {
            material.push(',');
        }
        material.push_str(&valuation.to_string());
    }
    material.push_str(";histogram=");
    for (valuation, count) in histogram {
        material.push_str(&format!("{valuation}:{count},"));
    }
    format!("padic-bucket-fnv64:{:016x}", fnv1a64(material.as_bytes()))
}

fn fnv1a64(bytes: &[u8]) -> u64 {
    let mut hash = 0xcbf29ce484222325u64;
    for byte in bytes {
        hash ^= *byte as u64;
        hash = hash.wrapping_mul(0x100000001b3);
    }
    hash
}

fn is_prime(value: u64) -> bool {
    if value < 2 {
        return false;
    }
    if value == 2 {
        return true;
    }
    if value % 2 == 0 {
        return false;
    }
    let mut divisor = 3;
    while divisor * divisor <= value {
        if value % divisor == 0 {
            return false;
        }
        divisor += 2;
    }
    true
}

fn mod_inverse(value: i128, modulus: i128) -> Option<i128> {
    let (gcd, x, _) = extended_gcd(value.rem_euclid(modulus), modulus);
    if gcd != 1 {
        None
    } else {
        Some(x.rem_euclid(modulus))
    }
}

fn extended_gcd(a: i128, b: i128) -> (i128, i128, i128) {
    if b == 0 {
        (a, 1, 0)
    } else {
        let (gcd, x, y) = extended_gcd(b, a % b);
        (gcd, y, x - (a / b) * y)
    }
}

impl Padic {
    /// Lowest digit (residue mod prime) of this p-adic element, as a `u128`
    /// in the range `[0, prime)`. This is a pure function of `residue` and
    /// `meta.prime` and never panics.
    pub fn lowest_digit(&self) -> u128 {
        if self.meta.prime == 0 {
            return self.residue;
        }
        self.residue % (self.meta.prime as u128)
    }

    /// True when `self` and `other` share the same `prime` and `precision`
    /// metadata. The residue values are ignored. Useful as a precondition
    /// for arithmetic that requires matching domain parameters.
    pub fn precision_consistent_with(&self, other: &Padic) -> bool {
        self.meta.prime == other.meta.prime && self.meta.precision == other.meta.precision
    }

    /// True when `self` and `other` share the same `prime` only (precision
    /// may differ). Useful when comparing low-precision values against
    /// higher-precision containers.
    pub fn same_prime_as(&self, other: &Padic) -> bool {
        self.meta.prime == other.meta.prime
    }
}

impl Padic {
    /// Base-`p` digit expansion of the residue, lowest digit first, up to
    /// `meta.precision` digits. The residue is interpreted modulo
    /// `p^precision` and any digits beyond `precision` are dropped.
    pub fn digits_base_p(&self) -> Vec<u64> {
        let prime = self.meta.prime as u128;
        let precision = self.meta.precision as usize;
        if prime == 0 || precision == 0 {
            return Vec::new();
        }
        let mut digits = Vec::with_capacity(precision);
        let mut remaining = self.residue;
        for _ in 0..precision {
            let d = (remaining % prime) as u64;
            digits.push(d);
            remaining /= prime;
        }
        digits
    }
}

impl Padic {
    /// Sum of all base-`p` digits of this element. Returns 0 when precision
    /// is 0 or the prime is 0.
    pub fn digit_sum(&self) -> u64 {
        self.digits_base_p().iter().sum()
    }
}

impl Padic {
    /// Position of the first non-zero base-`p` digit, lowest-digit-first
    /// indexing. Returns `None` when the residue is zero (all digits are
    /// zero) or when precision is 0. For non-zero residues, the returned
    /// index is the p-adic valuation of the residue (the number of
    /// trailing zeros in its base-p expansion).
    pub fn first_nonzero_digit(&self) -> Option<usize> {
        let digits = self.digits_base_p();
        digits.iter().position(|&d| d != 0)
    }
}