傅里叶变换与乘法
基本上所有的理工科都学过傅里叶变换，然而想深入理解傅里叶变换可不容易。就我自己而言，我一直只是把傅里叶变换及其逆变换当作信号在时域与频域之间等效变换的数学原理。后来学习的拉普拉斯变换也只被我当成信号在时域与复频域之间的等效变换的数学原理。 可能是因为要学的东西太多了（其实是自己笨和懒），我从没有去认真消化那些精巧的定理证明，也没有仔细琢磨那些定理到底有什么非凡意义。 一次偶然的机会让我知道原来傅里叶变换还被应用在大数乘法的算法中，最神奇的是基于傅里叶变换的FFT算法还是计算大数乘法最快的算法。这里的大数是指长度远远大于计算机硬件直接支持的字长的数。比如64位计算机可以直接支持2^64 以内的整数的直接相乘，但是对于更长的数（比如10^500以内）就不能直接计算而需要软件工程师使用专门的算法了。 基于FFT的乘法计算的复杂度为：$$O(n*\log(n))$$，而通常采用的Karatsuba算法的复杂度为：$$O(n^{1.585})$$这个发现让我意识到自己的知识局限性，并解开了一个纠缠多年的心结。
那是大一下学期一个平凡的周五晚上的最后一节课，有人在心中计划着周末的情人坡幽会而无心听课，有人认为课堂内容过于简单而把对老师拖堂和啰里啰唆的不满写在脸上，也有学霸坐着第一排埋头苦读看着复变函数。看着一群趾高气昂不可一世的学生，突然老师提起了计算机学院的online judge平台。她轻蔑的话语：“就算是最简单的问题也可以把你们难得屁滚尿流”并没有打击到大部分同学的自信心。术业有专攻怎么显然的道理大家都懂，这么明显的激将之言就是想让人去计算机的主场踢馆。按理说，情商上线的人是不会上当的。可怜我们寝室的几个单身狗，虽然记得高数老师的那句不要拿自己的兴趣爱好去挑战人家吃饭的本事，但是还是觉得就算技不如人，也能坚持几招而不败，怎么可能被人秒杀呢。于是那个周末，扯着高手在民间的口号，便去了人家的主场卖弄那点三脚猫功夫。
SOJ的第一关是求两个大数的和。多大数呢？10的500次方以内。看完题目，寝室的四条单身狗相视一笑：“呵呵，这有什么难的？”。从低位到高位，每位逐一相加在加上低位的进数，结果就出来了。动手码砖然后提交，顺利通关。四条单身狗瞬间自信心爆棚，“灭绝师太也太小看我们了吧，最简单的也不过小学三年级的水平”，没有注意到屏幕右下方上出现了一行小字“热身完毕，你的代码运行速度超过了4%的同学”，直接点击了下一关。
SOJ的第二关是求两个大数的积。多大数呢？还是10的500次方以内。同样的剧情不同的结局：超时。大家马上想到了不直接用字符串进行计算，而是转成整数数组。大家修改后再次提交，等待良久后，才看到结果是Accepted。可是这次屏幕右下方的小字“你的代码运行速度超过了0.3%的同学，请继续努力”总算引起了大家的注意。好奇的点击过后发现，我们的代码居然用了接近三分钟？第一名居然只用4毫秒?什么大部分挑战者的代码也只要300到500毫秒？而且第一名是一个学长在10年前提交的答案。那个下午大家激烈的讨论起哪些地方可以优化？不断尝试修改，可是运行时间还是要三分钟左右。那个夜晚的卧谈会总算换了个主题，一帮单身狗难得的没有探讨中日友好交流等其他类似话题， 在根本不知道什么是时间复杂度的时候瞎讨论如何优化一个算法。
后来我们通过查找资料了解到了Karatsuba算法，也初步了解了时间复杂度。我们采用Karatsuba算法，把时间缩短到几百毫秒的尺度，却怎么也达不到毫秒级的速度。这次偶然的发现，让我感觉得第一名的学长应该是采用的FFT算法。可惜我早已毕业，不能访问学校的Online judge 平台了。
理解FFT算法计算大数乘法需要一些数学基础。傅里叶变换是理工类本科的基础知识，他们的形式化定以、相关定理的推导和证明是现代科技的重要基础之一，是人类文明的科技树上的重要枝干。我建议找一本高数和信号与系统的教科书熟悉一下相关知识。我在这里记录一些基本公式作为自己的备忘。最初的傅里叶变换是连续时域到连续频域的变换：
$$F(\omega) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty f(t) e^{- i\omega t}\,dt$$
该变换是把一组函数映射为另一组函数的线性算子。注意这两组函数的定义域都是连续且无穷的。基于采用定理，可以推出离散且有限的离散傅里叶变换（DFT）。
$$\hat{x}[k]=\sum_{n=0}^{N-1} e^{-i\frac{2\pi}{N}nk}x[n] \qquad k = 0,1,\ldots,N-1.$$
直接计算DFT的复制度为O(n^2)，后来有人发明了FFT算法。FFT算法能够以O(n*log(n))的时间复制度来计算DFT。FFT的细节也很容易找到资料。
上述的数学基础我也曾走马观花的学过一遍，后来在工作过程中也经常复习以解决自己的一些疑惑。可是我怎么也想不到傅里叶变换还可以用来计算大数乘法。感兴趣的人可以在这里看看怎么样用FFT来计算大数乘法。
FFT计算大数乘法的过程是把两个大数看成两个向量分别与权重向量的点积，先把两个向量进行傅里叶变换，得到变换后的两个新向量，然后把两个新向量的元素直接相乘得到频域的乘积向量，再然后把乘积向量进行逆变换得到时域的卷积，最后把时域的卷积与权重向量进行点积运算得到最终结果。依据卷积定理，时域卷积对偶于频域乘积，这里的时域是傅里叶变换前的函数的定义域，频域傅里叶变换后的函数的定义域。
比如123456789 是[1,2,3,4,5,7,8,9] 与[10^8,10^7,10^6,10^5,10^4,10^3,10^2,10^1,10^0]的点积。把计算123456789×987654321当作先计算x:[1,2,3,4,5,6,7,8,9] 与y:[9,8,7,6,5,4,3,2,1]的卷积：z:[9, 26, 50, 80, 115, 154, 196, 240, 285, 240, 196, 154, 115,80, 50, 26, 9],最后再计算z与[10^16,10^15,10^14,10^13,10^12,10^11,10^10,10^9,10^8,10^7,10^6,10^5,10^4,10^3,10^2,10^1,10^0]的点积，即可得到结果121932631112635269。注意在实际应用中n必须是2的幂，并且傅里叶变换会把整数转换成浮点数，要考虑到浮点数运算过程中的精度损失。
import numpy as np
import random
N=2**13
ws=[10**i for i in range(2*N,0,-1)]
xs=[random.randint(0,9) if i>N else 0 for i in range(0,2*N)]
ys=[random.randint(0,9) if i>N else 0 for i in range(0,2*N)]
x1s=np.fft.rfft(xs)
y1s=np.fft.rfft(ys)
z1s=x1s*y1s
zs=np.fft.irfft(z1s)
zis=map(lambda x:int(x), np.rint(zs).tolist())
x=sum(map(lambda i,j:i*j,xs,ws))
y=sum(map(lambda i,j:i*j,ys,ws))
zi=sum(map(lambda i,j:i*j,zis,ws))
print(f"x={x},\ny={y},\nz={x*y}")
print(f"FFT multiply result:\nz={zi}")上面的代码在python 3.8 中运行测试过，自己用肉眼检查了一下，在10的八千一百九十二位以内的时候没发现问题（通过与Python 内带的大数计算进行对比）。
