| 12
 3
 4
 5
 6
 7
 8
 9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 
 | #include <bits/stdc++.h>using namespace std;
 #define pb push_back
 #define mp make_pair
 #define F first
 #define S second
 typedef long long LL;
 typedef pair<int, int> pii;
 const int N = 153;
 const double PI = acos(-1.0);
 double D[80][80][80];
 namespace FFT {
 static const int W = 1 << 22;
 static const double PI = acos(-1.0);
 struct Complex {
 double x, y;
 Complex(double _x = 0, double _y = 0) : x(_x), y(_y) {}
 Complex operator + (const Complex &t) const {return Complex(x + t.x, y + t.y);}
 Complex& operator += (const Complex &t) {x += t.x, y += t.y;return *this;}
 Complex operator - (const Complex &t) const {return Complex(x-t.x,y-t.y);}
 Complex& operator -= (const Complex &t) {x -= t.x, y -= t.y;return *this;}
 Complex operator * (const Complex &t) const {return Complex(x * t.x - y * t.y, x * t.y + y * t.x);}
 Complex operator / (const double &t) const {return Complex(x / t,y / t);}
 Complex& operator /= (const double &t) {x /= t, y /= t;return *this;}
 double real() {return x;}
 double imag() {return y;}
 };
 void fft(Complex a[], int n, int rev) {
 for (int i = 1, j = 0, k; i < n; ++i) {
 for (k = n >> 1; k > (j ^= k); k >>= 1);
 if (i < j) swap(a[i], a[j]);
 }
 for (int s = 1, ss = 2; s < n; s <<= 1,ss <<= 1) {
 Complex wn(cos(2 * PI * rev / ss), sin(2 * PI * rev / ss)), w;
 for (int i = 0, j ; i < n; i += ss) {
 for (j = i, w = 1; j < i + s; ++j, w = w * wn) {
 Complex t = w * a[j + s];
 a[j + s] = a[j] - t;
 a[j] = a[j] + t;
 }
 }
 }
 if (rev == -1) for (int i = 0; i < n; ++i) a[i] /= n;
 }
 }
 FFT::Complex A[FFT::W], B[FFT::W];
 inline double sqrr(double x) {
 return x * x * x * x;
 }
 int main() {
 int n, q;
 scanf("%d%d", &n, &q);
 for (int i = 1, x, y, z; i <= n; ++i) {
 scanf("%d%d%d", &x, &y, &z);
 --x, --y, --z;
 int t1 = N * N * x + N * y + z, t2 = N * N * (76 - x) + N * (76 - y) + (76 - z);
 A[t1].x += 1, B[t2].x += 1;
 }
 for (int i = 0; i <= 77; ++i)
 for (int j = 0; j <= 77; ++j)
 for (int k = 0; k <= 77; ++k)
 D[i][j][k] = sqrt(sqrr(i) + sqrr(j) + sqrr(k));
 FFT::fft(A, FFT::W, 1);
 FFT::fft(B, FFT::W, 1);
 for (int i = 0; i < FFT::W; ++i)    A[i] = A[i] * B[i];
 FFT::fft(A, FFT::W, -1);
 while (q--) {
 double ans = 0;
 int a, b, c, d;
 scanf("%d%d%d%d", &a, &b, &c, &d);
 for (int i = 0; i < FFT::W; ++i) {
 int cnt = int(A[i].real() + 0.5);
 if (!cnt)   continue;
 int x = i / (N * N) - 76, y = i / N % N - 76, z = i % N - 76;
 if (!x && !y && !z) continue;
 ans += fabs(a * x + b * y + c * z + d) * cnt / D[abs(x)][abs(y)][abs(z)];
 }
 ans /= n, ans /= n - 1;
 printf("%.10lf\n", ans);
 }
 return 0;
 }
 
 |