clock_curve_math/extensible/

fft.rs

//! Fast Fourier Transform (FFT) implementation for polynomial multiplication.
//!
//! This module provides O(n log n) multiplication of large integers using
//! the Fast Fourier Transform algorithm. It serves as a high-performance
//! alternative to Karatsuba and Toom-Cook multiplication for very large operands.
//!
//! Note: Full FFT operations require the "std" feature for trigonometric functions
//! and memory allocation. Without std, it falls back to schoolbook multiplication.

use crate::bigint::BigInt;
use crate::error::MathError;

/// Complex number representation for FFT
///
/// This struct provides basic complex arithmetic operations needed
/// for Fast Fourier Transform algorithms. It includes addition,
/// multiplication, and scaling operations optimized for FFT usage.
#[derive(Clone, Copy, Debug)]
pub struct Complex {
    /// The real part of the complex number
    pub real: f64,
    /// The imaginary part of the complex number
    pub imag: f64,
}

impl Complex {
    /// Create a new complex number
    ///
    /// # Arguments
    /// * `real` - The real part
    /// * `imag` - The imaginary part
    ///
    /// # Returns
    /// A new Complex number with the specified real and imaginary parts
    pub const fn new(real: f64, imag: f64) -> Self {
        Self { real, imag }
    }

    /// Create a complex number from a real number
    ///
    /// # Arguments
    /// * `real` - The real value to convert
    ///
    /// # Returns
    /// A complex number with the real value and zero imaginary part
    pub const fn from_real(real: f64) -> Self {
        Self::new(real, 0.0)
    }

    /// Complex addition
    ///
    /// # Arguments
    /// * `other` - The complex number to add
    ///
    /// # Returns
    /// The sum of self and other
    pub fn add(self, other: Self) -> Self {
        Self::new(self.real + other.real, self.imag + other.imag)
    }

    /// Complex subtraction
    ///
    /// # Arguments
    /// * `other` - The complex number to subtract
    ///
    /// # Returns
    /// The difference of self and other
    pub fn sub(self, other: Self) -> Self {
        Self::new(self.real - other.real, self.imag - other.imag)
    }

    /// Complex multiplication
    ///
    /// # Arguments
    /// * `other` - The complex number to multiply
    ///
    /// # Returns
    /// The product of self and other
    pub fn mul(self, other: Self) -> Self {
        Self::new(
            self.real * other.real - self.imag * other.imag,
            self.real * other.imag + self.imag * other.real,
        )
    }

    /// Scalar multiplication
    ///
    /// # Arguments
    /// * `scalar` - The real number to multiply by
    ///
    /// # Returns
    /// The complex number scaled by the scalar
    pub fn scale(self, scalar: f64) -> Self {
        Self::new(self.real * scalar, self.imag * scalar)
    }

    /// Complex exponential (Euler's formula)
    ///
    /// Computes e^(i*angle) = cos(angle) + i*sin(angle)
    ///
    /// # Arguments
    /// * `angle` - The angle in radians
    ///
    /// # Returns
    /// The complex exponential of the angle
    #[cfg(feature = "std")]
    pub fn exp(angle: f64) -> Self {
        Self::new(angle.cos(), angle.sin())
    }

    /// Fallback for no_std - requires std feature
    ///
    /// This function is not available in no_std environments and will panic.
    /// Use the std feature to enable FFT functionality.
    #[cfg(not(feature = "std"))]
    pub fn exp(_angle: f64) -> Self {
        panic!("FFT requires std feature for trigonometric functions")
    }
}

/// FFT multiplier for large integer multiplication
///
/// This struct provides Fast Fourier Transform-based multiplication
/// for very large integers. It offers O(n log n) complexity compared
/// to O(n²) for schoolbook multiplication, making it ideal for
/// cryptographic applications dealing with large numbers.
///
/// The FFT implementation requires the "std" feature for trigonometric
/// functions and memory allocation. Without std, it falls back to
/// schoolbook multiplication.
#[cfg(feature = "std")]
pub struct FFTMultiplier {
    // No fields needed for basic implementation
}

