一个人的数论

莫比乌斯反演神题

换言说，就算你推出柿子还是不会那个神奇的伯努利数。。。。

那么这样，我们先推一波柿子：

$f_d(n)
\\=\sum_{i=1}^{n}[gcd(i,n)=1]i^d
\\ =\sum_{i=1}^{n}i^d\sum_{t|gcd(i,n)}\mu (t)
\\ =\sum_{i=1}^{n}i^d\sum_{t|i \textit{ } and \textit{ } t|n}\mu(t)
\\ =\sum_{t|n}\mu(t)\sum_{i=1}^{\frac{n}{t}} (it)^d
\\ =\sum_{t|n}\mu(t)t^d\sum_{i=1}^{\frac{n}{t}}i^d$

发现柿子里后面的$\sum_{i=1}^{\frac{n}{t}}i^d$极为类似伯努利数的递推公式：

$S^m(n)=\frac{1}{m+1}\sum_{i=0}^{m}C_{m+1}^{i}B_in^{m+1-i}$

其中$S^m(n)=\sum_{i=0}^{n-1}i^d$

那么我们可以知道题目中的关于这一部分的答案。

这也就是伯努利数的应用，将一个看似不可循环的东西（比如说这题里面的$n$）变成了可以循环的东西（这道题里面用伯努利数把循环复杂度降到$d$）

然后题目就会变得可做，详细证明和代码实现可参考$oi-wiki$

我们先把后面的柿子写成多项式的形式代回到原式（里面的$a_i$就是多项式系数：$\frac{1}{m+1}\sum_{i=0}^{m}C_{m+1}^{i}B_i$）：

$\\=\sum_{t|n}\mu(t)t^d\sum_{i=0}^{d+1}a^i(\frac{n}{t})^i
\\=\sum_{i=0}^{d+1}a_i\sum_{t|n}\mu(t)t^d(\frac{n}{t})^i
\\=\sum_{i=0}^{d+1}a_i\sum_{t|n}\mu(t)N^it^{d-i}
\\=\sum_{i=0}^{d+1}a_iN^i\sum_{t|n}\mu(t)t^{d-i}$

这时我们发现后面的一部分$\sum_{t|n}\mu(t)t^{d-i}$是一个积性函数

设$g_n=\sum_{t|n}\mu(t)t^{d-i}=\sum_{t|n}\mu(t)id_{d-i}(t)$，而$\mu$和$id$都是积性函数，那么他们相乘也是积性函数

我看$oi-wiki$上面只说了积性函数卷积的积性，但没说这个东西所以以为自己错了，其实很显然，用乘法结合律和分配律就可以简单证明

就这样我们发现$g_n$是积性函数，那么$g_n=\prod_{j=1}^{\omega} g({p_j}^{c_j})$

题目里面给的是质因数分解，这个不用我们搞了，那么考虑贡献

拿出一个$\prod$里面的项来说事：

$g(p^c)=\sum_{t|p^c}\mu(t)id_{d-i}(t)$

发现莫比乌斯函数很棒，只有在$c$为$0/1$的时候整个柿子的答案不是$0$，那么

$\\g(p^c)=\sum_{t|p^c}\mu(t)id_{d-i}(t)
\\=\mu(1)id_{d-i}(1)+\mu(p)id_{d-i}(p)
\\=1-p^{d-i}$

这样，整个柿子的化简工作就大功告成了！

$f_d(n)=\sum_{i=0}^{d+1}a_iN^i\prod_{j=1}^{\omega}(1-p_j^{d-i})$

$a$可以直接筛出，后面的也可以直接搞，都是超级好实现的那种，

数据范围还贼小，比较舒服的$O(d^2+d\omega)$

#include<bits/stdc++.h>
#define int long long
using namespace std;
namespace AE86{
    inline int read(){
        int x=0,f=1;char ch=getchar();
        while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
        while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}return x*f;
    }inline void write(int x,char opt='\n'){
        char ch[20];int len=0;if(x<0)x=~x+1,putchar('-');
        do{ch[len++]=x%10+(1<<5)+(1<<4);x/=10;}while(x);
        for(int i=len-1;i>=0;--i)putchar(ch[i]);putchar(opt);}
}using namespace AE86;
const int mod=1e9+7;
int d,w,p[1001],t[1001],N=1,ans;
int C[150][150],inv[150],B[150],xi[150];
inline int qmo(int a,int b){int ans=1,c=mod;a%=c;while(b){if(b&1)ans=ans*a%c;b>>=1;a=a*a%c;}return ans;}
inline void getB(){
    for(int i=0;i<150;i++){C[i][0]=C[i][i]=1;for(int j=1;j<i;j++)C[i][j]=(C[i-1][j]+C[i-1][j-1])%mod;}
    inv[1]=1;for(int i=2;i<150;i++)inv[i]=(mod-mod/i)*inv[mod%i]%mod;
    B[0]=1;
    for(int i=1;i<150;i++){
        int ans=0;if(i==150-1) break;
        for(int j=0;j<i;j++) (ans+=C[i+1][j]*B[j]%mod)%=mod;
        ans=(ans*(-inv[i+1])%mod+mod)%mod; B[i]=ans;
    } for(int i=0;i<=d;i++) xi[d+1-i]=qmo(d+1,mod-2)*C[d+1][i]%mod*B[i]%mod;
}
namespace WSN{
    inline short main(){
        d=read(); w=read(); getB();
        for(int i=1;i<=w;i++) p[i]=read(),t[i]=read(),N=N*qmo(p[i],t[i])%mod;
        for(int i=1;i<=d+1;i++){
            int now=xi[i]*qmo(N,i)%mod,tmp=1;
            for(int j=1;j<=w;j++) tmp=tmp*(1-qmo(p[j],d-i+mod-1)+mod)%mod;
            (ans+=now*tmp%mod)%=mod;
        } write(ans);
        return 0;
    }
}
signed main(){return WSN::main();}