快速傅里叶变换 FFT （二）

快速傅里叶变换 FFT （二）

Chirp-Z Transform

FFT 可以求单位根为公比的等比数列在多项式函数上的点值。

Chirp-Z Transform 是一种可以求任意数为公比的任意长度的等比数列在多项式函数上的点值。

设 nn 次多项式 P(x)=∑i=0n−1aixiP(x)=\displaystyle\sum_{i=0}^{n-1}a_ix^i 和整数 c,mc,m ，则 Chirp-Z Transform 可以快速求出 P(1),P(c),P(c2)⋯P(cm−1)P(1),P(c),P(c^2)\cdots P(c^{m-1}) 。

引理 12 ：

ij=Ci+j2−Ci2−Cj2ij=C_{i+j}^2-C_i^2-C_j^2

证明显然。

设 ansi=P(ci)ans_i=P(c^i) ，则有：

ansi=∑j=0n−1ajcijans_i=\sum_{j=0}^{n-1} a_jc^{ij}

ansi=c−Ci2∑j=0n−1ajc−Cj2cCi+j2ans_i=c^{-C_i^2}\sum_{j=0}^{n-1} a_jc^{-C_j^2}c^{C_{i+j}^2}

这是减法卷积形式，可以直接求解。

题目： Chirp-Z Transform 。

主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int main(){
    cin>>n>>c>>m;
    n--;m--;
    for(int i=0;i<=n;i++)cin>>a[i],a[i]=1ll*a[i]*kpow(c,P-1-1ll*i*(i-1)/2%(P-1))%P;
    len=1,l=0;
    while(len<=n+m)len<<=1,l++;
    for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=0;i<=m+n;i++)b[n+m-i]=kpow(c,1ll*i*(i-1)/2%(P-1));
    NTT(a,1);NTT(b,1);for(int i=0;i<len;i++)a[i]=1ll*a[i]*b[i]%P;NTT(a,-1);
    int inv=kpow(len,P-2);
    reverse(a,a+n+m+1);
    for(int i=0;i<=m;i++)cout<<1ll*a[i]*inv%P*kpow(c,P-1-1ll*i*(i-1)/2%(P-1))%P<<" ";
    return 0;
}

前置内容

中国剩余定理（ CRT ）

给定 nn 个同余方程，满足 bib_i 均为质数：

