使用快速傅里叶变换 (FFT) 乘以大十进制数

介绍

将大十进制数相乘可能在计算上具有挑战性，尤其是在处理具有许多位数或多个小数位的数字时。传统的乘法方法对于极大的数字来说效率很低。这就是快速傅里叶变换 (FFT) 发挥作用的地方，它提供了一种强大而高效的算法，可以以惊人的速度进行大数相乘。

乘法中的应用

• FFT 通过将数字变换到频域、执行逐点乘法，然后应用逆 FFT 来实现多项式或大整数的快速乘法。

大数乘法的挑战

传统乘法方法的时间复杂度为 O(n²)，其中 n 是位数。对于非常大的数字，这在计算上变得昂贵。基于 FFT 的乘法算法将这种复杂性降低到 O(n log n)，使得处理大数的速度显着加快。

Cooley-Tukey FFT 的证明大纲

1. 离散傅立叶变换 (DFT) 的分解:

   • DFT 定义为：
     $$X_k=\underset{n=0}{\overset{N-1}{\Sigma }}{x}_n\cdot e^{-2\pi i\cdot kn\mathrm{/}N},$$Xk ​= n=0ΣN−1​ xn​⋅ e−2πi⋅kn /N,
     在哪里 $N$N 是输入信号的大小。
   • Cooley-Tukey FFT 将计算分解为更小的 DFT $N\mathrm{/}2$N/2 通过分离偶数索引和奇数索引项：
     $$X_k=\underset{n=0}{\overset{N\mathrm{/}2-1}{\Sigma }}{x}_{2n}\cdot e^{-2\pi i\cdot (2n)k\mathrm{/}N}\underset{n=0}{\overset{N\mathrm{/}2-1}{\Sigma }}{x}_{2n1}\cdot e^{-2\pi i\cdot (2n1)k\mathrm{/}N}\mathrm{。}$$Xk ​= n=0ΣN/2−1​ x2n​ ⋅e −2πi⋅(2n)k/N n=0ΣN/ 2−1​x2n 1​⋅ e−2πi⋅(2n 1)k/N.
   • 这减少为：
     $$X_k=\text{偶数项的 DFT}{W}_k\cdot \text{奇数项的 DFT},$$Xk​ =偶数项的 DFT Wk​⋅奇数项的 DFT,
     在哪里 $W_k=e^{-2\pi i\cdot k\mathrm{/}N}$Wk​ =e−2πi ⋅k/N .

2. 递归结构:

   • 每个大小的 DFT $N$N 被分成两个大小的 DFT $N\mathrm{/}2$N/2 ，导致递归结构。
   • 这种递归除法一直持续到大小的基本情况 $N=1$N=1 ，此时 DFT 只是输入值。

3. 蝴蝶操作:

   • 该算法使用蝶形运算合并较小 DFT 的结果：
     $$a^{\mathrm{\prime }}=u{W}_k\cdot v,b^{\mathrm{\prime }}=u-W_k\cdot v,$$一个′=u Wk​⋅v,b′ =u−Wk​⋅v,
     在哪里 $u$u 和 $v$v 是较小 DFT 的结果， $W_k$Wk​ 代表团结的根源。

4. 位反转排列:

   • 输入数组根据索引的二进制表示重新排序，以实现就地计算。

5. 时间复杂度:

   • 在每一层递归中，都有 $N$N 涉及单位根的计算，递归的深度是 ${\mathrm{日志}}_2(N)$log2 ​(N) .
   • 这会产生时间复杂度 $O(N\log N）$O(NlogN) .

逆FFT

• 逆 FFT 类似，但使用 $e^{2\pi i\cdot kn\mathrm{/}N}$e 2πi⋅kn/N 作为基础并对结果进行缩放 $1\mathrm{/}N$1/N .

了解 FFT 乘法算法

FFT 乘法算法的工作原理有几个关键步骤：

1. 预处理数字

   • 将输入数字转换为数字数组
   • 同时处理整数和小数部分
   • 将数组填充到最接近的 2 次幂以进行 FFT 计算

2. 快速傅立叶变换

   • 使用 FFT 将数字数组转换为频域
   • 这将乘法问题转化为频域中更简单的逐点乘法

3. 频域乘法

   • 对转换后的数组执行逐元素乘法
   • 利用复数运算进行高效计算

4. 逆FFT和结果处理

   • 将相乘后的数组变换回时域
   • 处理数字进位
   • 重构最终的十进制数

实施的关键组成部分

复数表示

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) { /* ... */ }
}

Complex 类对于执行 FFT 运算至关重要，它使我们能够操作实域和虚域中的数字。

快速傅立叶变换函数

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

FFT函数是算法的核心，有效地在时域和频域之间进行数字转换。

处理小数

该实现包括处理十进制数的复杂逻辑：

• 分隔整数和小数部分
• 跟踪小数点总数
• 使用正确的小数点位置重建结果

示例用例

// 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")

性能优势

• 时间复杂度：O(n log n) 与传统方法的 O(n²) 相比
• 精度：处理具有多个小数位的极大数字
• 效率：大量乘法运算速度明显加快

限制和注意事项

• 需要额外的内存来表示复数
• 浮点运算会影响精度
• 与传统乘法相比，实现更复杂

结论

FFT 乘法算法代表了一种有效地乘以大数的强大方法。通过利用频域变换，我们可以以惊人的速度和精度执行复杂的数学运算。

实际应用

• 科学计算
• 财务计算
• 密码学
• 大规模数值模拟

进一步阅读

• Cooley-Tukey FFT 算法
• 数论
• 计算数学

代码

完整的实现如下，为使用快速傅里叶变换方法乘以大十进制数提供了一个强大的解决方案。

/**
 * 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+$/, "")
  );
}

输出

// 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

以上是使用快速傅里叶变换 (FFT) 乘以大十进制数的详细内容。更多信息请关注PHP中文网其他相关文章！