clock_curve_math/extensible/

variable_precision.rs

//! Variable-precision arithmetic support.
//!
//! This module provides support for arithmetic with different limb counts
//! and field sizes, preparing for future cryptographic primitives.

#[cfg(feature = "alloc")]
extern crate alloc;
use super::multiplication::MultiplicationAlgorithm;
use crate::bigint::BigInt;
use crate::error::MathError;
#[cfg(feature = "alloc")]
use alloc::vec::Vec;

/// Variable-precision arithmetic support for cryptographic computations.
///
/// This struct provides a unified interface for modular arithmetic with different
/// precision levels and reduction techniques. It supports multiple limb counts
/// (64-bit words) and automatically precomputes constants for various reduction
/// algorithms to optimize performance for different modulus sizes.
///
/// # Supported Reduction Techniques
/// - **Montgomery Reduction**: Efficient for general moduli, provides constant-time
///   multiplication and reduction operations
/// - **Barrett Reduction**: Fast reduction using precomputed constants, avoids
///   expensive division operations
/// - **Solinas Reduction**: Specialized for moduli of the form 2^a ± 2^b ± 1,
///   used by curves like Curve25519 (p = 2^255 - 19)
///
/// # Performance Optimization
/// Precomputes all necessary constants during construction to minimize per-operation
/// overhead. Supports algorithm selection based on operand sizes and available
/// computational resources.
///
/// # Use Cases
/// - Elliptic curve cryptography with different field sizes
/// - Multi-precision arithmetic for cryptographic protocols
/// - Performance-critical modular computations
///
/// # Memory Layout
/// Values are stored as multi-precision integers using the specified number of
/// limbs. The Montgomery form uses R = 2^(limb_count × 64) as the Montgomery radix.
#[derive(Clone, Debug)]
#[allow(dead_code)]
pub struct VariablePrecision {
    /// The number of limbs (64-bit words) used for multi-precision arithmetic
    limb_count: usize,
    /// The modulus defining the finite field for all operations
    modulus: BigInt,
    /// Montgomery radix R = 2^(limb_count × 64)
    montgomery_r: BigInt,
    /// Montgomery constant R² mod modulus for conversion to Montgomery form
    montgomery_r2: BigInt,
    /// Montgomery inverse N' = -modulus⁻¹ mod R for Montgomery reduction
    montgomery_n_prime: BigInt,
    /// Barrett reduction constant μ = floor(2^(2×k) / modulus) where k = bit_length
    barrett_mu: Option<BigInt>,
    /// Solinas reduction constants for special moduli of form 2^a ± 2^b ± 1
    #[cfg(feature = "alloc")]
    solinas_constants: Option<Vec<BigInt>>,
}

