Address
https://www.luogu.org/problemnew/show/P4152
https://www.lydsy.com/JudgeOnline/problem.php?id=3434
http://uoj.ac/problem/54
Solution
考虑枚举第一个点的坐标
和
最后一个点的坐标
,
那么容易得出,这两个点连成的线段上的整点(不包括端点)个数为:
对答案的贡献:
考虑枚举 维向量 满足:
那么 就有 种合法取值。
得到答案为:
这是一个显然的 Mobius 反演。考虑把 解决掉:
是一个 次多项式,可以用 的时间求出,设 次项为 。
预处理 后:
可以得出,由于多项式 的系数只有对 整除(下取整)的结果,故 的取值不超过 种。
对下界分块之后使用 的前缀和计算。
复杂度
~~辣鸡 BZOJ 评测机卡常数~~
Code
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define For(i, a, b) for (i = a; i <= b; i++)
#define Mobius(i, n) for (i = 1; i <= n;)
using namespace std;
inline int read()
{
int res = 0; bool bo = 0; char c;
while (((c = getchar()) < '0' || c > '9') && c != '-');
if (c == '-') bo = 1; else res = c - 48;
while ((c = getchar()) >= '0' && c <= '9')
res = (res << 3) + (res << 1) + (c - 48);
return bo ? ~res + 1 : res;
}
const int MaxN = 1e5, E = 23, P = 14, N = MaxN + 5, MOD = 10007;
int tot, pri[N], miu[N], C[N][E], pw[N][P], f[N][E], sum[N][E][P],
n, c, m[P], a[P], b[P];
bool mark[N];
int Min(int a, int b) {return a < b ? a : b;}
void sieve()
{
int i, j, k;
For (i, 0, MaxN) C[i][0] = pw[i][0] = 1;
For (i, 1, MaxN)
{
For (j, 1, Min(18, i))
C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % MOD;
For (j, 1, 11) pw[i][j] = 1ll * pw[i][j - 1] * i % MOD;
}
mark[0] = mark[miu[1] = 1] = 1;
For (i, 2, MaxN)
{
if (!mark[i]) miu[pri[++tot] = i] = -1;
For (j, 1, tot)
{
if (1ll * i * pri[j] > MaxN) break;
mark[i * pri[j]] = 1;
if (i % pri[j] == 0) break;
else miu[i * pri[j]] = -miu[i];
}
}
For (k, 2, 20)
For (i, 1, MaxN)
{
int dl = MaxN / i;
For (j, 1, dl)
f[i * j][k] = (f[i * j][k] + MOD +
C[i - 1][k - 2] * miu[j]) % MOD;
For (j, 0, 11)
sum[i][k][j] = (sum[i - 1][k][j] +
f[i][k] * pw[i][j]) % MOD;
}
}
void prod(int x, int y, int n)
{
int i;
a[n + 1] = b[0] = 0;
For (i, 1, n + 1) b[i] = a[i - 1] * y % MOD;
For (i, 0, n + 1) a[i] = (a[i] * x - b[i] + MOD) % MOD;
}
void work()
{
int i, j, _min, ans = 0;
n = read(); c = read();
For (i, 1, n) m[i] = _min = read();
For (i, 1, n) _min = Min(_min, m[i]);
Mobius (i, _min)
{
int nxt = m[1] / (m[1] / i);
For (j, 2, n) nxt = Min(nxt, m[j] / (m[j] / i)), a[j] = 0;
a[a[0] = 1] = 0;
For (j, 1, n) prod(1ll * m[j] * (m[j] / i) % MOD,
(1ll * (m[j] / i) * (m[j] / i + 1) >> 1) % MOD, j - 1);
For (j, 0, n) ans = (ans + a[j] * (sum[nxt][c][j] -
sum[i - 1][c][j] + MOD)) % MOD;
i = nxt + 1;
}
printf("%d\n", ans);
}
int main()
{
int T = read();
sieve();
while (T--) work();
return 0;
}