Piyush Chauhan

Posted on Dec 14

Multiplying Large Decimal Numbers Using Fast Fourier Transform (FFT)

#javascript #algorithms #programming #multiplication

Introduction

Multiplying large decimal numbers can be computationally challenging, especially when dealing with numbers that have many digits or multiple decimal places. Traditional multiplication methods become inefficient for extremely large numbers. This is where the Fast Fourier Transform (FFT) comes to the rescue, providing a powerful and efficient algorithm for multiplying large numbers with remarkable speed.

Application in Multiplication

FFT enables fast multiplication of polynomials or large integers by transforming the numbers into the frequency domain, performing pointwise multiplication, and then applying the inverse FFT. Learn more about how this is implemented in actual production code by exploring the GMPLIB implementation here.

The Challenge of Large Number Multiplication

Traditional multiplication methods have a time complexity of O(n²), where n is the number of digits. For very large numbers, this becomes computationally expensive. The FFT-based multiplication algorithm reduces this complexity to O(n log n), making it significantly faster for large numbers.

Proof Outline for Cooley-Tukey FFT

Decomposition of the Discrete Fourier Transform (DFT):
- The DFT is defined as: $X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-2\pi i \cdot kn / N},$ where $N$ is the size of the input signal.
- The Cooley-Tukey FFT breaks the computation into smaller DFTs of size $N/2$ by separating even-indexed and odd-indexed terms: $X_k = \sum_{n=0}^{N/2-1} x_{2n} \cdot e^{-2\pi i \cdot (2n)k / N} + \sum_{n=0}^{N/2-1} x_{2n+1} \cdot e^{-2\pi i \cdot (2n+1)k / N}.$
- This reduces to: $X_k = \text{DFT of even terms} + W_k \cdot \text{DFT of odd terms},$ where $W_k = e^{-2\pi i \cdot k / N}$ .
Recursive Structure:
- Each DFT of size $N$ is split into two DFTs of size $N/2$ , leading to a recursive structure.
- This recursive division continues until the base case of size $N = 1$ , at which point the DFT is simply the input value.
Butterfly Operations:
- The algorithm merges results from smaller DFTs using the butterfly operations: $a' = u + W_k \cdot v, \quad b' = u - W_k \cdot v,$ where $u$ and $v$ are results from smaller DFTs and $W_k$ represents the roots of unity.
Bit-Reversal Permutation:
- The input array is reordered based on the binary representation of indices to enable in-place computation.
Time Complexity:
- At each level of recursion, there are $N$ computations involving roots of unity, and the depth of the recursion is $\log_2(N)$ .
- This yields a time complexity of $O(N \log N)$ .

Inverse FFT

The inverse FFT is similar but uses $e^{2\pi i \cdot kn / N}$ as the basis and scales the result by $1/N$ .

Understanding the FFT Multiplication Algorithm

The FFT multiplication algorithm works through several key steps:

Preprocessing the Numbers
- Convert the input numbers to arrays of digits
- Handle both integer and decimal parts
- Pad the arrays to the nearest power of 2 for FFT computation
Fast Fourier Transform
- Convert the number arrays into the frequency domain using FFT
- This transforms the multiplication problem into a simpler pointwise multiplication in the frequency domain
Frequency Domain Multiplication
- Perform element-wise multiplication of the transformed arrays
- Utilize complex number operations for efficient computation
Inverse FFT and Result Processing
- Transform the multiplied array back to the time domain
- Handle digit carries
- Reconstruct the final decimal number

Key Components of the Implementation

Complex Number Representation

class Complex {
  constructor(re = 0, im = 0) {
    this.re = re;  // Real part
    this.im = im;  // Imaginary part
  }

  // Static methods for complex number operations
  static add(a, b) { /* ... */ }
  static subtract(a, b) { /* ... */ }
  static multiply(a, b) { /* ... */ }
}

The Complex class is crucial for performing FFT operations, allowing us to manipulate numbers in both real and imaginary domains.

Fast Fourier Transform Function

