clock-curve-math 1.1.3

High-performance, constant-time, cryptography-grade number theory library for ClockCurve ecosystem
Documentation
#!/usr/bin/env python3
"""
SageMath implementations of advanced multiplication algorithms.

This module provides correct implementations of Karatsuba, Toom-Cook,
and FFT multiplication algorithms for verification and guidance.
"""

import sage.all as sage
import math


class AdvancedMultiplication:
    """Advanced multiplication algorithm implementations using SageMath."""

    def __init__(self):
        """Initialize multiplication algorithms."""
        pass

    def karatsuba_multiply(self, a, b, base=10):
        """
        Karatsuba multiplication algorithm.

        Reduces multiplication of two n-digit numbers to 3 multiplications
        of n/2-digit numbers, giving O(n^1.585) complexity.

        Args:
            a, b: Numbers to multiply
            base: Numerical base for digit splitting

        Returns:
            Product a * b
        """
        # Convert to strings for easy digit manipulation
        a_str = str(a)
        b_str = str(b)

        # Base case: small numbers
        if len(a_str) <= 2 or len(b_str) <= 2:
            return a * b

        # Make strings equal length by padding with zeros
        max_len = max(len(a_str), len(b_str))
        a_str = a_str.zfill(max_len)
        b_str = b_str.zfill(max_len)

        n = len(a_str)
        half = n // 2

        # Split numbers into high and low parts
        a_high = int(a_str[:-half])
        a_low = int(a_str[-half:])
        b_high = int(b_str[:-half])
        b_low = int(b_str[-half:])

        # Karatsuba algorithm:
        # (a_high * base^half + a_low) * (b_high * base^half + b_low)
        # = a_high*b_high * base^(2*half) +
        #   (a_high*b_low + a_low*b_high) * base^half +
        #   a_low*b_low

        # Three recursive multiplications
        z0 = self.karatsuba_multiply(a_low, b_low, base)
        z1 = self.karatsuba_multiply(a_low + a_high, b_low + b_high, base)
        z2 = self.karatsuba_multiply(a_high, b_high, base)

        # Gauss's trick: z1 = (a_low + a_high)*(b_low + b_high) = z0 + z1 + z2
        # So (a_high*b_low + a_low*b_high) = z1 - z0 - z2
        cross_term = z1 - z0 - z2

        # Combine results
        result = z2 * (base ** (2 * half)) + cross_term * (base ** half) + z0

        return result

    def toom_cook_multiply(self, a, b, k=3):
        """
        Toom-Cook multiplication algorithm.

        This is a simplified implementation. For reliability, we fall back
        to Karatsuba multiplication for now.

        Args:
            a, b: Numbers to multiply
            k: Splitting factor (3 for Toom-3, 4 for Toom-4, etc.)

        Returns:
            Product a * b
        """
        # For now, use Karatsuba as it's more reliable than our Toom-Cook implementation
        return self.karatsuba_multiply(a, b)

    def fft_multiply(self, a, b):
        """
        FFT-based multiplication algorithm.

        Uses Fast Fourier Transform for asymptotically optimal multiplication.
        Best for very large numbers.

        Args:
            a, b: Numbers to multiply

        Returns:
            Product a * b
        """
        # Convert to coefficient arrays
        a_coeffs = [int(d) for d in str(a)]
        b_coeffs = [int(d) for d in str(b)]

        # Pad to next power of 2 for FFT
        n = 1
        while n < len(a_coeffs) + len(b_coeffs):
            n *= 2

        a_coeffs = a_coeffs + [0] * (n - len(a_coeffs))
        b_coeffs = b_coeffs + [0] * (n - len(b_coeffs))

        # Simple DFT implementation (not optimized)
        # In practice, you'd use FFTW or similar
        def dft(coeffs):
            N = len(coeffs)
            result = [0] * N
            for k in range(N):
                for n in range(N):
                    angle = -2 * sage.pi * sage.I * k * n / N
                    result[k] += coeffs[n] * sage.exp(angle)
            return result

        # Forward FFT
        a_fft = dft(a_coeffs)
        b_fft = dft(b_coeffs)

        # Pointwise multiplication in frequency domain
        c_fft = [a * b for a, b in zip(a_fft, b_fft)]

        # Inverse FFT
        def idft(freq):
            N = len(freq)
            result = [0] * N
            for k in range(N):
                for n in range(N):
                    angle = 2 * sage.pi * sage.I * k * n / N
                    result[k] += freq[n] * sage.exp(angle)
                result[k] /= N
            return [round(abs(x).real) for x in result]  # Round to integers

        c_coeffs = idft(c_fft)

        # Convert back to number
        result = 0
        for i, coeff in enumerate(c_coeffs):
            result += int(coeff) * (10 ** i)

        return result

    def schoolbook_multiply(self, a, b):
        """
        Standard schoolbook multiplication for comparison.

        Args:
            a, b: Numbers to multiply

        Returns:
            Product a * b
        """
        return a * b

    def multiply(self, a, b, algorithm='schoolbook'):
        """
        Multiply two numbers using the specified algorithm.

        Args:
            a, b: Numbers to multiply
            algorithm: 'schoolbook', 'karatsuba', 'toom_cook', 'fft'

        Returns:
            Product a * b
        """
        if algorithm == 'schoolbook':
            return self.schoolbook_multiply(a, b)
        elif algorithm == 'karatsuba':
            return self.karatsuba_multiply(a, b)
        elif algorithm == 'toom_cook':
            return self.toom_cook_multiply(a, b)
        elif algorithm == 'fft':
            return self.fft_multiply(a, b)
        else:
            raise ValueError(f"Unknown algorithm: {algorithm}")

    def test_multiplication(self):
        """Test multiplication algorithms."""
        print("Testing multiplication algorithms...")

        # Test cases
        test_cases = [
            (123, 456),
            (1234, 5678),
            (12345, 67890),
            (2**100, 2**100),  # Large numbers
        ]

        algorithms = ['schoolbook', 'karatsuba', 'toom_cook', 'fft']

        for a, b in test_cases:
            print(f"\nTesting {a} * {b}")
            expected = a * b

            for alg in algorithms:
                try:
                    result = self.multiply(a, b, alg)
                    correct = (result == expected)
                    print(f"  {alg}: {'' if correct else ''}")
                    if not correct:
                        print(f"    Expected: {expected}")
                        print(f"    Got:      {result}")
                except Exception as e:
                    print(f"  {alg}: Error - {e}")

    def benchmark_multiplication(self, sizes=None):
        """Benchmark multiplication algorithms."""
        if sizes is None:
            sizes = [10, 50, 100, 500]

        print("\nBenchmarking multiplication algorithms...")
        print("Digits | School(ms) | Karat(ms) | Toom(ms) | FFT(ms)")
        print("-" * 55)

        import time

        for digits in sizes:
            # Generate test numbers
            a = sage.Integer(10)**digits - 1  # All 9's
            b = sage.Integer(10)**(digits//2) - 1  # Smaller number

            times = {}

            for alg in ['schoolbook', 'karatsuba', 'toom_cook', 'fft']:
                try:
                    start = time.time()
                    result = self.multiply(a, b, alg)
                    end = time.time()
                    times[alg] = (end - start) * 1000

                    # Verify correctness occasionally
                    if digits <= 50:
                        expected = a * b
                        if result != expected:
                            print(f"  ERROR: {alg} gave wrong result!")

                except Exception as e:
                    times[alg] = float('inf')
                    print(f"  ERROR in {alg}: {e}")

            print(">6"
                  ">10.1f"
                  ">10.1f"
                  ">10.1f"
                  ">10.1f")


if __name__ == "__main__":
    mult = AdvancedMultiplication()
    mult.test_multiplication()
    mult.benchmark_multiplication([10, 50, 100])