颂哈吉-施特拉森演算法

颂哈吉-施特拉森演算法（英语：Schönhage–Strassen algorithm）是渐近快速的大整数乘法算法。是由阿诺德·颂哈吉（英语：Arnold Schönhage）和沃尔克·施特拉森在1971年发明^[1]。若针对二个n位元的整数，其运行的位元复杂度（英语：bit complexity），若以大O符号表示，是${\displaystyle O(n\cdot \log n\cdot \log \log n)}$。演算法使用在有2ⁿ+1个元素的环上的迭代快速傅里叶变换，这是一种特别的数论转换。

颂哈吉-施特拉森演算法演算法是以整数乘法的FFT方法为基础。图中所列的是用单纯FFT法计算1234乘5678 = 7006652的过程。使用了模数337的数论转换，选择85为1的8次方根。为了说明方便，其中数字使用十进制表示，不用二进制表示。颂哈吉-施特拉森以此方式为基础，再加上negacyclic convolutions

颂哈吉（右）和施特拉森（左）在上沃尔法赫数学研究所下西洋棋，1979年

颂哈吉-施特拉森演算法是1971年至2007年之间，渐近最快的乘法演算法，2007年时有一个新的乘法演算法Fürer演算法（英语：Fürer's algorithm），其渐近复杂度较低^[2]，不过Fürer演算法只有在大到天文数字的程度时，其速度才会比颂哈吉-施特拉森演算法快，实务上不会使用Fürer演算法（参考银河式算法）。

在实务上，当数字大于2²¹⁵至2²¹⁷（十进制下的10,000到40,000位数）时，颂哈吉-施特拉森演算法的速度会比较早期的卡拉楚巴算法和图姆-库克算法要快^[3][4][5]。GNU多重精度运算库用这个演算法计算至少1728至7808个64位元字节的数值（十进制下的33,000至150,000位数），依硬体架构而定^[6]，有一个用Java实现的颂哈吉-施特拉森演算法，可以计算十进位超过74,000位数^[7]。

颂哈吉-施特拉森演算法的应用包括数学哲学（例如互联网梅森素数大搜索以及计算圆周率近似值），实务的应用包括克罗内克代入（英语：Kronecker substitution），其中将整系数多项式的乘法有效的简化为大数的乘法，GMP-ECM中用此算法来计算Lenstra椭圆曲线分解（英语：Lenstra elliptic curve factorization）^[8]。

细节

此段会说明颂哈吉-施特拉森演算法实现的细节，主要是依Crandall和Pomerance在《Prime Numbers: A Computational Perspective》书中的简介^[9]。其中说明的和颂哈吉原始的方法有些差异：其中用离散加权转换（discrete weighted transform）来快速计算负循环卷积（英语：negacyclic convolution）。另一个细节资讯的来源是高德纳的《计算机程序设计艺术》^[10]。此章节从建立blocks开始，最后会一步一步的说明演算法。

卷积

假设要将B基底的123和456用直式乘法相乘，B够大，因此计算中不会有进位的情形。结果会是：

		1	2	3
	×	4	5	6
		6	12	18
	5	10	15
4	8	12
4	13	28	27	18

数列(4, 13, 28, 27, 18)称为二个原始数列(1,2,3)和(4,5,6)的线性卷积（linear convolution）或非循环卷积（acyclic convolution）。只要有二个数列的线性卷积，要算出原数字在十进位以下的乘积就很简单了：只需要处理进位（例如，最右栏的数字，只保留8，将1进到上一栏的27）。此例中可以得到正确的乘积56088。

有另外几种也很有用的卷积。假设输入数列有n个元素（假设3个），线性卷积会有n+n−1个元素，若将最右边的n个元素加上最左边的 n−1个元素，可以产生循环卷积（cyclic convolution）：

	28	27	18
+		4	13
	28	31	31

若在循环卷积的结果再考虑其进位，所得结果等效于原来两数字的乘积mod Bⁿ − 1。在此例中，10³ − 1 = 999，将(28, 31, 31)计算其进位后得到 3141，而3141 ≡ 56088 (mod 999)。

相对的，若将最右边的n个元素减去最左边的n−1个元素，会产生负循环卷积（英语：negacyclic convolution）（negacyclic convolution）：

	28	27	18
−		4	13
	28	23	5

若在负循环卷积的结果再考虑其进位，所得结果等效于原来两数字的乘积mod Bⁿ + 1。此例中，10³ + 1 = 1001。将(28, 23, 5)计算进位后得到3035，而3035 ≡ 56088 (mod 1001)。负循环卷积可能会有负数，可以用借位的方式消除其负数。

卷积定理

颂哈吉-施特拉森演算法和其他以快速傅立叶变换为基础的乘法算法一样，在本质上都是以卷积定理为基础，可以快速的计算二个数列的循环卷积，卷积定理如下：

