• 多项式学习笔记(二) NTT


    前置芝士:

    原根:满足 (g^{varphi(p)-1} equiv 1 pmod p) , 且 (g^0,g^1,g^2...g^{varphi(p)-1}) 互不相同,就称 (g)(p) 的一个原根

    NTT

    在FFT的时候,我们一个比较大的问题就是精度处理。如果说题目要求系数对大质数取模的话,那么FFT可能在计算中溢出,这时候我们就需要另外一种算法 (NTT) (其实是 FFT 的变形)。

    在FFT 中,我们单位根满足以下性质:

    1. (w_{n}^{0} = w_{n}^{n} = 1)
    2. (w_{n}^{0},w_{n}^{1},w_{n}^{2}....w_{n}^{n}) 互不相同。
    3. (w_{n}^{i} = w_{2*n}^{2*i})
    4. (w_{n}^{k+{nover 2}} = -w_{n}^{k})
    5. (displaystylesum_{i=0}^{n-1} (w_{n}^{k})^{i} = 0)(k mid n)

    然后我们找到了单位根一个很好的替代品原根,我们设 (w_{n}^{i} = (g^{p-1over n})^i) 。带入之后可以发现它是满足以上性质的。

    性质1(w_{n}^{0} = w_{n}^{n} = 1)

    很简单, (w_{n}^{0} = (g^{p-1over n}) ^0 = 1) , (w_{n}^{n} = (g^{p-1over n})^n = g^{p-1} = 1)

    性质2: (w_{n}^{0},w_{n}^{1},w_{n}^{2}....w_{n}^{n}) 互不相同。

    根据原根的性质可得: (g^{0},g^1....g^{p-1}) 都互不相同,所以 性质2得证。

    性质3: (w_{2*n}^{2*i} = w_{n}^{i})

    (w_{2*n}^{2*i} = (g^{p-1over 2n}) ^{2i} = g^{{p-1over 2n} imes 2i} = (g^{p-1over n})^i)

    (w_{n}^{i} = (g^{p-1over n})^i) .

    所以 (w_{2*n}^{2*i} = w_{n}^{i}) ,性质三得证。

    性质4: (w_{n}^{k+{nover 2}} = -w_{n}^{k})

    (w_{n}^{k+{nover 2}} = (g^{p-1over n})^{k+{nover 2}} = (g^{p-1over n})^k (g^{p-1over n})^{nover 2} = (g^{p-1over n})^{k} g^{p-1over 2})

    有一个有意思的性质就是 (g^{p-1over 2} = -1)

    (ecause (g^{p-1over 2})^2 = 1)

    ( herefore g^{p-1over 2} = 1 或 g^{p-1over 2} = -1)

    (ecause g^{p-1over 2} eq g^{p-1})

    ( herefore g^{p-1over 2} = -1)

    所以 (w_{n}^{k+{nover 2}} = -(g^{p-1over n})^k = -w_{n}^{k})

    性质5: (displaystylesum_{i=0}^{n-1} (w_{n}^{k})^i = 0) (k mid n)

    根据等比数列的求和公式原式等价于 (1-(w_{n}^{k})^{n}over 1-w_{n}^{k})

    分母 (1-(w_{n}^{k})^n = 1-(g^{p-1over n})^{kn} = 1-(g^{p-1})^k = 0)

    性质五,得证。

    NTT 和FFT 的代码其实是差不多的。

    #include<iostream>
    #include<cstdio>
    #include<algorithm>
    #include<cmath>
    using namespace std;
    #define int long long
    const int N = 5e6+10;
    const int p = 998244353;
    int n,m,len,tim,a[N],b[N],Rev[N];
    inline int read()
    {
    	int s = 0,w = 1; char ch = getchar();
    	while(ch < '0' || ch > '9'){if(ch == '-')w = -1; ch = getchar();}
    	while(ch >= '0' && ch <= '9'){s = s * 10 + ch - '0'; ch = getchar();}
    	return s * w;
    }
    int ksm(int a,int b)
    {
    	int res = 1;
    	for(; b; b >>= 1)
    	{
    		if(b & 1) res = res * a % p;
    		a = a * a % p;
    	}
    	return res;
    }
    int inv(int n){return ksm(n,p-2);}
    void NTT(int *a,int flag)
    {
    	for(int i = 0; i < len; i++)
    	{
    		if(i < Rev[i]) swap(a[i],a[Rev[i]]);
    	}
    	for(int h = 1; h < len; h <<= 1)
    	{
    		int wn = ksm(3,(p-1)/(h<<1));//原根
    		if(flag == -1) wn = inv(wn);
    		for(int j = 0; j < len; j += (h<<1))
    		{
    			int w = 1;
    			for(int k = 0; k < h; k++)
    			{
    				int u = a[j+k];
    				int v = w * a[j+h+k] % p;
    				a[j+k] = (u + v) % p;
    				a[j+h+k] = (u-v+p) % p;
    				w = w * wn % p; //wni
    			}
    		}
    	}
    }
    signed main()
    {
    	n = read(); m = read();
    	for(int i = 0; i <= n; i++) a[i] = read();
    	for(int i = 0; i <= m; i++) b[i] = read();
    	len = 1, tim = 0;
    	while(len <= n+m) len <<= 1, tim++;
    	for(int i = 0; i < len; i++) Rev[i] = (Rev[i>>1]>>1) | ((i&1)<<(tim-1));
    	NTT(a,1); NTT(b,1);
    	for(int i = 0; i < len; i++) a[i] = a[i] * b[i] % p;
    	NTT(a,-1);
    	for(int i = 0; i <= n+m; i++) printf("%lld ",a[i]*inv(len) % p);
    	printf("
    ");
    	return 0;
    }
    
  • 相关阅读:
    【胡策篇】题目
    【学术篇】luogu3768 简单的数学题(纯口胡无代码)
    【学术篇】规律选手再次证明自己(舒老师的胡策题 T2 LX还在迷路)
    【模板篇】Link Cut Tree模板(指针)
    【学术篇】51nod 1238 最小公倍数之和
    【学术篇】2.28测试T2 线段 拓扑排序
    【学术篇】SPOJ FTOUR2 点分治
    【颓废篇】Py:从零开始的poj自动提交
    【学术篇】CF935E Fafa and Ancient Mathematics 树形dp
    安卓启动图去除顶部title和状态栏
  • 原文地址:https://www.cnblogs.com/genshy/p/14163037.html
Copyright © 2020-2023  润新知