#[cfg(feature = "std")]
impl FFTMultiplier {
    /// Create a new FFT multiplier
    ///
    /// # Returns
    /// A new FFTMultiplier instance ready for large integer multiplication
    pub fn new() -> Self {
        Self {}
    }

    /// Check if n is a power of 2
    ///
    /// # Arguments
    /// * `n` - The number to check
    ///
    /// # Returns
    /// true if n is a power of 2, false otherwise
    fn is_power_of_two(n: usize) -> bool {
        n > 0 && (n & (n - 1)) == 0
    }

    /// Find next power of 2 >= n
    ///
    /// This is used to pad input arrays to optimal FFT sizes.
    ///
    /// # Arguments
    /// * `n` - The minimum size required
    ///
    /// # Returns
    /// The smallest power of 2 that is >= n
    fn next_power_of_two(n: usize) -> usize {
        if n <= 1 {
            return 1;
        }
        // Use bit manipulation to find next power of 2
        let mut result = 1;
        while result < n {
            result <<= 1;
        }
        result
    }

    /// Perform bit-reversal permutation on complex array
    ///
    /// This reorders the array elements to match the Cooley-Tukey
    /// FFT algorithm's requirements, improving cache locality.
    ///
    /// # Arguments
    /// * `coeffs` - The complex coefficient array to permute in-place
    fn bit_reverse_permutation(coeffs: &mut [Complex]) {
        let n = coeffs.len();
        let log_n = (n as f64).log2() as usize;

        for i in 0..n {
            let mut j = 0;
            for k in 0..log_n {
                if (i & (1 << k)) != 0 {
                    j |= 1 << (log_n - 1 - k);
                }
            }
            if j > i {
                coeffs.swap(i, j);
            }
        }
    }

    /// Cooley-Tukey Fast Fourier Transform
    ///
    /// Implements the iterative Cooley-Tukey FFT algorithm for transforming
    /// a sequence of complex coefficients between time and frequency domains.
    ///
    /// # Arguments
    /// * `coeffs` - The complex coefficient array (modified in-place)
    /// * `inverse` - If true, compute inverse FFT; if false, forward FFT
    ///
    /// # Panics
    /// Panics if the length of coeffs is not a power of 2
    pub fn cooley_tukey_fft(&self, coeffs: &mut [Complex], inverse: bool) {
        let n = coeffs.len();

        // Bit-reversal permutation
        Self::bit_reverse_permutation(coeffs);

        // Iterative Cooley-Tukey algorithm
        let mut size = 2;
        while size <= n {
            let half_size = size / 2;
            let angle_step =
                if inverse { 2.0 } else { -2.0 } * std::f64::consts::PI / (size as f64);
            let wlen = Complex::exp(angle_step);

            let mut i = 0;
            while i < n {
                let mut w = Complex::new(1.0, 0.0);
                for j in 0..half_size {
                    let u = coeffs[i + j];
                    let v = coeffs[i + j + half_size].mul(w);
                    coeffs[i + j] = u.add(v);
                    coeffs[i + j + half_size] = u.sub(v);
                    w = w.mul(wlen);
                }
                i += size;
            }

            size *= 2;
        }

        // Scale for inverse transform
        if inverse {
            let scale = 1.0 / (n as f64);
            for coeff in coeffs.iter_mut() {
                *coeff = coeff.scale(scale);
            }
        }
    }

    /// Convert BigInt coefficients to complex Vec
    fn bigint_to_complex(coeffs: &[BigInt]) -> Vec<Complex> {
        coeffs
            .iter()
            .map(|c| {
                // Convert BigInt to f64 (this may lose precision for very large numbers)
                // For cryptographic use, we limit to reasonable sizes
                let real_val = c.to_u64().unwrap_or(0) as f64;
                Complex::from_real(real_val)
            })
            .collect()
    }

    /// Convert complex array back to rounded integer coefficients
    fn complex_to_coeffs(complex_coeffs: &[Complex]) -> Vec<i64> {
        complex_coeffs
            .iter()
            .map(|c| {
                // Round to nearest integer
                c.real.round() as i64
            })
            .collect()
    }