两个向量的循环卷积可以将这两个向量求其离散傅里叶变换（DFT），将所得向量依元素相乘，将所得结果再进行反离散傅里叶变换（IDFT）。

若依符号表示，可表示如下：

CyclicConvolution(X, Y) = IDFT(DFT(X) · DFT(Y))

若用快速傅里叶变换（FFT）来计算DFT及IDFT，再用递回的方式处理乘法演算法，来将DFT(X)和DFT(Y)相乘，就可以得到计算循环卷积的快速演算法。

在此演算法中，若计算负循环卷积会更有效：使用修改后的卷积定理版本（参考数论转换）也可以有相同效果。假设向量X和Y的长度是n，而a是2n阶的单位根（意思是a²ⁿ = 1，其他a的较小幂次都不等于1），可以定义第三个向量A，称为加权向量：

A = (a^j), 0 ≤ j < n

A⁻¹ = (a^−j), 0 ≤ j < n

负循环卷积可以表示为下式：

NegacyclicConvolution(X, Y) = A⁻¹ · IDFT(DFT(A · X) · DFT(A · Y))

这和循环卷积的公式类似，只是输入要先乘以A,输出要乘以A⁻¹。

代数环的选用

离散傅里叶变换是抽象的运算，因此可以在任意环的代数结构下运作。一般而言是在复数下进行，但实际上为了要有准确的结果，要处理复数运算到足够的精度，这种乘法不但慢，而且容易有错误。因此会用数论转换的方式进行，是在整数模N（N为特定整数）的底下进行。

就像在复数平面上，每一阶都可以找到原根一样，给定任何的阶数n，可以选择适当的N，使得b是在整数模N的系统下，n阶的原根（bⁿ ≡ 1 (mod N)，b其他较小的次幂，都不会是1 mod N）。

此演算法会花很多时间演算小数字的次幂，若用一个天真的方式来处理演算法，这会出现许多次。

1. 在FFT演算法中，原根b会一直的计算幂次，平方，并且和其他的数相乘。
2. 在计算原根a的次方，以计算加权向量A时，以及将A或A⁻¹和其他向量相乘的时候。
3. 在进行向量项对项的相乘时。

颂哈吉-施特拉森演算法的关键是选择模数N等于2ⁿ + 1，其中的n是要转换件数P=2^p的倍数。若标准系统，用二进位的形式表示大数值，有许多的好处：

• 所有的数字计算modulo 2ⁿ + 1，都可以只用移位和加法快速求得（这会在下一节解释。）
• 此转换会用到的所有原根都可以写成2^k，因此原根和其他数字的相乘只需要用到移位，原根的平方和幂次都只是在其指数上的变化。
• 转换向量的项对项递回乘法可以用逆循环卷积来计算，这比卷积要快，而且其结果会自动计算mod 2ⁿ + 1。

为了让递回乘法比较方便，会将颂哈吉-施特拉森演算法定位为二个数字的乘积mod 2ⁿ + 1（n为某整数），而不是一般的乘法。这不会失去演算法的一般性，因为可以选择够大的n，使得乘积mod 2ⁿ + 1就是乘积本身。

移位的优化

在此演算法中，有许多和二个次幂相乘或相除（包括所有的原根）的情形是可以用较小数字的移位和相加达成。原因是因为观察到以下的算式：

(2ⁿ)^k ≡ (−1)^k mod (2ⁿ + 1)

其中2ⁿ进制下的k位数，若用位值计数法表示，可以写成(d_k−1,...,d₁,d₀)，代表数值${\displaystyle \scriptstyle \sum _{i=0}^{k-1}d_{i}\cdot (2^{n})^{i}}$，针对每一个d_i可得0 ≤ d_i < 2ⁿ。

因此可以简化表示为二进制数字进行mod 2ⁿ + 1的计数：取最右边（最小位）的n位元，减去较高位的n位元，再加上再高位的n位元，一直计算到所有位元都已处理为止。若所得的数超过0到2ⁿ的范围，可以加或减模数2ⁿ + 1的倍数以进行正规化。例如，若n=3（因此模数是2³+1 = 9），要化简的数字是656，可以得到：

656 = 1010010000₂ ≡ 000₂ − 010₂ + 010₂ − 1₂ = 0 − 2 + 2 − 1 = −1 ≡ 8 (mod 2³ + 1).

甚至，针对很多位元的移位，还有简化的方式。假设数字A介于0到2ⁿ之间，希望乘以2^k。将k除以n可得到k = qn + r，其中r < n。因此可得：

A(2^k) = A(2^{qn + r}) = A[(2ⁿ)^q(2^r)] ≡ (−1)^q 左位移（英语：Shift_Left）(A, r) (mod 2ⁿ + 1).