impl VariablePrecision {
    /// Create a new variable precision arithmetic context.
    ///
    /// Initializes a context for modular arithmetic with the specified limb count
    /// and modulus. Precomputes all necessary constants for Montgomery, Barrett,
    /// and other reduction techniques to optimize subsequent operations.
    ///
    /// # Arguments
    /// * `limb_count` - Number of 64-bit limbs for multi-precision arithmetic (1-8)
    /// * `modulus` - The prime modulus defining the finite field
    ///
    /// # Returns
    /// A new `VariablePrecision` instance ready for modular arithmetic operations
    ///
    /// # Errors
    /// Returns `MathError::InvalidInput` if:
    /// - `limb_count` is 0
    /// - `limb_count` exceeds 8 (implementation limit)
    ///
    /// # Precomputation Details
    /// - **Montgomery Constants**: Computes R, R², and N' for Montgomery arithmetic
    /// - **Barrett Constant**: Computes μ for Barrett reduction
    /// - **Solinas Constants**: Placeholder for future Solinas prime support
    ///
    /// # Performance
    /// Construction involves modular exponentiation and inverse computation,
    /// which is expensive but done once. Subsequent operations are optimized.
    ///
    /// # Examples
    /// ```ignore
    /// use clock_curve_math::extensible::variable_precision::VariablePrecision;
    /// use clock_curve_math::BigInt;
    ///
    /// // Create context for 256-bit arithmetic (4 limbs × 64 bits)
    /// let modulus = BigInt::from_hex("fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f")?; // secp256k1 modulus
    /// let ctx = VariablePrecision::new(4, modulus)?;
    /// # Ok::<(), clock_curve_math::MathError>(())
    /// ```
    pub fn new(limb_count: usize, modulus: BigInt) -> Result<Self, MathError> {
        if limb_count == 0 || limb_count > 8 {
            return Err(MathError::InvalidInput);
        }

        // Compute Montgomery constants
        let r_limbs = limb_count as u32;
        let montgomery_r = BigInt::from_u64(1).shl(r_limbs * 64);

        // Compute R² = R^2 mod modulus where R = 2^(limb_count * 64)
        let montgomery_r2 = {
            let r_mod_n = super::utils::mod_pow(
                &BigInt::from_u64(2),
                &BigInt::from_u64((r_limbs * 64) as u64),
                &modulus,
            );
            super::utils::mod_pow(&r_mod_n, &BigInt::from_u64(2), &modulus)
        };

        // Compute N' = -N^(-1) mod R
        let montgomery_n_prime =
            if let Some(inv) = super::utils::mod_inverse(&modulus, &montgomery_r) {
                let mut result = BigInt::from_u64(0).sub(&inv);
                while result.ct_cmp(&BigInt::from_u64(0)) == core::cmp::Ordering::Less {
                    result = result.add(&montgomery_r);
                }
                result
            } else {
                // This shouldn't happen for prime moduli, but handle gracefully
                BigInt::from_u64(0)
            };

        // Compute Barrett reduction constant mu = floor(2^(2*k) / m) where k is bit length
        let bit_length = modulus.bit_length();
        let barrett_mu = if bit_length > 0 {
            BigInt::from_u64(1).shl(2 * bit_length).div_rem(&modulus).0
        } else {
            BigInt::from_u64(0)
        };

        Ok(Self {
            limb_count,
            modulus,
            montgomery_r,
            montgomery_r2,
            montgomery_n_prime,
            barrett_mu: Some(barrett_mu),
            solinas_constants: None, // Would be computed for Solinas primes
        })
    }

    /// Get the number of limbs used for multi-precision arithmetic.
    ///
    /// Returns the limb count specified during construction, which determines
    /// the maximum precision for arithmetic operations.
    ///
    /// # Returns
    /// The number of 64-bit limbs (1-8) used for multi-precision arithmetic
    pub fn limb_count(&self) -> usize {
        self.limb_count
    }

    /// Get the modulus defining the finite field.
    ///
    /// Returns a reference to the prime modulus used for all modular arithmetic
    /// operations in this context.
    ///
    /// # Returns
    /// A reference to the BigInt modulus
    pub fn modulus(&self) -> &BigInt {
        &self.modulus
    }

    /// Perform Montgomery multiplication with variable precision.
    ///
    /// Computes `(a * b * R^(-1)) mod modulus` where `R = 2^(limb_count * 64)`.
    /// This function implements the Montgomery multiplication algorithm for
    /// efficient modular arithmetic.
    ///
    /// # Arguments
    /// * `a` - First operand (must be in Montgomery form)
    /// * `b` - Second operand (must be in Montgomery form)
    ///
    /// # Returns
    /// The product `(a * b * R^(-1)) mod modulus` in Montgomery form
    ///
    /// # Algorithm
    /// Uses the standard Montgomery multiplication formula:
    /// 1. Compute t = a * b
    /// 2. Compute m = (t * N') mod R
    /// 3. Compute u = (t + m * N) / R
    /// 4. If u ≥ N, return u - N, else return u
    ///
    /// Where N is the modulus, R is the Montgomery radix, and N' is the
    /// precomputed Montgomery inverse.
    ///
    /// # Performance
    /// O(limb_count²) time complexity, constant-time for fixed limb counts.
    /// Significantly faster than standard modular multiplication for multiple operations.
    ///
    /// # Note
    /// Both operands must already be in Montgomery form. Use `to_montgomery()` to
    /// convert standard form values to Montgomery form before multiplication.
    pub fn montgomery_mul(&self, a: &BigInt, b: &BigInt) -> BigInt {
        // Montgomery multiplication algorithm
        let t = a.mul(b);
        let m = t
            .mul(&self.montgomery_n_prime)
            .div_rem(&self.montgomery_r)
            .1;
        let u = t.add(&m.mul(&self.modulus)).div_rem(&self.montgomery_r).1;

        // Final reduction
        if u.ct_cmp(&self.modulus) != core::cmp::Ordering::Less {
            u.sub(&self.modulus)
        } else {
            u
        }
    }