    /// Multiply two polynomials using FFT
    ///
    /// This implements the standard FFT-based polynomial multiplication algorithm
    /// with O(n log n) complexity instead of O(n²) for schoolbook multiplication.
    ///
    /// The algorithm:
    /// 1. Convert coefficients to complex numbers
    /// 2. Pad to next power of 2 length for optimal FFT performance
    /// 3. Compute forward FFT of both polynomials
    /// 4. Multiply pointwise in frequency domain
    /// 5. Compute inverse FFT to get result coefficients
    /// 6. Extract and round integer coefficients
    ///
    /// # Arguments
    /// * `a` - Coefficients of the first polynomial
    /// * `b` - Coefficients of the second polynomial
    ///
    /// # Returns
    /// Result coefficients of the product polynomial, or MathError on failure
    pub fn polynomial_multiply(
        &self,
        a: &[BigInt],
        b: &[BigInt],
    ) -> Result<Vec<BigInt>, MathError> {
        // Determine FFT size
        let n = Self::next_power_of_two(a.len() + b.len() - 1);

        // Convert to complex numbers and pad
        let mut a_complex = Self::bigint_to_complex(a);
        let mut b_complex = Self::bigint_to_complex(b);

        a_complex.resize(n, Complex::from_real(0.0));
        b_complex.resize(n, Complex::from_real(0.0));

        // Forward FFT
        self.cooley_tukey_fft(&mut a_complex, false);
        self.cooley_tukey_fft(&mut b_complex, false);

        // Pointwise multiplication in frequency domain
        let mut c_complex: Vec<Complex> = a_complex
            .iter()
            .zip(b_complex.iter())
            .map(|(x, y)| x.mul(*y))
            .collect();

        // Inverse FFT
        self.cooley_tukey_fft(&mut c_complex, true);

        // Convert back to integer coefficients
        let coeffs_i64 = Self::complex_to_coeffs(&c_complex);

        // Convert to BigInt and remove trailing zeros
        let mut result: Vec<BigInt> = coeffs_i64
            .iter()
            .map(|&x| {
                if x >= 0 {
                    BigInt::from_u64(x as u64)
                } else {
                    // Handle negative coefficients (shouldn't happen in multiplication)
                    BigInt::from_u64(0)
                }
            })
            .collect();

        // Remove trailing zeros
        while result.len() > 1 {
            if let Some(last) = result.last() {
                if last.is_zero() {
                    result.pop();
                } else {
                    break;
                }
            } else {
                break;
            }
        }

        Ok(result)
    }

    /// Multiply two large integers using FFT
    ///
    /// This method provides asymptotically optimal multiplication for very large integers
    /// by converting them to polynomial coefficients, multiplying using FFT, and converting
    /// back to integer representation.
    ///
    /// # Arguments
    /// * `a` - First large integer to multiply
    /// * `b` - Second large integer to multiply
    ///
    /// # Returns
    /// The product a * b as a BigInt, or MathError on failure
    ///
    /// # Performance
    /// Offers O(n log n) complexity for n-bit numbers vs O(n²) for schoolbook multiplication
    pub fn multiply_big_ints(&self, a: &BigInt, b: &BigInt) -> Result<BigInt, MathError> {
        // Use base 10^9 for coefficient conversion (good balance of precision and efficiency)
        const BASE: u64 = 1_000_000_000;

        // Convert to base-10^9 digits
        let a_coeffs = Self::big_int_to_base_digits(a, BASE);
        let b_coeffs = Self::big_int_to_base_digits(b, BASE);

        // Multiply polynomials
        let product_coeffs = self.polynomial_multiply(&a_coeffs, &b_coeffs)?;

        // Convert back to BigInt
        let result = Self::coeffs_to_big_int(&product_coeffs, BASE);

        Ok(result)
    }

