WC 2014 space.cpp

#include <cstdio>
#include <cstring>
#include <algorithm>
#define INF (1 << 30)
#define maxm 100200
#define mo 10007
#define M 100100
#define maxn 22
#define N 20
using namespace std;
int f[maxn][maxm],js[2][maxn],C[maxm][maxn];
int sum[maxm][maxn][12],gap[maxm],as[maxn];
int pri[maxm],mu[maxm],mark[maxm];
int i,j,k,l,m,n,a,b,mi,o,top,T,c,x,e,r,ans,just;

int Mo(int a) {
	return (a % 10007);
}

void init() {
	int i,j,d; mu[1] = 1;
	for (i = 2; i <= M; i++) {
		if (!mark[i]) {
			pri[++top] = i;
			mu[i] = -1;
		}
		for (j = 1; j <= top && i * pri[j] <= M; j++) {
			mark[i * pri[j]] = 1;
			if (i % pri[j]) mu[i * pri[j]] = -mu[i];
			else break;
		}
	}
	C[0][0] = 1;
	for (i = 1; i <= M; i++) {
		C[i][0] = 1;
		for (j = 1; j <= N; j++) C[i][j] = Mo(C[i - 1][j] + C[i - 1][j - 1]);
	}
	
	for (j = 2; j <= N; j++) {
		for (d = 1; d <= M; d++)
			for (i = 1; i * d <= M; i++)
				(f[j][i * d] += Mo(mu[i] * C[d - 1][j - 2])) %= 10007;
	}

	for (i = 1; i <= M; i++) {
		d = 1;
		for (j = 0; j <= 11; j++) {
			for (k = 2; k <= N; k++)
				sum[i][k][j] = Mo(f[k][i] * d);
			d = Mo(d * Mo(i));
		}
	}
	
	for (i = 1; i <= M; i++)
		for (j = 2; j <= N; j++)
			for (k = 0; k <= 11; k++) {
				sum[i][j][k] = sum[i][j][k] + sum[i - 1][j][k];
				if (sum[i][j][k] > mo) sum[i][j][k] -= mo;
			}
}

int getsum(int a) {
	return Mo(a * (a + 1) >> 1);
}

bool cmp(int a,int b) { return a > b; }

int main() {
	freopen("space.in","r",stdin);
	freopen("space.out","w",stdout);
	init();
	scanf("%d",&T);
	for (; T; T--) {
		scanf("%d %d",&n,&c); mi = INF; ans = 0;
		e = 0;
		for (i = 1; i <= n; i++) scanf("%d",&as[i]), mi = min(mi,as[i]);
		top = 0;
		for (i = 1; i <= n; i++) {
			for (a = 1; a <= as[i];) {
				gap[++top] = a;
				a = as[i] / (as[i] / a) + 1; 
			}
			gap[++top] = as[i] + 1;
		}
		sort(gap + 1,gap + top + 1);
		//memset(js,0,sizeof(js)); o = 0;
		
		for (i = 1; i <= top; i++) {
			if (gap[i] > mi) break;
			for (j = i; j <= top && gap[i] == gap[j]; j++);
			i = j - 1;
			l = gap[i]; r = gap[j] - 1;
			memset(js,0,sizeof(js));
			o = 0;
			js[1][0] = Mo(Mo(as[1]) * (Mo(as[1] / l)));
			js[1][1] = -getsum(Mo(as[1] / l));
			for (j = 2; j <= n; j++) {
				a = Mo(Mo(as[j]) * Mo(as[j] / l));
				b = -getsum(Mo(as[j] / l));
				js[o][0] = Mo(js[o ^ 1][0] * a);
				for (k = 1; k <= j; k++)
					js[o][k] = Mo(js[o ^ 1][k] * a + js[o ^ 1][k - 1] * b);
				o ^= 1;
			}
			for (j = 0; j <= n; j++)
				ans = Mo(ans + Mo(js[o ^ 1][j] * Mo(sum[r][c][j] - sum[l - 1][c][j])));
		}
		printf("%d\n",(ans % 10007 + 10007) % 10007);
 	}
	
	return 0;
}