{x≡a1(modb1)x≡a2(modb2)⋯x≡an(modbn)\begin{cases} x \equiv a_1 \pmod {b_1} \\ x \equiv a_2 \pmod {b_2} \\ \cdots \\ x \equiv a_n \pmod {b_n} \end{cases}

设 M=∏ibiM=\displaystyle\prod_i b_i ，求 xx 在 MM 意义下的解。

题目： 中国剩余定理 。

事实上，由于余数的可加性，设 xix_i 为满足如下方程的解：

{xi≡0(modb1)xi≡0(modb2)⋯xi≡ai(modbi)xi≡0(modbi+1)⋯xi≡0(modbn)\begin{cases} x_i \equiv 0 \pmod {b_1} \\ x_i \equiv 0 \pmod {b_2} \\ \cdots \\ x_i \equiv a_i \pmod{b_i}\\ x_i \equiv 0 \pmod{b_{i+1}}\\ \cdots\\ x_i \equiv 0 \pmod {b_n} \end{cases}

则有 x=∑ixix=\displaystyle\sum_i x_i 。

对于每个 ii ，由于余数具有可乘性，设 yiy_i 为满足如下方程的解：

{yi≡0(modb1)yi≡0(modb2)⋯yi≡1(modbi)yi≡0(modbi+1)⋯yi≡0(modbn)\begin{cases} y_i \equiv 0 \pmod {b_1} \\ y_i \equiv 0 \pmod {b_2} \\ \cdots \\ y_i \equiv 1 \pmod{b_i}\\ y_i \equiv 0 \pmod{b_{i+1}}\\ \cdots\\ y_i \equiv 0 \pmod {b_n} \end{cases}

则有 xi=aiyix_i=a_iy_i 。

不难发现，此时 yiy_i 必然是 ∏j≠ibj\displaystyle\prod_{j\ne i}b_j 的倍数。所以设 mi=Mbi,yi=mikim_i=\frac{M}{b_i},y_i=m_ik_i 。

由于第 ii 个同余方程也要满足，所以有 miki≡1(modbi)m_ik_i \equiv 1 \pmod{b_i} ，使用费马小定理可得 ki=mibi−2k_i=m_i^{b_i-2} 。

那么这样，最终就有：

x=∑i=1nmikiai mod Mx=\sum_{i=1}^n m_ik_ia_i \bmod M

主要代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
#include<bits/stdc++.h>
#define int __int128
using namespace std;
const int MAXN = 1e5 + 10;
int n;
int a[MAXN] , b[MAXN];
int s;

int exgcd(int a , int b , int &x , int &y) {
	if(!b) {
		x = 1;
		y = 0;
		return a;
	}
	int r = exgcd(b , a % b , x , y);
	int t = x;
	x = y;
	y = t - a / b * y;
	return r;
}

inline int read() {
	char c = getchar();
	int x = 0 , s = 1;
	while(c < '0' || c > '9') {
		if(c == '-')
			s = -1;
		c = getchar();
	}
	while(c >= '0' && c <= '9') {
		x = x * 10 + c - '0';
		c = getchar();
	}
	return x * s;
}

void write(int x){
	if(x > 9)
		write(x/10);
	putchar(x % 10 | 48);
}

signed main() {
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.tie(0);
	n = read();
	int ans = 1;
	for(register int i = 1;i <= n;i ++) {
		a[i] = read();
		b[i] = read();
		ans *= a[i];
	}
	for(register int i = 1;i <= n;i ++) {
		int k = ans / a[i];
		int x , y;
		exgcd(k , a[i] , x , y);
		s = s +  k * b[i] * x % ans;
	}
	write((s % ans + ans) % ans);
	return 0;
}

任意模数多项式乘法 MTT

MTT （ any Mod Theoretic Transforms ，任意模数快速数论变换）是一种支持在任意模数情况下的多项式乘法，可以做到精度与常数兼备。

题目： 任意模数多项式乘法 。

MTT 的实现分为两种：三模 NTT 和拆系数 FFT 。

三模 NTT

事实上，虽然模数不满足 NTT 使用条件，但是我们仍然可以选出一些满足 NTT 使用条件的模数，然后进行 CRT 合并答案即可。

我们选的三个模数必须满足乘积大于多项式乘法中所能达到的最大值，因为 CRT 只能求模数乘积意义下的解，多项式乘法最大值不能超过，不然会出错。

本题中，多项式乘法最大值会达到 109×109×105=102310^9 \times 10^9 \times 10^5=10^{23} ，可以选用的模数为 479×221+1,7×225+1,119×223+1479\times 2^{21}+1,7\times 2^{25}+1,119\times 2^{23}+1 。

然后就是模板了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#include <bits/stdc++.h>
using namespace std;
const int N=400010,P1=(1<<23)*119+1,P2=(1<<21)*479+1,P3=(1<<26)*7+1;
int n,m,p,a[N],b[N],r[N],len,l,x[N],y[N],z[N],tmpa[N],tmpb[N];
long long kpow(long long a,long long b,long long P){
    a%=P;b%=(P-1);
    long long t=1;
    while(b){
        if(b&1)t=t*a%P;
        a=a*a%P;
        b>>=1;
    }
    return t;
}
void NTT(int *f,int t,int P){
    for(int i=0;i<len;i++)if(i<r[i])swap(f[i],f[r[i]]);
    for(int i=2;i<=len;i<<=1){
        int G=kpow((t==1?3:kpow(3,P-2,P)),(P-1)/i,P);
        for(int j=0;j<len;j+=i){
            for(int k=j,g=1;k<(j+(i>>1));k++,g=1ll*g*G%P){
                int x=f[k],y=1ll*g*f[k+(i>>1)]%P;
                f[k]=(x+y)%P;
                f[k+(i>>1)]=(x-y+P)%P;
            }
        }
    }
}
void Mul(int *a,int *b,int *c,int P){
    for(int i=0;i<len;i++)tmpa[i]=a[i];
    for(int i=0;i<len;i++)tmpb[i]=b[i];
    NTT(tmpa,1,P);NTT(tmpb,1,P);
    for(int i=0;i<len;i++)tmpa[i]=1ll*tmpa[i]*tmpb[i]%P;
    NTT(tmpa,-1,P);
    int in=kpow(len,P-2,P);
    for(int i=0;i<len;i++)c[i]=tmpa[i]*1ll*in%P;
}
long long mul(long long a,long long b,long long P){
    a%=P;b%=P;
    long long t=0;
    while(b){
        if(b&1)t=(t+a)%P;
        a=(a+a)%P;
        b>>=1;
    }
    return t;
}
int main(){
    cin>>n>>m>>p;
    for(int i=0;i<=n;i++)cin>>a[i];
    for(int i=0;i<=m;i++)cin>>b[i];
    len=1,l=0;
    while(len<=n+m)len<<=1,l++;
    for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    Mul(a,b,x,P1);Mul(a,b,y,P2);Mul(a,b,z,P3);
    const long long I1=kpow(P2,P1-2,P1),I2=kpow(P1,P2-2,P2);
    const long long M=1ll*P1*P2;
    const long long B2=kpow(M,P3-2,P3);
    for(int i=0;i<=n+m;i++){
        long long c=(mul(x[i],mul(P2,I1,M),M)+mul(y[i],mul(P1,I2,M),M))%M;
        long long A=1ll*(z[i]%P3-c%P3+P3)%P3*B2%P3;
        cout<<(mul(A,M,p)+c%p)%p<<" ";
    }
    return 0;
}

可以发现，这种做法的问题主要在于常数大（不用 int128 甚至要用龟速乘），做了 99 次 FFT 。但是这种做法精度却很好， 100%100\% 正确。

拆系数 FFT

正常来讲， FFT 完全可以胜任 102310^{23} 量级下的运算，但是由于 double 精度实在太差，在模意义下几乎一定错。

于是就诞生了拆系数 FFT ，用来降低值域。

设相乘的两个多项式为 P(x),Q(x)P(x),Q(x) 且 P(x)Q(x)=R(x)P(x)Q(x)=R(x) ，设 W=215W=2^{15} ，且

P(x)=A(x)W+B(x),Q(x)=C(x)W+D(x)P(x)=A(x)W+B(x),Q(x)=C(x)W+D(x)

R(x)=P(x)Q(x)=(A(x)W+B(x))(C(x)W+D(x))=A(x)C(x)W2+(A(x)D(x)+B(x)C(x))W+B(x)D(x)R(x)=P(x)Q(x)=(A(x)W+B(x))(C(x)W+D(x))=A(x)C(x)W^2+(A(x)D(x)+B(x)C(x))W+B(x)D(x)

那么精度可以控制在一个可接受的范围内。

然后我们就要考虑如何求出 A(x)C(x),A(x)D(x),B(x)C(x),B(x)D(x)A(x)C(x),A(x)D(x),B(x)C(x),B(x)D(x) 。

12 次 FFT

这个只需要暴力就行。

这样的算法精度差，常数大，没啥用。

8 次 FFT

想要做到这个，有两种方案。

第一种是将 12 次 FFT 使用 FFT 的优化，三步变两步。

第二种是注意到如下性质：

FFT 过程中，是将两个序列正变换各一次，再将乘积逆变换一次。

那么我们可以先对 A,B,C,DA,B,C,D 经行 FFT ，在对 AC,AD,BC,BDAC,AD,BC,BD 经行逆变换，就可以达到 8 次了。

这种做法精度也差，常数也大，数据水可以过，数据强就没了。

第二种代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#include <bits/stdc++.h>
#define double long double
using namespace std;
const int N=300010;
#define ll long long
const double pi=acos(-1);
int n,m,len,r[N],l,p;
struct C{
    double x,y;
    C(){x=0;y=0;}
    C(double x,double y):x(x),y(y){}
}a[N],b[N],c[N],d[N],aa[N],bb[N],cc[N],dd[N],rja[N],rjb[N];
C operator+(C x,C y){return C(x.x+y.x,x.y+y.y);}
C operator-(C x,C y){return C(x.x-y.x,x.y-y.y);}
C operator*(C x,C y){return C(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
C operator/(C x,int y){return C(x.x/y,x.y/y);}
C operator~(C x){return C(x.x,-x.y);}
void FFT(C *f,int t){
    for(int i=0;i<len;i++)if(i<r[i])swap(f[i],f[r[i]]);
    for(int i=2;i<=len;i<<=1){
        C G(cos(pi/(i>>1)),t*sin(pi/(i>>1)));
        for(int j=0;j<len;j+=i){
            C g(1,0);
            for(int k=j;k<(j+(i>>1));k++,g=g*G){
                C x=f[k],y=g*f[k+(i>>1)];
                f[k]=x+y;
                f[k+(i>>1)]=x-y;
            }
        }
    }
    if(t==-1)
    for(int i=0;i<len;i++)f[i]=f[i]/len;
}
pair<int,int> split(int x,int c){return {x/c,x%c};}
ll Floor(C x){return (ll)(x.x+0.5);}
int main(){
    cin>>n>>m>>p;
    for(int i=0;i<=n;i++)cin>>rja[i].x;
    for(int i=0;i<=m;i++)cin>>rjb[i].x;
    l=0;len=1;
    while(len<=n+m)len<<=1,l++;
    for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)+((i&1)<<(l-1));
    const int ct=1<<15;
    for(int i=0;i<=m;i++)tie(d[i].x,c[i].x)=split(rjb[i].x,ct);
    for(int i=0;i<=n;i++)tie(b[i].x,a[i].x)=split(rja[i].x,ct);
    FFT(a,1);FFT(b,1);FFT(c,1);FFT(d,1);
    for(int i=0;i<len;i++)aa[i]=b[i]*d[i],bb[i]=a[i]*d[i],cc[i]=b[i]*c[i],dd[i]=a[i]*c[i];
    FFT(aa,-1);FFT(bb,-1);FFT(cc,-1);FFT(dd,-1);
    for(int i=0;i<=n+m;i++){
        ll ta=Floor(aa[i]),tb=Floor(bb[i]),tc=Floor(cc[i]),td=Floor(dd[i]);
        int ans1=ta%p*ct%p*ct%p;
        int ans2=(tb%p+tc%p)%p*ct%p;
        int ans3=td%p;
        cout<<((ans1+ans2)%p+ans3)%p<<" ";
    }   
    return 0;
}

6 次 FFT

注意到三步变两步的优化中，第二种方案是将两个实数序列合并到一起只进行一次 FFT 。

这种思想当然能用到 8 次 FFT 的第二种做法了。

将 A,BA,B 合并， C,DC,D 合并，可以做到 6 次 FFT 。

常数和精度到这一步就不错了，当然在某些题中依然不行。

5 次 FFT

在三步变两步中，核心是利用 FFT 可以同时计算复数上的多项式乘法，而我们只需要实数项这一特点进行优化的。

类似的，令 T(x)=C(x)+iD(x)T(x)=C(x)+iD(x) ，然后计算 A(x)T(x),B(x)T(x)A(x)T(x),B(x)T(x) ，就可以通过提取实部虚部得到答案了。

可以发现，全过程只需要对 A,B,TA,B,T 正变换，对 AT,BTAT,BT 逆变换，共只需 5 次 FFT !

常数已经很不错了，精度也还行（ FFT 的精度本来就没多好），这是一个比较好的算法。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#include <bits/stdc++.h>
#define double long double
using namespace std;
const int N=300010;
#define ll long long
const double pi=acos(-1);
int n,m,len,r[N],l,p;
struct C{
    double x,y;
    C(){x=0;y=0;}
    C(double x,double y):x(x),y(y){}
}a[N],b[N],aa[N],bb[N],rja[N],rjb[N],t[N];
C operator+(C x,C y){return C(x.x+y.x,x.y+y.y);}
C operator-(C x,C y){return C(x.x-y.x,x.y-y.y);}
C operator*(C x,C y){return C(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
C operator/(C x,int y){return C(x.x/y,x.y/y);}
C operator~(C x){return C(x.x,-x.y);}
void FFT(C *f){
    for(int i=0;i<len;i++)if(i<r[i])swap(f[i],f[r[i]]);
    for(int i=2;i<=len;i<<=1){
        C G(cos(pi/(i>>1)),sin(pi/(i>>1)));
        for(int j=0;j<len;j+=i){
            C g(1,0);
            for(int k=j;k<(j+(i>>1));k++,g=g*G){
                C x=f[k],y=g*f[k+(i>>1)];
                f[k]=x+y;
                f[k+(i>>1)]=x-y;
            }
        }
    }
}
void IFFT(C *f){
    FFT(f);reverse(f+1,f+len);
    for(int i=0;i<len;i++)f[i]=f[i]/len;
}
pair<int,int> split(int x,int c){return {x/c,x%c};}
ll Floor(double x){return (ll)(x+0.5);}
int main(){
    cin>>n>>m>>p;
    for(int i=0;i<=n;i++)cin>>rja[i].x;
    for(int i=0;i<=m;i++)cin>>rjb[i].x;
    l=0;len=1;
    while(len<=n+m)len<<=1,l++;
    for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)+((i&1)<<(l-1));
    const int ct=1<<15;
    for(int i=0;i<=m;i++)tie(t[i].y,t[i].x)=split(rjb[i].x,ct);
    for(int i=0;i<=n;i++)tie(b[i].x,a[i].x)=split(rja[i].x,ct);
    FFT(a);FFT(b);FFT(t);
    for(int i=0;i<len;i++)aa[i]=a[i]*t[i],bb[i]=b[i]*t[i];
    IFFT(aa);IFFT(bb);
    for(int i=0;i<=n+m;i++){
        ll ta=Floor(bb[i].y),tb=Floor(bb[i].x),tc=Floor(aa[i].y),td=Floor(aa[i].x);
        int ans1=ta%p*ct%p*ct%p;
        int ans2=(tb%p+tc%p)%p*ct%p;
        int ans3=td%p;
        cout<<((ans1+ans2)%p+ans3)%p<<" ";
    }   
    return 0;
}

4 次 FFT

上面的算法中，计算 A,BA,B 的 FFT 依旧可以用三步变两步优化成一次，这样可以一次做 A,BA,B ，一次 TT ，一次 ATAT ，一次 BTBT ，共四次即可。

常数非常好（比三模 NTT 好多了），精度和 FFT 基本相当，是一个十分优秀的算法。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#include <bits/stdc++.h>
#define double long double
using namespace std;
const int N=300010;
#define ll long long
const double pi=acos(-1);
int n,m,len,r[N],l,p;
struct C{
    double x,y;
    C(){x=0;y=0;}
    C(double x,double y):x(x),y(y){}
}a[N],b[N],aa[N],bb[N],rja[N],t[N];
const C I(0,1);
C operator+(C x,C y){return C(x.x+y.x,x.y+y.y);}
C operator-(C x,C y){return C(x.x-y.x,x.y-y.y);}
C operator*(C x,C y){return C(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
C operator/(C x,int y){return C(x.x/y,x.y/y);}
C operator~(C x){return C(x.x,-x.y);}
void FFT(C *f){
    for(int i=0;i<len;i++)if(i<r[i])swap(f[i],f[r[i]]);
    for(int i=2;i<=len;i<<=1){
        C G(cos(pi/(i>>1)),sin(pi/(i>>1)));
        for(int j=0;j<len;j+=i){
            C g(1,0);
            for(int k=j;k<(j+(i>>1));k++,g=g*G){
                C x=f[k],y=g*f[k+(i>>1)];
                f[k]=x+y;
                f[k+(i>>1)]=x-y;
            }
        }
    }
}
void IFFT(C *f){
    FFT(f);reverse(f+1,f+len);
    for(int i=0;i<len;i++)f[i]=f[i]/len;
}
pair<int,int> split(int x,int c){return {x/c,x%c};}
ll Floor(double x){return (ll)(x+0.5);}
int main(){
    cin>>n>>m>>p;
    for(int i=0;i<=n;i++)cin>>rja[i].x;
    for(int i=0;i<=m;i++)cin>>rja[i].y;
    l=0;len=1;
    while(len<=n+m)len<<=1,l++;
    for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)+((i&1)<<(l-1));
    const int ct=1<<15;
    for(int i=0;i<=m;i++)tie(t[i].y,t[i].x)=split(rja[i].y,ct);
    for(int i=0;i<=n;i++)tie(a[i].y,a[i].x)=split(rja[i].x,ct);
    FFT(a);FFT(t);a[len]=a[0];
    for(int i=0;i<len;i++)aa[i]=(a[i]+(~a[len-i]))/2*t[i],bb[i]=((~a[len-i])-a[i])/2*I*t[i];
    IFFT(aa);IFFT(bb);
    for(int i=0;i<=n+m;i++){
        ll ta=Floor(bb[i].y),tb=Floor(bb[i].x),tc=Floor(aa[i].y),td=Floor(aa[i].x);
        int ans1=ta%p*ct%p*ct%p;
        int ans2=(tb%p+tc%p)%p*ct%p;
        int ans3=td%p;
        cout<<((ans1+ans2)%p+ans3)%p<<" ";
    }   
    return 0;
}

可以发现，三模 NTT 和拆系数 FFT 各有千秋：三模 NTT 精度好，常数差；拆系数 FFT 精度差，常数好。

可以根据具体需求调整算法（一般都是拆系数 FFT ，好写又快）。

任意模数多项式乘法逆

题目： 任意模数多项式乘法逆 。

给定多项式 F(x)F(x) ，求一个多项式 G(x)G(x) ，使 F(x)G(x)≡1(modxn)F(x)G(x) \equiv 1 \pmod{x^n} 。

你只需要求出 G(x) mod 109+7G(x) \bmod 10^9+7 的结果。

这个直接套模板就行。最好写 4 次 FFT 。

不过由于牛顿迭代法需要多次调用 MTT ，直接写十分复杂，可以写封装。具体参见代码。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#include <bits/stdc++.h>
using namespace std;
const int N=600100,P=1000000007;
#define double long double
struct C{
    double x,y;
    C(){x=y=0;}
    C(double x,double y):x(x),y(y){}
    C operator+(C a){return C(x+a.x,y+a.y);}
    C operator-(C a){return C(x-a.x,y-a.y);}
    C operator*(C a){return C(x*a.x-y*a.y,x*a.y+y*a.x);}
    C operator~(){return C(x,-y);}
    C operator/(double a){return C(x/a,y/a);}
}tmp[N];
const double pi=acos(-1);
struct FFT_Structure{
    int len,l,r[N];
    C f[N];
    FFT_Structure(){len=l=0;}
    int resize(int n){
        len=1,l=0;
        while(len<n)len<<=1,l++;
        for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        return len;
    }
    void init(C a[]){for(int i=0;i<len;i++)f[i]=a[i];}
    void FFT(int t){
        for(int i=0;i<len;i++)if(i<r[i])swap(f[i],f[r[i]]);
        for(int i=2;i<=len;i<<=1){
            C G(cos(2*pi/i),t*sin(2*pi/i));
            for(int j=0;j<len;j+=i){
                C g(1,0);
                for(int k=j;k<(i>>1)+j;k++){
                    C x=f[k],y=g*f[k+(i>>1)];
                    f[k]=x+y;
                    f[k+(i>>1)]=x-y;
                    g=g*G;
                }
            }
        }
        if(t==-1)for(int i=0;i<len;i++)f[i]=f[i]/len;
    }
    C& operator[](int x){return f[x];}
}A,T,D;
#define ll long long
const int B=1<<15;
C split(int x){return C(x/B,x%B);}
//(Ax+B)(Cx+D)=(AC)x^2+(AD+BC)x+BD
//A(C+Di)=AC+ADi,B(C+Di)=BC+BDi
ll flr(double x){return (ll)(x+0.5);}
void MTT(int a[],int b[],int c[],int n){
    int len=A.resize(n<<1);T.resize(n<<1);D.resize(n<<1);
    for(int i=0;i<len;i++)tmp[i]=split(a[i]);A.init(tmp);
    for(int i=0;i<len;i++)tmp[i]=split(b[i]);T.init(tmp);
    A.FFT(1);T.FFT(1);
    A.f[len]=A.f[0];
    for(int i=0;i<len;i++){
        C ta=(A[i]+(~A[len-i]))/2;
        C tb=((~A[len-i])-A[i])*C(0,1)/2;
        D[i]=ta*T[i];
        T[i]=T[i]*tb;
    }
    D.FFT(-1);T.FFT(-1);
    for(int i=0;i<len;i++){
        ll ta=flr(D[i].x),tb=flr(D[i].y),tc=flr(T[i].x),td=flr(T[i].y);
        int ans1=ta%P*B%P*B%P;
        int ans2=(tb%P+tc%P)%P*B%P;
        int ans3=td%P;
        c[i]=((ans1+ans2)%P+ans3)%P;
    }
}
int a[N],n,b[N],c[N],ta[N],tb[N];
int kpow(int a,int b){
    int t=1;
    while(b){
        if(b&1)t=1ll*t*a%P;
        a=1ll*a*a%P;
        b>>=1;
    }
    return t;
}
int inv(int x){return kpow(x,P-2);}
int main(){
    cin>>n;
    for(int i=0;i<n;i++)cin>>a[i];
    b[0]=inv(a[0]);
    for(int tt=2;(tt>>1)<=n;tt<<=1){
        for(int i=0;i<tt;i++)ta[i]=a[i],c[i]=b[i];
        MTT(c,c,tb,tt);MTT(tb,ta,c,tt);
        for(int i=0;i<tt;i++)b[i]=(2ll*b[i]%P-c[i]+P)%P;
    }
    for(int i=0;i<n;i++)cout<<b[i]<<" ";
    return 0;
}

任意模数 Chirp-Z Transform

题目： 任意模数 Chirp-Z Transform 。

设 nn 次多项式 P(x)=∑i=0n−1aixiP(x)=\displaystyle\sum_{i=0}^{n-1}a_ix^i 和整数 c,mc,m ，求出 P(1),P(c),P(c2)⋯P(cm−1)P(1),P(c),P(c^2)\cdots P(c^{m-1}) 模 109+710^9+7 意义下的值。

这个也是直接写封装，套模板。但是这题数据范围达到 6×1056\times 10^5 ， FFT 的精度会有问题，三模 NTT 的常数巨大，无法通过。

事实上， FFT 中有一个 g=g*G 的操作，这个操作本来是要算 wnkw_n^k ，但是直接用三角函数很慢（因为三角函数会调用 nlog⁡nn\log n 次），所以用迭代代替。

但是，我们可以先预处理一个三角函数表，直接用三角函数算。然后 FFT 过程中直接查表，这样就可以提升精度，就可以通过了。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
#include <bits/stdc++.h>
using namespace std;
const int N=2100100,P=1000000007;
struct C{
    double x,y;
    C(){x=y=0;}
    C(double x,double y):x(x),y(y){}
    C operator+(C a){return C(x+a.x,y+a.y);}
    C operator-(C a){return C(x-a.x,y-a.y);}
    C operator*(C a){return C(x*a.x-y*a.y,x*a.y+y*a.x);}
    C operator~(){return C(x,-y);}
    C operator/(double a){return C(x/a,y/a);}
}tmp[N],qp[N];
const double pi=acos(-1);
struct FFT_Structure{
    int len,l,r[N];
    C f[N];
    FFT_Structure(){len=l=0;}
    int resize(int n){
        len=1,l=0;
        while(len<n)len<<=1,l++;
        for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        for(int i=2;i<=len;i<<=1){
            for(int j=0;j<(i>>1);j++){
                qp[2*len/i*j]=C(cos(2*j*pi/i),sin(2*j*pi/i));
            }
        }
        return len;
    }
    void init(C a[]){for(int i=0;i<len;i++)f[i]=a[i];}
    void FFT(int t){
        for(int i=0;i<len;i++)if(i<r[i])swap(f[i],f[r[i]]);
        for(int i=2,tt=1;i<=len;i<<=1,tt++){
        //    C G(cos(2*pi/i),t*sin(2*pi/i));
            for(int j=0;j<len;j+=i){
              //  C g(1,0);
                for(int k=j,tt=0;k<(i>>1)+j;k++,tt++){
                    C x=f[k],y=(t==1?qp[2*len/i*tt]:~qp[2*len/i*tt])*f[k+(i>>1)];
                    f[k]=x+y;
                    f[k+(i>>1)]=x-y;
                //    g=g*G;
                }
            }
        }
        if(t==-1)for(int i=0;i<len;i++)f[i]=f[i]/len;
    }
    C& operator[](int x){return f[x];}
}A,T,D;
#define ll long long
const int B=1<<15;
C split(int x){return C(x/B,x%B);}
//(Ax+B)(Cx+D)=(AC)x^2+(AD+BC)x+BD
//A(C+Di)=AC+ADi,B(C+Di)=BC+BDi
ll flr(double x){return (ll)(x+0.5);}
void MTT(int a[],int b[],int c[],int n){
    int len=A.resize(n);T.resize(n);D.resize(n);
    for(int i=0;i<len;i++)tmp[i]=split(a[i]);A.init(tmp);
    for(int i=0;i<len;i++)tmp[i]=split(b[i]);T.init(tmp);
    A.FFT(1);T.FFT(1);
    A.f[len]=A.f[0];
    for(int i=0;i<len;i++){
        C ta=(A[i]+(~A[len-i]))/2;
        C tb=((~A[len-i])-A[i])*C(0,1)/2;
        D[i]=ta*T[i];
        T[i]=T[i]*tb;
    }
    D.FFT(-1);T.FFT(-1);
    for(int i=0;i<len;i++){
        ll ta=flr(D[i].x),tb=flr(D[i].y),tc=flr(T[i].x),td=flr(T[i].y);
        int ans1=ta%P*B%P*B%P;
        int ans2=(tb%P+tc%P)%P*B%P;
        int ans3=td%P;
        c[i]=((ans1+ans2)%P+ans3)%P;
    }
}
const int BM=40000;
int n,m,a[N],b[N],c[N],k;
int bg[BM+10],sm[BM+10];
void init(int k){
    sm[0]=bg[0]=1;
    for(int i=1;i<=BM;i++)sm[i]=1ll*sm[i-1]*k%P;
    for(int i=1;i<=BM;i++)bg[i]=1ll*bg[i-1]*sm[BM]%P;
}
int kpow(int x){return 1ll*bg[x/BM]*sm[x%BM]%P;}
signed main(){
    cin>>n>>k>>m;init(k);
    if(k==0){
        return 1;
        int x,sum=0;
        cin>>x,sum+=x;
        for(int i=1,xx;i<n;i++)cin>>xx,(sum+=xx)%=P;
        cout<<sum<<' ';
        for(int i=1;i<m;i++)cout<<x<<" ";
        return 0;
    }
    for(int i=0,x;i<n;i++)cin>>x,a[n-i-1]=1ll*x*kpow((P-1-1ll*i*(i-1)/2%(P-1)))%P;
    for(int i=0;i<n+m;i++)b[i]=kpow((1ll*i*(i-1)/2%(P-1)));
    MTT(a,b,c,n+m);
    for(int i=0;i<m;i++)cout<<1ll*c[i+n-1]*kpow((P-1-1ll*i*(i-1)/2%(P-1)))%P<<" ";
    return 0;
}

完结撒花！