    /// Convert a value to Montgomery form for efficient modular arithmetic.
    ///
    /// Montgomery form represents the value x as `x * R mod modulus`, where
    /// R is a power of 2 larger than the modulus. This representation enables
    /// faster modular multiplication by avoiding expensive division operations.
    ///
    /// # Arguments
    /// * `value` - The value to convert to Montgomery form (0 ≤ value < modulus)
    ///
    /// # Returns
    /// The value converted to Montgomery form: `(value * R) mod modulus`
    ///
    /// # Mathematical Background
    /// In Montgomery arithmetic, numbers are represented in a transformed space
    /// where multiplication can be performed more efficiently. The conversion is:
    ///
    /// `montgomery_value = (value * R) mod modulus`
    ///
    /// Where R = 2^(limb_count × 64) is chosen to be larger than modulus.
    ///
    /// # Usage
    /// ```text
    /// let a_mont = ctx.to_montgomery(&a);
    /// let b_mont = ctx.to_montgomery(&b);
    /// let product_mont = ctx.montgomery_mul(&a_mont, &b_mont);
    /// let product = ctx.from_montgomery(&product_mont);
    /// ```
    ///
    /// # Performance
    /// Uses precomputed R² constant for efficient conversion via Montgomery multiplication.
    pub fn to_montgomery(&self, value: &BigInt) -> BigInt {
        // x * R mod N = montgomery_mul(x, R²)
        self.montgomery_mul(value, &self.montgomery_r2)
    }

    /// Convert a value from Montgomery form back to standard representation.
    ///
    /// Converts a Montgomery-form value back to its standard integer representation.
    /// This is the inverse operation of `to_montgomery()`.
    ///
    /// # Arguments
    /// * `value` - The Montgomery-form value to convert
    ///
    /// # Returns
    /// The value converted back to standard form: `(value * R^(-1)) mod modulus`
    ///
    /// # Mathematical Background
    /// Given a Montgomery value `mont_value = (x * R) mod modulus`, this function
    /// computes `x = (mont_value * R^(-1)) mod modulus` using Montgomery reduction.
    ///
    /// The algorithm performs Montgomery reduction to multiply by R^(-1) modulo modulus:
    /// 1. Compute m = (value * N') mod R
    /// 2. Compute u = (value + m * N) / R
    /// 3. Return u - N if u ≥ N, else return u
    ///
    /// # Usage
    /// ```text
    /// let a_mont = ctx.to_montgomery(&a);
    /// let b_mont = ctx.to_montgomery(&b);
    /// let product_mont = ctx.montgomery_mul(&a_mont, &b_mont);
    /// let product = ctx.from_montgomery(&product_mont); // Get final result
    /// ```
    ///
    /// # Performance
    /// O(limb_count²) time complexity using Montgomery reduction algorithm.
    pub fn from_montgomery(&self, value: &BigInt) -> BigInt {
        // Montgomery reduction to convert back to standard form
        let m = value
            .mul(&self.montgomery_n_prime)
            .div_rem(&self.montgomery_r)
            .1;
        let u = value
            .add(&m.mul(&self.modulus))
            .div_rem(&self.montgomery_r)
            .1;

        // Final reduction
        if u.ct_cmp(&self.modulus) != core::cmp::Ordering::Less {
            u.sub(&self.modulus)
        } else {
            u
        }
    }

