式子的推导网上有很多dalao写的都很好，就不再赘述了。
有两点是我没有想到的，再次写一下。
首先就是初始的式子。不知道直角三角形斜边上整数点个数=gcd(a,b)+1，更没有拓展到多维。
第二就是划出来60分的式子之后的优化，把后边的关于T的多项式的系数通过O(n2)时间算了出来$O(n^2*m)->O(n^3*\sqrt{m})$
十分巧妙。
1 #include<cstdio> 2 #include<cstring> 3 #include<iostream> 4 #include<algorithm> 5 #include<cmath> 6 #define N 100500 7 #define mod 10007 8 using namespace std; 9 int C[N][25],mu[N],h[N],nn,prime[N],tot,vis[N],g[25][N],pw[N][25],sum[25][25][N]; 10 void init(){ 11 C[0][0]=1; 12 for(int i=1;i<=100000;i++){ 13 C[i][0]=1; 14 for(int j=1;j<=20&&j<=i;j++) 15 C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod; 16 } 17 mu[1]=1;h[++nn]=1; 18 for(int i=2;i<=100000;i++){ 19 if(!vis[i]){ 20 prime[++tot]=i; 21 mu[i]=-1; 22 } 23 for(int j=1;j<=tot&&i*prime[j]<=100000;j++){ 24 vis[i*prime[j]]=1; 25 if(i%prime[j]==0){ 26 mu[i*prime[j]]=0; 27 break; 28 } 29 mu[i*prime[j]]=-mu[i]; 30 } 31 if(mu[i])h[++nn]=i; 32 } 33 for(int k=2;k<=20;k++){ 34 for(int i=1;i<=100000;i++) 35 for(int j=1;j<=nn&&h[j]*i<=100000;j++) 36 (g[k][i*h[j]]+=mu[h[j]]*C[i-1][k-2]+mod)%=mod; 37 } 38 for(int i=1;i<=100000;i++){ 39 pw[i][0]=1; 40 for(int j=1;j<=11;j++)pw[i][j]=pw[i][j-1]*i%mod; 41 } 42 for(int k=1;k<=20;k++) 43 for(int j=0;j<=11;j++) 44 for(int i=1;i<=100000;i++) 45 sum[k][j][i]=(sum[k][j][i-1]+g[k][i]*pw[i][j]%mod)%mod; 46 } 47 int T,n,c,m[25],ans,f[25]; 48 void work(int x){ 49 memset(f,0,sizeof f);f[0]=1; 50 for(int i=1;i<=n;i++){ 51 for(int j=i;j>0;j--) 52 f[j]=(f[j]*2%mod*m[i]%mod-f[j-1]*(m[i]/x+1)%mod+mod)%mod; 53 f[0]=f[0]*2*m[i]%mod; 54 } 55 } 56 int main(){ 57 init(); 58 scanf("%d",&T); 59 while(T--){ 60 scanf("%d%d",&n,&c);ans=0; 61 for(int i=1;i<=n;i++)scanf("%d",&m[i]); 62 sort(m+1,m+n+1); 63 for(int i=1,pos;i<=m[1];i=pos+1){ 64 pos=N; 65 for(int j=1;j<=n;j++)pos=min(pos,m[j]/(m[j]/i)); 66 int now=1; 67 for(int j=1;j<=n;j++) 68 now=(now*(m[j]/i)%mod*5004%mod); 69 work(i); 70 for(int j=0;j<=n;j++) 71 (ans+=now*f[j]%mod*(sum[c][j][pos]-sum[c][j][i-1]+mod)%mod)%=mod; 72 } 73 printf("%d\n",ans); 74 } 75 return 0; 76 }