算法竞赛如何训练数论这一块？ 第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
相关话题
怎样利用格理论，也就是 minkowski 基本定理来证明拉格朗日四平方和定理以及费马平方和定理?宁波大学的陈计教授是一个什么样的人？
为什么前N个自然数的最小公倍数约等于e^N?
1²+2²+…+n²求和公式的推导有哪些方法?
级数求积：是否有一般的收敛判别法？以及实例∏[p是素数] p/(p-1) 是否收敛?
下面这道数竞平面几何题求好的解题思路和方法？
对于 3 和 4 之间的整数 Bleem，你怎么看？
如何用数学知识解答「在进行社区大规模核酸检测时，分成几人一组进行混检效率最高」？
n! 和 n²，哪个更大呢？
黎曼猜想具体是如何推出素数定理的广义形式的？
下一个讨论
如果使 1÷0 有意义，那么应该等于多少？
相关的话题
所有质数的倒数的平方和的精确值是多少？【组合数学】这个魔术有什么策略吗？
到底有没有素数公式？素数公式的意义有多大？
数学家和物理学家思考方式是什么？
请问乘和乘以有区别吗？
诗歌《冰雹之路》，大家觉得如何？
有哪些看起来很简单但做起来很难的数学题？
有限域上为什么有x的m次方=e的解的个数不超过m？
如何才能在高考前证明哥德巴赫猜想？
无理数是否真的存在？
除了 3，4，5 以外是否还有别的三角形，它的三条边是连续自然数，它的面积也是自然数？
如何证明有理数加法群不是有限生成群？
拉马努金的那些壮观的公式，都是怎么发现的？
想问问数学大神们几个问题，学习数学的兴趣来源是什么？数学在生活中的具体应用？
阿里数学竞赛到底难不难？高考142分能不能去试试看？
对于 3 和 4 之间的整数 Bleem，你怎么看？
竞赛生和普通学生有哪些区别？
如何看待几何数论（geometry of numbers）这一数论分支？
全体自然数的发散级数和等于负十二分之一代表了什么？隐藏了一个天大的秘密吗？
请问下面这道题怎么解决？
圆周率π的这个用正切半角表示的无穷级数展开式怎么证明?
有哪些看起来很简单但做起来很难的数学题？
如果让一个高考数学理综满分总分700+的选手经过两年系统训练做奥林匹克竞赛题会比保送生做的差吗？
哪个高中竞赛性价比最高？
若两个正整数互质，如何证明它们的平方也互质？
这道函数问题怎么解决？
如何判断一个超级大的数是不是素数？
为什么IMO的国家队成员很少有女生？
如何证明一个无理数的整数倍数的小数部分在（0，1）上均匀分布？
如果打算证明黎曼猜想，请问从大一开始应该做什么数学基础准备？