算法竞赛如何训练数论这一块？ 第1页
1
duan-zhong-jie-66 网友的相关建议:
题主你好，我是浙江师范大学2016级数学专业的ACM集训队队员，曾在计蒜客实习，对ACM与OI竞赛都有了解，希望我的回答对你有所帮助。
如果你是大学生，可以考虑了解一下抽象代数，可以从一个奇妙的角度来理解数论。
下面是我整理的数论算法，全文非常长，在 markdown 格式下约 6w 字符，建议放到收藏夹吃灰。
文中某些代码仅为了帮助读者理解算法而写，建议理解后自己写模板。
文中不包含积性函数前缀和筛，例如min25筛，因为笔者认为这样的知识点太偏太难，考起来只是大佬们的自嗨。（划掉）（更新）听评论区说想要这些杂技内容，那我就把模板贴这了，没有原理解析，也不会有的。不过多项式求逆是没有的，因为我没遇到过，模板也就没写过。
另外，厚颜无耻地推广自己的bilibili账号，里面有我校自用的ACM课程，我要退役了，你们加油！
数论算法详解
1. 整数的取余运算
定义 带余除法 设 是整数，且 ，使得
且 ，称 为商， 为余数
显然带余除法中的商和余数都是唯一的，在下文中将商记为 ，将余数记为 ，“/”与“%”的运算优先级与乘除法相同，当然，在C语言中二者分别对应 a/b 与 a%b
1.1 整数的加减乘在 中的运算
1.1.1 定理 整数的加、减、乘对于取余的右分配律
注意，在计算减法时，通常需要加 ，防止变成负数
在计算乘法时，如果 较大（但不超过64位整数范围），可以使用快速乘法进行计算，原理与快速幂运算类似，复杂度为
ll fastMul(ll a,ll b,ll p){ a%=p; ll ans=0; while(b>0){ if(b&1)ans=(ans+a)%p; b>>=1; a=(a+a)%p; } return ans; } 如果需要更快的快速乘法，可以使用 long double 数据类型进行计算，复杂度为
ll modMul(ll a,ll b,ll p){ if (p<=1000000000) return a*b%p; else if (p<=1000000000000ll) return (((a*(b>>20)%p)<<20)+(a*(b&((1<<20)-1))))%p; else { ll d=(ll)floor(a*(long double)b/p+0.5); ll ret=(a*b-d*p)%p; if (ret<0) ret+=p; return ret; } } 1.2 整数的除法在 中的运算( 为素数)
虽然取余运算对于“+”、“-”、“×”不难，但通常情况下
如何计算 ？我们有时可以找一个乘法逆元 ，使得 ，那么就有
如果 是素数，有下面的定理
1.2.1 定理 费马小定理 设 是一个整数， 是一个素数，二者互质，那么
将上式改写一下，得到
因此取 即可，一般需要用快速幂计算，但要注意，与除数不能为0类似，要保证
ll inv(ll a,ll p){ return fpow(a,p-2); } 1.3 整数的除法在 中的运算
上面的方法给出了模数为素数的解决方案，如果模数不是素数，可以用下面的方法
1.3.1 定理 若 ，则
这样可以不使用逆元来求 中的除法
其他逆元相关内容将在下文中介绍
（如果你像我一样无聊，可以把整环 中的数写成这样一个类）
class mint{ static const ll mo=1e9+7; inline ll fpow(ll a,ll b){ll c=1;while(b){if(b&1)c=c*a%mo;b>>=1;a=a*a%mo;}return c;} inline ll inv(ll a){return fpow(a,mo-2);} inline ll norm(ll a){return a<0?(a%mo+mo):a%mo;} public: ll v; mint(){} mint(ll x):v(x){} mint &operator =(const ll b){v=norm(b);return *this;} mint &operator +=(const mint &b){v+=b.v;if(v>=mo)v-=mo;return *this;} mint &operator -=(const mint &b){v-=b.v;if(v<0)v+=mo;return *this;} mint &operator *=(const mint &b){v=v*b.v%mo;return *this;} mint &operator /=(const mint &b){v=v*inv(b.v)%mo;return *this;} mint operator + (const mint &b)const{mint a{*this};return a+=b;} mint operator - (const mint &b)const{mint a{*this};return a-=b;} mint operator * (const mint &b)const{mint a{*this};return a*=b;} mint operator / (const mint &b)const{mint a{*this};return a/=b;} mint operator ^ (const mint &b){return mint(fpow(v,b.v));} mint &operator++(){v++;if(v==mo)v-=mo;return *this;} mint &operator--(){v--;if(v==-1)v+=mo;return *this;} mint operator -()const{return mint(v?(mo-v):0);} bool operator ==(const mint &b)const{return v==b.v;} bool operator <(const mint &b)const{return v<b.v;} bool operator >(const mint &b)const{return v<b.v;} bool operator !=(const mint &b)const{return v!=b.v;} explicit operator ll()const{return v;} friend ostream &operator << (ostream &o,const mint & b){o<<b.v;return o;} friend istream &operator >> (istream &in,mint &a){in>>a.v;a.v%=mo;return in;} }; 2. 最大公因数与最小公倍数
2.1 最大公因数
顾名思义，最大公因数就是公因数中最大的那个，我们记 的最大公因数为 ，有如下性质
2.1.1 性质
2.1.2 性质
2.1.3 性质
2.1.4 性质
2.1.5 性质
2.1.1 性质是显然的，2.1.2 性质是辗转相减法的原理，2.1.3 性质可以视为2.1.2 性质的“一步到位”版本，2.1.4 性质指出多个数的最大公因数可以递推地进行求解，2.1.5 性质说明 对乘法有分配律
2.2 辗转相除法
根据2.1.3 性质，得到辗转相除法的参考代码模板
typedef long long ll; ll gcd(ll a, ll b){ return b?gcd(b,a%b):a; } 注意当 时，返回值为 gcd(b,a%b) 而不是 gcd(a%b,b) ，否则会不断递归导致栈溢出
2.3 最小公倍数
最小公倍数就是公倍数中最小的那个，我们记 的最小公倍数为 ，有如下性质
2.3.1 性质
下面是最小公倍数的参考代码模板
ll lcm(ll a, ll b){ return a/gcd(a,b)*b; } 注意是先除后乘，避免在中间过程中数据超出64位整数的范围
2.4 扩展欧几里得算法
考虑二元一次不定方程
其中 是已知的正整数，如何求出方程的解呢？
2.4.1 定理 上述方程有解的充要条件是 （ 是 的倍数）
可以理解为， 是 可以表示出的最小正整数
2.4.2 定理 方程 的所有解为
其中 是一组特解
2.4.3 定理 方程 的所有解为
其中 是方程 的一组特解
下面是参考代码模板：
ll ext_gcd(ll a,ll b,ll& x,ll& y){ ll d = a; if (!b){ x = 1;y = 0; }else{ d = ext_gcd(b,a%b,y,x); y -= a/b*x; } return d; } 在求逆元时，要找到 使得 ，实质上是求解方程 中的 ，因此可以用扩展欧几里得算法来求逆元，当然只有 时才有解，否则逆元不存在
2.4.4 推论 逆元的存在性 存在 的充要条件是
3. 同余方程组
有方程就有方程组
其中 是整数， 是正整数
3.1.1 中国剩余定理（孙子定理） 设上述方程组中 两两互质，则方程组的通解为
其中
下面是参考代码模板，需要调用前面的扩展欧几里得算法模板
ll Sunzi(ll *m,ll *a,int len){ ll lcm = 1; ll ans = 0; for (int i=0;i<len;i++){ ll k0,ki; ll d = ext_gcd(lcm,m[i],k0,ki); if ((a[i]-ans)%d!=0) return -1; else { ll t = m[i]/d; k0 = ( k0*(a[i]-ans)/d%t + t)%t; ans = k0*lcm + ans; lcm = lcm/d*m[i]; } } return ans; } 4. 素数
素数是只有1和本身两个因数的数，1不是素数
4.1 素数的判断
bool isPrime(ll n){ if(n==1)return false; for(ll i=2;i*i<=n;i++) if(n%i==0)return false; return true; } 这是最简单的素数的判断的参考代码模板，复杂度为
原理其实很简单，对于一个大于1的整数，如果 是它的一个大于 的因子，那么 是它的小于 的因子
在大多数情况下，这种判断方式的复杂度已经足够小了，如果要追求更高的效率，可以考虑 法
一个大于1的整数如果不是素数，那么一定有素因子，因此在枚举因子时只需要考虑可能为素数的因子即可。 法即枚举形如 的数，例如取 ，那么 都不可能为素数，只需要枚举形如 的数即可，这样复杂度降低了 。
下面的模板是 法 的版本
bool isPrime(ll n){ if(n==2||n==3||n==5)return 1; if(n%2==0||n%3==0||n%5==0||n==1)return 0; ll c=7,a[8]={4,2,4,2,4,6,2,6}; while(c*c<=n)for(auto i:a){if(n%c==0)return 0;c+=i;} return 1; } 如果 极大，可以使用素数测试算法，素数测试算法可以通过控制迭代次数来间接控制正确率，常用的有下面的Miller-Rabin方法
ll Rand(){ static ll x=(srand((int)time(0)),rand()); x+=1000003; if(x>1000000007)x-=1000000007; return x; } bool Witness(ll a,ll n){ ll t=0,u=n-1; while(!(u&1))u>>=1,t++; ll x=fpow(a,u,n),y; while(t--){ y=x*x%n; if(y==1 && x!=1 && x!=n-1)return true; x=y; } return x!=1; } bool MillerRabin(ll n,ll s){ if(n==2||n==3||n==5)return 1; if(n%2==0||n%3==0||n%5==0||n==1)return 0; while(s--){ if(Witness(Rand()%(n-1)+1,n))return false; } return true; } 当然， Rand() 怎么写可以自由发挥，这会影响其性能
4.2 素数筛
如果要求出不超过 的所有素数，素数筛是最好的选择，下面是一种朴素的筛法
void getPrime(bool p[],int n){ for(int i=1;i<=n;i++)p[i]=true; p[1]=false; for(int i=2;i<=n;i++){ if(p[i]){ for(int j=i+i;j<=n;j+=i)p[j]=false; } } } 这种方法的原理是从小到大将素数的倍数筛掉，复杂度为 ，注意到每个合数如果有多个素因子，那么就会被重复筛掉，造成复杂度的浪费，因此，用下面的方法可以保证每个合数只被它最小的素因子筛掉一遍，以 的复杂度解决上述问题
ll getPrime(ll n,bool vis[],ll prime[]){ ll tot=0; for(ll i=1;i<=n;i++)vis[i]=0; for(ll i=2;i<=n;i++){ if(!vis[i])prime[tot++]=i; for(ll j=0;j<tot;j++){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0)break; } } return tot; } 5. 组合数取余
5.1 组合数
5.1.1 定义 组合数 在 个不同元素中选取 个元素，不同的取法记为
5.2 杨辉三角
组合数与杨辉三角中的数字是一一对应的
杨辉三角的自然数形式
杨辉三角的组合数形式
按照上面的写法，杨辉三角的第 行第 列即为
注意到上图中每个数等于其左边的数与上边的数（如果有的话）之和，这就是杨辉恒等式
5.3 组合数取余的求法
在ACM竞赛中，我们常常需要计算 ，可以参考下面几种方法
- 如果 很小（不超过50），可以用C++的库函数
double tgamma(double x)，这是一个欧拉积分
在整数点处的取值满足
因此代码可以这么写
ll C(ll n,ll m){ return (ll)round(tgamma(n+1)/tgamma(m+1)/tgamma(n-m+1)); } 效率并不高，但是对于追求手速来说足够了
- 如果 不大，可以开 的空间，可以利用杨辉恒等式来预处理组合数表
const ll mo=1e9+7; ll C[1005][1005]; void getC(int n){ for(int i=0;i<=n;i++){ for(int j=0;j<=i;j++){ if(j==0 || j==i) C[i][j]=1; else C[i][j]=(C[i-1][j-1]+C[i-1][j])%mo; } } } - 如果 比较大，可以开 的空间，可以利用前文所述的逆元来求解，当然，要保证 是素数
const ll mo=1e9+7; ll C(ll n,ll m){ static ll M=0,inv[N],mul[N],invMul[N]; while(M<=n){ if(M){ inv[M]=M==1?1:(mo-mo/M)*inv[mo%M]%mo; mul[M]=mul[M-1]*M%mo; invMul[M]=invMul[M-1]*inv[M]%mo; } else mul[M]=1,invMul[M]=1; M++; } return mul[n]*invMul[m]%mo*invMul[n-m]%mo; } 上面的代码中用 的复杂度处理了 的逆元，处理 次 的询问的总复杂度为
- 如果 更大， 是素数，可以用Lucas定理来求解
5.3.1 Lucas定理
若 是素数，则
其中
即将 表示成 进制形式
5.3.2 推论
ll Lucas(ll n,ll m,ll p){ ll ans=1; while(n|m)ans=ans*C(n%P,m%P)%P,n/=P,m/=P; return ans; } 6. 与 的结构
6.1 的结构
6.1.1 算数基本定理 任何一个大于1的整数 ，都可以唯一地表示成素数乘积的形式
其中 是素数
对于一个较大的数，有用来分解其因数的 Pollard Rho 算法
typedef long long ll; typedef pair<int,int> pii; typedef pair<ll,ll> pll; namespace pollard_rho{ const int C=2307; const int S=10; typedef pair<ll,int> pli; mt19937 rd(time(0)); vector<ll>ve; ll gcd(ll a,ll b){return b?gcd(b,a%b):a;} ll mul(ll a,ll b,ll mod){return (__int128)a*b%mod;} ll power(ll a,ll b,ll mod){ ll res=1;a%=mod; while(b){ if(b&1)res=mul(res,a,mod); a=mul(a,a,mod); b>>=1; } return res; } bool check(ll a,ll n){ ll m=n-1,x,y; int j=0; while(!(m&1))m>>=1,j++; x=power(a,m,n); for(int i=1;i<=j;x=y,i++){ y=mul(x,x,n); if(y==1&&x!=1&&x!=n-1)return 1; } return y!=1; } bool miller_rabin(ll n){ ll a; if(n==1)return 0; if(n==2)return 1; if(!(n&1))return 0; for(int i=0;i<S;i++)if(check(rd()%(n-1)+1,n))return 0; return 1; } ll pollard_rho(ll n,int c){ ll i=1,k=2,x=rd()%n,y=x,d; while(1){ i++;x=(mul(x,x,n)+c)%n,d=gcd(y-x,n); if(d>1&&d<n)return d; if(y==x)return n; if(i==k)y=x,k<<=1; } } void findfac(ll n,int c){ if(n==1)return ; if(miller_rabin(n)){ ve.push_back(n); return ; } ll m=n; while(m==n)m=pollard_rho(n,c--); findfac(m,c); findfac(n/m,c); } vector<pli> solve(ll n){ vector<pli>res; ve.clear(); findfac(n,C); sort(ve.begin(),ve.end()); for(auto x:ve){ if(res.empty()||res.back().fi!=x)res.push_back({x,1}); else res.back().se++; } return res; } } 6.1.2 推论 若
则
6.1.3 定理 的阶乘中 的幂次为
6.1.4 定理
6.2
的结构
定理6.2.1 是循环群，即存在 ，使得
这样的 称为 的原根
素数一定有原根，原根不唯一，部分合数也有原根
原根一般不大，暴力枚举即可
1000000007的原根为5，998244353的原根为3
7. 离散对数与BSGS算法
考虑求解方程
这样的 称为离散对数，可以写为
Baby step giant step 算法
设 （ 为某常正整数），则原方程可以写成
将 存入表 (table，C++中可以用unordered_map) 中，然后枚举 ，在表中查找 即可，复杂度为 ，取 ，那么复杂度为
ll bsgs(ll a,ll b,ll p){ static unordered_map<ll,ll> tab; tab.clear(); ll u=(ll)sqrt(p)+1; ll now=1,step; rep(i,0,u-1){ ll tmp=b*inv(now,p)%p; if(!tab.count(tmp))tab[tmp]=i; (now*=a)%=p; } step=now; now=1; for(ll i=0;i<p;i+=u){ if(tab.count(now))return i+tab[now]; (now*=step)%=p; } return -1; } 8. 高次同余方程
考虑求解方程
先求 的原根 ，设 ， ，用BSGS算法求出 ，方程可写成
进而有
用扩展欧几里得算法求出 ，也就求出了
ll ModEquationSolve(ll a,ll y,ll p){ a%=p-1; ll g=primitiveRoot(p),t=bsgs(g,y,p),z,z_,d=ext_gcd(a,p-1,z,z_); if(t%d!=0)return -1; ll tmp=(p-1)/d; z=(t/d*z%tmp+tmp)%tmp; return fpow(g,z,p); } 9. 积性函数
9.1 积性函数
9.1.1 定义 积性函数 映射 如果满足：任意 且 ，有
那么称 为积性函数
最简单的积性函数 与
其中 为示性函数，当条件 成立时取1，否则取0
9.1.2 定理 积性函数必满足
9.2 欧拉函数
9.2.1 定义 欧拉函数
欧拉函数是积性函数，且有
vector<ll> factors(ll x){ vector<ll> fac; ll y=x; for(ll i=2;i*i<=x;i++){ if(y%i==0){ fac.push_back(i); while(y%i==0)y/=i; if(y==1)return fac; } } if(y!=1)fac.push_back(y); return fac; } ll Euler(ll n){ vector<ll> fac=factors(n); ll ans=n; for(auto p:fac)ans=ans/p*(p-1); return ans; } 9.2.2 欧拉定理 若 是大于1的整数， ，则
这又是一个求逆元的方法，其结果与扩展欧几里得算法相同
9.3 莫比乌斯函数
9.3.1 定义 莫比乌斯函数
ll Mobius(ll n){ vector<ll> fac=factors(n); for(auto p:fac)n/=p; return n>1?0:(fac.size()&1)?-1:1; } 9.3.2 定理 莫比乌斯反演
9.3.2' 定理 莫比乌斯反演的另一种形式
9.4 积性函数筛
积性函数可以通过前文中的线性素数筛在线性时间复杂度内求出，某些取值依赖于自变量的素因子的数论函数也能通过其求出
#define rep(i,a,b) for(ll i=a;i<=b;i++) const ll N=10005; //积性数论函数筛 void Number_Theory_Function_Solve(vector<ll> &a,ll n,string function_name){ static bool vis[N]; static ll prime[N]; a.resize(n+1); ll tot=0; rep(i,1,n)vis[i]=0; a[1]=1; if(function_name=="IsPrime"){ /* { 1, if x is a prime f(x) = { { 0, otherwise 素数判断函数不是积性函数，但可以用积性函数筛求解 */// rep(i,2,n){ if(!vis[i])prime[tot++]=i,a[i]=1; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ a[i*prime[j]]=0; break; } a[i*prime[j]]=0; } } a[1]=0; }else if(function_name=="E"){ /* { 1, x=1 f(x) = I(x=1) = { { 0, otherwise 积性函数群中的单位元 */// rep(i,2,n)a[i]=0; }else if(function_name=="Euler"){ /* n f(x) = sum I(gcd(i,x)=1) = x product (1-1/p) i=1 p|x, p is Prime 欧拉函数 */// rep(i,2,n){ if(!vis[i])prime[tot++]=i,a[i]=i-1; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ a[i*prime[j]]=a[i]*prime[j]; break; } a[i*prime[j]]=a[i]*a[prime[j]]; } } }else if(function_name=="Mobius"){ /* { 1, x=1 f(x) = { (-1)^k, n=p1 p2 ... pk, pi is prime { 0, otherwise 莫比乌斯函数 */// rep(i,2,n){ if(!vis[i])prime[tot++]=i,a[i]=-1; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ a[i*prime[j]]=0; break; } a[i*prime[j]]=a[i]*a[prime[j]]; } } }else if(function_name=="FactorCounter"){ /* n f(x) = sum I(i|x) i=1 因子个数 */// static ll t[N]; rep(i,2,n){ if(!vis[i])prime[tot++]=i,a[i]=2,t[i]=1; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ a[i*prime[j]]=a[i]/(t[i]+1)*(t[i]+2); t[i*prime[j]]=t[i]+1; break; } a[i*prime[j]]=a[i]*a[prime[j]]; t[i*prime[j]]=1; } } }else if(function_name=="FactorSum"){ /* n f(x) = sum I(i|x)*i i=1 因子和 */// static ll pw[N]; rep(i,2,n){ if(!vis[i])prime[tot++]=i,a[i]=i+1,pw[i]=i; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ a[i*prime[j]]=(pw[i]*prime[j]*prime[j]-1)/(prime[j]-1)*a[i/pw[i]]; pw[i*prime[j]]=pw[i]*prime[j]; break; } a[i*prime[j]]=a[i]*a[prime[j]]; pw[i*prime[j]]=prime[j]; } } }else if(function_name=="MinFactor"){ rep(i,2,n){ if(!vis[i])prime[tot++]=i,a[i]=i; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ a[i*prime[j]]=prime[j]; break; } a[i*prime[j]]=prime[j]; } } }else{ assert(false); } } 9.5 数论函数的狄利克雷乘积
9.5.1 定义 数论函数的狄利克雷乘积
9.5.2 定理 所有积性数论函数构成一个群，该群有单位元
任意积性数论函数 ，都存在唯一的积性数论函数 ，使得
且
称 为 的狄利克雷逆函数
9.5.3 推论 数论函数间的关系 设 为欧拉函数， 为莫比乌斯函数， 为单位数论函数， 为常数函数（ ）， 为恒等函数（ ），那么
bool mulFunCheck(const vector<ll> &val){ ll n=val.size()-1; static bool vis[N]; static ll prime[N],pw[N],t[N],tot=0; if(val[1]!=1)return false; rep(i,1,n)vis[i]=0; rep(i,2,n){ if(!vis[i])prime[tot++]=i,pw[i]=i,t[i]=1; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ if(val[i*prime[j]]!=val[i/pw[i]]*val[pw[i]*prime[j]])return false; pw[i*prime[j]]=pw[i]*prime[j],t[i*prime[j]]=t[i]+1; break; } if(val[i*prime[j]]!=val[i]*val[prime[j]])return false; pw[i*prime[j]]=prime[j],t[i*prime[j]]=1; } } return true; } //狄利克雷乘积 h(n)=sum_{d|n}f(d)g(n/d) vector<ll> DirichletProduct(const vector<ll> &f,const vector<ll> &g){ ll n=f.size()-1; vector<ll> h(n+1,0); rep(i,1,n)h[i]=0; rep(i,1,n)rep(j,1,n/i)h[i*j]+=f[i]*g[j]; return h; } //狄利克雷逆函数 f*g=g*f=[n==1] vector<ll> DirichletInversion(const vector<ll> &f){ ll n=f.size()-1; vector<ll> g(n+1,0); g[1]=1; rep(i,1,n)rep(j,2,n/i)g[i*j]-=f[j]*g[i]; return g; } //Mobius反演 f(n)=sum(g(d),d|n) g(n)=sum(f(d)mu(n/d),d|n) vector<ll> MobiusInversion(const vector<ll> &f){ ll n=f.size()-1; static vector<ll> mu; Number_Theory_Function_Solve(mu,n,"Mobius"); return DirichletProduct(mu,f); } 10. 多项式乘积算法
试想这样一个问题，求两个多项式
的乘积
使用传统的方法至少需要 的复杂度，如何加速上述过程？
10.1 快速傅里叶变换
首先考虑如何用其他方式表示多项式 .
任取n个不同的数(可以是整数、实数，甚至是复数)
将其代入 中，就得到一个线性方程组
只要 足够大，就能够唯一地确定一个多项式，换言之，上述方程组可以表示一个多项式，将这两种多项式的表示方法分别称为系数表示和点值表示.
利用快速傅里叶变换来求多项式乘积的总体思路是
- 选取合适的 个不同的数
- 将多项式 与 转化为点值表示(称为离散傅里叶变换，简称 )
- 计算 的点值表示
- 将 转化为系数表示(称为逆离散傅里叶变换，简称 )
1. 选取合适的 个不同的数
我们选取复数域上 的 个不同的值(或称 个 次单位复根)作为 的值，即
至于指数形式的复数 ，用大家所熟知的欧拉公式即可求得其代数形式
经过简单计算可知
当 为偶数时
其中 为 除以 的余数，上述两等式将在后文中使用.
我们为什么要费尽周折选取如此复杂的点呢?是为了使用快速傅里叶变换.
2. 将多项式 与 转化为点值表示
考虑多项式 ，当 时(当不满足该条件时，向 补充系数为0的高次项来扩大 使其满足该条件)，将其化为两个多项式
则有
进而
也就是说，要求 在 个不同点处的值，只需要求 与 在 个不同点处的值，由于 ，可对 与 重复进行上述过程，最终经过 步后得到 个函数
之后回推得到 的点值表示，上述过程就是快速傅里叶变换的过程，复杂度为 即 .
当然，还需要对 进行同样的变换.
3. 计算 的点值表示
点值表示的优点是可以快速地求出两个选取了相同点值的多项式的乘积，例如多项式
与多项式
的乘积
只需要 的复杂度即可求得.
4. 将 转化为系数表示
下面以 为例，讲解如何将多项式从点值表示转化为系数表示，此过程又称多项式的插值.
将 的点值表示写成矩阵形式 即
此处矩阵 中的1可视为 .
现在我们已知的是 和 ，要求的是 ， 是一范德蒙德矩阵，可求得其逆矩阵
因此 即
也就是说，只需将 与 对换，将 换成 ，再乘上系数 ，进行类似步骤2的变换，即可进行逆快速傅里叶变换，算法复杂度同样为 .
按照上述方法将 转化为系数表示，本题得解.
#include<bits/stdc++.h> using namespace std; const int maxn=1<<18; //typedef complex<double> C; struct C { double a,b; C(){} C(double a,double b):a(a),b(b){} C operator = (int a){*this=C(a*1.0,0);return *this;} C operator + (const C &t){return C(a+t.a,b+t.b);} C operator - (const C &t){return C(a-t.a,b-t.b);} C operator * (const C &t){return C(a*t.a-b*t.b,a*t.b+b*t.a);} }; C wn(int n,int f) { return C(cos(acos(-1.0)/n),f*sin((acos(-1.0))/n)); } C inv(int n) { return C(1.0/n,0); } C a[maxn],b[maxn],c[maxn]; int g[maxn]; void FFT(C *a,int n,int f) { for(int i=0;i<n;i++)if(i>g[i])swap(a[i],a[g[i]]); for(int i=1;i<n;i<<=1) { C w=wn(i,f),x,y; for(int j=0;j<n;j+=i+i) { C e;e=1; for(int k=0;k<i;e=e*w,k++) { x=a[j+k]; y=a[j+k+i]*e; a[j+k]=x+y; a[j+k+i]=x-y; } } } if(f==-1) { C Inv=inv(n); for(int i=0;i<n;i++)a[i]=a[i]*Inv; } } void conv(C *a,int n,C *b,int m,C *c) { int k=0,s=2; while((1<<k)<max(n,m)+1)k++,s<<=1; for(int i=1;i<s;i++)g[i]=(g[i/2]/2)|((i&1)<<k); FFT(a,s,1); FFT(b,s,1); for(int i=0;i<s;i++)c[i]=a[i]*b[i]; FFT(c,s,-1); } int main() { int n,m; scanf("%d%d",&n,&m); for(int i=0;i<=n;i++)scanf("%lf",&a[i]); for(int i=0;i<=m;i++)scanf("%lf",&b[i]); conv(a,n,b,m,c); for(int i=0;i<=n+m;i++)printf("%d ",(int)(c[i].a+0.5)); return 0; } 10.2 快速数论变换
快速傅里叶变换可以用于计算两个实系数多项式的乘积，但毕竟有精度问题，对于整系数多项式和模 意义下的乘法，是否有无损精度的算法呢？
FFT中精度损失的关键出在第一步—— 的选取，算法的实现依赖一个基本的等式
这个等式导致， 有 个不同的取值，但它们的平方只有 个不同的取值，这 个值的平方只有 个不同的取值，以此类推
若 满足某些条件时，是可以找到整数 的，例如当 ， 时，取
此时恰好满足类似的性质，后续的算法步骤类似，但这种变换就不是离散傅里叶变换（DFT）了，而是数论变换（NTT）
const ll mo=998244353; ll fpow(ll a, ll b){ ll ans=1; while(b>0){if(b&1)ans=ans*a%mo;b>>=1;a=a*a%mo;} return ans; } ll D(ll x) { ((x>=mo) && (x-=mo)) || ((x<0) && (x+=mo)); return x; } void NTT(ll a[],ll n,ll op) { for(ll i=1,j=n>>1;i<n-1;++i) { if(i<j) swap(a[i],a[j]); ll k=n>>1; while(k<=j) { j-=k; k>>=1; } j+=k; } for(ll len=2;len<=n;len<<=1) { ll rt=fpow(3,(mo-1)/len); for(ll i=0;i<n;i+=len) { ll w=1; for(ll j=i;j<i+len/2;++j) { ll u=a[j],t=1LL*a[j+len/2]*w%mo; a[j]=D(u+t),a[j+len/2]=D(u-t); w=1LL*w*rt%mo; } } } if(op==-1) { reverse(a+1,a+n); ll in=fpow(n,mo-2); for(ll i=0;i<n;++i) a[i]=1LL*a[i]*in%mo; } } vector<ll> Conv(vector<ll> const &A,vector<ll> const &B,ll N) { static ll a[2000005],b[2000005]; auto Make2=[](ll x)->ll { return 1<<((32-__builtin_clz(x))+((x&(-x))!=x)); }; ll n=Make2(A.size()+B.size()-1); for(ll i=0;i<n;++i) { a[i]=i<A.size()?A[i]:0; b[i]=i<B.size()?B[i]:0; } NTT(a,n,1);NTT(b,n,1); for(ll i=0;i<n;++i) a[i]=1LL*a[i]*b[i]%mo; NTT(a,n,-1); vector<ll> C(N); for (ll i=0;i<N;i++) C[i]=a[i]; return C; } 杂技
1. 积性函数前缀和
1.1 min25筛
用于求积性函数 的前 项和，要求 在素数处为一多项式，在素数的整数幂处可以直接求值。
#define rep(i,a,b) for(ll i=a;i<=b;i++) typedef long long ll; class MultiplicativeFunction{ private: static const ll SQRTN=1100000; ll prime[SQRTN],pcnt,cnt; bool vis[SQRTN]; ll d,n; ll val[SQRTN*2],id1[SQRTN],id2[SQRTN]; ll G[SQRTN*2],fpval[SQRTN*2]; ll g[SQRTN*2]; function<ll(ll,ll)> f_prime_pow; ll sum_pow_n_k(ll n,ll k){ if(k==0)return n; if(k==1)return n*(n+1)/2; if(k==2)return n*(n+1)*(2*n+1)/6; return -1; } void getPrime(ll n){ pcnt=0; fill(vis+1,vis+1+d,0); rep(i,2,n){ if(!vis[i])prime[++pcnt]=i; rep(j,1,pcnt){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0)break; } } } ll powll(ll x,ll k){ ll ans=1; while(k--)ans*=x; return ans; } ll &index(ll x){return x<=d?id1[x]:id2[n/x];} void sieve(ll x,ll k){ rep(i,1,pcnt)fpval[i]=powll(prime[i],k)+fpval[i-1]; rep(i,1,cnt)g[i]=sum_pow_n_k(val[i],k)-1; rep(i,1,pcnt){ if(prime[i]*prime[i]>n)break; rep(j,1,cnt){ if(prime[i]*prime[i]>val[j])break; g[j]=g[j]-powll(prime[i],k)*(g[index(val[j]/prime[i])]-fpval[i-1]); } } } ll S(ll x,ll now){ if(x<=1 || x<prime[now])return 0; ll ans=G[index(x)]-fpval[now-1]; rep(i,now,pcnt){ if(prime[i]*prime[i]>x)break; ll prod=prime[i]; for(ll k=1;prod*prime[i]<=x;k++,prod*=prime[i]){ ans+=f_prime_pow(prime[i],k)*S(x/prod,i+1)+f_prime_pow(prime[i],k+1); } } return ans; } public: ll solve(ll x,const vector<ll> &c,function<ll(ll,ll)> fp){ f_prime_pow=fp; n=x;d=sqrt(n);cnt=0; getPrime(d); for(ll i=1;i<=n;i=n/(n/i)+1){ val[++cnt]=n/i; index(val[cnt])=cnt; } fill(G+1,G+1+cnt,0); for(int k=0;k<(int)c.size();k++){ if(c[k]!=0){ sieve(x,k); rep(i,1,cnt)G[i]+=c[k]*g[i]; } } rep(i,1,pcnt){ fpval[i]=0; for(ll k=(ll)c.size()-1;k>=0;k--){ fpval[i]=fpval[i]*prime[i]+c[k]; } fpval[i]+=fpval[i-1]; } return S(n,1)+1; } }; 使用说明：调用 solve 函数，c 是函数 在素数处的解析式（多项式）的系数，fp 是 在素数的整数幂处的取值计算函数。
例：
MultiplicativeFunction f; //Euler函数 cout<<f.solve(n,{-1,1},[](ll p,ll k){ll a=1;while(k--)a*=p;return a/p*(p-1);})<<endl; //Mobius函数 cout<<f.solve(n,{-1},[](ll p,ll k){return k==1?-1:0;})<<endl; 1.2 另一种筛
用于求积性函数 的前 项和，要求找一个函数 ，使得 与 的狄利克雷乘积、 可以快速求前 项和。
class sieve{ private: static const ll SQRTN=100005; ll n,d,cnt; ll val[SQRTN*2],id1[SQRTN],id2[SQRTN],s[SQRTN*2]; vector<ll> sum; bool vis[SQRTN*2]; ll &index(ll x){return x<=d?id1[x]:id2[n/x];} function<ll(ll)> sum_g,sum_f_g; ll S(ll x){ if(x<sum.size())return sum[x]; if(x==1)return 1; ll p=index(x); if(vis[p])return val[p]; val[p]=sum_f_g(x); for(ll i=2;i<=x;i=x/(x/i)+1){ val[p]-=(sum_g(x/(x/i))-sum_g(i-1))*S(x/i); } vis[p]=1; return val[p]; } public: void init(function<ll(ll)> sg,function<ll(ll)> sfg,const vector<ll> &v={}){sum_g=sg;sum_f_g=sfg;sum=v;} ll operator () (ll x){ n=x; d=sqrt(x); cnt=0; for(ll i=1;i<=x;i=x/(x/i)+1){ val[++cnt]=x/i; vis[cnt]=0; index(val[cnt])=cnt; } return S(n); } }; 使用说明：调用 init 函数进行初始化，sg 为 的前 项和，sfg 为 的前 项和，如果需要，可以预处理前若干项的和放到 v 中进行剪枝；直接使用 () 进行求和。
例：
sieve f; //Mobius函数 mu(x) g=1 f*g=E f.init([](ll x){return x;},[](ll x){return x>=1;}); cout<<f(n)<<endl; //Euler函数 phi(x) g=1 f*g=id f.init([](ll x){return x;},[](ll x){return x*(x+1)/2;}); cout<<f(n)<<endl; //x^k phi(x) g=id^k f*g=id^(k+1) f.init([](ll x){return x*(x+1)/2;},[](ll x){return x*(x+1)*(2*x+1)/6;}); cout<<f(n)<<endl; //x^k mu(x) g=id^k f*g=E f.init([](ll x){return x*(x+1)/2;},[](ll x){return x>=1;}); cout<<f(n)<<endl; 一个简化版本，使用方式相同：
class sieve2{ private: vector<ll> sum; unordered_map<ll,ll> mp; function<ll(ll)> sum_g,sum_f_g; ll S(ll x){ if(x<sum.size())return sum[x]; if(mp.find(x)!=mp.end())return mp[x]; ll ans=sum_f_g(x); for(ll i=2;i<=x;i=x/(x/i)+1){ ans-=(sum_g(x/(x/i))-sum_g(i-1))*S(x/i); } mp[x]=ans; return ans; } public: void init(function<ll(ll)> sg,function<ll(ll)> sfg,const vector<ll> &v={0,1}){ mp.clear(); sum_g=sg;sum_f_g=sfg;sum=v; } ll operator () (ll x){ return S(x); } }; 2. 二次剩余
Shanks 算法，求解二次同余方程 ， 是素数。
ll fpow(ll a,ll b,ll c){a%=c;ll ans=1;while(b>0){if(b&1)ans=ans*a%c;b>>=1;a=a*a%c;}return ans;} bool IsQuadraticResidue(ll a,ll p){//判断a是否为p的二次剩余 p是素数 return fpow(a,(p-1)>>1,p)==1; } ll Shanks(ll a,ll p){//求解二次同余方程x^2=a(mod p) p是素数 if(a==0)return 0; ll q=p-1,e=0; while(!(q&1))q>>=1,e++; static mt19937_64 rd(time(0)); ll n=rd()%(p-1)+1;//随机选取p的一个非二次剩余，若p为定值，n也可为定值 while(IsQuadraticResidue(n,p))n=rd()%(p-1)+1; ll z=fpow(n,q,p),y=z,r=e,x=fpow(a,(q-1)>>1,p),b=a*x%p*x%p; x=a*x%p; while(b!=1){ ll temp=b*b%p,m=1; while(temp!=1)(temp*=temp)%=p,m++; if(m==r)return -1; ll t=y; rep(i,1,r-m-1)(t*=t)%=p; y=t*t%p,r=m%p,x=x*t%p,b=b*y%p; } return x; } 3. 组合数
在 Lucas 定理的基础上修改，可以求 除以 的余数， 可以非素数。
class Combine{ private: ll fpow(ll a,ll b,ll c){a%=c;ll ans=1;while(b>0){if(b&1)ans=ans*a%c;b>>=1;a=a*a%c;}return ans;} void factorization(ll x,vector<ll> &p,vector<ll> &w,vector<ll> &pw){ p.clear(),w.clear(),pw.clear(); ll y=x; for(ll i=2;i*i<=x;i++)if(y%i==0){ p.push_back(i),w.push_back(0),pw.push_back(1); while(y%i==0)y/=i,w.back()++,pw.back()*=i; } if(y!=1)p.push_back(y),w.push_back(1),pw.push_back(y); } ll powll(ll a,ll b){ll ans=1;while(b)ans*=a,b--;return ans;} ll ext_gcd(ll a,ll b,ll& x,ll& y){ ll d=a; if(!b)x=1,y=0; else d=ext_gcd(b,a%b,y,x),y-=a/b*x; return d; } ll inv(ll a,ll b){ ll x,y,d=ext_gcd(a,b,x,y); assert(d==1); return (x%b+b)%b; } ll crt(const vector<ll> &m,const vector<ll> &b){ ll ans=0,M=1; for(auto i:m)M*=i; for(int i=0;i<(int)m.size();i++){ (ans+=M/m[i]*inv(M/m[i],m[i])%M*b[i])%=M; } return ans; } public: ll mod; vector<ll> fac[105],p,w,pw,res; void init(ll mo){ mod=mo; factorization(mod,p,w,pw); res.resize(p.size()); for(int i=0;i<(int)p.size();i++){ fac[i].resize(pw[i]); fac[i][0]=1; for(int j=1;j<pw[i];j++){ fac[i][j]=(j%p[i]==0)?fac[i][j-1]:(fac[i][j-1]*j%pw[i]); } } } Combine(){} Combine(ll mod):mod(mod){init(mod);} ll factorial(ll n,ll i){ if(n<p[i])return fac[i][n]; else return fpow(fac[i][pw[i]-1],n/pw[i],pw[i])*fac[i][n%pw[i]]%pw[i]*factorial(n/p[i],i)%pw[i]; } ll C_pw(ll n,ll m,ll i){ ll c=0; for(ll j=p[i];j<=n;j*=p[i])c+=n/j; for(ll j=p[i];j<=m;j*=p[i])c-=m/j; for(ll j=p[i];j<=n-m;j*=p[i])c-=(n-m)/j; return factorial(n,i)*inv(factorial(m,i),pw[i])%pw[i]*inv(factorial(n-m,i),pw[i])%pw[i]*fpow(p[i],c,pw[i])%pw[i]; } ll operator () (ll n,ll m){ assert(m<=n && m>=0); for(int i=0;i<(int)p.size();i++)res[i]=C_pw(n,m,i); return crt(pw,res); } }C; 4. 一个求和技巧（常被称为类欧几里得算法）
求 ， 表示取整。
ll sum_pow(ll n,ll k){ if(k==0)return n; else if(k==1)return n*(n+1)/2; else if(k==2)return n*(n+1)*(2*n+1)/6; else if(k==3)return n*n*(n+1)*(n+1)/4; else if(k==4)return n*(2*n+1)*(n+1)*(3*n*n+3*n-1)/30; else assert(false); } ll EuclidLike1(ll a,ll b,ll c,ll n){ if(a==0)return b/c*(n+1); else if(a>=c || b>=c)return (a/c)*sum_pow(n,1)+(b/c)*(n+1)+EuclidLike1(a%c,b%c,c,n); else return (a*n+b)/c*n-EuclidLike1(c,c-b-1,a,(a*n+b)/c-1); } 5. 模数任意的 FFT
模数任意（可取 甚至非素数）的 FFT，你可以忘掉 NTT 了。
模板需要 define 相应的关键字（FFT_conv,FFT_conv_mod,FFT_conv_mod_fast）来启用相关功能。
namespace FFT{ // FFT_conv:普通多项式乘积 // FFT_conv_mod:多项式乘积取模 // FFT_conv_mod_fast:优化的多项式乘积取模 #define FFT_conv_mod_fast const int maxn=1<<20; //typedef complex<double> cp; struct cp{ long double a,b; cp(){} cp(long double a,long double b):a(a),b(b){} cp(ll a):a(a*1.0),b(0){} cp operator + (const cp &t)const{return cp(a+t.a,b+t.b);} cp operator - (const cp &t)const{return cp(a-t.a,b-t.b);} cp operator * (const cp &t)const{return cp(a*t.a-b*t.b,a*t.b+b*t.a);} cp conj()const{return cp(a,-b);} }; cp wn(int n,int f){ const long double pi=acos(-1.0); return cp(cos(pi/n),f*sin(pi/n)); } int g[maxn]; void DFT(cp *a,int n,int f){ for(int i=0;i<n;i++)if(i>g[i])swap(a[i],a[g[i]]); for(int i=1;i<n;i<<=1){ cp w=wn(i,f),x,y; for(int j=0;j<n;j+=i+i){ cp e(1,0); for(int k=0;k<i;e=e*w,k++){ x=a[j+k]; y=a[j+k+i]*e; a[j+k]=x+y; a[j+k+i]=x-y; } } } if(f==-1){ cp Inv(1.0/n,0); for(int i=0;i<n;i++)a[i]=a[i]*Inv; } } #ifdef FFT_conv cp a[maxn],b[maxn],c[maxn]; vector<ll> conv(const vector<ll> &u,const vector<ll> &v,ll sz){ int k=0,s=1,n=(int)u.size()-1,m=(int)v.size()-1; while(s<n+m+1)k++,s<<=1; for(int i=1;i<s;i++)g[i]=(g[i/2]/2)|((i&1)<<(k-1)); for(int i=0;i<s;i++)a[i]=i<=n?u[i]:0; for(int i=0;i<s;i++)b[i]=i<=m?v[i]:0; DFT(a,s,1); DFT(b,s,1); for(int i=0;i<s;i++)c[i]=a[i]*b[i]; DFT(c,s,-1); vector<ll> ans(sz+1,0); for(int i=0;i<=sz;i++)ans[i]=i<s?ll(round(c[i].a)):0; return ans; } #endif #ifdef FFT_conv_mod cp a1[maxn],a2[maxn],b1[maxn],b2[maxn],c[maxn]; vector<ll> conv_mod(const vector<ll> &u,const vector<ll> &v,int sz,ll mod){ int k=0,s=1,n=(int)u.size()-1,m=(int)v.size()-1,M=sqrt(mod)+1; while(s<n+m+1)k++,s<<=1; for(int i=1;i<s;i++)g[i]=(g[i/2]/2)|((i&1)<<(k-1)); for(int i=0;i<s;i++){ if(i<=n)a1[i]=u[i]%mod/M,b1[i]=u[i]%mod%M; else a1[i]=b1[i]=0; if(i<=m)a2[i]=v[i]%mod/M,b2[i]=v[i]%mod%M; else a2[i]=b2[i]=0; } DFT(a1,s,1); DFT(b1,s,1); DFT(a2,s,1); DFT(b2,s,1); vector<ll> ans(sz+1,0); for(int i=0;i<s;i++)c[i]=a1[i]*a2[i]; DFT(c,s,-1); for(int i=0;i<min(sz+1,s);i++)(ans[i]+=ll(round(c[i].a))%mod*M%mod*M)%=mod; for(int i=0;i<s;i++)c[i]=a1[i]*b2[i]+a2[i]*b1[i]; DFT(c,s,-1); for(int i=0;i<min(sz+1,s);i++)(ans[i]+=ll(round(c[i].a))%mod*M)%=mod; for(int i=0;i<s;i++)c[i]=b1[i]*b2[i]; DFT(c,s,-1); for(int i=0;i<min(sz+1,s);i++)(ans[i]+=ll(round(c[i].a))%mod)%=mod; return ans; } #endif #ifdef FFT_conv_mod_fast cp a[maxn],b[maxn],Aa[maxn],Ab[maxn],Ba[maxn],Bb[maxn]; vector<ll> conv_mod_fast(const vector<ll> &u,const vector<ll> &v,int sz,ll mod){ int k=0,s=1,n=(int)u.size()-1,m=(int)v.size()-1,M=sqrt(mod)+1; while(s<n+m+1)k++,s<<=1; for(int i=1;i<s;i++)g[i]=(g[i/2]/2)|((i&1)<<(k-1)); for(int i=0;i<s;i++){ if(i<=n)a[i].a=u[i]%mod%M,a[i].b=u[i]%mod/M; else a[i]=0; if(i<=m)b[i].a=v[i]%mod%M,b[i].b=v[i]%mod/M; else b[i]=0; } DFT(a,s,1); DFT(b,s,1); for(int i=0;i<s;i++){ int j=(s-i)%s; cp t1=(a[i]+a[j].conj())*cp(0.5,0); cp t2=(a[i]-a[j].conj())*cp(0,-0.5); cp t3=(b[i]+b[j].conj())*cp(0.5,0); cp t4=(b[i]-b[j].conj())*cp(0,-0.5); Aa[i]=t1*t3,Ab[i]=t1*t4,Ba[i]=t2*t3,Bb[i]=t2*t4; } for(int i=0;i<s;i++){ a[i]=Aa[i]+Ab[i]*cp(0,1); b[i]=Ba[i]+Bb[i]*cp(0,1); } DFT(a,s,-1); DFT(b,s,-1); vector<ll> ans(sz+1,0); for(int i=0;i<min(s,sz+1);i++){ ll t1=(ll)round(a[i].a)%mod; ll t2=(ll)round(a[i].b)%mod; ll t3=(ll)round(b[i].a)%mod; ll t4=(ll)round(b[i].b)%mod; (ans[i]+=t1+(t2+t3)*M%mod+t4*M*M)%=mod; } return ans; } #endif } 6. 多项式插值
6.1 Lagrange 插值
#define rep(i,a,b) for(ll i=a;i<=b;i++) const ll mo=1e9+7; ll fpow(ll a,ll b){ll ans=1;while(b>0){if(b&1)ans=ans*a%mo;b>>=1;a=a*a%mo;}return ans;} ll lagrange(vector<ll> x,vector<ll> y,ll X){ auto p=y.begin(); ll ans=0; for(auto k:x){ ll a=*p++%mo,b=1; for(auto j:x)if(j!=k)a=(X-j)%mo*a%mo,b=(k-j)%mo*b%mo; ans=(ans+mo+a*fpow(b,mo-2)%mo)%mo; } return ans; } 插值点连续的情况下， ，时间复杂度优化到 ， 是插值点的个数。
#define rep(i,a,b) for(ll i=a;i<=b;i++) ll Lagrange(const vector<ll> &y,ll x,ll mod){ ll n=y.size()-1; static const ll N=100005; static ll lx[N],rx[N],inv[N],mul[N]; rep(i,1,n)inv[i]=(i==1)?1:((mod-mod/i)*inv[mod%i]%mod); rep(i,0,n)lx[i]=rx[i]=(x-1-i)%mod; rep(i,1,n)(lx[i]*=lx[i-1])%=mod; for(ll i=n-1;i>=0;i--)(rx[i]*=rx[i+1])%=mod; rep(i,0,n)mul[i]=(i==0)?1:(mul[i-1]*inv[i])%mod; ll ans=0; rep(i,0,n){ ll now=1; if(i>0)(now*=lx[i-1]*mul[i]%mod)%=mod; if(i<n)(now*=rx[i+1]*mul[n-i]%mod)%=mod; if((n-i)&1)now=-now; (ans+=y[i]*now)%=mod; } return ans+(ans<0)*mod; } 二维的 Lagrange 插值
#define rep(i,a,b) for(ll i=a;i<=b;i++) type lagrange2(vector<type> x,vector<type> y,vector<vector<type> > z,type X,type Y){ int M=x.size()-1,N=y.size()-1; type ans=0; rep(m,0,M)rep(n,0,N){ type t=z[m][n]; rep(i,0,M)if(i!=m)t*=(X-x[i])/(x[m]-x[i]); rep(i,0,N)if(i!=n)t*=(Y-y[i])/(y[n]-y[i]); ans+=t; } return ans; } 6.2 牛顿插值
#define rep(i,a,b) for(ll i=a;i<=b;i++) ll fpow(ll a,ll b){ll ans=1;while(b>0){if(b&1)ans=ans*a%mo;b>>=1;a=a*a%mo;}return ans;} class NewtonPoly{ public: ll f[105],d[105],x[105],n=0; void add(ll X,ll Y){ x[n]=X,f[n]=Y%mo; rep(i,1,n)f[n-i]=(f[n-i+1]-f[n-i])%mo*fpow((x[n]-x[n-i])%mo,mo-2)%mo; d[n++]=f[0]; } ll cal(ll X){ ll ans=0,t=1; rep(i,0,n-1)ans=(ans+d[i]*t)%mo,t=(X-x[i])%mo*t%mo; return ans+mo*(ans<0); } }; 6.3 利用插值求 前 项和
#define rep(i,a,b) for(ll i=a;i<=b;i++) const ll mo=1e9+7; ll sum_pow_n_k(ll n,ll k){ static ll inv[1005],d[1005][1005],init_flag=1; if(init_flag){ init_flag=0; const ll all=1002; inv[1]=1; rep(i,2,all)inv[i]=(mo-mo/i)*inv[mo%i]%mo; d[1][1]=1; rep(i,2,all)rep(j,1,i)d[i][j]=(j*j*d[i-1][j]+(j-1)*d[i-1][j-1])%mo*inv[j]%mo; } n%=mo; ll ans=0,now=n; rep(i,1,k+1){ (ans+=d[k+1][i]*now)%=mo; (now*=(n-i))%=mo; } return ans+(ans<0)*mo; } 7. 线性递推数列
7.1 BM 算法
根据数列前若干项求递推公式。
例如斐波那契数列，传入 {1,1,2,3,5,8,13}，传出 {1,-1,-1}（ ）。
#define rep(i,a,b) for(ll i=a;i<=b;i++) const ll N=25; const ll mo=1e9+7; ll fpow(ll a,ll b){ll ans=1;while(b>0){if(b&1)ans=ans*a%mo;b>>=1;a=a*a%mo;}return ans;} //BM算法 求线性递推数列的递推公式 vector<ll> BM(const vector<ll> &s){ vector<ll> C={1},B={1},T; ll L=0,m=1,b=1; rep(n,0,s.size()-1){ ll d=0; rep(i,0,L)d=(d+s[n-i]%mo*C[i])%mo; if(d==0)m++; else{ T=C; ll t=mo-fpow(b,mo-2)*d%mo; while(C.size()<B.size()+m)C.push_back(0); rep(i,0,B.size()-1)C[i+m]=(C[i+m]+t*B[i])%mo; if(2*L>n)m++; else L=n+1-L,B=T,b=d,m=1; } } return C; } 7.2 线性递推数列求单项值
利用多项式求 ，复杂度为 ， 为递推公式长度。
s 为前若干项，C 为递推公式系数， 从 开始。
#define rep(i,a,b) for(ll i=a;i<=b;i++) ll polySolve(const vector<ll> &s,const vector<ll> &C,ll n){ if(n<s.size())return s[n]; static ll g[N],f[N],d[N]; ll k=(ll)C.size()-1,w=1; rep(i,0,k)f[i]=i==1,d[i]=i==k?1:C[k-i]; while((w<<1)<=n)w<<=1; while(w>>=1){ rep(i,0,k+k-2)g[i]=0; rep(i,0,k-1)if(f[i])rep(j,0,k-1)(g[i+j]+=f[i]*f[j])%=mo; for(ll i=k+k-2;i>=k;i--)if(g[i])rep(j,1,k)(g[i-j]-=g[i]*d[k-j])%=mo; rep(i,0,k-1)f[i]=g[i]; if(w&n)for(ll i=k;i>=0;i--)f[i]=i==k?f[i-1]:(i==0?-f[k]*d[i]:(f[i-1]-f[k]*d[i]))%mo; } ll ans=0; rep(i,0,k-1)(ans+=f[i]*s[i])%=mo; return ans+(ans<0)*mo; } 7.3 RS 算法
当模数非素数时，使用 RS 算法替代 BM 算法。
模板是网上搜刮来的，它终结了线性递推题目在ACM中的历史，模板作者留了一个 FFT 的优化接口，供下一个毒瘤出题人使用。
struct LinearSequence { typedef vector<long long> vec; static void extand(vec &a, size_t d, ll value = 0) { if (d <= a.size()) return; a.resize(d, value); } static vec BerlekampMassey(const vec &s, ll mod) { function<ll(ll)> inverse = [&](ll a) { return a == 1 ? 1 : (ll) (mod - mod / a) * inverse(mod % a) % mod; }; vec A = {1}, B = {1}; ll b = s[0]; for (size_t i = 1, m = 1; i < s.size(); ++i, m++) { ll d = 0; for (size_t j = 0; j < A.size(); ++j) { d += A[j] * s[i - j] % mod; } if (!(d %= mod)) continue; if (2 * (A.size() - 1) <= i) { auto temp = A; extand(A, B.size() + m); ll coef = d * inverse(b) % mod; for (size_t j = 0; j < B.size(); ++j) { A[j + m] -= coef * B[j] % mod; if (A[j + m] < 0) A[j + m] += mod; } B = temp, b = d, m = 0; } else { extand(A, B.size() + m); ll coef = d * inverse(b) % mod; for (size_t j = 0; j < B.size(); ++j) { A[j + m] -= coef * B[j] % mod; if (A[j + m] < 0) A[j + m] += mod; } } } return A; } static void exgcd(ll a, ll b, ll &g, ll &x, ll &y) { if (!b) x = 1, y = 0, g = a; else { exgcd(b, a % b, g, y, x); y -= x * (a / b); } } static ll crt(const vec &c, const vec &m) { int n = c.size(); ll M = 1, ans = 0; for (int i = 0; i < n; ++i) M *= m[i]; for (int i = 0; i < n; ++i) { ll x, y, g, tm = M / m[i]; exgcd(tm, m[i], g, x, y); ans = (ans + tm * x * c[i] % M) % M; } return (ans + M) % M; } static vec ReedsSloane(const vec &s, ll mod) { auto inverse = [](ll a, ll m) { ll d, x, y; exgcd(a, m, d, x, y); return d == 1 ? (x % m + m) % m : -1; }; auto L = [](const vec & a, const vec & b) { int da = (a.size() > 1 || (a.size() == 1 && a[0])) ? a.size() - 1 : -1000; int db = (b.size() > 1 || (b.size() == 1 && b[0])) ? b.size() - 1 : -1000; return max(da, db + 1); }; auto prime_power = [&](const vec & s, ll mod, ll p, ll e) { // linear feedback shift register mod p^e, p is prime vector<vec> a(e), b(e), an(e), bn(e), ao(e), bo(e); vec t(e), u(e), r(e), to(e, 1), uo(e), pw(e + 1);; pw[0] = 1; for (int i = pw[0] = 1; i <= e; ++i) pw[i] = pw[i - 1] * p; for (ll i = 0; i < e; ++i) { a[i] = {pw[i]}, an[i] = {pw[i]}; b[i] = {0}, bn[i] = {s[0] * pw[i] % mod}; t[i] = s[0] * pw[i] % mod; if (t[i] == 0) { t[i] = 1, u[i] = e; } else { for (u[i] = 0; t[i] % p == 0; t[i] /= p, ++u[i]); } } for (size_t k = 1; k < s.size(); ++k) { for (int g = 0; g < e; ++g) { if (L(an[g], bn[g]) > L(a[g], b[g])) { ao[g] = a[e - 1 - u[g]]; bo[g] = b[e - 1 - u[g]]; to[g] = t[e - 1 - u[g]]; uo[g] = u[e - 1 - u[g]]; r[g] = k - 1; } } a = an, b = bn; for (int o = 0; o < e; ++o) { ll d = 0; for (size_t i = 0; i < a[o].size() && i <= k; ++i) { d = (d + a[o][i] * s[k - i]) % mod; } if (d == 0) { t[o] = 1, u[o] = e; } else { for (u[o] = 0, t[o] = d; t[o] % p == 0; t[o] /= p, ++u[o]); int g = e - 1 - u[o]; if (L(a[g], b[g]) == 0) { extand(bn[o], k + 1); bn[o][k] = (bn[o][k] + d) % mod; } else { ll coef = t[o] * inverse(to[g], mod) % mod * pw[u[o] - uo[g]] % mod; int m = k - r[g]; extand(an[o], ao[g].size() + m); extand(bn[o], bo[g].size() + m); for (size_t i = 0; i < ao[g].size(); ++i) { an[o][i + m] -= coef * ao[g][i] % mod; if (an[o][i + m] < 0) an[o][i + m] += mod; } while (an[o].size() && an[o].back() == 0) an[o].pop_back(); for (size_t i = 0; i < bo[g].size(); ++i) { bn[o][i + m] -= coef * bo[g][i] % mod; if (bn[o][i + m] < 0) bn[o][i + m] -= mod; } while (bn[o].size() && bn[o].back() == 0) bn[o].pop_back(); } } } } return make_pair(an[0], bn[0]); }; vector<tuple<ll, ll, int>> fac; for (ll i = 2; i * i <= mod; ++i) if (mod % i == 0) { ll cnt = 0, pw = 1; while (mod % i == 0) mod /= i, ++cnt, pw *= i; fac.emplace_back(pw, i, cnt); } if (mod > 1) fac.emplace_back(mod, mod, 1); vector<vec> as; size_t n = 0; for (auto && x : fac) { ll mod, p, e; vec a, b; tie(mod, p, e) = x; auto ss = s; for (auto && x : ss) x %= mod; tie(a, b) = prime_power(ss, mod, p, e); as.emplace_back(a); n = max(n, a.size()); } vec a(n), c(as.size()), m(as.size()); for (size_t i = 0; i < n; ++i) { for (size_t j = 0; j < as.size(); ++j) { m[j] = get<0>(fac[j]); c[j] = i < as[j].size() ? as[j][i] : 0; } a[i] = crt(c, m); } return a; } LinearSequence(const vec &s, const vec &c, ll mod) : init(s), trans(c), mod(mod), m(s.size()) {} LinearSequence(const vec &s, ll mod, bool is_prime = true) : mod(mod) { vec A; if (is_prime) A = BerlekampMassey(s, mod); else A = ReedsSloane(s, mod); if (A.empty()) A = {0}; m = A.size() - 1; trans.resize(m); for (int i = 0; i < m; ++i) { trans[i] = (mod - A[i + 1]) % mod; } reverse(trans.begin(), trans.end()); init = {s.begin(), s.begin() + m}; } ll calc(ll n) { if (mod == 1) return 0; if (n < m) return init[n]; vec v(m), u(m << 1); int msk = !!n; for (ll m = n; m > 1; m >>= 1) msk <<= 1; v[0] = 1 % mod; for (int x = 0; msk; msk >>= 1, x <<= 1) { fill_n(u.begin(), m * 2, 0); x |= !!(n & msk); if (x < m) u[x] = 1 % mod; else {// can be optimized by fft/ntt for (int i = 0; i < m; ++i) { for (int j = 0, t = i + (x & 1); j < m; ++j, ++t) { u[t] = (u[t] + v[i] * v[j]) % mod; } } for (int i = m * 2 - 1; i >= m; --i) { for (int j = 0, t = i - m; j < m; ++j, ++t) { u[t] = (u[t] + trans[j] * u[i]) % mod; } } } v = {u.begin(), u.begin() + m}; } ll ret = 0; for (int i = 0; i < m; ++i) { ret = (ret + v[i] * init[i]) % mod; } return ret; } vec init, trans; ll mod; int m; }; 8. 数值积分
8.1 复化 Simpson 公式
#define rep(i,a,b) for(ll i=a;i<=b;i++) //复化Simpson公式 O(n) template<typename type,typename fun> type Simpson(fun f,type a,type b,int n){ type h=(b-a)/n,s=0; rep(i,0,n-1)s+=h/6.0*(f(a+h*i)+4.0*f(a+h*(i+0.5))+f(a+h*(i+1.0))); return s; } 8.2 Romberg 算法
#define rep(i,a,b) for(ll i=a;i<=b;i++) //Romberg算法 O(2^n) template<typename type,typename fun> type Romberg(fun f,type a,type b,int n){ assert(n>=3); type h=b-a,T[n+2],S[n+2],C[n+2],R[n+2]; T[1]=h/2.0*(f(a)+f(b)); rep(t,0,n-1){ type s=0; rep(i,0,(1<<t)-1)s+=f(a+h*(i+0.5)); T[t+2]=T[t+1]/2.0+h*s/2.0; h/=2.0; } rep(i,1,n)S[i]=T[i+1]*4.0/3.0-T[i]/3.0; rep(i,1,n-1)C[i]=S[i+1]*16.0/15.0-S[i]/15.0; rep(i,1,n-2)R[i]=C[i+1]*64.0/63.0-C[i]/63.0; return R[n-2]; } 1
相关话题
为何常用偶数进制却少见奇数进制？参加今年（2020）的丘成桐大学生数学竞赛是怎样一种体验？
读美本还是清北？优劣势是什么？
阿贝尔变换强大在哪里？
为什么俄罗斯人在通常被认为考验智力的领域大有作为？
证明了黎曼猜想就能马上得到素数公式吗？
黎曼 ζ 函数为什么要那么解析延拓？
这道函数问题怎么解决？
为什么俄罗斯人在通常被认为考验智力的领域大有作为？
在五大基础学科竞赛中你有哪些笑话和梗？
下一个讨论
如果使 1÷0 有意义，那么应该等于多少？
相关的话题
数学史上有哪些比较著名的猜想因为有反例的存在而没有成为定理？我想证明自然数有穷可行吗？
阿贝尔变换强大在哪里？
如何用初等数论知识证明26是唯一夹在一个平方数和立方数间的正整数？
有没有什么数字的某个幂次方等于0？
如果有一天上帝给了数学家素数的通项公式，这会对数学界有什么影响？
关于p进数域？
该n元不等式如何解？（张端阳的题）?
二分查找有几种写法？它们的区别是什么？
蝴蝶定理有多少种证法？
有没有学习信息竞赛的网课或一对一？
如何严格证明下面这个不等式？
为什么同样走了高考外的路，竞赛生被尊敬而艺术生不受待见?
这个结论是对的吗？能否初等证明？
有没有一个数可以既是完美数又是完全平方数？
中国的数物信奥赛一等奖得主总之来自什么样的社会阶层？
法国的数学水平那么强，为什么在 IMO 上的成绩却很一般？
如何看待O(n log n)时间的整数乘法算法?
有人在p-adic数域Qp上研究过类似球堆积这样的几何数论问题吗？
竞赛组合题的成绩可以通过训练得到显著提高吗？
无理数是否真的存在？
【组合数学】这个魔术有什么策略吗？
求证：关于菲尔兹奖得主舒尔茨的这个非常特殊的说法，是否属实？
有没有一个数可以既是完美数又是完全平方数？
为什么数学教材里，学生首先学习的就是算术，却不学习作为基础的集合与逻辑？
能否看看竞赛生们的书架？
如何用组合数学证明 (n²)! 能被 (n!)^(n＋1) 整除？
根号素数的有限组合是否一定是无理数?
哪些数学定理在直觉上是对的，但证明起来很困难？
有没有哪个素数可以以多种方式写成两个正整数的平方和？