整系数Discrete Fourier transform技巧

Discrete Fourier transform

向量的DFT 定义为:

其中是primitive root of unity of order ,复数时通常取

Inverse transform

性质:

因此

Number theoretic transform

DFT可以用于复数以外的其他ring,常用于

使用128 bits模数需要高效的u64*u64%u64,其中模数是常数。

硬件除法指令(32 bits、64 bits)

DIV指令性能很低。

1
2
3
4
5
6
extern inline u64 mul_mod(u64 a, u64 b, u64 m)
{
u64 r;
asm("mulq %2\n\tdivq %3" : "=&d"(r), "+a"(a) : "rm"(b), "rm"(m) : "cc");
return r;
}

到AVX-512也没有提供把两个64 bits乘数的积放在一个128 bits寄存器的指令,GCC没有提供用乘法、移位等模拟的u128除以u64的常量除法。

64位mantissa浮点数(32 bits、64 bits)

当模数时可以用64位mantissa浮点数计算u64*u64%u64

由等式

两边模,得

即用u64乘法计算的低64位,减去的低64位。其中,可以用64位mantissa浮点数(Intel x87 80-bit double-extended precision)计算,再round成u64

round时若向上取整了,减数会大于被减数。若,可以根据差的符号位判断。

1
2
3
4
5
u64 mul_mod(u64 a, u64 b, u64 P)
{
u64 x = a*b, r = x - P*u64((long double)a*b/P+0.5);
return i64(r) < 0 ? r + P : r;
}

存储的倒数,用(long double)a*b*Q代替(long double)a*b/P能快些。此时会引入额外的误差,Matters Computational说适用于,原因不明。

编译器生成的常量除法(32 bits)

对于固定的模,GCC/llvm可以生成u64%u32的高效代码。llvm的lib/Transforms/Utils/IntegerDivision.cpp

Montgomery modular multiplication+Barret reduction(32 bits)

Faster arithmetic for number-theoretic transforms,用乘法和移位代替除法。

快速计算u64%u32,用乘法和移位代替除法。设为大于等于的整数,为负整数。

的估计值,若大于等于则减去的倍数。

因此

,则,即估计值最多小2,最多两个conditional move指令(if (r >= P) r -= P;)即可修正余数。

,则估计值最多小1。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
extern inline u64 barrett_30_2(u64 a, u64 P, u64 M)
{ // 2^29 < P < 2^30, M = floor(2^61/P)
u64 r = a-((a>>28)*M>>33)*P;
if (r >= P) r -= P;
return r;
}

extern inline u64 barrett_30_1(u64 a, u64 P, u64 M)
{ // 2^29 < P < 2^30, M = floor(2^60/P)
u64 r = a-((a>>29)*M>>31)*P;
if (r >= P) r -= P;
if (r >= P) r -= P;
return r;
}

当模数为常量时,该算法不如编译器生成的常量除法。若模数不固定时可以考虑使用。

Cyclic convolution

如果要计算linear convolution,比如多项式乘法,可以把长度补到,求cyclic convolution。

下面讨论用于计算整系数向量卷积的discrete Fourier transform。

精度

使用complex<double>计算convolution,需要保证结果每一项的实数部分在(Number.MIN_SAFE_VALUENumber.MAX_SAFE_INTEGER)范围内,是double能精确表示的最大整数。采取round half to even规则,均表示为,无法区分。

设每项系数的绝对值小于等于,那么convolution结果每一项绝对值小于等于,若则可放心使用complex<double>

complex<double>还要受到浮点运算误差影响。根据Roundoff Error Analysis of the Fast Fourier Transform,没仔细看,relative error均值为log2(n)*浮点运算精度*变换前系数最大值,对于结果,这个量达到就很可能出错,增长速度可以看作是,不如。因此通常不必担心浮点运算的误差。

对于模的number theoretic transform,,若则可放心使用。

998244353, 897581057, 880803841小于,两倍不超过INT_MAX,且可表示为,其中为2的幂,适合用作number theoretic transform的模。

系数取自的uniform distribution,则系数均值为,方差为。若把系数平移至,则系数均值,方差为。若其中之一independent and identically distributed,则方差会很小。可以用Chebyshev's inequality等估计系数绝对值,上界可以减小。即使不是independent and identically distributed,也可以用来计算,是independent and identically distributed uniform

方案0:sqrt decomposition(FFT, NTT)

为接近的整数,分解,则:

适当选择可以使的系数小于等于,convolution结果系数最大值为,比原来的小。

求出后,计算等式右边四个convolution,带权相加即得到原convolution。

需要4次长为的DFT、1次长为的inverse DFT。

可以使用Toom-2 (Karatsuba)计算,减少为3次DFT、1次inverse DFT。

容易扩展到cube root decomposition等。对于number theoretic transform,分成份需要个DFT和个IDFT,不如用Chinese remainder theorem。

优化0:imaginary unit用于进位(FFT)

接近且尽可能小。

分解

𝕒𝟘𝕒𝟙𝕓𝟘𝕓𝟙

每一项绝对值变为倍了。右边同余于,提取虚部与实部即可得到:

需要2次长为的DFT、1次长为的inverse DFT。

优化1:正交计算两个实系数向量DFT(FFT)

的共轭的DFT可由的DFT求出:

分解后,根据上面的公式,用计算,同法计算。然后用计算出,同法计算出

需要2次长为的DFT、2次长为的inverse DFT。

奇偶项优化(FFT)

该优化可以和其他方式叠加。

看作多项式,同样地,看作多项式

奇次项, 偶次项,同样地,定义的系数为A1(x)n$,高位用填充。

用正交计算两个实系数向量DFT的方式,用2次长度为(之前都是 )的DFT计算, 循环右移1位的DFT的第项等于,因此根据的DFT的系数可以得到的DFT的系数。

构造长为的向量

的实部为结果的偶次项系数,虚部为结果的奇次项系数。

需要2次长为的DFT、1次长为的inverse DFT。

方案1:Chinese remainder theorem(NTT)

个可用于number theoretic transform的质数,使,计算个NTT,之后用Chinese remainder theorem合并。

有至少两种算法。

经典算法(Gauss's algorithm)

Gauss之前也有很多人提出。

对于每个用Blankinship's algorithm计算