这题确实在能力范围之外，现在学反演的瓶颈找到了，就是除了套路之外啥也不敢想，然后应该自己设计的线性筛总是设计不了
（这题想了也没法做……这求多项式的手段真的让我大开眼界了……）
Description
给定一个大整数 (n) 的质因数分解形式和一个整数 (d)
求不大于 (n) 且和 (n) 互质的数的 (d) 次方和
(dle 1000)
Solution
上手推式子
[f_d(n)=sum_{p|n} p^dmu(p) sum_{i=1}^{lfloor frac n p
floor} i^d
]
然后我们发现这东西并不可做，在 (n) 这么大的情况下开始毫无头绪
以下参考了这篇题解link
设(S_d(n)=sumlimits_{i=1}^n i^d)
那么则有
[f_d(n)=sum_{p|n} p^dmu(p) S_d(lfloor frac n p
floor)
]
猜 (S_d(n)) 是一个多项式，那么有如下表示（好像有证明？link)
这一部分等过一段时间填坑吧……（先背个结论）
[S_d(n) =sum_{i=0}^{d+1} a_i n^i
]
我们这里可以用任意插值 (+) 高斯消元得到这个多项式
接着推式子，得到
[f_d(n)=sum_{i=0}^{d+1} a_i imes (sum_{p|n} p^d (frac n p)^i mu (p) )
]
然后看看后面的部分就是两个函数 (u(x)=x^d mu (x)) 和 (v(x)=x^i) 的狄利克雷卷积
令 (h=u * v)
(u,v) 都是积性函数，所以它们的狄利克雷卷积也是积性函数（可以简单证明一下，最后发现推出来要证 (mu) 是积性函数，所以就证毕了），所以要线性筛
因为它是个积性函数，所以可以把 (n) 分解质因数然后乘起来
(h_i(p^a)) 把 (mu) 的值为 (0) 的部分消掉之后剩下了：
[p^{ai}(1-p^{d-i})
]
然后是代码实现……（怎么感觉都是知识点门槛很高但是没啥水准）
[inom n m
]
我数据范围开小了卡了 (3h)
Code
#include<bits/stdc++.h>
using namespace std;
#define int long long
namespace yspm{
inline int read()
{
int res=0,f=1; char k;
while(!isdigit(k=getchar())) if(k=='-') f=-1;
while(isdigit(k)) res=res*10+k-'0',k=getchar();
return res*f;
}
const int mod=1e9+7,N=1010;
inline int ksm(int x,int y)
{
if(y<0) return ksm(ksm(x,-y),mod-2);
int res=1; for(;y;y>>=1,(x*=x)%=mod) if(y&1) res=res*x%mod;
return res;
}
int ans,n,d,p[N][2],m[N][N],res[N],sum[N];
inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
signed main()
{
d=read(); n=read();
for(int i=1;i<=n;++i) p[i][0]=read(),p[i][1]=read();
for(int i=0;i<=d+2;++i)
{
int res=0; for(int j=1;j<=i+1;++j) res=add(res,ksm(j,d));
m[i][d+2]=res; m[i][0]=1;
for(int j=1;j<=d+1;++j) m[i][j]=m[i][j-1]*(i+1)%mod;
}
for(int i=0;i<=d+2;++i)
{
int t=i;
for(int j=i;j<=d+2;++j) if(m[j][i]){t=j; break;}
for(int j=i;j<=d+2;++j) swap(m[i][j],m[t][j]);
for(int j=i+1;j<=d+2;++j)
{
int x=m[j][i]*ksm(m[i][i],mod-2)%mod;
for(int k=i;k<=2+d;++k) m[j][k]=(m[j][k]-m[i][k]*x%mod+mod)%mod;
}
}
for(int i=d+2;i>=0;--i)
{
m[i][d+2]=m[i][d+2]*ksm(m[i][i],mod-2)%mod;
for(int j=i-1;j>=0;--j) m[j][d+2]-=m[i][d+2]*m[j][i]%mod,m[j][d+2]=(m[j][d+2]+mod)%mod;
}
for(int i=0;i<=d+1;++i) res[i]=m[i][d+2];
for(int i=0;i<=d+1;++i)
{
int tmp=1;
for(int j=1;j<=n;++j)
{
tmp=tmp*ksm(p[j][0],p[j][1]*i)%mod*((1-ksm(p[j][0],d-i)+mod)%mod)%mod;
}
ans=add(ans,res[i]*tmp%mod);
} cout<<ans<<endl;
return 0;
}
}
signed main(){return yspm::main();}