因为A≤ 2ⁿ，且r < n。左移r位元后最多有2n−1位元，因此只需要一次位移和减法（之后可能要正规化）。

最终，若要除以2^k，可将此段的第一个式子平方，可得：

2²ⁿ ≡ 1 (mod 2ⁿ + 1)

因此

A/2^k = A(2^−k) ≡ A(2^{2n − k}) = Shift_Left(A, 2n − k) (mod 2ⁿ + 1).

演算法

此演算法有分割、评估（前向FFT）、各项相乘、内插（反FFT）以及合并的阶段，类似Karatsuba算法和Toom-Cook算法。

假设输入数字是x和y，以及整数N，以下演算法可以计算xy mod 2^N + 1的结果。假设N够大，使得乘积取模数后的结果仍是乘积本身。

1. 将输入的数字分割为向量X和Y，各有2^k个元素，其中2^k是N的因素（例如将12345678 → (12, 34, 56, 78)）。
2. 为了运算方便，需要用较小的N以进行递回乘法。因此选择n是大于2N/2^k + k，而且可被2^k整除的最小数字。
3. 利用负循环卷积计算XY mod 2ⁿ + 1：
   1. 将X和Y和加权向量A相乘，作法是使用移位（将第j项左移jn/2^k）
   2. 用数论FFT计算X和Y的DFT（将所有乘法用移位表示：针对第2^k次原根，使用2^2n/2k）。
   3. 递回的应用此演算法，将转换后的X和Y对应项相乘。
   4. 计算所得向量的IDFT，得到结果向量C（用移位代替乘法），这对应内插阶段。
   5. 将结果向量C和A⁻¹相乘，用移位来进行。.
   6. 调整符号：结果的一些项可能是负的。可以计算C的第j项的最大可能正值(j + 1)2^2N/2k，若有超过，再减去模数2ⁿ + 1。
4. 最后，处理进位mod 2^N + 1，得到最终的结果。

输入数字最佳的分割组数和${\displaystyle {\sqrt {N}}}$成比例，其中N是输入数字的位元数，此设法会让执行时间的复杂度为O(N log N log log N),^[1][9]，因此需依规则设定参数k。实务上，会依输入数字大小以及演算法架构而定，一般会介于4到16之间^[8]。

在步骤2，已观察到以下的现象：

• 输入向量的每一项最多只有n/2^k个位元。
• 二个输入向量项的乘积最多有2n/2^k个位元。
• 卷积的每一个元素是最多2^k个乘积的和，因此不会超过 2n/2^k + k个位元。
• n一定要可被2^k整除，以确保在步骤1的递回呼叫中2^k 整除N"的条件会成立。

参考资料

1. A. Schönhage and V. Strassen, "Schnelle Multiplikation großer Zahlen （页面存档备份，存于互联网档案馆）", Computing 7 (1971), pp. 281–292.
2. Martin Fürer, "Faster integer multiplication 互联网档案馆的存档，存档日期2013-04-25.", STOC 2007 Proceedings, pp. 57–66.
3. Van Meter, Rodney; Itoh, Kohei M. Fast Quantum Modular Exponentiation. Physical Review. A. 2005, 71 (5): 052320. S2CID 14983569. arXiv:quant-ph/0408006 . doi:10.1103/PhysRevA.71.052320.
4. Overview of Magma V2.9 Features, arithmetic section 互联网档案馆的存档，存档日期2006-08-20.: Discusses practical crossover points between various algorithms.
5. Luis Carlos Coronado García, "Can Schönhage multiplication speed up the RSA encryption or decryption? Archived", University of Technology, Darmstadt (2005)
6. MUL_FFT_THRESHOLD. GMP developers' corner. [3 November 2011]. （原始内容存档于24 November 2010）.
7. An improved BigInteger class which uses efficient algorithms, including Schönhage–Strassen. Oracle. [2014-01-10].
8. Pierrick Gaudry, Alexander Kruppa, and Paul Zimmermann. A GMP-based Implementation of Schönhage–Strassen’s Large Integer Multiplication AlgorithmArchived. Proceedings of the 2007 International Symposium on Symbolic and Algebraic Computation, pp.167–174.
9. R. Crandall & C. Pomerance. Prime Numbers – A Computational Perspective. Second Edition, Springer, 2005. Section 9.5.6: Schönhage method, p. 502. ISBN 0-387-94777-9
10. Donald E. Knuth, The Art of Computer Programming（计算机程序设计艺术）, Volume 2: Seminumerical Algorithms (3rd Edition), 1997. Addison-Wesley Professional, ISBN 0-201-89684-2. Section 4.3.3.C: Discrete Fourier transforms, pg.305.