bzoj3434 [Wc2014]时空穿梭

式子的推导网上有很多dalao写的都很好，就不再赘述了。

有两点是我没有想到的，再次写一下。

首先就是初始的式子。不知道直角三角形斜边上整数点个数=gcd(a,b)+1，更没有拓展到多维。

第二就是划出来60分的式子之后的优化，把后边的关于T的多项式的系数通过O(n2)时间算了出来$O(n^2*m)->O(n^3*\sqrt{m})$

十分巧妙。

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
#define N 100500
#define mod 10007
using namespace std;
int C[N][25],mu[N],h[N],nn,prime[N],tot,vis[N],g[25][N],pw[N][25],sum[25][25][N];
void init(){
    C[0][0]=1;
    for(int i=1;i<=100000;i++){
        C[i][0]=1;
        for(int j=1;j<=20&&j<=i;j++)
            C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
    }
    mu[1]=1;h[++nn]=1;
    for(int i=2;i<=100000;i++){
        if(!vis[i]){
            prime[++tot]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=tot&&i*prime[j]<=100000;j++){
            vis[i*prime[j]]=1;
            if(i%prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
        if(mu[i])h[++nn]=i;
    }
    for(int k=2;k<=20;k++){
        for(int i=1;i<=100000;i++)
            for(int j=1;j<=nn&&h[j]*i<=100000;j++)
                (g[k][i*h[j]]+=mu[h[j]]*C[i-1][k-2]+mod)%=mod;
    }
    for(int i=1;i<=100000;i++){
        pw[i][0]=1;
        for(int j=1;j<=11;j++)pw[i][j]=pw[i][j-1]*i%mod;
    }
    for(int k=1;k<=20;k++)
        for(int j=0;j<=11;j++)
            for(int i=1;i<=100000;i++)
                sum[k][j][i]=(sum[k][j][i-1]+g[k][i]*pw[i][j]%mod)%mod;
}
int T,n,c,m[25],ans,f[25];
void work(int x){
    memset(f,0,sizeof f);f[0]=1;
    for(int i=1;i<=n;i++){
        for(int j=i;j>0;j--)
            f[j]=(f[j]*2%mod*m[i]%mod-f[j-1]*(m[i]/x+1)%mod+mod)%mod;
        f[0]=f[0]*2*m[i]%mod;
    }
}
int main(){
    init();
    scanf("%d",&T);
    while(T--){
        scanf("%d%d",&n,&c);ans=0;
        for(int i=1;i<=n;i++)scanf("%d",&m[i]);
        sort(m+1,m+n+1);
        for(int i=1,pos;i<=m[1];i=pos+1){
            pos=N;
            for(int j=1;j<=n;j++)pos=min(pos,m[j]/(m[j]/i));
            int now=1;
            for(int j=1;j<=n;j++)
                now=(now*(m[j]/i)%mod*5004%mod);
            work(i);
            for(int j=0;j<=n;j++)
                (ans+=now*f[j]%mod*(sum[c][j][pos]-sum[c][j][i-1]+mod)%mod)%=mod;
        }
        printf("%d\n",ans);
    }
    return 0;
}