    /// Convert BigInt to base digits
    ///
    /// Decomposes a large integer into digits in the specified base.
    /// This is used to convert integers to polynomial coefficients for FFT multiplication.
    ///
    /// # Arguments
    /// * `n` - The BigInt to convert
    /// * `base` - The numerical base for digit decomposition
    ///
    /// # Returns
    /// Vector of digits in the specified base (least significant first)
    fn big_int_to_base_digits(n: &BigInt, base: u64) -> Vec<BigInt> {
        if n.is_zero() {
            return vec![BigInt::from_u64(0)];
        }

        let mut digits = Vec::new();
        let mut current = n.clone();
        let base_big = BigInt::from_u64(base);

        while !current.is_zero() {
            let (quotient, remainder) = current.div_rem(&base_big);
            digits.push(remainder);
            current = quotient;
        }

        digits
    }

    /// Convert coefficient array back to BigInt
    ///
    /// Reconstructs a BigInt from its base-digit representation.
    ///
    /// # Arguments
    /// * `coeffs` - The digit coefficients (least significant first)
    /// * `base` - The numerical base used for the digits
    ///
    /// # Returns
    /// The reconstructed BigInt value
    fn coeffs_to_big_int(coeffs: &[BigInt], base: u64) -> BigInt {
        let mut result = BigInt::from_u64(0);
        let base_big = BigInt::from_u64(base);

        for (i, coeff) in coeffs.iter().enumerate() {
            if i == 0 {
                result = coeff.clone();
            } else {
                // result += coeff * (base^i)
                // For efficiency, we use iterative multiplication instead of pow
                let mut power = BigInt::from_u64(1);
                for _ in 0..i {
                    power = power.mul(&base_big);
                }
                let term = coeff.mul(&power);
                result = result.add(&term);
            }
        }

        result
    }
}

/// Fallback FFT multiplier for no_std environments
///
/// This implementation provides API compatibility when the std feature is not available.
/// It falls back to basic schoolbook multiplication since FFT requires trigonometric
/// functions and complex arithmetic which are not available in no_std.
#[cfg(not(feature = "std"))]
pub struct FFTMultiplier;

#[cfg(not(feature = "std"))]
impl FFTMultiplier {
    /// Create a new FFT multiplier (no_std fallback)
    ///
    /// # Returns
    /// A new FFTMultiplier instance that uses schoolbook multiplication
    pub const fn new() -> Self {
        Self
    }

    /// Multiply using schoolbook multiplication in no_std
    ///
    /// Since FFT requires std features for trigonometric functions,
    /// this falls back to basic BigInt multiplication.
    ///
    /// # Arguments
    /// * `a` - First integer to multiply
    /// * `b` - Second integer to multiply
    ///
    /// # Returns
    /// The product a * b
    pub fn multiply_big_ints(&self, a: &BigInt, b: &BigInt) -> Result<BigInt, MathError> {
        // In no_std, fall back to schoolbook multiplication
        // FFT requires std feature for trigonometric functions
        Ok(mul_schoolbook(a, b))
    }
}

/// Schoolbook multiplication for fallback
///
/// Basic polynomial multiplication using BigInt's built-in multiplication.
/// This serves as the fallback when FFT is not available.
fn mul_schoolbook(a: &BigInt, b: &BigInt) -> BigInt {
    a.mul(b)
}

