@codechef - MGCH3D@ 3D Queries
@description@
在三维空间内有 N 个不同的点，请计算下面式子的值 Q 次：
其中 A、B、C 和 D 在每次计算中都会被重新指定。
输入格式
输入数据第一行包含两个整数 N 和 Q，分别表示点的数量和询问的次数。
接下来的 N 行，每行包含三个整数 Xi、Yi 和 Zi 表示点坐标。
接下来的 Q 行，每行包含 4 个整数 A、B、C 和 D，表示一次询问。
输出格式
对于每次询问，输出一行包含相应的答案，精确到 10^(−6) 以上。
数据范围
• 2 ≤ N ≤ 777777
• 1 ≤ Q ≤ 77
• 1 ≤ Xi,Yi,Zi ≤ 77
• 1 ≤ Ai,Bi,Ci ≤ 77
• 1 ≤ Di ≤ 7777
• 愿数字 7 给你带来好运
样例数据
输入
10 5
45 70 41
9 1 43
1 68 8
70 76 7
1 19 33
71 70 53
42 54 71
11 13 30
16 63 25
30 24 34
56 61 29 7328
63 32 18 365
37 41 11 2332
36 19 43 7432
68 55 46 6338
输出
6.692386875130186
1.323651551014940
2.269817185835997
6.783038317971530
5.816449269601737
@solution@
还以为是什么神仙推式子数学题。。。结果没想到就是道 sb 题 -_-。。。
考虑题目给定的式子，实际上与 xi, xj 这些变量无关而只与 (xi - xj), (yi - yj), (zi - zj) 这三个量有关。
于是我们只需要统计各个三元组分别的出现次数，就可以在 (77*2)^3 时间内完成一次询问。
考虑在一维的时候，实际上就是给你 n 个不同的数 x[1...n]，对于每个非零的 a 求 \((\sum_{x[i]-x[j]=a}1)\)。
woc 这么像一个卷积，那我们就构造一个卷积得了。
令 \(f(p) = \sum_{i=1}^{n}p^{x_i}, g(p) = \sum_{i=1}^{n}p^{-x_i}\)，把 f 和 g 作卷积就好了。
负指数幂加一个常数就变成正指数了。
三维的就用高维卷积，实际上就是把它转为 2*77 进制数，最低位为 x，次第位为 y，以此类推。
因为点互不相同，所以一个点对一个三元组的贡献最多为 1，那么答案不超过 N。所以 ntt 不会产生问题。
时间复杂度 O((77*2)^3*q + (77*2)3*log((77*2)3))
@accepted code@
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
const int G = 3;
const int MOD = 998244353;
const int N = 77;
const int M = 4194304;
int pow_mod(int b, int p) {
int ret = 1;
while( p ) {
if( p & 1 ) ret = 1LL*ret*b%MOD;
b = 1LL*b*b%MOD;
p >>= 1;
}
return ret;
}
int pw[25], ipw[25];
void ntt(int *A, int len, int type) {
for(int i=0,j=0;i<len;i++) {
if( i < j ) swap(A[i], A[j]);
for(int k=(len>>1);(j^=k)<k;k>>=1);
}
for(int i=1;(1<<i)<=len;i++) {
int s = (1<<i), t = (s>>1);
int u = (type == 1 ? pw[i] : ipw[i]);
for(int j=0;j<len;j+=s) {
for(int k=0,p=1;k<t;k++,p=1LL*p*u%MOD) {
int x = A[j+k], y = 1LL*A[j+k+t]*p%MOD;
A[j+k] = (x + y)%MOD, A[j+k+t] = (x + MOD - y)%MOD;
}
}
}
if( type == -1 ) {
int inv = pow_mod(len, MOD-2);
for(int i=0;i<len;i++)
A[i] = 1LL*A[i]*inv%MOD;
}
}
void init() {
for(int i=0;i<25;i++)
pw[i] = pow_mod(G, (MOD - 1)/(1<<i)), ipw[i] = pow_mod(pw[i], MOD - 2);
}
int n, q;
double a, b, c, d;
double pw4(double x) {
return x*x*x*x;
}
double func(double x, double y, double z) {
return fabs(a*x + b*y + c*z + d)/sqrt(pw4(x) + pw4(y) + pw4(z))/n/(n-1);
}
int A[2*N + 5][2*N + 5][2*N + 5], B[2*N + 5][2*N + 5][2*N + 5], C[2*N + 5][2*N + 5][2*N + 5];
int f[M], g[M];
int main() {
init(); scanf("%d%d", &n, &q);
for(int i=1;i<=n;i++) {
int X, Y, Z; scanf("%d%d%d", &X, &Y, &Z);
X--, Y--, Z--; A[X][Y][Z]++;
}
for(int i=0;i<N;i++)
for(int j=0;j<N;j++)
for(int k=0;k<N;k++)
B[i][j][k] = A[N-i-1][N-j-1][N-k-1];
for(int i=0;i<2*N-1;i++)
for(int j=0;j<2*N-1;j++)
for(int k=0;k<2*N-1;k++)
f[(i*(2*N-1) + j)*(2*N-1) + k] = A[i][j][k], g[(i*(2*N-1) + j)*(2*N-1) + k] = B[i][j][k];
ntt(f, M, 1), ntt(g, M, 1);
for(int i=0;i<M;i++)
f[i] = 1LL*f[i]*g[i]%MOD;
ntt(f, M, -1);
for(int i=0;i<2*N-1;i++)
for(int j=0;j<2*N-1;j++)
for(int k=0;k<2*N-1;k++)
C[i][j][k] = f[(i*(2*N-1) + j)*(2*N-1) + k];
/*
for(int i=0;i<N;i++)
for(int j=0;j<N;j++)
for(int k=0;k<N;k++)
for(int p=0;p<N;p++)
for(int q=0;q<N;q++)
for(int r=0;r<N;r++)
C[i+p][j+q][k+r] = (C[i+p][j+q][k+r] + 1LL*A[i][j][k]*B[p][q][r]%MOD)%MOD;
*/
for(int i=1;i<=q;i++) {
scanf("%lf%lf%lf%lf", &a, &b, &c, &d);
double ans = 0;
for(int i=-(N-1);i<=N-1;i++)
for(int j=-(N-1);j<=N-1;j++)
for(int k=-(N-1);k<=N-1;k++)
if( (i || j || k) && C[N-1+i][N-1+j][N-1+k] )
ans += C[N-1+i][N-1+j][N-1+k]*func(i, j, k);
printf("%.10lf\n", ans);
}
}
@details@
其实高维卷积还可以通过一维一维的 ntt 来搞，这样子 log 会小些，但是边长会由 77*2 变为 256，反而会慢些。
看来以后想题不能直接想复杂的，要坚定地认为“既然出题人出得出来那它一定可以做”。