    /// Perform Barrett reduction for fast modular reduction.
    ///
    /// Barrett reduction uses precomputed constants to avoid expensive divisions.
    /// This is faster than standard modular reduction for large moduli.
    ///
    pub fn barrett_reduce(&self, value: &BigInt) -> BigInt {
        if let Some(mu) = &self.barrett_mu {
            // Barrett reduction: x mod m ≈ x - floor(x * mu / 2^k) * m
            let k = self.modulus.bit_length();
            let q = value.mul(mu).shr(k); // Approximate quotient
            let r = value.sub(&q.mul(&self.modulus)); // Approximate remainder

            // Correct the result if necessary
            if r.ct_cmp(&self.modulus) != core::cmp::Ordering::Less {
                r.sub(&self.modulus)
            } else {
                r
            }
        } else {
            // Fallback to standard reduction
            value.div_rem(&self.modulus).1
        }
    }

    /// Perform Solinas reduction for special moduli of cryptographic interest.
    ///
    /// Solinas reduction is a specialized algorithm for moduli of the form
    /// `2^a ± 2^b ± 1`. These moduli are common in elliptic curve cryptography,
    /// including Curve25519 (p = 2^255 - 19) and other standardized curves.
    ///
    /// # Arguments
    /// * `value` - The value to reduce modulo the special modulus
    ///
    /// # Returns
    /// The reduced value: `value mod modulus` (0 ≤ result < modulus)
    ///
    /// # Algorithm
    /// For a modulus p = 2^a ± 2^b ± 1, Solinas reduction exploits the special
    /// form to perform reduction efficiently:
    ///
    /// 1. Compute q = floor(x / 2^a)
    /// 2. Compute correction = q * (2^a ± 2^b ± 1)
    /// 3. Compute result = x - correction
    /// 4. Final adjustment to ensure result ∈ [0, p-1]
    ///
    /// This avoids expensive division by using bit shifts and the special
    /// structure of the modulus.
    ///
    /// # Supported Moduli
    /// - **Curve25519**: p = 2^255 - 19
    /// - **Ed25519**: Same as Curve25519
    /// - **NIST P-256**: p = 2^256 - 2^224 + 2^192 + 2^96 - 1 (Solinas-like)
    /// - Other moduli of the form 2^a ± 2^b ± 1
    ///
    /// # Performance
    /// O(1) time complexity for fixed parameter sizes, using only bit operations
    /// and additions. Extremely fast for supported moduli.
    ///
    /// # Configuration
    /// Requires precomputed Solinas constants specifying the parameters (a, b)
    /// and signs for the modulus form. Currently not implemented - falls back
    /// to standard modular reduction.
    ///
    /// # Note
    /// This implementation is currently a placeholder. Full Solinas reduction
    /// requires additional constants and parameter detection logic.
    pub fn solinas_reduce(&self, value: &BigInt) -> BigInt {
        #[cfg(feature = "alloc")]
        if let Some(constants) = &self.solinas_constants {
            // Solinas reduction for moduli of the form 2^a ± 2^b ± 1
            // The constants contain the parameters (a, b, sign1, sign2)

            // For a modulus p = 2^a ± 2^b ± 1, we can reduce x mod p using:
            // x mod p = x - floor(x / 2^a) * (2^a ± 2^b ± 1) + adjustments

            // Extract parameters from constants
            // constants[0] = a, constants[1] = b
            // constants[2] = sign for 2^b term, constants[3] = sign for constant term
            if constants.len() >= 4 {
                let a = constants[0].to_u64().unwrap_or(0) as u32;
                let b = constants[1].to_u64().unwrap_or(0) as u32;
                let sign_b = constants[2].to_u64().unwrap_or(0) as i32;
                let sign_1 = constants[3].to_u64().unwrap_or(0) as i32;

                // Compute x div 2^a and x mod 2^a
                let x_shift_a = value.shr(a);
                let mask = BigInt::from_u64(1).shl(a).sub(&BigInt::from_u64(1));
                let _x_low_a = value.and(&mask);

                // Compute the correction term: floor(x / 2^a) * (2^a ± 2^b ± 1)
                let mut correction = x_shift_a.shl(a); // floor(x / 2^a) * 2^a

                if sign_b != 0 {
                    let b_term = x_shift_a.shl(b);
                    if sign_b > 0 {
                        correction = correction.add(&b_term);
                    } else {
                        correction = correction.sub(&b_term);
                    }
                }

                if sign_1 != 0 {
                    if sign_1 > 0 {
                        correction = correction.add(&x_shift_a);
                    } else {
                        correction = correction.sub(&x_shift_a);
                    }
                }

                // Compute result = x - correction
                let mut result = value.sub(&correction);

                // Final reduction: ensure result is in [0, p-1]
                if result.ct_cmp(&BigInt::from_u64(0)) == core::cmp::Ordering::Less {
                    result = result.add(&self.modulus);
                }

                while result.ct_cmp(&self.modulus) != core::cmp::Ordering::Less {
                    result = result.sub(&self.modulus);
                }

                result
            } else {
                // Invalid constants, fall back to standard reduction
                value.div_rem(&self.modulus).1
            }
        } else {
            // No Solinas constants available, fall back to standard reduction
            value.div_rem(&self.modulus).1
        }
    }

