/**//* 题意:问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
搜索
最新评论
|
|