 /**//*
题意:问n+1个n维点的球心,即其他点到该点距离相等
设距离为R,很容里列出n+1个方程,但是有2次项,只要下一行减上一行即可
得到一个n个方程n个未知数,用高斯消元
矩阵为
A[i][j] = (V[i+1][j] - V[i][j])*2
A[i][n+1] = V[i+1][j]^2 - V[i][j]^2 (只是增广的那一项)

会超long long 但用java高精度会超时!
1 <= N <= 50, |xi| <= 10^17 而且答案是整数解且 |Xi| <= 10^17
只要运算过程中取模一个大于 2*10^17的素数即可
Ax = B
=> Ax %p = B %p
=> A%p x = B %p
即原方程的解等价于系数为A%p的方程的解

注意:
由于答案可以是负数,%p的话结果是[0,p),不对
所以一开始对所有坐标都平移一下,都加上10^17,最后答案减去10^17即正确了

看http://www.cppblog.com/AmazingCaddy/archive/2010/08/25/124608.html

我自己写的高斯消元超时,太慢了
他这种是对row这一行都除以A[row][row]使得A[row][row]为1,很快
由于取余,所以可以直接除(乘以逆元),不怕整不整除问题!!
*/
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>

using namespace std;


const long long MOD = 1000000000000000003LL;
const long long p = 1000000000000000003LL;
const long long INF = 100000000000000000LL;
const int MAXN = 60;

typedef long long ll;

long long v[MAXN][MAXN];
long long A[MAXN][MAXN];

inline long long imod(long long a)
  {
return (a%MOD + MOD)%MOD;
}

inline long long iabs(long long a)
  {
return a>0?a:-a;
}

inline long long imul(long long a,long long b)
  {
long long ans = 0 , n = iabs(b);
while(n)
 {
if(n&1)
 {
ans += a;
if(ans >= MOD) ans -= MOD;
}
a <<= 1;
if(a >= MOD ) a -= MOD;
n>>=1;
}
if(b<0)ans= imod(-ans);
return ans;
}

inline long long ipow(long long a,long long b)
  {
if(b==1)return a%MOD;
long long ans = ipow(a,b>>1);
ans = imul(ans,ans);
if(b&1)ans = imul(ans,a);
return ans;
}

inline long long inv(long long a)
  {
return ipow(a,MOD-2);
}

inline long long gcd(long long a,long long b)
  {
return b?gcd(b,a%b):a;
}

void printA(int n)
  {
for(int i = 1 ; i<= n; i++)
 {
for(int j =1; j<= n+1 ; j++)
printf("%I64d ",A[i][j]);
puts("");
}
}

void getA(int n)
  {
for(int i = 1; i<= n; i++)
 {
A[i][n+1] = 0;
for(int j = 1; j <= n ; j++)
 {
A[i][j] = imod((v[i+1][j] - v[i][j])*2);
A[i][n+1] = imod(A[i][n+1] + imul(v[i+1][j],v[i+1][j]) - imul(v[i][j],v[i][j]));
}
}
//printA(n);
}

void guass(int n)
  {
for(int row = 1; row <= n ; row ++)
 {
int k = row;
long long Max = A[row][row];
for(int i = row+1; i<=n; i++)
if(A[i][row]>Max)Max=A[k=i][row];
if(k!=row)
 {
for(int j = row ; j <= n+1 ; j++)
swap(A[k][j] , A[row][j]);
}
long long div = inv(A[row][row]);
//使得A[row][row] = 1;
for(int j = row + 1 ; j <= n+1 ; j++)
 {
A[row][j] = imul(A[row][j],div);//由于取余,所以可以直接除(乘以逆元),不怕整不整除问题!!
for(int i = row+1; i <= n ; i++)
 {
A[i][j] = imod(A[i][j] - imul(A[row][j],A[i][row]));
}
}
}
for(int row = n ; row >= 1; row --)
 {
for(int j = n ; j> row ; j--)
 {
A[row][n+1] = imod(A[row][n+1] - imul(A[j][j],A[row][j]) );
}
A[row][row] = A[row][n+1];
}
}

void print(int n)
  {
for(int i = 1; i <= n; i++)
 {
if(i>1)putchar(' ');
printf("%I64d",A[i][i]-INF);
}
puts("");
}

int main()
  {
#ifdef ONLINE_JUDGE
#else
freopen("in","r",stdin);
#endif

int T , t = 1;
for(scanf("%d",&T);T--;)
 {
int n;
scanf("%d",&n);
for(int i = 1; i <= n+1; i ++)
for(int j= 1; j <= n; j++)
 {
scanf("%I64d",&v[i][j]);
v[i][j] += INF;//平移一下,这样答案才是正确的
}
getA(n);
guass(n);
printf("Case %d:\n",t++);
print(n);
}
return 0;
}
|
|
常用链接
随笔分类
Links
搜索
最新评论

|
|