题意:问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,很快
由于取余,所以可以直接除(乘以逆元),不怕整不整除问题!!
设距离为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,很快
由于取余,所以可以直接除(乘以逆元),不怕整不整除问题!!
ps: hdu 3571 高斯消元 取余 不能用高精度
求逆元是用了费马小定理,很强悍啊。。。
即:若p是素数,a与p互素,则
a^(p-1)≡1(mod p)。。。。。。
*/
Accepted | 3571 | 390MS | 308K | 2758 B |
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
using namespace std;
const long long mod = 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 ll imod(ll a)
{
return (a%mod+mod)%mod;
}
inline ll iabs(ll a)
{
return a>0?a:-a;
}
inline ll imul(ll a,ll b)
{
ll ans=0,n=iabs(b);
while(n)
{
if(n&1)
{
ans+=a;
if(ans>=mod)这里优化很重要,如果直接写成ans=imod(ans);直接TLE。。。
ans-=mod;
}
a<<=1;
if(a>=mod)
a-=a;
n>>=1;
}
if(b<0)
ans=imod(-ans);
return ans;
}
inline ll ipow(ll a,ll b)
{
if(b==1)return a%mod;
ll ans=ipow(a,b>>1);
ans=imul(ans,ans);
if(b&1)
ans=imul(ans,a);
return ans;
}
inline ll inv(ll a)
{
return ipow(a,mod-2);
}
void printA(int n)
{
for(int i=0;i<n;i++)
{
for(int j=0;j<=n;j++)
printf("%I64d ",a[i][j]);
printf("\n");
}
}
void getA(int n)
{
for(int i=0;i<n;i++)
{
a[i][n]=0;
for(int j=0;j<n;j++)
{
a[i][j]=imod((v[i+1][j]-v[i][j])*2);
a[i][n]=imod(a[i][n]+imul(v[i+1][j],v[i+1][j])-imul(v[i][j],v[i][j]));
//printf("a[%d][j]=%I64d\n",i,a[i][j]);
}
}
//printA(n);
}
void gauss(int n)
{
int k,col,max,j,i,id;
for(k=0;k<n;k++)
{
max=a[k][k];id=k;
for(i=k+1;i<n;i++)
{
if(max<a[i][k])
{
max=a[i][k];
id=i;
}
}
if(id!=k)
{
for(i=k;i<=n;i++)
swap(a[k][i],a[id][i]);
}
ll div=inv(a[k][k]);//使得A[row][row] = 1;利用费马小定理
for(j=k+1;j<=n;j++)
{
a[k][j]=imul(a[k][j],div);//由于取余,所以可以直接除(乘以逆元),不怕整不整除问题!!
for(i=k+1;i<n;i++)
{
a[i][j]=imod(a[i][j]-imul(a[k][j],a[i][k]));
}
}
}
for(i=n-1;i>=0;i--)
{
for(j=n-1;j>i;j--)
{
a[i][n]=imod(a[i][n]-imul(a[i][j],a[j][j]));
}
a[i][i]=a[i][n];
}
}
void print(int n)
{
printf("%I64d",a[0][0]-INF);
for(int i=1;i<n;i++)
{
printf(" %I64d",a[i][i]-INF);
}
printf("\n");
}
int main()
{
int t,cas=1,i,j,n;
scanf("%d",&t);
while(t--)
{
scanf("%d",&n);
for(i=0;i<=n;i++)
for(j=0;j<n;j++)
{
scanf("%I64d",&v[i][j]);
v[i][j]+=INF;//平移一下,这样答案才是正确的
}
getA(n);
gauss(n);
printf("Case %d:\n",cas++);
print(n);
}
return 0;
}