组合数学求和问题解析-优快云博客

题目大意

求 $\sum_{i=1}^nC_{(i,n)}^m$

做法

随手化一下式子变成
$\sum_{d|n}C_{d}^m*\phi(\frac{n}{d})$
发现很想狄利克雷卷积的形式，不妨尝试凑出另一个积性函数。
$\frac{1}{m!}\sum_{d|n}\sum_{i=0}^ms(m,i)*(-1)^{m-i}*d^i*\phi(\frac{n}{d})$
$\frac{1}{m!}\sum_{i=0}^ms(m,i)*(-1)^{m-i}\sum_{d|n}d^i*\phi(\frac{n}{d})$
现在就好了！
两个积性函数狄利克雷卷积还是积性函数。
设 $g_i(n)=\sum_{d|n}d^i*\phi(\frac{n}{d})$
答案是 $\frac{1}{m!}\sum_{i=0}^ms(m,i)*(-1)^{m-i}\Pi_{p^k|n}g_i(p^k)$
考虑计算质数的次幂。
$g_i(p^k)=\sum_{j=0}^kp^{ij}*\phi(p^{k-j})$
$g_i(p^k)=(p-1)\sum_{j=0}^{k-1}p^{k-1+(i-1)j}+p^{ik}$
分解质因数用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);
}