HackerRank-Coprime Power Sum-类min25筛的DP

link：https://www.hackerrank.com/challenges/coprime-power-sum/problem
题意：给一个长度为 \(n(1\leq n\leq 50)\) 的序列 \(s\)，以及 \(0\leq k\leq 10\)，\(1\leq m\leq 10^{12}\)，对每个 \(i:1\leq i\leq m\)，如果不是任何一个 \(s_j\) 的倍数，就加上 \(i^k\) 的贡献，求最后的贡献是多少，保证 \(s\) 两两互质

---

两两互质可以类似当素数来考虑：因为这样如果某个数 \(x\) 可以表示成 \(s\) 中若干数字的乘积，那么表示一定是唯一的。

类似min25筛的DP，设 \(G_i(n)\) 表示用前 \(i\) 个 \(s\) 筛完 \(1,\dots,n\) 中的数字后的答案，显然 \(f(n)=n^k\) 是个完全积性函数，那么 \(G_i(n)=G_{i-1}(n)-\)这一次会被筛掉的部分的贡献，也就是 \(s_i^k\times G_{i-1}(\lfloor\frac{n}{s_i}\rfloor)\) ，这样就足够了，不重不漏，因为那些公倍数一定在之前就被筛过了。（这里和min25不一样，min25这部分的DP还带有一个“最小质因子”的限制）

然后就是自然数幂和的部分，Lagrange插值处理一下

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
#define fastio ios_base::sync_with_stdio(0);cin.tie(0);cout.tie(0)
using namespace std;
typedef long long ll;
constexpr int MOD=1e9+7;
template <unsigned M_>struct ModInt{static constexpr unsigned M=M_;unsigned x;constexpr ModInt():x(0U){}constexpr ModInt(int x_):x(((x_%=static_cast<int>(M))<0)?(x_+static_cast<int>(M)):x_){}constexpr ModInt(ll x_):x(((x_%=MOD)<0)?(x_+M):x_){}ModInt & operator +=(const ModInt &a){x=((x+=a.x)>=M)?(x-M):x;return *this;}ModInt & operator -=(const ModInt &a){x=((x-=a.x)>=M)?(x+M):x;return *this;}ModInt & operator *=(const ModInt &a){x=(static_cast<unsigned long long>(x)*a.x)%MOD;return *this;}ModInt & operator /=(const ModInt &a){return (*this*=a.inv());}ModInt pow(ll e)const{if(e<0)return inv().pow(-e);ModInt a=*this,b=1;for(;e;e>>=1){if(e&1)b*=a;a*=a;}return b;}ModInt inv()const{unsigned a=M,b=x;int y=0,z=1;for(;b;){const unsigned q=a/b;const unsigned c=a-q*b;a=b;b=c;const int w=y-static_cast<int>(q)*z;y=z;z=w;}assert(a==1U);return ModInt(y);}ModInt operator+()const {return *this;}ModInt operator-()const {ModInt a; a.x = x ? (M - x) : 0U; return a; }ModInt operator+(const ModInt &a)const { return (ModInt(*this) += a); }ModInt operator-(const ModInt &a)const { return (ModInt(*this) -= a); }ModInt operator*(const ModInt &a)const { return (ModInt(*this) *= a); }ModInt operator/(const ModInt &a)const { return (ModInt(*this) /= a); }template <class T> friend ModInt operator+(T a, const ModInt &b) { return (ModInt(a) += b); }template <class T> friend ModInt operator-(T a, const ModInt &b) { return (ModInt(a) -= b); }template <class T> friend ModInt operator*(T a, const ModInt &b) { return (ModInt(a) *= b); }template <class T> friend ModInt operator/(T a, const ModInt &b) { return (ModInt(a) /= b); }friend ostream &operator<<(ostream &os, const ModInt &a) { return os << a.x; }friend istream &operator>>(istream &is, ModInt &a) { return is >> a.x; }
};
typedef ModInt<MOD> Mint;
constexpr int N=105;
Mint fact[N],inv_fact[N];
vector<vector<Mint>> y;
Mint Lagrange(const long long & n,const int &k,const vector<Mint> &y){//static Mint pre[15],suf[15];// static vector<Mint> x(k+1);Mint ans=0;    // for(int i=0;i<=k;i++)x[i]=i;//pre[0]=n;suf[k+1]=1;for(int i=1;i<=k;i++)pre[i]=pre[i-1]*(n-i);for(int i=k;i>=0;i--)suf[i]=suf[i+1]*(n-i);for(int i=0;i<=k;i++){Mint up=(i==0?1:pre[i-1])*suf[i+1];Mint down=inv_fact[i]*inv_fact[k-i]*((k-i&1)?MOD-1:1);ans=(ans+y[i]*up*down) ;}return ans;
}
void init(){fact[0]=1;for(int i=1;i<N;i++)fact[i]=fact[i-1]*i;inv_fact[N-1]=fact[N-1].inv();for(int i=N-1;i>=1;i--)inv_fact[i-1]=inv_fact[i]*i;y.resize(11);for(int k=0;k<=10;k++){y[k].resize(k+10);y[k][0]=0;for(int i=1;i<=k+5;i++)y[k][i]=y[k][i-1]+Mint(i).pow(k);}
}
namespace Min25{constexpr int maxs=2e6;//2sqrt(n)long long global_n,lis[maxs + 1];int lim,cnt;int le[maxs+1],ge[maxs+1];Mint G[maxs+1];#define idx(v) (v<=lim?le[v]:ge[global_n/(v)])Mint init(const ll &n,const int &k,const vector<long long> &s){//calc Fprimeglobal_n=n;cnt=0;for(long long i=1,j,v;i<=n;i=n/j+1){j=n/i;v=j;lis[++cnt]=j;(j<=lim?le[j]:ge[n/j])=cnt;G[cnt]=Lagrange(j,k+1,y[k]);}for(int j=1;j<=s.size()-1;j++){auto p=s[j];auto pk=Mint(p).pow(k);for(int i=1;lis[i]>=p;i++){const int id=idx(lis[i]/p);G[i]-=pk*(G[id]);}}return G[idx(n)];}
}
void solve(){int n,k;long long m;cin>>n>>k>>m;vector<long long> s(n+1);rep(i,1,n)cin>>s[i];Min25::lim=sqrt(m);cout<<Min25::init(m,k,s)<<endl;
}
int main(){fastio;init();int tc;cin>>tc;while(tc--)solve();return 0;
}