http://acm.hdu.edu.cn/showproblem.php?pid=3571比赛前一天晚上刚跟haozi研究了最小包围球,其中要根据四个点求出球心坐标,然后我提出了这么拓展到N维的情况。今天比赛的时候居然就出了这么一题,虽然知道可以用高斯消元,但是这题的数据好大,有10
17那么大,而且要是整数解,用double精度根本不够。我想利用大整数,要队友写了一个java的,每次消元的时候取最小公倍数,结果TLE了。。。比赛的时候只有haozi他们队和电子科大的一个队伍过了,Orz!!!
赛后,问了一下haozi,因为最后的答案在[ -10
17 , 10
17]
的范围内,可以取一个大于2×10
17的素数p,计算的过程中都 mod p,由于坐标有正有负,可以先将球心坐标都加上10
17 这么一个偏移量,最后求出答案之后在剪掉。
PS: 高斯消元改自浙大模板
1 #include <cstdio>
2 #include <iostream>
3 #include <map>
4 #include <queue>
5 #include <complex>
6 #include <algorithm>
7 #include <cmath>
8 #include <cstring>
9 using namespace std;
10 typedef __int64 ll;
11
12 const ll p = 1000000000000000003LL;
13 const int maxn = 60;
14 const ll inf = 100000000000000000LL;
15
16 ll mod( ll a, ll n ) { return ( a % n + n ) % n; }
17
18 ll mul_mod ( ll a, ll b )
19 {
20 ll ret = 0, tmp = a % p;
21 while( b )
22 {
23 if( b & 1 )
24 if( ( ret += tmp ) >= p )
25 ret -= p;
26 if( ( tmp <<= 1 ) >= p ) tmp -= p;
27 b >>= 1;
28 }
29 return ret;
30 }
31
32 void gcd( ll a, ll b, ll & d, ll & x, ll & y )
33 {
34 if( ! b ) { d = a; x = 1; y = 0; }
35 else
36 {
37 gcd( b, a % b, d, y, x );
38 y -= x * ( a / b );
39 }
40 }
41
42 ll inv( ll a, ll n )
43 {
44 ll d, x, y;
45 gcd( a, p, d, x, y );
46 return d == 1 ? mod( x, p ) : -1;
47 }
48
49 ll ABS( ll a ) { return ( a < 0 ? -a : a ); }
50
51 int Gauss( int n, ll a[][maxn], ll b[] )
52 {
53 int i, j, k, row;
54 ll maxp, t;
55 for( k = 0; k < n; k++ )
56 {
57 for( maxp = 0, i = k; i < n; i++ )
58 if( ABS( a[i][k] ) > ABS( maxp ) )
59 maxp = a[row=i][k];
60 if( maxp == 0 ) return 0;
61 if( row != k )
62 {
63 for( j = k; j < n; j++ )
64 t = a[k][j], a[k][j] = a[row][j], a[row][j] = t;
65 t = b[k], b[k] = b[row], b[row] = t;
66 }
67 ll INV = inv( maxp, p );
68 for( j = k + 1; j < n; j++ )
69 {
70 a[k][j] = mul_mod( a[k][j], INV );
71 for( i = k + 1; i < n; i++ )
72 a[i][j] = mod( a[i][j] - mul_mod( a[i][k], a[k][j] ), p );
73 }
74 b[k] = mul_mod( b[k], INV );
75 for( i = k + 1; i < n; i++ )
76 b[i] = mod( b[i] - mul_mod( b[k], a[i][k] ), p );
77 }
78 for( i = n - 1; i >= 0; i-- )
79 for( j = i + 1; j < n; j++ )
80 b[i] = mod( b[i] - mul_mod( a[i][j], b[j] ), p );
81 return 1;
82 }
83
84 int main(int argc, char *argv[])
85 {
86 int cas, n;
87 ll a[maxn][maxn], b[maxn], c[maxn][maxn], d[maxn];
88 scanf("%d",&cas);
89 for( int t = 1; t <= cas; t++ )
90 {
91 scanf("%d",&n);
92 for( int i = 0; i <= n; i++ )
93 {
94 b[i] = 0;
95 for( int j = 0; j < n ; j++ )
96 {
97 scanf("%I64d",&a[i][j]);
98 a[i][j] += inf;
99 b[i] = ( b[i] + mul_mod( a[i][j], a[i][j] ) ) % p;
100 a[i][j] = ( a[i][j] + a[i][j] ) % p;
101 }
102 }
103 for( int i = 0; i < n; i++ )
104 {
105 for( int j = 0; j < n; j++ )
106 c[i][j] = mod( a[i][j] - a[n][j], p );
107 d[i] = mod( b[i] - b[n], p );
108 }
109 Gauss( n, c, d );
110 //gauss_cpivot( n, c, d );
111 printf("Case %d:\n",t);
112 printf("%I64d",d[0]-inf);
113 for( int i = 1; i < n; i++ )
114 printf(" %I64d",d[i]-inf);
115 printf("\n");
116 }
117 return 0;
118 }