组合和

题目大意

求∑ni=1Cm(i,n)

做法

随手化一下式子变成
∑d|nCmd∗ϕ(nd)
发现很想狄利克雷卷积的形式，不妨尝试凑出另一个积性函数。
1m!∑d|n∑mi=0s(m,i)∗(−1)m−i∗di∗ϕ(nd)
1m!∑mi=0s(m,i)∗(−1)m−i∑d|ndi∗ϕ(nd)
现在就好了！
两个积性函数狄利克雷卷积还是积性函数。
设gi(n)=∑d|ndi∗ϕ(nd)
答案是1m!∑mi=0s(m,i)∗(−1)m−iΠpk|ngi(pk)
考虑计算质数的次幂。
gi(pk)=∑kj=0pij∗ϕ(pk−j)
gi(pk)=(p−1)∑k−1j=0pk−1+(i−1)j+pik
分解质因数用pollard_rho。

#include<cstdio>#include<algorithm>#include<cmath>#define fo(i,a,b) for(i=a;i<=b;i++)using namespace std;typedef long long ll;const int maxn=3000+10,mo=1998585857;ll ss[maxn][maxn];ll p[maxn],c[maxn];ll i,j,k,l,t,n,m,tot,top,ans,num;ll qsm(ll x,ll y){    if (!y) return 1;    ll t=qsm(x,y/2);    t=t*t%mo;    if (y%2) t=t*x%mo;    return t;}ll calcg(ll i,ll p,ll k){    ll j,t=0;    fo(j,0,k-1) (t+=qsm(p%mo,k-1+(i-1)*j))%=mo;    t=t*((p-1)%mo)%mo;    (t+=qsm(p%mo,i*k))%=mo;    return t;}ll random(ll x){    ll t=rand()%10000;    t=t*10000+rand()%10000;    t=t*10000+rand()%10000;    t=t*10000+rand()%10000;    return t%x;}ll mul(ll a,ll b,ll p){    ll tmp=(a*b-(ll)((long double)a/p*b+1e-8)*p);    return tmp<0?tmp+p:tmp;}ll quicksortmi(ll x,ll y,ll mo){    if (!y) return 1;    ll t=quicksortmi(x,y/2,mo);    t=mul(t,t,mo);    if (y%2) t=mul(t,x,mo);    return t;}bool Miller_Rabin(ll n){    if (n==1) return 0;    int s=5,t=0,i;    ll a,p,k=n-1;    while (k%2==0) k/=2,t++;    while (s--){        a=random(n-1)+1;        p=a=quicksortmi(a,k,n);        fo(i,1,t){            a=mul(a,a,n);            if (a==1&&p!=1&&p!=n-1) return 0;            p=a;        }        if (a!=1) return 0;    }    return 1;}ll gcd(ll a,ll b){    return b?gcd(b,a%b):a;}ll Pollard_Rho(ll n){    ll k,x,y,c,d,i=1;    while (1){        c=random(n-1);        k=2;y=x=random(n);        i=1;        while (1){            y=(mul(y,y,n)+c)%n;            d=gcd(abs(x-y),n);            if (i==k) x=y,k<<=1;            i++;            if (d!=1) break;        }        if (d!=n) return d;    }}void work(ll n){    if (Miller_Rabin(n)){        //ans=max(ans,n);        p[++top]=n;        return;    }    ll p=Pollard_Rho(n);    ll q=n/p;    work(p);work(q);}int main(){    freopen("sum.in","r",stdin);freopen("sum.out","w",stdout);    scanf("%lld%lld",&n,&m);    ss[0][0]=1;    fo(i,1,m){        fo(j,1,i) ss[i][j]=(ss[i-1][j-1]+ss[i-1][j]*(i-1)%mo)%mo;    }    /*t=floor(sqrt(n));    k=n;    fo(i,2,t)        while (k%i==0){            p[++top]=i;            k/=i;        }    if (k>1) p[++top]=k;*/    work(n);    sort(p+1,p+top+1);    tot=top;    top=0;    fo(i,1,tot)        if (p[i]!=p[i-1]){            p[++top]=p[i];            c[top]=1;        }        else c[top]++;    fo(i,0,m){        num=1;        fo(j,1,top) num=num*calcg(i,p[j],c[j])%mo;        num=num*ss[m][i]%mo;        if ((m-i)%2) (ans-=num)%=mo;else (ans+=num)%=mo;    }    fo(i,1,m) ans=ans*qsm(i,mo-2)%mo;    (ans+=mo)%=mo;    printf("%lld\n",ans);}