题目描述
定义：
[F(0)=1,F(n)=n
\
A_{i}^j=frac{F(i)}{F(j)}
]
求
[sum_{i=1}^nsum_{j=1}^isum_{k=1}^{j}gcd(A_{i-j}^j imes F(j),A_{j-k}^k imes F(k))
]
(mathcal{Solution})
式子看起来挺花里胡哨的，把里面的东西简单展开一下其实就是求
[sum_{i=1}^nsum_{j=1}^isum_{k=1}^{j}gcd((i-j)!,(j-k)!)
\
sum_{i=1}^nsum_{j=1}^isum_{k=1}^{j}min(i-j,j-k)!
]
遇到多个求和号在一起有很多突破的方向，这里选择考虑贡献。
这个 (min) 似乎很难搞，但是终究是对若干个阶乘求和，则考虑对于一个 (d) ，计算 (d!) 的贡献，容易发现 (d) 在 (left[1,left lfloor frac{n-1}{2} ight floor ight ])内。
上面那个式子放在数轴上考虑，其实就是取出两个连在一起的区间 ([k,j],[j,i])，答案加上这两个区间较短区间的长度的阶乘。（这里的区间长度为右端点减左端点，即所求式子中的形式，后同）
把这两个区间合并成一个区间 ([k,i])，因为较短的区间长度为 (d)，则 ([k,i]) 长度至少为 (2d)。
当 ([k,i]) 的长度 (len=2d) 时没有区别，当 (len>2d) 时较短区间在左侧或右侧都可以。
那么就可以写出答案的式子：
[sum_{d=1}^{left lfloor frac{n-1}{2}
ight
floor}d!left(sum_{len=2d+1}^{n}sum_{i=1}^{n-len}2+sum_{i=1}^{n-2d}1
ight)
]
(d,len) 均为上文定义，(i) 是枚举的大区间的左端点，也就是上文中 ([k,i]) 的 (k)。
后面的括号的左半部分是 (len>2d) 时的贡献次数，右半部分是 (len=2d) 时的贡献次数。
化简一下：
[sum_{d=1}^{left lfloor frac{n-1}{2}
ight
floor}d!left(2sum_{len=2d+1}^{n}(n-len)+(n-2d)
ight)
\
sum_{d=1}^{left lfloor frac{n-1}{2}
ight
floor}d!left(sum_{len=0}^{n-2d}len+(n-2d)
ight)
\
sum_{d=1}^{left lfloor frac{n-1}{2}
ight
floor}d!left(2 imesfrac{(n-2d)(n-2d-1)}{2}+(n-2d)
ight)
\
sum_{d=1}^{left lfloor frac{n-1}{2}
ight
floor}d!left(n-2d
ight)^2
]
最后这个东西长的十分好看，事实上也很好求，把里面的完全平方式展开一下：
[n^2sum_{d=1}^{left lfloor frac{n-1}{2}
ight
floor}d!-nsum_{d=1}^{left lfloor frac{n-1}{2}
ight
floor}d!d+sum_{d=1}^{left lfloor frac{n-1}{2}
ight
floor}d!4d^2
]
发现里面三个求和号都能在 (mathcal{O}(n)) 的复杂度内预处理，预处理后直接计算即可。
时间复杂度 (mathcal{O}(max{n}+T))
#include<iostream>
#include<cstdio>
#define ll long long
const ll mod = 10086001;
template <typename T> T Max(T x, T y) { return x > y ? x : y; }
template <typename T> T Min(T x, T y) { return x < y ? x : y; }
template <typename T> T Add(T x, T y) { return (x + y > mod) ? (x + y - mod) : (x + y); }
template <typename T> T Mod(T x) { return (x >= mod) ? (x - mod) : x; }
template <typename T>
T& read(T& r) {
r = 0; bool w = 0; char ch = getchar();
while(ch < '0' || ch > '9') w = ch == '-' ? 1 : 0, ch = getchar();
while(ch >= '0' && ch <= '9') r = r * 10 + (ch ^ 48), ch = getchar();
return r = w ? -r : r;
}
const int N = 1000100;
ll fac[N], f1[N], f2[N], f3[N];
int T, n[N], mx;
signed main() {
read(T);
for(int i = 1; i <= T; ++i) mx = Max(mx, read(n[i]));
fac[0] = f1[0] = 1;
for(int i = 1; i <= mx; ++i)
fac[i] = fac[i-1] * i % mod,
f1[i] = Add(f1[i-1], fac[i]),
f2[i] = Add(f2[i-1], fac[i] * 4 * i % mod),
f3[i] = Add(f3[i-1], fac[i] * 4 * i % mod * i % mod);
for(int i = 1; i <= T; ++i) {
int m = (n[i]-1)/2;
printf("%lld
", Add(Mod(1ll * n[i] * n[i] % mod * f1[m] % mod - n[i] * f2[m] % mod + mod), f3[m]));
}
return 0;
}