炒冷饭，一道16年前的北美OUG SQL挑战赛回顾 - l1t1/note GitHub Wiki
最近研究高精度计算，其中频繁用到了快速傅里叶转换(FFT), 想起16年前我参与讨论的一个帖子，利用ITPUB论坛查询功能检索到了：首届NoCOUG国际SQL挑战赛。
冠军得主用的就是FFT，当时不理解，现在理解得稍微多点了。
题目如下：
有一张表记录了骰子的六个面的面值极其出现概率：
CREATE TABLE die (
face_id NUMBER(2) NOT NULL CHECK (face_id > 0) PRIMARY KEY,
face_value NUMBER(2) NOT NULL CHECK (face_value > 0),
probability NUMBER(10,10) NOT NULL CHECK (probability >=0 AND probability <= 1)
);
INSERT INTO die VALUES (1, 1, 1/6 + 1/12);
INSERT INTO die VALUES (2, 3, 1/6 + 1/12);
INSERT INTO die VALUES (3, 4, 1/6 + 1/12);
INSERT INTO die VALUES (4, 5, 1/6 - 1/12);
INSERT INTO die VALUES (5, 6, 1/6 - 1/12);
INSERT INTO die VALUES (6, 8, 1/6 - 1/12);
注意和传统的1-6骰子不同，概率也不是平均出现的。
现在要用SQL求出在掷N次后各个总面值及其概率。
这题可以用简单地两表交叉连接（笛卡尔积）表示投掷2次，以此类推。所以难点在于怎么用固定SQL编写适应动态N的算法以及当N很大时怎么提高效率。当时还有一个额外的挑战，出题时（2029年3月）Oracle最高版本11.1还不支持递归CTE（https://blog.itpub.net/17086096/viewspace-1967966/）。2010年newkid在书里还专门写了一章介绍递归CTE(Oracle叫子查询因子化)。
帖主newkid的解法是这样的，注释是我加的。
VAR N NUMBER;
EXEC :N :=2; /*绑定变量N表示投掷次数*/
WITH d AS ( /*将die表自连接N层模拟投N次, 将第N层的各次投掷结果组合用字符串拼接起来，分隔符是~ ，1行表示一种组合*/
SELECT ROWNUM AS id
,SYS_CONNECT_BY_PATH(face_value,'~')||'~' AS v
,SYS_CONNECT_BY_PATH(probability,'~')||'~' AS p
FROM die
WHERE LEVEL=:N
CONNECT BY LEVEL<=:N
)
,d2 AS ( /*用一个1-N的表去拆分字符串，把表示面值的字符串拆分后加总，表示概率的字符串拆分后先求自然对数再加总再用它算e的次方，相当于相乘*/
SELECT id
,SUM(SUBSTR(v,INSTR(V,'~',1,l)+1,INSTR(V,'~',1,l+1)-INSTR(V,'~',1,l)-1)) AS v
,EXP(SUM(LN(SUBSTR(p,INSTR(p,'~',1,l)+1,INSTR(p,'~',1,l+1)-INSTR(p,'~',1,l)-1)))) AS p
FROM d
,(SELECT LEVEL as l FROM DUAL CONNECT BY LEVEL<=:N)
GROUP BY id
)
SELECT v /* 按总面值分组汇总 */
,SUM(p) AS p
FROM d2
GROUP BY v
ORDER BY 1;
观察题目，次数相加，系数相乘是否和多项式乘法相似，同一个多项式多次与自己相乘就是乘方。把面值作为次数，概率作为系数，原题就变成多项式乘方。
现在是2025年5月，numpy有多项式乘方计算函数polypow，第一个参数是对应次数的系数列表，从小到大排列，第一项代表常数项。结果是个数组，其元素是积的系数列表，同样是第一项代表常数项，翻译过来就是面值2（第3项）系数为0.0625, 以此类推。
>>> from numpy.polynomial import polynomial as P
>>> d=[0, 1/4, 0,1/4,1/4,1/12,1/12,0,1/12]
>>> P.polypow(d, 2)
array([0. , 0. , 0.0625 , 0. , 0.125 ,
0.125 , 0.10416667, 0.16666667, 0.10416667, 0.125 ,
0.04861111, 0.05555556, 0.04861111, 0.01388889, 0.01388889,
0. , 0.00694444])
多项式两次方的结果与如下sql的计算结果一致：
SELECT face_value, SUM (probability) AS probability
FROM
(
SELECT
d1.face_value + d2.face_value AS face_value,
d1.probability * d2.probability AS probability
FROM die d1 CROSS JOIN die d2
)
GROUP BY face_value
ORDER BY 1;
经过测试，polypow计算上述多项式300次方是0.001秒，3000次方是0.03秒。
现在我理解了为什么有人想到FFT来计算, 因为FFT一个用途就是把多项式计算复杂度从O(N^2)降低到O(N*LOGN)，用于高精度计算（本质上10进制数是各位的10^N多项式之和）。
上述函数的源代码见这里，注意最后用卷积来提高效率。
# note: this is more efficient than `pu._pow(polymul, c1, c2)`, as it
# avoids calling `as_series` repeatedly
return pu._pow(np.convolve, c, pow, maxpower)
而卷积运算正好与FFT相关np.convolve。现在人工智能模型运算中卷积也用得很多。
See Also
--------
scipy.signal.fftconvolve : Convolve two arrays using the Fast Fourier
Transform.
不用FFT，用快速幂算法也能提高一些效率，这就是原帖中出题者的解法，矩阵相乘，7次=4次+2次+1次，把连接次数按二进制位分解。
以下是Qwen 3按照我的提示给出的结果
import numpy as np
def poly_pow(poly, power):
"""
使用快速幂算法计算多项式的N次幂。
参数:
poly (array-like): 输入多项式，按升幂排列的系数列表或数组。
power (int): 非负整数，要计算的幂次。
返回:
result (np.ndarray): 结果多项式系数数组。
"""
result = np.array([1]) # 初始化为常数项 1（即多项式 1）
current = np.array(poly, dtype=float)
while power > 0:
if power % 2 == 1:
result = np.polymul(result, current) # result *= current
current = np.polymul(current, current) # current ^= 2
power //= 2
return result
# 示例用法：
if __name__ == "__main__":
poly = [1, 2] # 表示 1 + 2x
n = 5
res = poly_pow(poly, n)
print("结果多项式系数:", res)
时代在变，原帖中很多获奖选手的解题帖早就打不开了，而数学思想永流传, AI在这方面可以做很多事。
参考资料：