    /// Perform multi-limb multiplication using a specified algorithm.
    ///
    /// This method provides fine-grained control over the multiplication algorithm
    /// used for computing `a * b`. Different algorithms offer various trade-offs
    /// between performance, memory usage, and computational complexity.
    ///
    /// # Arguments
    /// * `a` - First operand
    /// * `b` - Second operand
    /// * `algorithm` - The multiplication algorithm to use
    ///
    /// # Returns
    /// The product `a * b` computed using the specified algorithm
    ///
    /// # Supported Algorithms
    /// - **`Standard`/`Schoolbook`**: Traditional long multiplication, O(n²)
    /// - **`Karatsuba`**: Divide-and-conquer algorithm, O(n^1.585)
    /// - **`ToomCook`**: Polynomial evaluation/interpolation, O(n^1.465)
    /// - **`FFT`**: Fast Fourier Transform based, O(n log n)
    ///
    /// # Algorithm Selection Guidelines
    /// - **Small operands** (< 1000 bits): Use `Standard` - overhead not worth it
    /// - **Medium operands** (1000-10000 bits): Use `Karatsuba` - good balance
    /// - **Large operands** (10000-50000 bits): Use `ToomCook` - better asymptotic
    /// - **Very large operands** (>50000 bits): Use `FFT` - optimal asymptotic
    ///
    /// # Performance Notes
    /// - Algorithm selection affects both time and space complexity
    /// - Some algorithms may have higher constant factors
    /// - FFT requires additional feature flags and may not be available
    ///
    /// # Examples
    /// ```ignore
    /// use crate::extensible::multiplication::MultiplicationAlgorithm;
    ///
    /// // Use Karatsuba for medium-sized numbers
    /// let product = ctx.mul_with_algorithm(&a, &b, MultiplicationAlgorithm::Karatsuba);
    ///
    /// // Use FFT for very large numbers (if available)
    /// let big_product = ctx.mul_with_algorithm(&x, &y, MultiplicationAlgorithm::FFT);
    /// ```
    pub fn mul_with_algorithm(
        &self,
        a: &BigInt,
        b: &BigInt,
        algorithm: MultiplicationAlgorithm,
    ) -> BigInt {
        match algorithm {
            MultiplicationAlgorithm::Standard | MultiplicationAlgorithm::Schoolbook => a.mul(b),
            MultiplicationAlgorithm::Karatsuba => {
                super::multiplication::mul_karatsuba_standalone(a, b)
            }
            MultiplicationAlgorithm::ToomCook => {
                super::multiplication::mul_toom_cook_standalone(a, b)
            }
            MultiplicationAlgorithm::FFT => super::multiplication::mul_fft_standalone(a, b),
        }
    }
}