算法竞赛如何训练数论这一块？第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页