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