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 Fprime
        global_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;
}