【HDU】3571 N-dimensional Sphere

题意：给出n维的n+1个坐标，求一个坐标，使得n+1到该坐标的距离相等。距离定义为{(x1, x2 ... xN)| ∑(xi-Xi)^2 = R^2 (i=1,2,...,N) }

首先，列出n+1个方程，显然二次项系数都为1，所以可以先把二次项消去，得到n个n元一次方程。

然后，开始高斯消元。由于java在n=50时，很久都出不了答案，c++又会爆longlong，所以我们可以选择模P。

ax=b的解等于ax=b(mod P)的解（x<P），这个很好证明。

所以在模P下，加减乘除都很容易搞定。除法用乘以逆元，乘法用二分加法，拓展欧几里德部分的乘法是不会溢出的。

为了方便处理，把输入所有的值加上一个偏移量。

  1 #include<cstdio>
  2 #include<cstring>
  3 #include<algorithm>
  4 typedef long long LL;
  5 #define MAXN 60
  6 #define P 200000000000000003LL
  7 #define S 100000000000000000LL
  8 using namespace std;
  9 LL x[MAXN], g[MAXN][MAXN], a[MAXN][MAXN], b[MAXN][MAXN];
 10 int n;
 11 inline LL Mod(LL x) {
 12     if (x >= P)
 13         return x - P;
 14     return x;
 15 }
 16 LL MulMod(LL a, LL b) {
 17     LL res;
 18     for (res = 0; b; b >>= 1) {
 19         if (b & 1)
 20             res = Mod(res + a);
 21         a = Mod(a + a);
 22     }
 23     return res;
 24 }
 25 LL ExtGcd(LL a, LL b, LL &x, LL &y) {
 26     if (b == 0) {
 27         x = 1;
 28         y = 0;
 29         return a;
 30     }
 31     LL t, d;
 32     d = ExtGcd(b, a % b, x, y);
 33     t = x;
 34     x = y;
 35     y = t - a / b * y;
 36     return d;
 37 }
 38 LL InvMod(LL a, LL n) {
 39     LL x, y;
 40     ExtGcd(a, n, x, y);
 41     return (x % n + n) % n;
 42 }
 43 void Gauss() {
 44     int i, j, k;
 45     LL inv, tmp;
 46     for (i = 0; i < n; i++) {
 47         for (j = i; j < n; j++) {
 48             if (g[j][i])
 49                 break;
 50         }
 51         if (i != j) {
 52             for (k = i; k <= n; k++)
 53                 swap(g[i][k], g[j][k]);
 54         }
 55         inv = InvMod(g[i][i], P);
 56         for (j = i + 1; j < n; j++) {
 57             if (g[j][i]) {
 58                 tmp = MulMod(g[j][i], inv);
 59                 for (k = i; k <= n; k++) {
 60                     g[j][k] -= MulMod(tmp, g[i][k]);
 61                     g[j][k] = (g[j][k] % P + P) % P;
 62                 }
 63             }
 64         }
 65     }
 66     for (i = n - 1; i >= 0; i--) {
 67         tmp = 0;
 68         for (j = i + 1; j < n; j++) {
 69             tmp += MulMod(x[j], g[i][j]);
 70             if (tmp >= P)
 71                 tmp -= P;
 72         }
 73         tmp = g[i][n] - tmp;
 74         tmp = (tmp % P + P) % P;
 75         x[i] = MulMod(tmp, InvMod(g[i][i], P));
 76     }
 77 }
 78 int main() {
 79     int c, ca = 1;
 80     int i, j;
 81     LL tmp;
 82     scanf("%d", &c);
 83     while (c--) {
 84         scanf("%d", &n);
 85         memset(g, 0, sizeof(g));
 86         memset(b, 0, sizeof(b));
 87         for (i = 0; i <= n; i++) {
 88             for (j = 0; j < n; j++) {
 89                 scanf("%I64d", &a[i][j]);
 90                 a[i][j] += S;
 91                 b[i][n] += MulMod(a[i][j], a[i][j]);
 92                 if (b[i][n] >= P)
 93                     b[i][n] -= P;
 94             }
 95         }
 96         for (i = 0; i < n; i++) {
 97             for (j = 0; j < n; j++) {
 98                 tmp = a[i + 1][j] - a[i][j];
 99                 tmp = (tmp % P + P) % P;
100                 g[i][j] = MulMod(tmp, 2);
101             }
102             g[i][n] = b[i + 1][n] - b[i][n];
103             g[i][n] = (g[i][n] % P + P) % P;
104         }
105         Gauss();
106         printf("Case %d:\n", ca++);
107         printf("%I64d", x[0] - S);
108         for (i = 1; i < n; i++)
109             printf(" %I64d", x[i] - S);
110         putchar('\n');
111     }
112     return 0;
113 }
posted on 2012-09-04 15:18  DrunBee  阅读(259)  评论(0)  编辑  收藏  举报