BZOJ4314 倍數？倍數！

好神仙啊....

---

題意

在$ [0,n) $中選$ k$個不一樣的數使和爲$ n$的倍數

求方案數

$ n \leq 10^9, \ k \leq 10^3$

---

題解

k能夠放大到1e6的

先不考慮$ k$的限制

對答案構建多項式$ f(x)=\prod\limits_{i=0}^{n-1}(x^i+1)$

答案就是這個多項式全部次數爲$ n$的倍數的項的係數和

考慮單位根反演

$$ans=\frac{1}{n}\sum_{i=0}^{n-1}\prod_{j=0}^{n-1}(w_n^{ij}+1)$$

設$ d=\gcd(n,i),t=\frac{n}{d}$

$$ans=\frac{1}{n}\sum_{d|n}\sum_{i=0}^{t-1}(\prod_{j=0}^{t-1}(w_t^{ij}+1))^d[\gcd(t,i)=1]$$

因爲$\gcd(t,i)=1$，能夠去掉單位根指數上的$ i$

$$ans=\frac{1}{n}\sum_{d|n}\sum_{i=0}^{t-1}(\prod_{j=0}^{t-1}(w_t^{j}+1))^d[\gcd(t,i)=1]$$

考慮$ \prod\limits_{j=0}^{t-1}(w_t^{j}+1)$是什麼

根據定義可知$ w_t^{0..t-1}$是$ x^t-1=0$的$ n$個根

所以有$ x^t-1=\prod\limits_{i=0}^{t-1}(x-w_t^i)$

討論$ n$的奇偶性可得$ \prod\limits_{j=0}^{t-1}(w_t^{j}+1)=1-(-1)^t$

再用歐拉函數進行化簡得$$ans=\frac{1}{n}\sum_{d|n}\phi(t)(1-(-1)^t)^d$$

 

而後考慮有$ k$這個限制怎麼作

咱們再添加一個新變量$ y$，以$ y$爲主元構建多項式$ f(y)=\prod\limits_{i=0}^{n-1}(yx^i+1)$

咱們要求的就是這個多項式$ y^k$的係數

用跟上面相同的方法能夠化簡得最後的答案多項式爲$$ans=\frac{1}{n}\sum_{d|n}\phi(t)(1-(-y)^t)^d$$

因爲只須要知道$y^k$的係數，直接展開就行了

跑的飛快

---

代碼

#include<ctime>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#include<vector>
#define p 1000000007
#define rt register int
#define ll long long
using namespace std;
inline ll read(){
    ll x=0;char zf=1;char ch=getchar();
    while(ch!='-'&&!isdigit(ch))ch=getchar();
    if(ch=='-')zf=-1,ch=getchar();
    while(isdigit(ch))x=x*10+ch-'0',ch=getchar();return x*zf;
}
void write(ll y){if(y<0)putchar('-'),y=-y;if(y>9)write(y/10);putchar(y%10+48);}
void writeln(const ll y){write(y);putchar('\n');}
int k,m,n,x,y,z,cnt,ans;
int phi[1010],ss[1010];bool pri[1010];
int njc[1010],inv[1010];
int ksm(int x,int y=p-2){
    int ans=1;
    for(;y;y>>=1,x=1ll*x*x%p)if(y&1)ans=1ll*ans*x%p;
    return ans;
}
int C(int x,int y){
    int ans=1;
    for(rt i=x;i>=x-y+1;i--)ans=1ll*ans*i%p;
    return 1ll*ans*njc[y]%p;
}
int main(){
    n=read();k=read();phi[1]=1;
    for(rt i=0;i<=1;i++)njc[i]=inv[i]=1;
    for(rt i=2;i<=k;i++){
        inv[i]=1ll*inv[p%i]*(p-p/i)%p;
        njc[i]=1ll*njc[i-1]*inv[i]%p;
    }
    for(rt i=2;i<=k;i++){
        if(!pri[i])ss[++cnt]=i,phi[i]=i-1;
        for(rt j=1;j<=cnt&&i*ss[j]<=k;j++){
            phi[i*ss[j]]=phi[i]*phi[ss[j]];
            pri[i*ss[j]]=1;
            if(i%ss[j]==0){
                phi[i*ss[j]]=phi[i]*ss[j];
                break;
            } 
        }
    }
    int ans=0,invn=ksm(n);
    for(rt d=1;d<=k;d++)if(n%d==0&&k%d==0){
        const int v=k/d;
        int tag=1;
        if((v&1)&&(d&1^1))tag=-tag;
        (ans+=1ll*tag*phi[d]%p*invn%p*C(n/d,k/d)%p)%=p;
    }
    cout<<(ans+p)%p;
    return 0;
}