/// Standalone FFT multiplication function
///
/// This is the main entry point for FFT-based multiplication in the library.
/// It automatically chooses between FFT and fallback algorithms based on input size
/// and available features.
///
/// # Algorithm Selection:
/// - Uses FFT for very large numbers (>10,000 bits) when std feature is enabled
/// - Falls back to Toom-Cook multiplication for medium-sized numbers
/// - Uses schoolbook multiplication for small numbers or when std is unavailable
///
/// # Arguments
/// * `a` - First large integer to multiply
/// * `b` - Second large integer to multiply
///
/// # Returns
/// The product a * b as a BigInt
///
/// # Performance
/// Provides O(n log n) complexity for large integers where FFT is beneficial
pub fn mul_fft_standalone(a: &BigInt, b: &BigInt) -> BigInt {
    // FFT multiplication for extremely large operands
    // FFT provides O(n log n) complexity vs O(n^2) for schoolbook multiplication

    let a_bits = a.bit_length();
    let b_bits = b.bit_length();
    let max_bits = core::cmp::max(a_bits, b_bits);

    // Use FFT threshold - in practice this would be much lower with full implementation
    const FFT_THRESHOLD_BITS: u32 = 10000; // Use FFT for numbers > 10,000 bits

    if max_bits < FFT_THRESHOLD_BITS {
        return super::multiplication::mul_toom_cook_standalone(a, b);
    }

    // Use FFT for very large numbers (requires std feature)
    #[cfg(feature = "std")]
    {
        let multiplier = FFTMultiplier::new();
        match multiplier.multiply_big_ints(a, b) {
            Ok(result) => result,
            Err(_) => {
                // Fall back to Toom-Cook if FFT fails
                super::multiplication::mul_toom_cook_standalone(a, b)
            }
        }
    }

    #[cfg(not(feature = "std"))]
    {
        // Fall back to Toom-Cook in no_std environments
        super::multiplication::mul_toom_cook_standalone(a, b)
    }
}

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

    /// Test basic complex number arithmetic operations
    #[test]
    fn test_complex_arithmetic() {
        let z1 = Complex::new(1.0, 2.0);
        let z2 = Complex::new(3.0, 4.0);

        let sum = z1.add(z2);
        assert!((sum.real - 4.0).abs() < 1e-10);
        assert!((sum.imag - 6.0).abs() < 1e-10);

        let product = z1.mul(z2);
        assert!((product.real - -5.0).abs() < 1e-10);
        assert!((product.imag - 10.0).abs() < 1e-10);
    }

    /// Test power-of-two detection utility function
    #[cfg(feature = "std")]
    #[test]
    fn test_power_of_two() {
        assert!(FFTMultiplier::is_power_of_two(1));
        assert!(FFTMultiplier::is_power_of_two(2));
        assert!(FFTMultiplier::is_power_of_two(4));
        assert!(FFTMultiplier::is_power_of_two(8));
        assert!(!FFTMultiplier::is_power_of_two(3));
        assert!(!FFTMultiplier::is_power_of_two(6));
    }

    /// Test next power-of-two calculation
    #[cfg(feature = "std")]
    #[test]
    fn test_next_power_of_two() {
        assert_eq!(FFTMultiplier::next_power_of_two(1), 1);
        assert_eq!(FFTMultiplier::next_power_of_two(2), 2);
        assert_eq!(FFTMultiplier::next_power_of_two(3), 4);
        assert_eq!(FFTMultiplier::next_power_of_two(5), 8);
        assert_eq!(FFTMultiplier::next_power_of_two(9), 16);
    }

    /// Test polynomial multiplication using FFT
    #[cfg(feature = "std")]
    #[test]
    fn test_fft_polynomial_multiply() {
        let multiplier = FFTMultiplier::new();

        // (1 + 2x + 3x²) * (4 + 5x) = 4 + 13x + 22x² + 15x³
        let a = vec![
            BigInt::from_u64(1),
            BigInt::from_u64(2),
            BigInt::from_u64(3),
        ];
        let b = vec![BigInt::from_u64(4), BigInt::from_u64(5)];

        let result = multiplier.polynomial_multiply(&a, &b).unwrap();

        assert_eq!(result.len(), 4);
        assert_eq!(result[0].to_u64().unwrap(), 4);
        assert_eq!(result[1].to_u64().unwrap(), 13);
        assert_eq!(result[2].to_u64().unwrap(), 22);
        assert_eq!(result[3].to_u64().unwrap(), 15);
    }

    /// Test big integer multiplication using FFT
    #[test]
    fn test_big_int_multiply() {
        #[cfg(feature = "std")]
        let multiplier = FFTMultiplier::new();

        #[cfg(not(feature = "std"))]
        let multiplier = FFTMultiplier;

        let a = BigInt::from_u64(123);
        let b = BigInt::from_u64(456);
        let expected = BigInt::from_u64(56088);

        let result = multiplier.multiply_big_ints(&a, &b).unwrap();
        assert_eq!(result, expected);
    }

    /// Test standalone FFT multiplication function
    #[test]
    fn test_fft_standalone_fallback() {
        let a = BigInt::from_u64(123);
        let b = BigInt::from_u64(456);
        let expected = BigInt::from_u64(56088);

        let result = mul_fft_standalone(&a, &b);
        assert_eq!(result, expected);
    }
}