function fft(a, invert = false) {
  // Bit reversal preprocessing
  // Butterfly operations in frequency domain
  // Optional inverse transformation
}

The FFT function is the core of the algorithm, transforming numbers between time and frequency domains efficiently.

Handling Decimal Numbers

The implementation includes sophisticated logic for handling decimal numbers:

Separating integer and decimal parts
Tracking total decimal places
Reconstructing the result with the correct decimal point placement

Example Use Cases

// Multiplying large integers
fftMultiply("12345678901234567890", "98765432109876543210")

// Multiplying very large different size integers
fftMultiply("12345678901234567890786238746872364872364987293795843790587345", "9876543210987654321087634875782369487239874023894")

// Multiplying decimal numbers
fftMultiply("123.456", "987.654")

// Handling different decimal places
fftMultiply("1.23", "45.6789")

// Handling different decimal places with large numbers
fftMultiply("1234567890123456789078623874687236487236498.7293795843790587345", "98765432109876543210876348757823694.87239874023894")

Performance Advantages

Time Complexity: O(n log n) compared to O(n²) of traditional methods
Precision: Handles extremely large numbers with multiple decimal places
Efficiency: Significantly faster for large number multiplications

Limitations and Considerations

Requires additional memory for complex number representations
Precision can be affected by floating-point arithmetic
More complex implementation compared to traditional multiplication

Conclusion

The FFT multiplication algorithm represents a powerful approach to multiplying large numbers efficiently. By leveraging frequency domain transformations, we can perform complex mathematical operations with remarkable speed and precision.

Practical Applications

Scientific computing
Financial calculations
Cryptography
Large-scale numerical simulations

Code

The complete implementation is following, providing a robust solution for multiplying large decimal numbers using the Fast Fourier Transform approach.

/**
 * Fast Fourier Transform (FFT) implementation for decimal multiplication
 * @param {number[]} a - Input array of real numbers
 * @param {boolean} invert - Whether to perform inverse FFT
 * @returns {Complex[]} - Transformed array of complex numbers
 */
class Complex {
  constructor(re = 0, im = 0) {
    this.re = re;
    this.im = im;
  }

  static add(a, b) {
    return new Complex(a.re + b.re, a.im + b.im);
  }

  static subtract(a, b) {
    return new Complex(a.re - b.re, a.im - b.im);
  }

  static multiply(a, b) {
    return new Complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re);
  }
}

function fft(a, invert = false) {
  let n = 1;
  while (n < a.length) n <<= 1;
  a = a.slice(0);
  a.length = n;

  const angle = ((2 * Math.PI) / n) * (invert ? -1 : 1);
  const roots = new Array(n);
  for (let i = 0; i < n; i++) {
    roots[i] = new Complex(Math.cos(angle * i), Math.sin(angle * i));
  }

  // Bit reversal
  for (let i = 1, j = 0; i < n; i++) {
    let bit = n >> 1;
    for (; j & bit; bit >>= 1) {
      j ^= bit;
    }
    j ^= bit;
    if (i < j) {
      [a[i], a[j]] = [a[j], a[i]];
    }
  }

  // Butterfly operations
  for (let len = 2; len <= n; len <<= 1) {
    const halfLen = len >> 1;
    for (let i = 0; i < n; i += len) {
      for (let j = 0; j < halfLen; j++) {
        const u = a[i + j];
        const v = Complex.multiply(a[i + j + halfLen], roots[(n / len) * j]);
        a[i + j] = Complex.add(u, v);
        a[i + j + halfLen] = Complex.subtract(u, v);
      }
    }
  }

  if (invert) {
    for (let i = 0; i < n; i++) {
      a[i].re /= n;
      a[i].im /= n;
    }
  }

  return a;
}

/**
 * Multiply two decimal numbers using FFT
 * @param {string} num1 - First number as a string
 * @param {string} num2 - Second number as a string
 * @returns {string} - Product of the two numbers
 */
