BZOJ3434 WC2014时空穿梭（莫比乌斯反演）

　　考虑枚举相邻点距离差的比例。显然应使比例值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……Σa_c-2 [gcd(a₁,a₂,……,a_c-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<'0'||c>'9')) c=getchar();return c;}
int gcd(int n,int m){return m==0?n:gcd(m,n%m);}
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
int T,n,c,a[30],b[30],mobius[N],prime[N],g[30][N],f[15][30][N],p[30],C[N][30],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[1]=1;
    for (int i=2;i<=N-10;i++)
    {
        if (!flag[i]) prime[++cnt]=i,mobius[i]=-1;
        for (int j=1;j<=cnt&&prime[j]*i<=N-10;j++)
        {
            flag[prime[j]*i]=1;
            if (i%prime[j]==0) break;
            mobius[prime[j]*i]=-mobius[i];
        }
    }
    C[0][0]=1;
    for (int i=1;i<=N-10;i++)
        for (int j=0;j<=min(20,i);j++)
        C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
    for (int c=2;c<=20;c++)
        for (int i=1;i<=N-10;i++)
        if (mobius[i])
            for (int j=i;j<=N-10;j+=i)
            g[c][j]=(g[c][j]+mobius[i]*(c-2>j/i-1?0:C[j/i-1][c-2])+P)%P;
    int s=1;
    for (int c=2;c<=20;c++)
        for (int i=1;i<=N-10;i++)
        {
            int s=1;
            for (int k=0;k<=11;k++)
            {
                f[k][c][i]=(f[k][c][i-1]+g[c][i]*s)%P;
                s=s*i%P;
            }
        }
    T=read();
    while (T--)
    {
        n=read(),c=read();int m=N,ans=0;
        for (int i=1;i<=n;i++) m=min(m,a[i]=read());
        for (int i=1;i<=m;i++)
        {
            int t=m;p[0]=1;for (int j=1;j<=n;j++) t=min(t,a[j]/(b[j]=a[j]/i)),p[j]=0;
            for (int j=1;j<=n;j++) 
            {
                for (int k=j;k>=1;k--)
                p[k]=(1ll*p[k]*b[j]*a[j]-1ll*p[k-1]*(b[j]+1)*b[j]/2)%P;
                p[0]=1ll*p[0]*b[j]*a[j]%P;
            }
            for (int j=0;j<=n;j++)
            ans=(ans+p[j]*(f[j][c][t]-f[j][c][i-1])%P+P)%P;
            i=t;
        }
        cout<<ans<<endl;
    }
    return 0;
}