「BZOJ4314」倍数？倍数！

题目链接： https://www.lydsy.com/JudgeOnline/problem.php?id=4314。

神仙题。。感觉每一步都好巧妙啊。。

答案多项式为：

$$F(x)=[y^k]\prod_{i=0}^{n-1}(x^iy+1)$$

这个生成函数有两个变量，但是应该很好理解，就是说如果选当前数字就是$x^i$，否则就是$1$，一共只能选$k$个$x^i$所以拿一个$y$来限制一下，最后取$y^k$这一项。

那么答案就是这个多项式指数为$n$的倍数的系数和。

拿单位根反演搞一下，答案就是：

$$\frac{1}{n}\sum_{i=0}^{n-1}F(\omega_n^i)=[y^k]\frac{1}{n}\sum_{i=0}^{n-1}\prod_{j=0}^{n-1}(\omega_n^{ij}y+1)$$

枚举$n$的约数$d$，令$d=n/t$，式子可以写成：

$$[y^k]\frac{1}{n}\sum_{d|n}\sum_{i=0}^{t-1}\left(\prod_{j=0}^{n-1}(\omega_n^{idj}y+1)\right)[\gcd(i,t)=1]$$

仔细想想可以知道每个$id$只会被算一次。

那么根据单位根的性质，上下可以同除$d$，然后那个连乘就会有一个长度为$n/d$的循环节，所以：

$$[y^k]\frac{1}{n}\sum_{d|n}\sum_{i=0}^{t-1}\left(\prod_{j=0}^{t-1}(\omega_t^{ij}y+1)\right)^d[\gcd(i,t)=1]$$

由于$\gcd(i,t)=1​$，所以单位根指数上的$i​$可以去掉，因为在$\bmod t​$意义下每个数都能被取到，所以：

$$[y^k]\frac{1}{n}\sum_{d|n}\sum_{i=0}^{t-1}\left(\prod_{j=0}^{t-1}(\omega_t^{j}y+1)\right)^d[\gcd(i,t)=1]$$

考虑单位根的本质，$\omega_n^{0,1,\cdots,n-1}$是$x^n-1=0$的$n$个解，那么可以得到等式：

$$x^n-1=\prod_{i=0}^{n-1}(x-\omega_n^i)$$

如果我们代入$x=-y$可以得到：

$$y^n-(-1)^n=\prod_{i=0}^{n-1}(\omega_n^i+y)$$

所以如果我们改一下最开始的式子，改成：

$$\prod_{i=0}^{n-1}(x^i+y)$$

那么取$y^{n-k}$项，意义是一样的，所以式子可以变成：

$$[y^k]\frac{1}{n}\sum_{d|n}\sum_{i=0}^{t-1}\left(\prod_{j=0}^{t-1}(\omega_t^{j}+y)\right)^d[\gcd(i,t)=1]$$ $$[y^{n-k}]\frac{1}{n}\sum_{d|n}\sum_{i=0}^{t-1}\left(y^t-(-1)^t\right)^d[\gcd(i,t)=1]$$

注意到中间一块和$i$没啥关系，所以继续变一下：

$$[y^{n-k}]\frac{1}{n}\sum_{d|n}\varphi(t)\left(y^t-(-1)^t\right)^d$$

由于我们只需要知道第$n-k$项，暴力算二项式展开就好了。

复杂度大概是$O(\sqrt n)$。

#include<bits/stdc++.h>
using namespace std;

void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}

void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');}

#define lf double
#define ll long long 

#define pii pair<int,int >
#define vec vector<int >

#define pb push_back
#define mp make_pair
#define fr first
#define sc second

#define FOR(i,l,r) for(int i=l,i##_r=r;i<=i##_r;i++)

const int maxn = 1e6+10;
const int inf = 1e9;
const lf eps = 1e-8;
const int mod = 1e9+7;

int n,k,pri[maxn],vis[maxn],phi[maxn],tot,fac[maxn],ifac[maxn],inv[maxn];

int qpow(int a,int x) {
    int res=1;
    for(;x;x>>=1,a=1ll*a*a%mod) if(x&1) res=1ll*res*a%mod;
    return res;
}

void gen() {
    phi[1]=1;
    for(int i=2;i<=k;i++) {
        if(!vis[i]) pri[++tot]=i,phi[i]=i-1;
        for(int j=1;j<=tot&&i*pri[j]<=k;j++) {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) {phi[i*pri[j]]=phi[i]*pri[j];break;}
            phi[i*pri[j]]=phi[i]*phi[pri[j]];
        }
    }
    fac[0]=ifac[0]=inv[0]=inv[1]=1;
    // for(int i=1;i<=k;i++) fac[i]=1ll*fac[i-1]*i%mod;
    // ifac[k]=qpow(fac[k],mod-2);
    // for(int i=k-1;i;i--) ifac[i]=1ll*ifac[i+1]*(i+1)%mod;
    for(int i=2;i<=k;i++) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
}

int ans;

// int c(int a,int b) {return 1ll*fac[a]*ifac[b]%mod*ifac[a-b]%mod;}

int c(int a,int b) {
    int res=1;
    for(int i=1;i<=b;i++) res=1ll*res*(a-i+1)%mod*inv[i]%mod;
    return res;
}

void calc(int d) {
    int t=n/d;
    if((n-k)%t!=0) return ;
    int p=(n-k)/t,b=t&1?1:mod-1;
    ans=(ans+1ll*phi[t]*c(d,d-p)%mod*qpow(b,d-p)%mod)%mod;
}

int main() {
    read(n),read(k),gen();
    for(int i=1;i*i<=n;i++) if(n%i==0) calc(i),i*i==n?void():calc(n/i);
    write(1ll*ans*qpow(n,mod-2)%mod);
    return 0;
}