考虑枚举相邻点距离差的比例。显然应使比例值gcd为1以保证不重复统计。确定比例之后,各维坐标的方案数就可以分开考虑。设比例之和为k,则若坐标上限为m,该维坐标取值方案数即为Σm-ki (i=1~⌊m/k⌋),也即⌊m/k⌋·m-k·(⌊m/k⌋+1)·⌊m/k⌋/2,设其为f(m,k)。总方案数即将各维方案数相乘,设为F(k)。
于是得到答案即为ΣkΣa1Σa2……Σac-2 [gcd(a1,a2,……,ac-2,k)=1]·F(k)。套路一波,得到Σk F(k)·(Σd μ(d)·g(k/d)) (d|k),其中g(n)为将n划成c-1份的方案数,也即C(n-1,c-2)。对每个询问暴力一遍,复杂度即为O(Tnm)。注意这里的组合数只能递推,因为值域比模数还大。
卡卡常就过了毕竟正解的复杂度也优不到哪里去。
考虑优化。容易想到整除分块。固定⌊m/k⌋后,要求和的部分是一个关于k的n次多项式。分治NTT暴力求出多项式系数对各项求个前缀和即可。
#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
#define N 100010
#define P 10007
char getc(){char c=getchar();while ((c<'A'||c>'Z')&&(c<'a'||c>'z')&&(c<''||c>'')) c=getchar();return c;}
int gcd(int n,int m){return m==?n:gcd(m,n%m);}
int read()
{
int x=,f=;char c=getchar();
while (c<''||c>'') {if (c=='-') f=-;c=getchar();}
while (c>=''&&c<='') x=(x<<)+(x<<)+(c^),c=getchar();
return x*f;
}
int T,n,c,a[],b[],mobius[N],prime[N],g[][N],f[][][N],p[],C[N][],fac[N],inv[N],cnt;
bool flag[N];
int main()
{
#ifndef ONLINE_JUDGE
freopen("bzoj3434.in","r",stdin);
freopen("bzoj3434.out","w",stdout);
const char LL[]="%I64d\n";
#else
const char LL[]="%lld\n";
#endif
mobius[]=;
for (int i=;i<=N-;i++)
{
if (!flag[i]) prime[++cnt]=i,mobius[i]=-;
for (int j=;j<=cnt&&prime[j]*i<=N-;j++)
{
flag[prime[j]*i]=;
if (i%prime[j]==) break;
mobius[prime[j]*i]=-mobius[i];
}
}
C[][]=;
for (int i=;i<=N-;i++)
for (int j=;j<=min(,i);j++)
C[i][j]=(C[i-][j-]+C[i-][j])%P;
for (int c=;c<=;c++)
for (int i=;i<=N-;i++)
if (mobius[i])
for (int j=i;j<=N-;j+=i)
g[c][j]=(g[c][j]+mobius[i]*(c->j/i-?:C[j/i-][c-])+P)%P;
int s=;
for (int c=;c<=;c++)
for (int i=;i<=N-;i++)
{
int s=;
for (int k=;k<=;k++)
{
f[k][c][i]=(f[k][c][i-]+g[c][i]*s)%P;
s=s*i%P;
}
}
T=read();
while (T--)
{
n=read(),c=read();int m=N,ans=;
for (int i=;i<=n;i++) m=min(m,a[i]=read());
for (int i=;i<=m;i++)
{
int t=m;p[]=;for (int j=;j<=n;j++) t=min(t,a[j]/(b[j]=a[j]/i)),p[j]=;
for (int j=;j<=n;j++)
{
for (int k=j;k>=;k--)
p[k]=(1ll*p[k]*b[j]*a[j]-1ll*p[k-]*(b[j]+)*b[j]/)%P;
p[]=1ll*p[]*b[j]*a[j]%P;
}
for (int j=;j<=n;j++)
ans=(ans+p[j]*(f[j][c][t]-f[j][c][i-])%P+P)%P;
i=t;
}
cout<<ans<<endl;
}
return ;
}