数学记录
数学记录(无码版)
多项式全家桶
多项式
定义
形如 \(\sum a_nx^n\) 的求和式若为有限项相加，那么称作多项式 记作\(f(x)=\sum_{n=0}^m a_nx^n\)
对于如 \(\sum_{n=0}^{\infty}a_nx^n\) 有一个非负整数次幂乘以一个常数，称作幂级数
运算
加减运算
两个多项式 \(F(x)=\sum_{n\geq0}a_nx^n\) 与 \(G(x)=\sum_{n\geq0}b_nx^n\) 有
\[F(x)\pm G(x)=\sum_{n\geq0}(a_n\pm b_n)x^n
\]
卷积（乘法）
两个多项式 \(F(x)=\sum_{n\geq0}a_nx^n\) 与 \(G(x)=\sum_{n\geq0}b_nx^n\) 有
\[F(x)\times G(x)=\sum_{i\geq0}a_ix^i \sum_{j\geq0}b_jx^j\\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =\sum_{n\geq0}x^n\sum_{i\geq0}a_ib_{n-i}
\]
FFT
单位根
对于 \(n\) 次单位根为 \(\omega_n\) 表示在复平面上以原点为圆心的半径为1的园，并从实部正半轴开始逆时针将园 \(n\) 等分后的原点与第一个 \(n\) 等分点构成的复数。
在弧度制下，该角弧度为 \(\frac{2\pi}{n}\) ，那么
\[\omega_n^k=\cos (k\times \frac{2\pi}{n}) +i\sin (k\times \frac{2\pi}{n})
\]
那么可以推导出：
\[(\omega_n^n)=(\omega_n^0)=1\\
\omega_{2n}^{2k}=\cos(2k\times\frac{2\pi}{2n})+i\sin(2k\times\frac{2\pi}{2n})=\omega_n^k\\
\omega_n^{\frac{n}{2}}=\cos(\frac{n}{2}\times\frac{2\pi}{n})+i\sin(\frac{n}{2}\times\frac{2\pi}{n})=\cos(\pi)+i\sin(\pi)=-1\\
\omega_n^{k+\frac n2}可以理解为在选择\frac{2k\pi}n 的基础上继续旋转\frac{2\pi \frac{n}2}{n}=\pi,即多旋转180^\circ,\\从向量的角度看，可以感性理解为其相乘，则\omega_n^{k+\frac n2}=-\omega_n^{k}
\]
\[(\omega_n^k)^2=(\cos (k\times \frac{2\pi}{n}) +i\sin (k\times \frac{2\pi}{n}))\times(\cos (k\times \frac{2\pi}{n}) +i\sin (k\times \frac{2\pi}{n}))\\
=\cos^2(k\times \frac{2\pi}{n})-\sin^2(k\times \frac{2\pi}{n}) +i2\cos(k\times \frac{2\pi}{n})\sin(k\times \frac{2\pi}{n})\\
\\\because \begin{bmatrix}
\cos \alpha&-\sin\alpha\\
\sin \alpha&\cos\alpha\\
\end{bmatrix}
\begin{bmatrix}
\cos \alpha\\
\sin \alpha
\end{bmatrix}=
\begin{bmatrix}
\cos^2\alpha-\sin^2\alpha\\
2\cos\alpha\sin\alpha\\
\end{bmatrix}\\
\therefore \cos 2\alpha=\cos^2\alpha-\sin^2\alpha\\
\sin 2\alpha=2\cos\alpha\sin\alpha\\
\therefore(\omega_n^k)^2=\cos(2k\times\frac{2\pi}{n})+i\sin(2k\times\frac{2\pi}{n})=\omega_n^{2k}
\]
\[可以感性理解一下(\omega_n^i)(\omega_n^j)=(\omega_n^{i+j})，因为其可以看作向量旋转
\]
DFT
考虑将这 \(n\) 个 \(n\) 次单位根带入一个 \(n\) 次多项式，获得其点值表示法，复杂度肯定是 \(O(n^2)\) ，但是如果对下标进行奇偶性分离：
\[A(x)=\sum_{i=0}^{n-1} a_ix^i=a_0+a_1x^1+a_2x^2+\dots +a_{n-2}x^{n-2}+a_{n-1}x^{n-1}\\
=a_0+a_2x^2+\dots+a_{n-2}x^{n-2}
\\+a_1x^1+a_3x^3+\dots+a_{n-1}x^{n-1}\\
令A_1(x)=a_0+a_2x^1+\dots+a_{n-2}x^{\frac n2-1}\\A_2(x)=a_1+a_3x^1+\dots+a_{n-1}x^{\frac n2-1}\\
那么A(x)=A_1(x^2)+xA_2(x^2)\\
将其中一个单位根 \omega_n^k \ (k<\frac n2) 带入：\\
A(\omega_n^k)=A_1((\omega_n^k)^2)+\omega_n^kA_2((\omega_n^k)^2)\\
=A_1(\omega_n^2k)+\omega_n^kA_2(\omega_n^2k)\\
若将\omega_n^{k+\frac n2}带入：\\
A`(-\omega_n^k)=A_1(\omega_n^2k)-\omega_n^kA_2(\omega_n^2k)\\
\]
我们可以发现，\(A\) 与 \(A`\) 只有一个常数项不同，所以在求 \(A\) 时可以 \(O(1)\) 算出第二个，所以 \(k\) 的取值范围时 \([0,\frac n 2-1]\) ，然后发现 \(A_1,A_2,A\) 的性质相同，所以可以递归转化，复杂度为 \(O(n\log_2 n)\)。
IDFT
现在已知 \((y_0,y_1,y_2,\dots,y_{n-1})\) 是 \((a_0,a_1,a_2,\dots,a_{n-1})\) 的点值表示法。
假设有个向量 \((c_0,c_1,c_2,\dots,c_{n-1})\) ，即多项式 \(B(x)=\sum\limits_{i=0}^{n-1}y_ix^i\) 在 \((\omega_n^0,\omega_n^{-1},\omega_n^{-2},\dots,\omega_n^{-(n-1)})\) 上的值。
那么：
\[\begin{aligned}
c_k&=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i\\
&=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i\\
&=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j(\omega_n^{-k})^i)\\
&=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^i)^j(\omega_n^{-k})^i\\
&=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^j)^i(\omega_n^{-k})^i\\
&=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i\\
&=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\
\end{aligned}
\]
\[设S(x,n)=\sum_{i=0}^{n-1}x^i\\
当x=\omega_n^k时：\\
S(\omega_n^k,n)=1+(\omega_n^k)+(\omega_n^k)^2+\dots+(\omega_n^k)^{n-1}\\
(\omega_n^k)S(\omega_n^k,n)=(\omega_n^k)+(\omega_n^k)^2+(\omega_n^k)^3+\dots+(\omega_n^k)^n\\
两式相减：\\
(\omega_n^k-1)S(\omega_n^k,n)=(\omega_n^k)^n-1\\
S(\omega_n^k,n)=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}\\
\therefore S(\omega_n^k,n)=\frac{1-1}{\omega_n^k-1}\\
但是在k=0时S(\omega_n^k,n)=n\\
所以在原式c_k=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\
当i\neq k时，后一个\sum中值为0\\
当i= k时，后一个\sum中值为n\\
故c_k=na_k\\
a_k=\frac{c_k}{n}
\]
总结
那么两个多项式相乘可以先将其分别转化为点值表示法，然后对于位置相乘得到结果的点值表示法，然后在跑一遍那个将多项式转化为点值的函数，将点值转化为多项式，最后依据 \(IDFT\) 除以即可。
#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=3e6+10;
const int mod=998244353;
const int inf=1e9+10;
const double PI=acos(-1.0);
inline int read(){
int res=0,v=1;
char c=getchar();
while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
return res*v;
}
struct comple{
double a,b;
comple(double a=0,double b=0):
a(a),b(b){;}
comple operator +(const comple &x)const{return comple(a+x.a,b+x.b);}
comple operator -(const comple &x)const{return comple(a-x.a,b-x.b);}
comple operator *(const comple &x)const{return comple(a*x.a-b*x.b,a*x.b+b*x.a);}
};
void Fast(int now,comple *a,int tag){
if(now==1)return;
comple a1[(now>>1)+10],a2[(now>>1)+10];
for(int i=0;i<(now>>1);++i){
a1[i]=a[i<<1];
a2[i]=a[i<<1|1];
}
Fast(now>>1,a1,tag);
Fast(now>>1,a2,tag);
comple W=comple(cos(2.0*PI/now),tag*(sin(2.0*PI/now))),w=comple(1,0);
for(int i=0;i<(now>>1);++i,w=w*W){
comple z=w*a2[i];
a[i]=a1[i]+z;
a[i+(now>>1)]=a1[i]-z;
}
}
comple a[Max],b[Max];
bool Med;
signed main(){
int n=read();
int m=read();
for(int i=0;i<=n;++i)a[i].a=raed();
for(int i=0;i<=m;++i)b[i].a=read();
int pos=1;
while(pos<=n+m)pos<<=1;
Fast(pos,a,1);
Fast(pos,b,1);
for(int i=0;i<=pos;++i)a[i]=a[i]*b[i];
Fast(pos,a,-1);
for(int i=0;i<=n+m;++i)cout << (int)(a[i].a/pos+0.5)<< ' ';
cerr<< "Time: "<<clock()/1000.0 << "s\n";
cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
return 0;
}
实际上对于上述算法还可以进行优化，就是不递归的写法。
因为我们发现一个数的位置是其元下标二进制翻转后的位置，所有可以预处理然后枚举长度合并。
#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=3e6+10;
const int mod=998244353;
const int inf=1e9+10;
const double PI=acos(-1.0);
inline int read(){
int res=0,v=1;
char c=getchar();
while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
return res*v;
}
struct comple {
double a,b;
comple (double a=0,double b=0):
a(a),b(b){;}
comple operator +(const comple &x)const{return comple(a+x.a,b+x.b);}
comple operator -(const comple &x)const{return comple(a-x.a,b-x.b);}
comple operator *(const comple &x)const{return comple(a*x.a-b*x.b,a*x.b+b*x.a);}
};
int rev[Max];
int limit;
void Fast(comple *a,int tag){
for(int i=0;i<limit;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int mid=1;mid<limit;mid<<=1){
comple W=comple(cos(PI/mid),tag*sin(PI/mid));
for(int R=mid<<1,j=0;j<limit;j+=R){
comple w=comple(1,0);
for(int l=0;l<mid;++l,w=w*W){
comple x=a[l+j],y=a[l+j+mid]*w;
a[l+j]=x+y;
a[l+j+mid]=x-y;
}
}
}
}
comple a[Max],b[Max];
bool Med;
signed main(){
int n=read();
int m=read();
for(int i=0;i<=n;++i)a[i].a=read();
for(int i=0;i<=m;++i)b[i].a=read();
limit=1;
int len=0;
while(limit<=n+m)limit<<=1,++len;
for(int i=0;i<=limit;++i)rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
Fast(a,1);
Fast(b,1);
for(int i=0;i<=limit;++i)a[i]=a[i]*b[i];
Fast(a,-1);
for(int i=0;i<=n+m;++i)cout << (int)(a[i].a/limit+0.5) << ' ';
cout << "\n";
cerr<< "Time: "<<clock()/1000.0 << "s\n";
cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
return 0;
}
NTT
原根
将最小的满足 \(a^l\equiv1\pmod p\) 的 \(l\) 称作 \(a\) 在模 \(p\) 意义下的阶。
如果 \(g\) 在模 \(p\) 意义下的阶是 \(\varphi(p)\)，那么称 \(g\) 是 \(p\) 的原根。
求原根可以这样做：如果 \(\gcd(g,p)=1\)，设 \(p_1,p_2 \dots p_k\) 为 \(p\) 的所有质因子，如果 \(g\) 是 \(p\) 的原根就满足对于任意 \(i\in[1,k],g^{\frac{\varphi(p)}{p_i}}\not\equiv1\pmod p\)
可以发现原根有一个非常优秀的性质就是对于在 \(1\) 到 \(\varphi(p)\) 中的任意一个幂次都是不一样的，而且是以 \(\varphi(p)\) 为长度的一个循环节，这与 \(\omega\) 的性质很像。而且当 \(p\) 为质数时，如果 \(\varphi(p)\) 有很多 \(2\) 这个因子就十分优秀，像 \(998244353=2^{23}\times 7\times 17+1\) 它的原根是 \(3\) 。
那么现在我们可以用 \(g^{\frac{\varphi(p)}{n}}\) 代替 \(n\) 次单位根就行。
#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=3e6+10;
const int mod=998244353;
const int inf=1e9+10;
const int g=3;
const int invg=332748118;
template <int mod>
struct modint{
int val;
static int norm(const int &x){return x<0?x+mod:x;}
modint inv()const{
int a=val,b=mod,u=1,v=0,t;
while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
return modint(u);
}
modint():val(0){}
modint(const int &m):val(norm(m)){}
modint(const ll &m):val(norm(m%mod)){}
modint operator -()const{return modint(norm(-val));}
bool operator ==(const modint &x){return val==x.val;}
bool operator !=(const modint &x){return val!=x.val;}
bool operator <=(const modint &x){return val<=x.val;}
bool operator >=(const modint &x){return val>=x.val;}
bool operator >(const modint &x){return val>x.val;}
bool operator <(const modint &x){return val<x.val;}
modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
modint& operator >>=(const modint &x){return val>>=x.val,*this;}
modint& operator <<=(const modint &x){return val<<=x.val,*this;}
modint& operator ^=(const modint &x){return val^=x.val,*this;}
modint operator >>(const modint &x){return modint(*this)>>=x;}
modint operator <<(const modint &x){return modint(*this)<<=x;}
modint& operator /=(const modint &x){return *this*=x.inv();}
modint operator +(const modint &x){return modint(*this)+=x;}
modint operator -(const modint &x){return modint(*this)-=x;}
modint operator *(const modint &x){return modint(*this)*=x;}
modint operator /(const modint &x){return modint(*this)/=x;}
modint operator ^(const modint &x){return modint(*this)^=x;}
friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};
typedef modint<998244353>m98;
inline int read(){
int res=0,v=1;
char c=getchar();
while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
return res*v;
}
m98 ksm(m98 a,int b){
m98 ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a;
a=a*a;
}
return ans;
}
m98 a[Max],b[Max];
int rev[Max],limit;
void Fast(m98 *a,int tag){
for(int i=0;i<limit;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int mid=1;mid<limit;mid<<=1){
m98 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
for(int R=mid<<1,j=0;j<limit;j+=R){
m98 w=1;
for(int p=0;p<mid;++p,w=w*W){
m98 x=a[p+j],y=w*a[p+j+mid];
a[p+j]=x+y;
a[p+j+mid]=x-y;
}
}
}
if(tag==-1){
m98 inv=m98(limit).inv();
for(int i=0;i<limit;++i)a[i]=a[i]*inv;
}
}
bool Med;
signed main(){
int n=read();
int m=raed();
for(int i=0;i<=n;++i)a[i].val=read();
for(int i=0;i<=m;++i)b[i].val=raed();
limit=1;
int len=0;
while(limit<=n+m)limit<<=1,++len;
for(int i=0;i<=limit;++i)rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
Fast(a,1);
Fast(b,1);
for(int i=0;i<=limit;++i)a[i]=a[i]*b[i];
Fast(a,-1);
for(int i=0;i<=n+m;++i)cout << a[i] << ' ';
cout << "\n";
cerr<< "Time: "<<clock()/1000.0 << "s\n";
cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
return 0;
}
/*
*/
多项式牛顿迭代
给定多项式 \(g(x)\) 求在模 \(x^n\) 下的 \(f(x)\) 使其满足：
\[g(f(x))\equiv0\pmod {x^n}
\]
通过倍增和泰勒展开(不会)得到：
\[f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}\pmod{x^n}
\]
多项式求逆
设给定的多项式是 \(h(x)\)，那么有：
\[g(f(x))=\frac{1}{f(x)}-h(x)\equiv 0\pmod{x^n}
\]
于是有：
\[f(x)\equiv f_0(x)-\frac{\frac{1}{f_0(x)}-h(x)}{-\frac{1}{f_0^2(x)}}\pmod{x^n}\\
\equiv 2f_0(x)-f_0^2(x)h(x)\pmod{x^n}\\
\]
而当 \(n=1\) 时，\([x^0]f_0(x)=[x^0]h(x)^{-1}\)
多项式开根
设给定的多项式是 \(h(x)\)，那么有：
\[g(f(x))\equiv f^2(x)-h(x)\pmod{x^n}\\
\]
于是有：
\[f(x)\equiv f_0(x)-\frac{f^2_0(x)-h(x)}{2f_0(x)}\pmod{x^n}\\
\equiv\frac{f_0^2(x)+h(x)}{2f_0(x)}
\]
多项式对数函数
假设我们要求 \(G(x)\equiv\ln(F(x))\pmod{x^n}\) ，设函数 \(f(x)=\ln(x)\)，那么我们要求的是：
\[G(x)\equiv f(F(x))\pmod{x^n}\\
G'(x)=f(F(x))'\\
G'(x)=f'(F(x))F'(x)\\
\]
注意到 \(\ln(x)\) 的导数是 \(\frac{1}{x}\)，所以有：
\[G'\equiv \frac{F'(x)}{F(x)} \pmod{x^n}
\]
因为我们只关心余数，所以可以求逆就行。
多项式指数函数
设给定的多项式是 \(h(x)\)，那么有：
\[g(f(x))\equiv\ln(f(x))-h(x)\pmod{x^n}
\]
于是有：
\[f(x)\equiv f_0(x)-\frac{\ln(f_0(x))-h(x)}{\frac{1}{f_0(x)}}\pmod{x^n}\\
\equiv(f_0(x))(1-\ln(f_0(x))+h(x))\pmod{x^n}
\]
多项式快速幂
设给定的多项式是 \(h(x)\)，那么求：
\[G(x)\equiv h(x)^k \pmod{x^n}\\
\ln(G(x))\equiv k\ln(h(x))\pmod{x^n}\\
G(x)\equiv e^{k\ln(h(x))}\pmod{x^n}
\]
Code
#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=1e6+10;
const int mod=998244353;
const int inf=1e9+10;
template <int mod>
struct modint{
int val;
static int norm(const int &x){return x<0?x+mod:x;}
modint inv()const{
int a=val,b=mod,u=1,v=0,t;
while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
return modint(u);
}
modint():val(0){}
modint(const int &m):val(norm(m)){}
modint(const ll &m):val(norm(m%mod)){}
modint operator -()const{return modint(norm(-val));}
bool operator ==(const modint &x){return val==x.val;}
bool operator !=(const modint &x){return val!=x.val;}
bool operator <=(const modint &x){return val<=x.val;}
bool operator >=(const modint &x){return val>=x.val;}
bool operator >(const modint &x){return val>x.val;}
bool operator <(const modint &x){return val<x.val;}
modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
modint& operator >>=(const modint &x){return val>>=x.val,*this;}
modint& operator ^=(const modint &x){return val^=x.val,*this;}
modint operator >>(const modint &x){return modint(*this)>>=x;}
modint operator <<(const modint &x){return modint(*this)<<=x;}
modint& operator /=(const modint &x){return *this*=x.inv();}
modint operator +(const modint &x){return modint(*this)+=x;}
modint operator -(const modint &x){return modint(*this)-=x;}
modint operator *(const modint &x){return modint(*this)*=x;}
modint operator /(const modint &x){return modint(*this)/=x;}
modint operator ^(const modint &x){return modint(*this)^=x;}
friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};
typedef modint<998244353>m98;
inline int read(){
int res=0,v=1;
char c=getchar();
while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
return res*v;
}
m98 ksm(m98 a,int b){
m98 ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a;
a=a*a;
}
return ans;
}
m98 inv[Max],frac[Max];
void init(int len) {
frac[0]=inv[0]=1;
for(int i=1;i<=len;++i)frac[i]=frac[i-1]*i;
inv[len]=frac[len].inv();for(int i=len-1;i;--i)inv[i]=inv[i+1]*(i+1);
for(int i=1;i<=len;++i)inv[i]*=frac[i-1];
}
int rev[Max];
void init(int lim,int len){for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));}
struct Poly{ //导和积都到n，否则到n-1;
const int g=3;
const int invg=332748118;
m98 a[Max];
void out(int n,int x=0){for(int i=0;i<n;++i)cout << a[i].val << " ";if(!x)cout << "\n";}
void inte(int n){for(int i=n;i>=0;--i)a[i]=a[i-1]*inv[i];a[0]=0;}
void ReadIn(int n){for(int i=0;i<n;++i)a[i].val=read();}
void der(int n){for(int i=0;i<=n;++i)a[i]=a[i+1]*(i+1);}
m98& operator [](const int x){return a[x];}
void NTT(int tag,int lim){
for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int mid=1;mid<lim;mid<<=1){
m98 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
for(int R=mid<<1,j=0;j<lim;j+=R){
m98 w=1;
for(int p=0;p<mid;++p,w=w*W){
m98 x=a[j+p],y=a[j+p+mid]*w;
a[j+p]=x+y;
a[j+p+mid]=x-y;
}
}
}
if(tag==-1){
m98 inv=m98(lim).inv();
for(int i=0;i<lim;++i)a[i]*=inv;
}
}
}a;
void GetInv(Poly &y,Poly &x,int n){ //第二个为参数，第一个为返回值
static Poly b[2],c;
int lim=2,len=1,bas=1;
int pos=0;
b[pos][0]=x[0].inv();
while(bas<=(n<<1)){
init(lim,len);
pos^=1;
for(int i=0;i<bas;++i)b[pos][i]=b[pos^1][i]<<1;
for(int i=0;i<=lim;++i)c[i]=x[i];
for(int i=lim>>1;i<=lim;++i)c[i]=0;
b[pos^1].NTT(1,lim);c.NTT(1,lim);
for(int i=0;i<=lim;++i)b[pos^1][i]=b[pos^1][i]*b[pos^1][i]*c[i];
b[pos^1].NTT(-1,lim);
for(int i=0;i<bas;++i)b[pos][i]-=b[pos^1][i],b[pos^1][i]=0;
bas<<=1,lim<<=1,++len;
}
for(int i=0;i<n;++i)y[i]=b[pos][i];
for(int i=0;i<bas;++i)b[0][i]=b[1][i]=c[i]=0;
}
void GetSqrt(Poly &y,Poly &x,int n){
static Poly t1,t2,s;
s[0]=1;
int lim=2,len=1,bas=1;
while(bas<=(n<<1)){
for(int i=0;i<bas;++i)t1[i]=x[i];
for(int i=0;i<bas;++i)t2[i]=s[i]<<1;
GetInv(t2,t2,bas);init(lim,len);
t2.NTT(1,lim);t1.NTT(1,lim);
for(int i=0;i<=lim;++i)t1[i]=t1[i]*t2[i];
t1.NTT(-1,lim);
for(int i=0;i<bas;++i)s[i]=s[i]/2+t1[i];
bas<<=1,lim<<=1,++len;
}
for(int i=0;i<n;++i)y[i]=s[i];
for(int i=0;i<bas;++i)t1[i]=t2[i]=s[i]=0;
}
void GetLn(Poly &y,Poly&x,int n){
static Poly t1,t2;
for(int i=0;i<n;++i)t1[i]=t2[i]=x[i];
t1.der(n);GetInv(t2,t2,n);
int lim=1,len=0;
while(lim<=n<<1)lim<<=1,++len;
init(lim,len);
t1.NTT(1,lim),t2.NTT(1,lim) ;
for(int i=0;i<lim;++i)t1[i]*=t2[i];
t1.NTT(-1,lim);t1.inte(n) ;
for(int i=0;i<n;++i)y[i]=t1[i];
for(int i=0;i<lim;++i)t1[i]=t2[i]=0;
}
void GetExp(Poly &y,Poly &x,int n){
static Poly s,tmp;s[0]=1;
int bas=1,lim=2,len=1;
while(bas<=n<<1){
for(int i=0;i<bas;++i)tmp[i]=s[i];
GetLn(tmp,tmp,bas);
for(int i=0;i<bas;++i)tmp[i]=x[i]-tmp[i];
tmp[0]+=1;init(lim,len);tmp.NTT(1,lim),s.NTT(1,lim) ;
for(int i=0;i<lim;++i)s[i]=s[i]*tmp[i];
s.NTT(-1,lim) ;
for(int i=bas;i<lim;++i)s[i]=0;
lim<<=1,bas<<=1,++len;
}
for(int i=0;i<n;++i)y[i]=s[i];
for(itn i=0;i<bas;++i)s[i]=tmp[i]=0;
}
void expow(Poly &y,Poly &x,int n,m98 k){
static Poly s;
for(itn i=0;i<n;++i)s[i]=x[i];GetLn(s,s,n);
for(int i=0;i<n;++i)s[i]*=k;GetExp(s,s,n);
for(int i=0;i<n;++i)y[i]=s[i];
for(int i=0;i<n;++i)s[i]=0;
}
Poly b;
m98 k;
void In(){
k=0;
char c=getchar();
while(c<'0'||c>'9')c=getchar();
while(c>='0'&&c<='9'){k=(k<<3)+(k<<1)+(c^48),c=getchar();}
}
bool Med;
signed main(){
init(200000);
cerr<< "Time: "<<clock()/1000.0 << "s\n";
cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
return 0;
}
/*
*/
第二类斯特林数
定义
第二类斯特林数 \(\begin{Bmatrix}n\\m\end{Bmatrix}\) 也记作 \(S(n,m)\)，表示将 \(n\) 个区分的小球放入 \(m\) 个相同的盒子且不允许有空盒子的方案数。
递推公式
\[\begin{Bmatrix}n\\m\end{Bmatrix}=\begin{Bmatrix}n-1\\m-1\end{Bmatrix}+m\begin{Bmatrix}n-1\\m\end{Bmatrix}
\]
通过组合意义来解释：
当加入一个球时可以将其单独放入一个空盒子中，方案数为 \(\begin{Bmatrix}n-1\\m-1\end{Bmatrix}\) ，
也可以将其放入一个非空的盒子中，方案数为 \(m\begin{Bmatrix}n-1\\m\end{Bmatrix}\)。
通过加法原理得到递推公式。
通向公式
\[\begin{Bmatrix}n\\m\end{Bmatrix}=\sum_{i=0}^m\frac{(-1)^{m-i}i^n}{i!(m-i)!}
\]
可以使用容斥证明该式子：
- 将 \(n\) 个区分的球放入 \(m\) 个区分的盒子中且允许有空盒子的方案为 \(A_m\)，
- 将 \(n\) 个区分的球放入 \(m\) 个区分的盒子但不能有空盒子的方案是 \(B_m\)。
那么有 :
\[A_m=m^n=\sum_{i=0}^m\binom{m}{i}B_i
\]
通过二项式反演可得：
\[\begin{aligned}
B_m&=\sum_{i=0}^m(-1)^{m-i}\binom{m}{i}A_i\\&=\sum_{i=0}^m (-1)^{m-i}\frac{m!i^n}{i!(m-i)!}
\end{aligned}
\]
因为 \(B_m\) 正好 \(\begin{Bmatrix}n\\m\end{Bmatrix}\) 的 \(i!\) 倍，所以得到通向公式：
\[\begin{Bmatrix}n\\m\end{Bmatrix}=\sum_{i=0}^m\frac{(-1)^{m-i}i^n}{i!(m-i)!}
\]
第二类斯特林数·行
发现其通向公式
\[\begin{Bmatrix}n\\i\end{Bmatrix}=\sum\limits_{j=0}^i\frac{(-1)^{i-j}j^n}{j!(i-j)!}\\
=\sum_{j=0}^i\frac{j^n}{j!}.\frac{(-1)^{i-j}}{(i-j)!}
\]
这是一个卷积的形式，所以可以 NTT 求。
#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=1e6+10;
const int mod=167772161;
const int inf=1e9+10;
template <int mod>
struct modint{
int val;
static int norm(const int &x){return x<0?x+mod:x;}
modint inv()const{
int a=val,b=mod,u=1,v=0,t;
while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
return modint(u);
}
modint():val(0){}
modint(const int &m):val(norm(m)){}
modint(const ll &m):val(norm(m%mod)){}
modint operator -()const{return modint(norm(-val));}
bool operator ==(const modint &x){return val==x.val;}
bool operator !=(const modint &x){return val!=x.val;}
bool operator <=(const modint &x){return val<=x.val;}
bool operator >=(const modint &x){return val>=x.val;}
bool operator >(const modint &x){return val>x.val;}
bool operator <(const modint &x){return val<x.val;}
modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
modint& operator >>=(const modint &x){return val>>=x.val,*this;}
modint& operator ^=(const modint &x){return val^=x.val,*this;}
modint operator >>(const modint &x){return modint(*this)>>=x;}
modint operator <<(const modint &x){return modint(*this)<<=x;}
modint& operator /=(const modint &x){return *this*=x.inv();}
modint operator +(const modint &x){return modint(*this)+=x;}
modint operator -(const modint &x){return modint(*this)-=x;}
modint operator *(const modint &x){return modint(*this)*=x;}
modint operator /(const modint &x){return modint(*this)/=x;}
modint operator ^(const modint &x){return modint(*this)^=x;}
friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};
typedef modint<167772161>m16;
inline int read(){
int res=0,v=1;
char c=getchar();
while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
return res*v;
}
int rev[Max];
void init(int lim,int len){
for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
}
m16 ksm(m16 a,int b){
m16 ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a;
a=a*a;
}
return ans;
}
struct Poly{
const m16 g=3;
const m16 invg=(mod+1)/3;
m16 a[Max];
m16& operator [](const int x){return a[x];}
void NTT(int tag,int lim){
for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int mid=1;mid<lim;mid<<=1){
m16 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
for(int R=mid<<1,j=0;j<lim;j+=R){
m16 w=1;
for(int p=0;p<mid;++p,w=w*W){
m16 x=a[p+j],y=a[p+j+mid]*w;
a[p+j]=x+y;
a[p+j+mid]=x-y;
}
}
}
if(tag==-1){
m16 inv=m16(lim).inv();
for(int i=0;i<lim;++i)a[i]*=inv;
}
}
};
Poly a,b;
m16 frac[Max],inv[Max];
bool Med;
signed main(){
int n=read();
inv[0]=frac[0]=1;
for(int i=1;i<=n;++i)frac[i]=frac[i-1]*i;
inv[n]=frac[n].inv();
for(int i=n-1;i>=1;--i)inv[i]=inv[i+1]*(i+1);
for(int i=0;i<=n;++i){
a[i]=ksm(i,n)*inv[i];
b[i]=inv[i]*(i%2==0?1:-1);
}
int lim=1,len=0;
while(lim<=n+n)lim<<=1,++len;
init(lim,len);
a.NTT(1,lim),b.NTT(1,lim) ;
for(int i=0;i<lim;++i)a[i]*=b[i];
a.NTT(-1,lim) ;
for(int i=0;i<=n;++i)cotu << a[i] << ' ';
cotu << "\n";
cerr<< "Time: "<<clock()/1000.0 << "s\n";
cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
return 0;
}
第二类斯特林数·列
先假设每个盒子也要区分，考虑对于一个盒子计算，其生成函数是 \(G(x)=\sum\limits_{i\geq 1}\frac{x^i}{i!}=(e^x-1)\)，其第 \(i\) 项表示在这个盒子放入了 \(i\) 个球。 这样对于 \(m\) 个盒子，\(F(x)=G(x)^m\)，那么如果有 \(n\) 个球，答案就是 \([x^n]G(x)^m\) 由于目前的盒子是区分的，所以要乘以一个 \(m!\)，故有 \(n\) 个球的答案是 \(m![x^n]G(x)^m\)。这样可以多项式快速幂求，都是发现多项式 \(G(x)\) 的首项是 \(0\) 所以要多项式平移。这里使用的是 EGF 其卷积形式如下 \(\sum\limits_{i\geq 0}x^i\sum\limits_{j=0}^i\frac{1}{i!(i-j)!}a_ib_j=\sum\limits_{i\geq 0} \frac{x^i}{i!}\sum\limits_{j=0}^i\binom{i}{j}a_ib_j\) 就是考虑了 \(a\),\(b\) 合并的顺序，而最后要乘一个 \(i!\) 去消掉 EGF 的 \(i!\) 这个系数。
#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=2e6+10;
const int mod=167772161;
const int inf=1e9+10;
template <int mod>
struct modint{
int val;
static int norm(const int &x){return x<0?x+mod:x;}
modint inv()const{
int a=val,b=mod,u=1,v=0,t;
while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
return modint(u);
}
modint():val(0){}
modint(const int &m):val(norm(m)){}
modint(const ll &m):val(norm(m%mod)){}
modint operator -()const{return modint(norm(-val));}
bool operator ==(const modint &x){return val==x.val;}
bool operator !=(const modint &x){return val!=x.val;}
bool operator <=(const modint &x){return val<=x.val;}
bool operator >=(const modint &x){return val>=x.val;}
bool operator >(const modint &x){return val>x.val;}
bool operator <(const modint &x){return val<x.val;}
modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
modint& operator >>=(const modint &x){return val>>=x.val,*this;}
modint& operator ^=(const modint &x){return val^=x.val,*this;}
modint operator >>(const modint &x){return modint(*this)>>=x;}
modint operator <<(const modint &x){return modint(*this)<<=x;}
modint& operator /=(const modint &x){return *this*=x.inv();}
modint operator +(const modint &x){return modint(*this)+=x;}
modint operator -(const modint &x){return modint(*this)-=x;}
modint operator *(const modint &x){return modint(*this)*=x;}
modint operator /(const modint &x){return modint(*this)/=x;}
modint operator ^(const modint &x){return modint(*this)^=x;}
friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};
typedef modint<167772161>m16;
inline int read(){
int res=0,v=1;
char c=getchar();
while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
return res*v;
}
m16 ksm(m16 a,int b){
m16 ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a;
a=a*a;
}
return ans;
}
m16 inv[Max],frac[Max],Inv[Max];
void init(int len) {
frac[0]=inv[0]=1;
for(int i=1;i<=len;++i)frac[i]=frac[i-1]*i;
Inv[len]=frac[len].inv();for(int i=len-1;i;--i)Inv[i]=Inv[i+1]*(i+1);
for(int i=1;i<=len;++i)inv[i]=Inv[i]*frac[i-1];
}
int rev[Max];
void init(int lim,int len){for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));}
struct Poly{ //导和积都到n，否则到n-1;
const m16 g=3;
const m16 invg=(mod+1)/3;
m16 a[Max];
void out(int n,int x=0){for(int i=0;i<n;++i)cout << a[i].val << " ";if(!x)cout << "\n";}
void inte(int n){for(int i=n;i>=0;--i)a[i]=a[i-1]*inv[i];a[0]=0;}
void ReadIn(int n){for(int i=0;i<n;++i)a[i].val=read();}
void der(int n){for(int i=0;i<=n;++i)a[i]=a[i+1]*(i+1);}
m16& operator [](const int x){return a[x];}
void NTT(int tag,int lim){
for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int mid=1;mid<lim;mid<<=1){
m16 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
for(int R=mid<<1,j=0;j<lim;j+=R){
m16 w=1;
for(int p=0;p<mid;++p,w=w*W){
m16 x=a[j+p],y=a[j+p+mid]*w;
a[j+p]=x+y;
a[j+p+mid]=x-y;
}
}
}
if(tag==-1){
m16 inv=m16(lim).inv();
for(int i=0;i<lim;++i)a[i]*=inv;
}
}
}a;
void GetInv(Poly &y,Poly &x,int n){ //第二个为参数，第一个为返回值
static Poly b[2],c;
int lim=2,len=1,bas=1;
int pos=0;
b[pos][0]=x[0].inv();
while(bas<=(n<<1)){
init(lim,len);
pos^=1;
for(int i=0;i<bas;++i)b[pos][i]=b[pos^1][i]<<1;
for(int i=0;i<=lim;++i)c[i]=x[i];
for(int i=lim>>1;i<=lim;++i)c[i]=0;
b[pos^1].NTT(1,lim);c.NTT(1,lim);
for(int i=0;i<=lim;++i)b[pos^1][i]=b[pos^1][i]*b[pos^1][i]*c[i];
b[pos^1].NTT(-1,lim);
for(int i=0;i<bas;++i)b[pos][i]-=b[pos^1][i],b[pos^1][i]=0;
bas<<=1,lim<<=1,++len;
}
for(int i=0;i<n;++i)y[i]=b[pos][i];
for(int i=0;i<bas;++i)b[0][i]=b[1][i]=c[i]=0;
}
void GetLn(Poly &y,Poly&x,int n){
static Poly t1,t2;
for(int i=0;i<n;++i)t1[i]=t2[i]=x[i];
t1.der(n);GetInv(t2,t2,n);
int lim=1,len=0;
while(lim<=n<<1)lim<<=1,++len;
init(lim,len);
t1.NTT(1,lim),t2.NTT(1,lim) ;
for(int i=0;i<lim;++i)t1[i]*=t2[i];
t1.NTT(-1,lim);t1.inte(n) ;
for(int i=0;i<n;++i)y[i]=t1[i];
for(int i=0;i<lim;++i)t1[i]=t2[i]=0;
}
void GetExp(Poly &y,Poly &x,int n){
static Poly s,tmp;s[0]=1;
int bas=1,lim=2,len=1;
while(bas<=n<<1){
for(int i=0;i<bas;++i)tmp[i]=s[i];
GetLn(tmp,tmp,bas);
for(int i=0;i<bas;++i)tmp[i]=x[i]-tmp[i];
tmp[0]+=1;init(lim,len);tmp.NTT(1,lim),s.NTT(1,lim) ;
for(int i=0;i<lim;++i)s[i]=s[i]*tmp[i];
s.NTT(-1,lim) ;
for(int i=bas;i<lim;++i)s[i]=0;
lim<<=1,bas<<=1,++len;
}
for(int i=0;i<n;++i)y[i]=s[i];
for(itn i=0;i<bas;++i)s[i]=tmp[i]=0;
}
void expow(Poly &y,Poly &x,int n,int k){
static Poly s;
for(itn i=0;i<=n;++i)s[i]=x[i];GetLn(s,s,n);
for(int i=0;i<=n;++i)s[i]*=k;GetExp(s,s,n);
for(int i=0;i<=n;++i)y[i]=s[i];
for(int i=0;i<=n;++i)s[i]=0;
}
Poly b;
bool Med;
signed main(){
int n=read()+1;
int k=read();
init(n<<3);
for(int i=1;i<n;++i)b[i-1]=Inv[i];
expow(b,b,n,k);
for(int i=0;i<n;++i){
if(i<k)cout << 0 << " ";
else cout << (b[i-k]*Inv[k]*frac[i])<< ' ';
}
cerr<< "Time: "<<clock()/1000.0 << "s\n";
cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
return 0;
}
/*
P5401 [CTS2019] 珍珠
设 \(cnt_i\) 表示颜色为 \(i\) 的珍珠个数。那么对于一个合法序列满足：
\[\begin{aligned}
\sum_{i=1}^{D} \lfloor\frac{cnt_i}{2}\rfloor&\geq m\\
\sum_{i=1}^D\frac{cnt_i-(cnt_i\mod 2)}{2}&\geq m\\
\sum_{i=1}^Dcnt_i-(cnt_i\mod2)&\geq2m\\
n-\sum_{i=1}^Dcnt_i\mod2&\geq 2m\\
\sum_{i=1}^D cnt_i \mod2&\leq n-2m
\end{aligned}
\]
然后可以特判两种情况，分别是 \(n<2m\) 和 \(2m\leq n-D\) ，答案分别为 0 和 \(D^n\) 。
对于其他情况，我们发现答案只与每种珍珠出现的奇偶性有关，那么我们设 \(odd_i\) 表示恰好有 \(i\) 个奇数的方案数，那么 \(ans=\sum\limits_{i=0}^{n-2m} odd_i\) 。
看到恰好可以考虑容斥，定义 \(f_i\) 表示钦定 \(i\) 个是奇数，其余的不进行限制的方案数，那么 \(f_i=\sum\limits_{j=i}^{D} \binom{j}{i} odd_j\)
由二项式反演可以得到：
\[\begin{aligned}
odd_i&=\sum_{j=i}^D (-1)^{j-i}\binom{j}{i}f(j)\\
&=\sum_{j=i}^D(-1)^{j-i}\frac{j!}{i!(j-i)!} f(j)\\
&=\frac{1}{i!} \sum_{j=i}^D \frac{(-1)^{j-i}}{(j-i)!}j!f(j)\\
\end{aligned}
\]
通过EGF的卷积定义：
\[\begin{aligned}
(F*G)[k]&=\sum_{i+j=k} \binom{k}{i}F[i]G[j]\\
\frac{(F*G)[k]}{k!}&=\sum_{i+j=k}\frac{F[i]}{i!}\frac{G[j]}{j!}
\end{aligned}
\]
我们可以发现上面式子中：
\[\sum_{j=i}^D \frac{(-1)^{j-i}}{(j-i)!}j!f(j)
\]
是个卷积形式。
考虑先求 \(f(i)\) ，因为奇数位是1的EGF是：
\[\begin{aligned}
\sum\limits_{i=0}\frac{[i\%2==1]}{i!}x^i&=\frac{e^x-e^{-x}}{2}\\
\therefore f_k&=\binom{D}{k}n^k(e^x)^{D-k}\\
&=\binom{D}{k}n![x^n]\frac{1}{2^k}(e^x-e^{-x})^k(e^x)^{D-k}\\
&=\binom{D}{k}\frac{n!}{2^k}[x^n](e^x-e^{-x})^k(e^x)^{D-k}\\
\end{aligned}
\]
由二项式定理 \((x+y)^k=\sum\limits_{i=0}^k\binom{k}{i}x^iy^{k-i}\) 可知:
\[\begin{aligned}
f_k&=\binom Dk \frac{n!}{2^k}[x^n]\sum_{i=0}^k\binom kie^{ix}(-e^x)^{k-i}(e^x)^{D-k}\\
&=\binom Dk \frac{n!}{2^k}\sum_{i=0}^k\binom ki(-1)^{k-i}e^{ix-(k-i)x+(D-k)x}[x^n]\\
&=\binom Dk \frac{n!}{2^k}\sum_{i=0}^k\binom ki(-1)^{k-i}[x^n]e^{(D-2(k-i))x}\\
\because [x^n]e^{ax}&=\frac{a^n}{n!}\\
\therefore f_k&=\binom Dk\frac{n!}{2^k}\sum_{i=0}^k\binom ki(-1)^{k-i}\frac{(D-2(k-j))^n}{n!}\\
&=\frac{D!}{k!(D-k)!}\frac{n!}{2^k}\sum_{i=0}^k\frac{k!}{i!(i-k)!}(-1)^{k-i}\frac{(D-2(k-j))^n}{n!}\\
令 j&=k-i\\
\therefore f_k&=\frac{D!}{k!(D-k)!}\frac{n!}{2^k}\sum_{j=0}^k\frac{k!}{j!(k-j)!}(-1)^j\frac{(D-2j)^n}{n!}\\
&=\frac{D!}{(D-k)!2^k}\sum_{i=0}^k\frac{(-1)^i(D-2i)^n}{i!}\frac{1}{(k-i)!}
\end{aligned}
\]
显然后面求和号中是一个EGF的卷积形式，于是我们可以求出 \(f_k\) 。
#include <bits/stdc++.h>
#define int long long
#define reaD read
#define raed read
#define haed head
#define cotu cout
using namespace std;
const int Max=3e6+10;
const int mod=998244353;
const int g=3;
const int gi=332748118;
inline int read(){
int res=0,v=1;
char c=getchar();
while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
return res*v;
}
int a[Max],b[Max];
namespace Ntt{
int ksm(int a,int b){
int ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a%mod;
a=a*a%mod;
}
return ans;
}
int p[Max],lim,pos;
void rever(int lim,int pos){
for(int i=1;i<lim;++i){
p[i]=(p[i>>1]>>1)|((i&1)<<(pos-1));
}
}
void NTT(int *a,int tag){
rever(lim,pos);
for(int i=0;i<lim;++i){
if(i<p[i])swap(a[i],a[p[i]]);
}
for(int mid=1;mid<lim;mid<<=1){
int W=ksm(tag==1?g:gi,(mod-1)/(mid<<1))%mod;
for(int r=mid<<1,j=0;j<lim;j+=r){
int w=1;
for(int k=0;k<mid;++k,w=w*W%mod){
int q=w*a[j+k+mid]%mod;
int l=a[j+k];
a[j+k]=(l+q)%mod;
a[j+k+mid]=(l-q+mod)%mod;
}
}
}
if(tag==-1){
int inv=ksm(lim,mod-2) ;
for(int i=0;i<=lim;++i)a[i]=a[i]*inv%mod;
}
}
void Solve(int *a,int *b){
NTT(a,1);
NTT(b,1);
for(int i=0;i<=lim;++i)a[i]=a[i]*b[i]%mod;
NTT(a,-1);
}
}
using namespace Ntt;
int inv[Max],frac[Max];
void init(int n){
frac[0]=1;
for(int i=1;i<=n;++i){
frac[i]=frac[i-1]*i%mod;
inv[i]=ksm(frac[i],mod-2);
}
inv[0]=1;
}
signed main(){
int d,n,m;
d=read(),n=read(),m=read();
init(200000);
if(n<2*m){
cout << 0 << "\n";
return 0;
}
if(2*m<=n-d){
cout << ksm(d,n) %mod << "\n";
return 0;
}
//处理f
for(int i=0;i<=d;++i){
a[i]=((i%2==0?1:-1)*ksm(d-2*i,n)%mod*inv[i]%mod+mod)%mod;;
}
for(int i=0;i<=d;++i){
b[i]=inv[i];
}
lim=1;
pos=0;
while(lim<=2*d)lim<<=1,++pos;
rever(lim,pos);
Solve(a,b);
for(int i=0;i<=d;++i)a[i]=a[i]*frac[d]%mod*inv[d-i]%mod*ksm(ksm(2,i)%mod,mod-2)%mod*frac[i]%mod;
for(int i=d+1;i<=lim;++i)a[i]=b[i]=0;
//处理odd
for(int i=0;i<=d;++i){
b[d-i]=(((i%2==0)?inv[i]:-inv[i])%mod+mod)%mod;
}
Solve(a,b);
int ans=0;
for(int i=0;i<=n-2*m;++i){
ans=(ans+a[i+d]*inv[i]%mod)%mod;
}
printf("%lld\n",ans);
return 0;
}
P4980 【模板】Polya 定理
\[\begin{aligned}
ans&=\frac1n\sum_{i=1}^nn^{\gcd(n,i)}\\
&=\frac1n\sum_{d|n}n^d\sum_{i=1}^{n}[\gcd(n,i)==d]\\
&=\frac1n\sum_{d|n}n^d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}[\gcd(\frac{n}{d},i)==1]\\
&=\frac1n\sum_{d|n}n^d\varphi(\frac{n}{d})
\end{aligned}
\]
#include <bits/stdc++.h>
#define int long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
const int Max=1e5+10;
const int mod=1e9+7;
const int inf=1e9+10;
inline int read(){
int res=0,v=1;
char c=getchar();
while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
return res*v;
}
int ksm(int a,int b){
int ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a%mod;
a=a*a%mod;
}
return ans;
}
int phi(int x){
int ans=x;
for(int i=2;i*i<=x;++i){
if(x%i==0){
ans-=ans/i;
while(x%i==0)x/=i;
}
}
if(x>1)ans-=ans/x;
return ans;
}
int GetAns(int n,int d){
return ksm(n,d)*phi(n/d)%mod;
}
int Get(){
int n=read();
int ans=0;
for(int i=1;i*i<=n;++i){
if(n%i==0){
ans=(ans+GetAns(n,i))%mod;
if(i*i!=n)ans=(ans+GetAns(n,n/i))%mod;
}
}
return ans*ksm(n,mod-2)%mod;
}
signed main(){
int T=read();
while(T--)cout << Get() << "\n";
return 0;
}
P10322 高洁（Purity）
\[\sum_{i=1}^ni^{v(i)+1}\\
d=\prod_{i=1}p_i^{t_i}\\
设t_{max}为\max(t_1,t_2\dots)
\\
令：f(x)=\sum_{i=1}^ni^{x+1}[v(i)==x]\\
\therefore ans=\sum_{i=1}^{t_{max}}f(i)\\
设d_t=\prod_{i=1}p_{i}^{\lceil\frac{t_i}{k}\rceil}\\
\therefore d|i^k等价于d_k|i\\
并且[v(i)==k] 等价于[d_{k-1}\nmid i]\&\&[d_k|i]\\
\therefore f(k)=
\begin{cases}
\sum\limits_{i=1}^ni[d_{t_{max}}\nmid i]&&&k=0\\
\sum\limits_{i=1}^ni^2[d_1|i]&&&k=1\\
\sum\limits_{i=1}^ni^{k+1}[d_{k-1}\nmid i][d_k|i]&&&k>1
\end{cases}\\
\]
当 \(k=0\) 时：
\[\begin{aligned}
f(k)&=\sum_{i=1}^ni[1-d_{t_{max}}|i]\\
&=\sum_{i=1}^ni-\sum_{i=1}^ni[d_{t_{max}}|i]\\
&=\sum_{i=1}^ni-d_{t_{max}}\sum_{i=1}^{\lfloor\frac{n}{d_{t_{max}}}\rfloor}i
\end{aligned}
\]
当 \(k=1\) 时：
转枚举 \(i\) 变为枚举 \(i\) 的因子
\[\begin{aligned}
f(k)&=\sum_{i=1}^{\lfloor\frac{n}{d_1}\rfloor}(id_1)^2\\
&=d_1^2\sum_{i=1}^{\lfloor\frac{n}{d_1}\rfloor}(i)^2\\
\end{aligned}
\]
当 \(k>1\) 时：
\[\begin{aligned}
f(k)&=\sum_{i=1}^{\lfloor\frac{n}{d_k}\rfloor}(id_k)^{k+1}[\frac{d_{k-1}}{d_k}\nmid i]\\
&=d_k^{k+1}\sum_{i=1}^{\lfloor\frac{n}{d_k}\rfloor}i^{k+1}[1-d_{k-1}|i]\\
&=d_{k}^{k+1}(\sum_{i=1}^{\lfloor\frac{n}{d_k}\rfloor}i^{k+1}-\sum_{i=1}^{\lfloor\frac{n}{d_k}\rfloor}i^{k+1}[\frac{d_{k-1}}{d_k}|i])\\
&=d_{k}^{k+1}(\sum_{i=1}^{\lfloor\frac{n}{d_k}\rfloor}i^{k+1}-\sum_{i=1}^{\lfloor\frac{n}{d_{k-1}}\rfloor}(i\frac{d_{k-1}}{d_k})^{k+1})\\
&=d_k^{k+1}\sum_{i=1}^{\lfloor\frac{n}{d_k}\rfloor}i^{k+1}-d_{k-1}^{k+1}\sum_{i=1}^{\lfloor\frac{n}{d_{k-1}}\rfloor}i^{k+1}
\end{aligned}
\]
设 \(sum(n,k)=\sum\limits_{i=1}^ni^k\)
\[f(k)=
\begin{cases}
sum(n,1)-d_{t_{max}}sum(\lfloor\frac{n}{d_{t_{max}}}\rfloor,1)&&k=0\\
d_1^2sum(\lfloor\frac{n}{d_1}\rfloor,2)&&k=1\\
d_k^{k+1}sum(\lfloor\frac{n}{d_k}\rfloor,k+1)-d_{k-1}^{k+1}sum(\lfloor\frac{n}{d_{k-1}}\rfloor,k+1)&&k>1
\end{cases}
\]
#include <bits/stdc++.h>
#define int __int128
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
const int Max=1e5+10;
const int mod=998244353;
const int inf=1e9+10;
namespace Fread {
const int SIZE=1<<21;
char buf[SIZE],*S,*T;
inline char getchar() {
if(S==T){
T=(S=buf)+fread(buf,1,SIZE,stdin);
if(S==T)return '\n';
}
return *S++;
}
}
namespace Fwrite {
const int SIZE=1<<21;
char buf[SIZE],*S=buf,*T=buf+SIZE;
inline void flush(){
fwrite(buf,1,S-buf,stdout);
S=buf;
}
inline void putchar(char c){
*S++=c;
if(S==T)flush();
}
struct POPOSSIBLE{
~POPOSSIBLE(){flush();}
}ztr;
}
#define getchar Fread :: getchar
#define putchar Fwrite :: putchar
namespace Fastio{
struct Reader{
template<typename T>
Reader& operator >> (T& x) {
char c=getchar();
T f=1;
while(c<'0'||c>'9'){
if(c=='-')f=-1;
c=getchar();
}
x=0;
while(c>='0'&&c<='9'){
x=x*10+(c-'0');
c=getchar();
}
x*=f;
return *this;
}
Reader& operator >> (char& c) {
c=getchar();
while(c==' '||c=='\n')c=getchar();
return *this;
}
Reader& operator >> (char* str) {
int len=0;
char c=getchar();
while(c==' '||c=='\n')c=getchar();
while(c!=' '&&c!='\n'&&c!='\r'){
str[len++]=c;
c=getchar();
}
str[len]='\0';
return *this;
}
Reader(){}
}cin;
struct Writer{
template<typename T>
Writer& operator << (T x) {
if(x==0){putchar('0');return *this;}
if(x<0){putchar('-');x=-x;}
static int sta[45];
int top=0;
while(x){sta[++top]=x%10;x/=10;}
while(top){putchar(sta[top]+'0');--top;}
return *this;
}
Writer& operator << (char c) {
putchar(c);
return *this;
}
Writer& operator << (char* str) {
int cur=0;
while(str[cur])putchar(str[cur++]);
return *this;
}
Writer& operator << (const char* str) {
int cur=0;
while(str[cur])putchar(str[cur++]);
return *this;
}
Writer(){}
}cout;
}
#define cin Fastio :: cin
#define cout Fastio :: cout
#define endl Fastio :: endl
inline int read(){
int x;
cin>> x;
return x;
}
int ksm(int a,int b){
int ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a%mod;
a=a*a%mod;
}
return ans%mod;
}
int P[Max],q[Max],s[Max];
int n,d,t,a[Max];
int sum(int n,int k){
P[0]=1;
q[0]=n%mod;
for(int i=1;i<=k+1;++i){
P[i]=P[i-1]*i%mod;
q[i]=q[i-1]*(n-i)%mod;
}
s[k+2]=1;
for(int i=k+1;i;--i){
s[i]=s[i+1]*(n-i)%mod;
}
int ans=0;
int res=0;
for(int i=1;i<=k+1;++i){
res=(res+ksm(i,k)%mod)%mod;
ans=(ans+res*q[i-1]%mod*s[i+1]%mod*ksm(P[i]*((k-i+1)%2==0?1:-1)*P[k+1-i]%mod,mod-2)%mod)%mod;
}
return (ans%mod+mod)%mod;
}
int f(int k){
if(k==0)return ((sum(n,1)-a[t]*sum(n/a[t],1)%mod)%mod+mod)%mod;
if(k==1)return a[1]*a[1]%mod*sum(n/a[1],2)%mod;
return ((ksm(a[k],k+1)*sum(n/a[k],k+1)%mod-ksm(a[k-1],k+1)*sum(n/a[k-1],k+1)%mod)%mod+mod)%mod;
}
vector<int>p;
vector<int>b;
void work(){
n=read();
d=read();
if(d==1){
cout << sum(n,2) << "\n";
return;
}
t=0;
p.clear();
b.clear();
for(int i=2;i*i<=d;++i){
if(d%i==0){
p.push_back(i); //质因数分解
b.push_back(1);
int z=b.size()-1;
d/=i;
while(d%i==0){
b[z]++;
d/=i;
}
t=max(t,b[z]);
}
}
if(d!=1){
p.push_back(d);
b.push_back(1);
int z=b.size()-1;
t=max(t,b[z]);
}
int z=p.size();
for(int i=1;i<=t;++i){
a[i]=1;
for(int j=0;j<z;++j){
a[i]=a[i]*(ksm(p[j],ceil(1.0*b[j]/i)))%mod; //计算d_k
}
}
int ans=0;
for(int i=0;i<=t;++i){
ans=(ans+f(i)%mod)%mod;
}
cout << ans << "\n";
}
signed main(){
int t=read();
while(t--)work();
return 0;
}
不知道为什么要开 __int128 才过，不知道哪里爆 long long 了。
P6271 [湖北省队互测2014] 一个人的数论
设 \(F(x)=\sum\limits_{i=1}^ni^k\)，\(G(x)=\sum\limits_{i=1}^ni^k[\gcd(i,n)==1]\)
\[\begin{aligned}
\because F(x)&=\sum_{d|x}d^kG(\frac{x}{d})\\
\end{aligned}
\]
由莫反可知：
\[\begin{aligned}
G(x)&=\sum_{d|x}d^k\mu(d)F(\frac{x}{d})\\
\end{aligned}
\]
因为F(x)是应该x+1次多项式，即:
\[\begin{aligned}
\\F(x)&=\sum_{i=0}^xf_ix^i\\
\therefore G(x)&=\sum_{d|x}d^k\mu(d)\sum_{i=0}^{k+1}f_i(\frac{x}{d})^i\\
&=\sum_{i=0}^{k+1}f_ix^i\sum_{d|x}d^{k-i}\mu(d)\\
\end{aligned}
\]
由于 \(\sum_{d|x}d^{k-i}\mu(d)\) 为积性函数
\[\begin{aligned}
\because n&=\prod_{p_i|n}p_i^{a_i}
\\
\therefore G(x)&=\sum_{i=0}^{k+1}f_ix^i\prod_{p_j|x}\sum_{z=0}^{a_i}\mu(p_j^z)p_j^{z(k-i)}\\
\therefore \mu(p^k)&=
\begin{cases}
1&&&k=0\\
-1&&&k=1\\
0&&&k>1
\end{cases}\\
\therefore G(x)&=\sum_{i=0}^{k+1}f_ix^i\prod_{p_j|x}(1-p_j^{k-i})\\
\therefore G(n)&=\sum_{i=0}^{k+1}f_in^i\prod_{p_j|n}(1-p_j^{k-i})\\
\end{aligned}
\]
那么只需要插出 \(F(x)\) 的多项式就行。
#include <bits/stdc++.h>
#define int long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
const int Max=1e5+10;
const int mod=1e9+7;
const int inf=1e9+10;
inline int read(){
int res=0,v=1;
char c=getchar();
while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
return res*v;
}
int ksm(int a,int b){
int ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a%mod;
a=a*a%mod;
}
return ans;
}
int a[Max],f[Max];
int len;
void Mul(int c,int b){
++len;
for(int i=len;i>=1;--i){
a[i]=(a[i]*b+a[i-1]*c)%mod;
}
a[0]=a[0]*b%mod;
}
void mul(int v){
for(int i=0;i<=len;++i){
a[i]=a[i]*v%mod;
}
}
void add(){
for(int i=0;i<=len;++i){
f[i]=(f[i]+a[i])%mod;
a[i]=0;
}
len=0;
a[0]=1;
}
void init(int k){
a[0]=1;
len=0;
int val=0;
for(int i=1;i<=k+2;++i){
val=(val+ksm(i,k))%mod;
int m=ksm(val,mod-2);
for(int j=1;j<=k+2;++j){
if(i==j)continue;
Mul(1,mod-j);
m=m*(i+mod-j)%mod;
}
mul(ksm(m,mod-2));
add();
}
}
int p[Max];
signed main(){
int k=read(),w=read();
int n=1;
for(int i=1;i<=w;++i){
p[i]=read();
n=n*ksm(p[i],read())%mod;
}
init(k);
// for(int i=1;i<=k+2;++i)cout<< f[i] << ' ';
// cout << "\n";
int ans=0;
int in=n;
for(int i=1;i<=k+2;++i){
int pos=in*f[i]%mod;
in=in*n%mod;
for(int j=1;j<=w;++j){
pos=pos*(1+mod-ksm(p[j],k-i-1+mod)%mod)%mod;
}
ans=((ans+pos)%mod+mod)%mod;
}
cout << ans << "\n";
return 0;
}
P10269 实力派
第一个问题：
\[\begin{aligned}
Ans_1&=\sum_{T\subseteq S\&\&|T|=k}[\gcd_{x\in T}x=1]\\
&=\sum_{T\subseteq S\&\&|T|=k}\sum_{i|\gcd_{x\in T}x}\mu(i)\\
&=\sum_{i=1}^n\mu(i)C_{cnt_i}^k\\
\end{aligned}
\]
\(cnt_i\) 表示 \(i\) 的倍数的出现次数，\(n\) 为值域。
第二个问题：
\[考虑容斥，如果一个数作为\gcd那么其他选出来的数必定是其倍数，\\并且不是其倍数的倍数\\
例如：当i作为\gcd时，ans_i=i(C_{cnt_i}^k-C_{cnt_{2i}}^k-C_{cnt_{3i}}^k-C_{cnt_{5i}}^k+C_{cnt_{6i}}^k\dots)\\
可以发现ans_i=i\sum_{j=1}\mu(j)C_{cnt_j}^k\\
\begin{aligned}
\therefore Ans_2&=\sum_{i=1}^nans_i=\sum_{i=1}^ni\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\mu(j)C_{cnt_{ij}}^k\\
&=\sum_{i=1}^nC_{cnt_i}^k\sum_{j|i}\mu(j)\frac{i}{j}\\
&=\sum_{i=1}^nC_{cnt_i}^k\sum_{j|i}\mu(j)id(\frac ij)\\
\because \mu\ * id &=\varphi\\
\therefore Ans_2&=\sum_{i=1}^nC_{cnt_i}^k\varphi(i)\\\
\end{aligned}
\]
现在这样是 \(O(nm)\) 的暴力，肯定过不去，考虑优化。可以发现这两个式子只与 \(cnt_i\) 的值有关，那么我们可以将式子改写为：
\[\begin{aligned}
Ans_1&=\sum_{i=1}^vC_{i}^k\sum_{j=1}^n\mu(j)[cnt_j==i]\\
Ans_2&=\sum_{i=1}^vC_{i}^k\sum_{j=1}^n\varphi(j)[cnt_j==i]\\
v&=\max_{i=1}^n\ cnt_i
\end{aligned}
\]
可以发现两个式子后面一个求和号中的值是不随 \(k\) 变化而变化的，那么我们可以花费 \(O(n)\) 的代价算出来其和为 \(S1_i\),\(S2_i\)，故：
\[Ans_1=\sum_{i=1}^vC_i^kS1_i\\
Ans_2=\sum_{i=1}^vC_i^kS2_i\\
\]
由于 \(cnt_i\) 的种类不会过多，大约是 \(O(\sqrt n)\) 种，可以暴力了。
#include <bits/stdc++.h>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
const int Max=1e6+10;
int mod=998244353;
const int inf=1e9+10;
namespace fast_IO {
#define IOSIZE 1000000
char ibuf[IOSIZE], obuf[IOSIZE], *p1 = ibuf, *p2 = ibuf, *p3 = obuf;
#define getchar() ((p1==p2)and(p2=(p1=ibuf)+fread(ibuf,1,IOSIZE,stdin),p1==p2)?(EOF):(*p1++))
#define putchar(x) ((p3==obuf+IOSIZE)&&(fwrite(obuf,p3-obuf,1,stdout),p3=obuf),*p3++=x)
#define isdigit(ch) (ch>47&&ch<58)
#define isspace(ch) (ch<33)
template<typename T> inline T read() { T s = 0; int w = 1; char ch; while (ch = getchar(), !isdigit(ch) and (ch != EOF)) if (ch == '-') w = -1; if (ch == EOF) return false; while (isdigit(ch)) s = s * 10 + ch - 48, ch = getchar(); return s * w; }
template<typename T> inline bool read(T &s) { s = 0; int w = 1; char ch; while (ch = getchar(), !isdigit(ch) and (ch != EOF)) if (ch == '-') w = -1; if (ch == EOF) return false; while (isdigit(ch)) s = s * 10 + ch - 48, ch = getchar(); return s *= w, true; }
template<typename T> inline void print(T x) { if (x < 0) putchar('-'), x = -x; if (x > 9) print(x / 10); putchar(x % 10 + 48); }
inline bool read(char &s) { while (s = getchar(), isspace(s)); return true; }
inline bool read(char *s) { char ch; while (ch = getchar(), isspace(ch)); if (ch == EOF) return false; while (!isspace(ch)) *s++ = ch, ch = getchar(); *s = '\000'; return true; }
inline void print(char x) { putchar(x); }
inline void print(char *x) { while (*x) putchar(*x++); }
inline void print(const char *x) { for (int i = 0; x[i]; i++) putchar(x[i]); }
inline bool read(std::string& s) { s = ""; char ch; while (ch = getchar(), isspace(ch)); if (ch == EOF) return false; while (!isspace(ch)) s += ch, ch = getchar(); return true; }
inline void print(std::string x) { for (int i = 0, n = x.size(); i < n; i++) putchar(x[i]); }
inline bool read(bool &b) { char ch; while(ch=getchar(), isspace(ch)); b=ch^48; return true; }
inline void print(bool b) { putchar(b+48); }
template<typename T, typename... T1> inline int read(T& a, T1&... other) { return read(a) + read(other...); }
template<typename T, typename... T1> inline void print(T a, T1... other) { print(a), print(other...); }
struct Fast_IO { ~Fast_IO() { fwrite(obuf, p3 - obuf, 1, stdout); } } io;
template<typename T> Fast_IO& operator >> (Fast_IO &io, T &b) { return read(b), io; }
template<typename T> Fast_IO& operator << (Fast_IO &io, T b) { return print(b), io; }
#define cout io
#define cin io
#define endl '\n'
} using namespace fast_IO;
inline int read(){
int x;
cin>> x;
return x;
}
int n,m;
ll ksm(ll a,int b){
ll ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a%mod;
a=a*a%mod;
}
return ans;
}
int pri[Max],phi[Max],flag[Max],pos,mu[Max];
void pre(int n){
mu[1]=phi[1]=1;
for(int i=2;i<=n;++i){
if(!flag[i]){
pri[++pos]=i;
phi[i]=i-1;
mu[i]=-1;
}
for(int j=1;j<=pos&&pri[j]*i<=n;++j){
flag[i*pri[j]]=1;
if(i%pri[j]==0){
phi[i*pri[j]]=phi[i]*pri[j];
break;
}
mu[i*pri[j]]-=mu[i];
phi[i*pri[j]]=phi[i]*(pri[j]-1);
}
}
}
int a[Max];
ll p[Max];
ll inv[Max];
ll C(int n,int m){
return p[n]*inv[m]%mod*inv[n-m]%mod;//ksm(p[m],mod-2)%mod*ksm(p[n-m],mod-2)%mod;
}
ll val[Max];
ll ans1[Max],ans2[Max];
void add(ll &a,ll b){
a=a+b;
a+=(a<0?mod:0);
a-=((a>=mod)?mod:0);
}
signed main(){
n=read();
m=read();
mod=read();
int v=0;
p[0]=1;
for(int i=1;i<=n;++i){
a[i]=read();
v=max(v,a[i]);
p[i]=p[i-1]*i%mod;
val[a[i]]++;
}
pre(max(n,v));
for(int i=1;i<=v;++i){
for(int j=2*i;j<=v;j+=i) {
val[i]+=val[j];
}
}
inv[n]=ksm(p[n],mod-2);
for(int i=n-1;i>=0;--i){
inv[i]=inv[i+1]*(i+1)%mod;
}
for(int i=1;i<=max(n,v);++i)mu[i]=mu[i]+(mu[i]<0?mod:0);
for(int i=1;i<=v;++i){
for(int k=1;k<=val[i];++k){
add(ans2[k],C(val[i],k)*phi[i]%mod);
add(ans1[k],C(val[i],k)*mu[i]%mod);
}
}
for(int i=1;i<=m;++i){
int x=read();
cotu << ans1[x] << ' ' << ans2[x] << "\n";
}
return 0;
}
P7360 「JZOI-1」红包
\[\prod_{i_1=1}^n\prod_{i_2=1}^n\dots\prod_{i_k=i}^nlcm(i_1,i_2\dotsi_k)\mod 998244353\\
考虑对每个质数计算其贡献：\\
a_{p,j}表示在j质因数分解以后的质因子p的个数\\（后文可能会省略p，若省略则默认为p_i）
\\
ans=\prod_{p_i|n}p_i^{\sum\limits_{j_1=1}^n\sum\limits_{j_2=1}^n\dots\sum\limits_{j_k=1}^n\max(a_{p_i,j_1},a_{p_i,j_2},\dots a_{p_i,j_k})}\\
尝试对p_i计算\sum\limits_{j_1=1}^n\sum\limits_{j_2=1}^n\dots\sum\limits_{j_k=1}^n\max(a_{p_i,j_1},a_{p_i,j_2},\dots a_{p_i,j_k})\\
假设:t=\max(a_{j_1},a_{j_2},\dots a_{j_k})\\
那么满足\max(a_{j_1},a_{j_2},\dots a_{j_k})<t时的j的可供取值的个数是：
n-b_t\\等等，好像有锅\\
b_{p，j}表示在[1,n]中拥有只是j个质因子p的数的个数（省略同a\\
那么存在k^{(n-b_t)}种取值使得\max(a_{j_1},a_{j_2}\dots a_{j_k})<t\\
那么是不是可以对于每个p暴力枚举t，做差得到\max(a_{j_1},a_{j_2}\dots a_{j_k})=t的个数。\\
ans=\prod_{p_i|n}p_i^{\sum\limits_{t=1}^{p_i^t\leq n}t\times(k^{n-b_{t+1}}-k^{n-b_i})}
\\
问题关键在于求b_{p,j},这个我们可以花费 O(n\log n)的时间对每个数质因数分解，\\然后求前缀和就行\\
\]
现在问题在于暴力（或前缀和）的复杂度是多少（二者复杂度相同），可以列出计算式：
\[\sum_{p\in prime\&\&p<n}\log_pn
\]
通过打表发现在 \(1e6\) 时这个值大约是 \(75000\)
但是好像复杂度有问题（我也不会算），主要是还要加快速幂，先搁一下
P4128 [SHOI2006] 有色图
因为有 \(Burnside\) 引理：
\[|X/G|=\frac{1}{|G|}\sum_{g\in G}|X^g|\\
\]
\(|X/G|\) 是群 \(G\) 作用在集合 \(X\) 的轨道。
\(X^g\) 是 \(g\) 作用在 \(X\) 上的不变的集合。
那么我们要求的就是 \(|X/G|\)。
由于对顶点重排是一个置换群，每一种置换是有若干轮换构成的，我们可以发现有若干边是等价的，即他们的颜色必须相同，而对于每个等价类都有 \(m\) 种选择，那么如果有 \(k\) 个等价类，那么就会存在 \(m^k\) 种染色方式，使其在经过若干次该操作后会恢复成原状。（注意：一个等价类并不等同于一个轮换，因为是经过若干操作后变化原状，并不是一次操作。）
现在问题的关键在于求解 \(k\)，如果我们将一个置换拆解为 \(len\) 个轮换，第 \(i\) 个轮换长度是 \(b_i\) 。
先考虑在同一个轮换中的贡献：我们可以发现，在同一个轮换的同一个等价类中的边长度是一样的，因为这样才能在若干次操作后恢复原因。由于在一个轮换中我们选择的是点，那么他们之间两两连边也应该取出，那么就构成了一个完全图。而在一个有 \(a\) 个点的完全图中长度不同的边有 \(\lfloor\frac{a}{2}\rfloor\) 中，于是这 \(len\) 个轮换的贡献是 \(\sum\limits_{i=1}^{len}\lfloor\frac{b_i}{2}\rfloor\) 。
这是在同一个轮换中，但是在一个轮换中我们只取出其内部连边，然而两两轮换中的连边还没有分类。那么如何统计两两轮换间的贡献呢？
我们可以发现在只两个轮换 \(i\) 和 \(j\) 时 ，只需要经过 \(\operatorname{lcm}(b_i,b_j)\) 次操作就会恢复原样，那么其等价类长度为 \(\operatorname{lcm}(b_i,b_j)\)，由于有 \(b_ib_j\) 条边，那么会产生 \(\frac{b_ib_j}{\operatorname{lcm}(b_i,b_j)}=\gcd(b_i,b_j)\) 个等价类，那么两两轮换中的贡献是：\(\sum\limits_{i=1}^{len}\sum\limits_{j=1}^{i-1}\gcd(b_i,b_j)\)。
那么可以得到 \(k=\sum\limits_{i=1}^{len}\lfloor\frac{b_i}{2}\rfloor+\sum\limits_{i=1}^{len}\sum\limits_{j=1}^{i-1}\gcd(b_i,b_j)\)，但是新的问题又出现了，就是该群的置换个数即 \(|G|\) 是 \(n!\) 的，我们没法去直接枚举所有置换。但是我们发现答案好像只与 \(b\) 数组有关系，那么我们是不是只用去找的 \(b\) 数组的所有情况然后去重就行了呢？
说做就做：将这 \(n\) 个数随便放置的方案是 \(n!\) 种，但是对于一个长度为 \(b_i\) 的轮换会重复计算 \(b_i!\) 种，那么总共有 \(\frac{n!}{\prod\limits_ib_i!}\) ，这是考虑不同轮换之间的方案，而对于一个轮换内部，其实是一个园排列，那么对于一个长度是 \(b_i\) 的轮换会存在 \((b_i-1)!\) 种不同方案，那么总方案就是 \(\prod\limits_i(b_i-1)!\) 种。有一开始的式子我们发现答案与轮换的顺序没有关系，那么我们还要去除因为顺序所带来的重复，假设 \(b_i\) 出现了 \(cnt_i\) 次，那么重复次数为 \(cnt_i!\)，所以全部重复次数为 \(\prod\limits_icnt_i!\)。
那么我们将最后式子写出来：
\[\begin{aligned}
k&=\sum\limits_{i=1}^{len}\lfloor\frac{b_i}{2}\rfloor+\sum\limits_{i=1}^{len}\sum\limits_{j=1}^{i-1}\gcd(b_i,b_j)\\
|X/G|&=\frac{1}{|G|}\sum_{g\in G}\frac{m^{k}\times n!\times \prod(b_i-1)!}{\prod (b_i!cnt_i!)}\\
\end{aligned}
\]
由于 \(|G|\) 是 \(n!\)，那么整理这个式子得到：
\[\begin{aligned}
k&=\sum\limits_{i=1}^{len}\lfloor\frac{b_i}{2}\rfloor+\sum\limits_{i=1}^{len}\sum\limits_{j=1}^{i-1}\gcd(b_i,b_j)\\
|X/G|&=\sum_{g\in G}\frac{m^{k}}{\prod (b_icnt_i!)}\\
\end{aligned}
\]
我们发现 \(b\) 是 \(n\) 的一种分拆形式，并且 \(n\) 本身不大，那么我们可以暴力枚举，计算贡献。
#include <bits/stdc++.h>
#define int long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
const int Max=1e2+10;
int mod=998244353;
const int inf=1e9+10;
inline int read(){
int res=0,v=1;
char c=getchar();
while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
return res*v;
}
int frac[Max],b[Max];
int n,ans=0,m;
void pre(){
frac[0]=1;
for(int i=1;i<=n;++i){
frac[i]=frac[i-1]*i%mod;
}
}
int ksm(int a,int b){
int ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a%mod;
a=a*a%mod;
}
return ans;
}
int work(int len){
int pos=0;
int u=1;
for(int i=1;i<=len;++i){
pos=pos+(b[i]/2); //一个环
u=u*b[i]%mod;
for(int j=1;j<i;++j){
pos=pos+__gcd(b[i],b[j]);//环交贡献
}
}
for(int i=1;i<=len;){
int j=i+1;
while(b[j]==b[i]&&j<=len)++j; //去重
u=u*frac[j-i]%mod;
i=j;
}
return ksm(u%mod,mod-2)*ksm(m,pos)%mod;
}
//注意pos的值不能取模
void dfs(int now,int last,int sum){ //分拆
if(!sum){
ans=(ans+work(now-1))%mod;
return;
}
for(int i=last;i<=sum;++i){
b[now]=i;
dfs(now+1,i,sum-i);
}
}
signed main(){
n=read();
m=read();
mod=read();
pre();
dfs(1,1,n);
cout << ans << "\n";
return 0;
}
然后对于 P4727 [HNOI2009] 图的同构计数 可以将问题转换为每条边染不染色，即 \(m=2\) 的情况求解。
P4619 [SDOI2018] 旧试题
令 \((i,j)\) 表示 \([\gcd(i,j)==1]\)
\[\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^Cd(ijk)\\
\]
先考虑 \(C=1\) 时：
有结论 \(d(ij)=\sum\limits_{x|i}\sum\limits_{y|i}(x,y)\) 。因为我们可以对于每个质因子单独考虑，对于质因子 \(P\) 它在 \(i\) 中出现了 \(p_i\) 次，在 \(j\) 中出现 \(p_j\) 次，那么它会在 \(ij\) 中出现 \(p_i+p_j\) 次。我们知道 \(d(i)=\prod\limits_{P\in pri\&\&P|i} (p_i+1)\) 因为对于一个 \(P\) 你可以选择的次数是在 \([0,p_i]\) 中，根据乘法原理可知。那么质因子 \(P\) 在 \(d(i,j)\) 中的贡献也是 \(p_i+p_j+1\)。而当 \([\gcd(x,y)==1]\) 时，对于质因子 \(P\) 它要么不出现，要么只在 \(x\) 或 \(y\) 中出现（废话）。这样它的贡献是 \(p_i+p_j+1\) 现在考虑是不是乘法原理，因为对于每个 \(P\) 是相互独立的，假设在 \(x\) 中出现了 \(a_x\) 次，那么是不会影响另外质因子 \(P^{'}\) 。所以可以证明。
那么将这个结论推广到 \(C\neq 1\) 时，那么有 \(d(ijk)=\sum\limits_{x|i}\sum\limits_{y|i}\sum\limits_{z|k}(x,y)(y,z)(x,k)\) 。
那么：
\[\begin{aligned}
ans&=\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{x|i}\sum_{y|i}\sum_{z|k}(x,y)(y,z)(x,k)\\
&=\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C(x,y)(y,x)(x,k)\frac{A}{x}
\end{aligned}
\]
List Generation
考虑答案为从点 \((0,0)\) 到点 \((n,m)\) 的只向上或右走的路径长度之和，考虑枚举转折点，为了不重，将转折点分为两类：
- 从左边来然后向上走
- 从下边来然后向右走
于是考虑选 \(i\) 转折点 1 的方案有 \(\binom ni\binom mi\) ，然后路径上一共有 \(n+m-1\) 个点，由于已经钦定了 \(i\) 个点了，那么剩下 \(s=n+m-i-1\) 个点可选可不选，对于这 \(s\) 个点的贡献是：
\[\begin{aligned}
\sum_{j=0}^s\binom sj(i+j+2)&=2^s(i+2)+\sum_{j=0}^s\binom sj j\\
&=2^s(i+2)+\frac{\sum_{j=0}^s\binom sj s}{2}\\
&=2^s(i+2)+s2^{s-1}\\
\therefore ans&=\sum_{i=0}^{\min(n,m)}\binom ni \binom mi(2^s(i+2)+s2^{s-1})
\end{aligned}
\]
#include <bits/stdc++.h>
#define pii pair<int,int>
#define int long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=1e7+10;
const int mod=1e9+7;
const int inf=1e9+10;
inline int read(){
int res=0,v=1;
char c=getchar();
while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
return res*v;
}
int frc[Max],pr[Max],inv[Max];
int ksm(int a,int b){
int ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a%mod;
a=a*a%mod;
}
return ans;
}
void pre(){
pr[0]=1;
for(int i=1;i<=1e7;++i){
pr[i]=pr[i-1]*2%mod;
}
inv[0]=frc[0]=1;
for(int i=1;i<=5e6+2;++i){
frc[i]=frc[i-1]*i%mod;
}
int z=5e6+2;
inv[z]=ksm(frc[z],mod-2);
for(int i=z-1;i>=1;--i){
inv[i]=inv[i+1]*(i+1)%mod;
}
}
int C(int n,int m){ //n选m
if(n<0||m<0||n<m)return 0;
return frc[n]*inv[m]%mod*inv[n-m]%mod;
}
void work(){
int n=read();
int m=read();
int ans=0;
for(int i=0;i<=min(n,m);++i){
int s=n+m-i-1;
ans=(ans+(C(n,i)*C(m,i)%mod*(pr[s]*(i+2)%mod+(s?s*pr[s-1]%mod:0))%mod))%mod;
}
cout << ans << "\n";
}
bool Med;
signed main(){
pre();
int T=read();
while(T--)work();
cerr<< "Time: "<<clock()/1000.0 << "s\n";
cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
return 0;
}
P6031 CF1278F Cards 加强版
因为每次洗牌是独立的，所以一次抽中王牌的概率 \(p=\frac{1}{m}\)，所以有：
\[ans=\sum_{i=0}^n\binom nip^i(1-p)^{n-i}i^k
\]
用第二类斯特林数拆幂次：
因为 \(m^n\) 等于将 \(n\) 个区分的小球放入 \(m\) 个区分的盒子中，允许有空的方案数，用第二类斯特林数表示就是：
\[m^n=\sum_{i=0}^{\min(n,m)} \binom{m}{i}i!\begin{Bmatrix}n\\i\end{Bmatrix}
\]
代入原式可得：
\[\begin{aligned}
ans&=\sum_{i=0}^n\binom ni p^i(1-p)^{n-i}\sum_{j=0}^{k} \binom ij j!\begin{Bmatrix}k\\j\end{Bmatrix}\\
&=\sum_{j=0}^{k}j!\begin{Bmatrix}k\\j\end{Bmatrix}\sum_{i=0}^n\binom ni \binom ij p^i(1-p)^{n-i}\\
\end{aligned}
\]
考虑第二个求和号中的内容：
\[\begin{aligned}
&=\sum_{i=0}^n\binom ni\binom ijp^i(1-p)^{n-i}\\
&=\sum_{i=0}^n \frac{n!i!}{i!(n-i)!j!(i-j)!}p^i(1-p)^{n-i}\\
&=\frac{1}{j!}\sum_{i=0}^n\frac{n!}{(n-i)!(i-j)!}p^i(1-p)^{n-i}\\
&=\frac{1}{j!}\sum_{i=0}^n\frac{(n-j)!}{(n-i)!(i-j)!}\frac{n!}{(n-j)!}p^i(1-p)^{n-i}\\
&=\frac{n!}{(n-j)!j!}\sum_{i=0}^n\binom{n-j}{n-i}p^i(1-p)^{n-i}\\
&=\binom nj\sum_{i=0}^n\binom{n-j}{i-j}p^i(1-p)^{n-i}\\
&=p^j\binom nj \sum_{i=0}^{n-j}\binom{n-j}{i}p^i(1-p)^{n-i-j}\\
&=p^j\binom{n}{j}\\
\therefore ans&=\sum_{j=0}^kj!\begin{Bmatrix}k\\j\end{Bmatrix}p^j\binom nj
\end{aligned}
\]
\[\because m^n=\sum_{i=0}^m\binom mi i!\begin{Bmatrix}n\\i\end{Bmatrix}
\]
通过二项式反演得到：
\[\begin{Bmatrix}n\\m\end{Bmatrix}=\frac{1}{m!}\sum_{i=0}^m(-1)^{m-i}\binom mi i^n
\]
代回原式得到：
\[\begin{aligned}
ans&=\sum_{j=0}^k\frac{1}{j!}\sum_{i=0}^j(-1)^{i-j}\binom ji i^kj!\binom njp^j\\
&=\sum_{i=0}^k(-1)^ii^k\sum_{j=i}^k\binom nj\binom ji(-p)^j\\
&=\sum_{i=0}^k(-1)^ii^k\sum_{j=i}^k\frac{n!j!}{j!i!(n-j)!(j-i)!}(-p)^j\\
&=\sum_{i=0}^k(-1)^ii^k\binom ni\sum_{j=i}^k \binom {n-i}{j-i}(-p)^j\\
&=\sum_{i=0}^ki^kp^i\binom ni\sum_{j=0}^{k-i}\binom {n-i}{j}(-p)^j\\
\end{aligned}
\]
设 \(S(i)=\sum\limits_{j=0}^{k-i}\binom {n-i}j(-p)^j\)，裂开组合数：
\[\begin{aligned}
S(i)&=\sum_{j=0}^{k-i}\left( \binom{n-i-1}{j}+\binom{n-i-1}{j-1}\right)(-p)^j\\
&=S(i+1)+\binom{n-i-1}{k-i}(-p)^{k-i}+\sum_{j=0}^{k-i}\binom{n-i-1}{j-1}(-p)^j\\
\end{aligned}
\]
令 \(s=S(i+1)+\binom{n-i-1}{k-i}(-p)^{k-i}\),
\[\begin{aligned}
S(i)&=s+(-p)\sum_{j=0}^{k-i}\binom{n-i-1}{j-1}(-p)^{j-1}\\
&=s+(-p)\sum_{j=0}^{k-i-1}\binom{n-i-1}j(-p)^j\\
\therefore S(i)&=c+(-p)S(i+1)=(1-p)S(i+1)+\binom{n-i-1}{k-i}(-1)^{k-i}\\(S(k)&=1)\\
\therefore ans&=\sum_{i=0}^ki^k\binom ni p^i S(i)
\end{aligned}
\]
#include <bits/stdc++.h>
#define pii pair<int,int>
#define int long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=1e7+10;
const int mod=998244353;
const int inf=1e9+10;
inline int read(){
int res=0,v=1;
char c=getchar();
while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
return res*v;
}
int ksm(int a,int b) {
int ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a%mod;
a=a*a%mod;
}
return ans;
}
int pri[Max],idk[Max],cnt;
int inv[Max];
int k;
void pre(int n){
inv[1]=idk[1]=1;
for(int i=2;i<=n;++i){
if(!idk[i]){
pri[++cnt]=i;
idk[i]=ksm(i,k);
inv[i]=ksm(i,mod-2);
}
for(int j=2;j<=cnt&&i*pri[j]<=n;++j){
idk[i*pri[j]]=idk[i]*idk[pri[j]]%mod;
inv[i*pri[j]]=inv[i]*inv[pri[j]]%mod;
if(i%pri[j]==0)break;
}
}
}
int n,m;
int frac[Max],Inv[Max];
void init(){
frac[0]=Inv[0]=1;
for(int i=1;i<=n;++i)frac[i]=frac[i-1]*i%mod;
Inv[n]=ksm(frac[n],mod-2);
for(int i=n-1;i>=1;--i)Inv[i]=Inv[i+1]*(i+1)%mod;
}
int C(int n,int m){
return frac[n]*Inv[m]%mod*Inv[n-m]%mod;
}
void Work1(){ //n<=k
int ans=0;
init();
pre(n);
int p=ksm(m,mod-2);
int Q=((1-p)%mod+mod)%mod;
int Z=ksm(Q,n);
Q=ksm(Q,mod-2);
int z=1;
for(int i=0;i<=n;++i){
ans=(ans+C(n,i)*z%mod*Z%mod*idk[i]%mod)%mod;
z=z*p%mod;
Z=Z*Q%mod;
}
cout << ans << "\n";
}
int S[Max];
void Work2(){
pre(k);
int p=ksm(m,mod-2);
int c=1,x=1;
S[k]=1;
for(int i=k-1;i;--i){
c=c*(n-i-1)%mod*inv[k-i]%mod;
x=((-1ll*x*p)%mod+mod)%mod;
S[i]=(S[i+1]*((1-p)+mod)%mod+c*x)%mod;
}
c=x=1;
int ans=0;
for(int i=1;i<=k;++i){
x=x*p%mod;
c=c*(n-i+1)%mod*inv[i]%mod;
ans=(ans+S[i]*c%mod*idk[i]%mod*x%mod)%mod;
}
cout << ans << "\n";
}
bool Med;
signed main(){
n=read();
m=read();
k=read();
if(n<=k)Work1();
else Work2();
cerr<< "Time: "<<clock()/1000.0 << "s\n";
cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
return 0;
}