function fftMultiply(num1, num2) {
  // Handle zero cases
  if (num1 === "0" || num2 === "0") return "0";

  // Parse and separate integer and decimal parts
  const parseNumber = (numStr) => {
    const [intPart, decPart] = numStr.split(".");
    return {
      intPart: intPart || "0",
      decPart: decPart || "",
      totalDecimalPlaces: (decPart || "").length,
    };
  };

  const parsed1 = parseNumber(num1);
  const parsed2 = parseNumber(num2);

  // Combine numbers removing decimal point
  const combinedNum1 = parsed1.intPart + parsed1.decPart;
  const combinedNum2 = parsed2.intPart + parsed2.decPart;

  // Total decimal places
  const totalDecimalPlaces =
    parsed1.totalDecimalPlaces + parsed2.totalDecimalPlaces;

  // Convert to digit arrays (least significant first)
  const a = combinedNum1.split("").map(Number).reverse();
  const b = combinedNum2.split("").map(Number).reverse();

  // Determine result size and pad
  const resultSize = a.length + b.length;
  const fftSize = 1 << Math.ceil(Math.log2(resultSize));

  // Pad input arrays
  while (a.length < fftSize) a.push(0);
  while (b.length < fftSize) b.push(0);

  // Convert to complex arrays
  const complexA = a.map((x) => new Complex(x, 0));
  const complexB = b.map((x) => new Complex(x, 0));

  // Perform FFT
  const fftA = fft(complexA);
  const fftB = fft(complexB);

  // Pointwise multiplication in frequency domain
  const fftProduct = new Array(fftSize);
  for (let i = 0; i < fftSize; i++) {
    fftProduct[i] = Complex.multiply(fftA[i], fftB[i]);
  }

  // Inverse FFT
  const product = fft(fftProduct, true);

  // Convert back to integer representation
  const result = new Array(resultSize).fill(0);
  for (let i = 0; i < resultSize; i++) {
    result[i] = Math.round(product[i].re);
  }

  // Handle carries
  for (let i = 0; i < result.length - 1; i++) {
    if (result[i] >= 10) {
      result[i + 1] += Math.floor(result[i] / 10);
      result[i] %= 10;
    }
  }

  // Remove leading zeros and convert to string
  while (result.length > 1 && result[result.length - 1] === 0) {
    result.pop();
  }

  // Insert decimal point
  const resultStr = result.reverse().join("");
  if (totalDecimalPlaces === 0) {
    return resultStr;
  }

  // Handle case where result might be shorter than decimal places
  if (resultStr.length <= totalDecimalPlaces) {
    return "0." + "0".repeat(totalDecimalPlaces - resultStr.length) + resultStr;
  }

  // Insert decimal point
  return (
    resultStr.slice(0, -totalDecimalPlaces) +
    "." +
    resultStr.slice(-totalDecimalPlaces).replace(/0+$/, "")
  );
}

Output

// Example Usage - Self verify using Python
console.log(
  "Product of integers:",
  fftMultiply("12345678901234567890", "98765432109876543210")
);
console.log("Product of decimals:", fftMultiply("123.456", "987.654"));
console.log("Product of mixed decimals:", fftMultiply("12.34", "56.78"));
console.log(
  "Product with different decimal places:",
  fftMultiply("1.23", "45.6789")
);
console.log(
  "Product with large integers:",
  fftMultiply(
    "12345678901234567890786238746872364872364987293795843790587345",
    "9876543210987654321087634875782369487239874023894"
  )
);
const num1 = "1234567890123456789078623874687236487236498.7293795843790587345";
const num2 = "98765432109876543210876348757823694.87239874023894";
console.log("Product:", fftMultiply(num1, num2));

Product of integers: 1219326311370217952237463801111263526900
Product of decimals: 121931.812224
Product of mixed decimals: 700.6652
Product with different decimal places: 56.185047
Product with large integers: 121932631137021795232593613105722759976860134207381319681901040774443113318245930967231822167723255326824021430
Product: 121932631137021795232593613105722759976860134207381319681901040774443113318245.93096723182216772325532682402143

DEV Community