http://acm.hdu.edu.cn/showproblem.php?pid=3571

比賽前一天晚上剛跟haozi研究了最小包圍球,其中要根據(jù)四個點求出球心坐標(biāo),然后我提出了這么拓展到N維的情況。今天比賽的時候居然就出了這么一題,雖然知道可以用高斯消元,但是這題的數(shù)據(jù)好大,有1017那么大,而且要是整數(shù)解,用double精度根本不夠。我想利用大整數(shù),要隊友寫了一個java的,每次消元的時候取最小公倍數(shù),結(jié)果TLE了。。。比賽的時候只有haozi他們隊和電子科大的一個隊伍過了,Orz!!!

賽后,問了一下haozi,因為最后的答案在[ -1017  , 1017] 的范圍內(nèi),可以取一個大于2×1017的素數(shù)p,計算的過程中都 mod p,由于坐標(biāo)有正有負(fù),可以先將球心坐標(biāo)都加上1017  這么一個偏移量,最后求出答案之后在剪掉。
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     forint t = 1; t <= cas; t++ )
 90     {
 91         scanf("%d",&n);
 92         forint i = 0; i <= n; i++ )
 93         {
 94             b[i] = 0;
 95             forint 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         forint i = 0; i < n; i++ )
104         {
105             forint 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         forint i = 1; i < n; i++ )
114             printf(" %I64d",d[i]-inf);
115         printf("\n");
116     }
117     return 0;
118 }