1. 程式人生 > >BZOJ 1013: [JSOI2008]球形空間產生器sphere(高斯消元)

BZOJ 1013: [JSOI2008]球形空間產生器sphere(高斯消元)

題目:https://www.lydsy.com/JudgeOnline/problem.php?id=1013

思路:
存在二次項,考慮兩式相減可以把所有未知數的二次項消掉,
n+1 個等式用第一個與後面的做差,形成n個不等式,然後
高斯消元即可。

程式碼:

#include<cstdio>
#include<cmath>
#include<algorithm>
using std::swap;
#define MAXN 15
#define EPS 1e-8
typedef double Matrix[MAXN][MAXN];//係數矩陣
int Rank;//有效方程的行數,Rank之後的方程x係數為0
double X[MAXN];//解
bool Free[MAXN];//是否為自由變數
double aa[MAXN],bb[MAXN];
Matrix A;//係數矩陣
//浮點數比較
int fcmp(double a,double b)
{
    if((a-b)>EPS)
        return 1;
    else if((a-b)>-EPS)
        return 0;
    return -1;
}
//高斯·約當消元
void Gauss(int n,int m)
{
    int r,c,mxr;
    for(r=1,c=1;r<=n&&c<m;r++,c++)
    {
        //尋找絕對值最大(選主)
        mxr=r;
        for(int i=r+1;i<=n;i++)
            if(fcmp(fabs(A[i][c]),fabs(A[mxr][c]))>0)
                mxr=i;
        //第c項在第r個方程之後係數都為0
        if(fcmp(A[mxr][c],0.0)==0)
        {r--;continue;}//執行下一項
        //選好主,交換方程
        if(mxr!=r)swap(A[mxr],A[r]);
        //消元
        for(int i=1;i<=n;i++)
            if(i!=r&&fcmp(A[i][c],0.0)!=0)
                for(int j=m;j>=c;j--)
                    A[i][j]-=A[r][j]/A[r][c]*A[i][c];
    }
    Rank=r-1;
}
//判斷是否有解
bool check(int n,int m)
{
    //判斷是否無解
    for(int i=Rank+1;i<=n;i++)
        if(fcmp(A[i][m],0.0)!=0)//在x係數為0時,結果不為0
            return 0;
    //初始所有都是未知變數
    for(int i=1;i<m;i++)
        Free[i]=1;
    //計算解
    for(int i=Rank,cnt=0,pos;i>0;i--)
    {
        cnt=0;//記數當前方程未知變數個數
        for(int j=1;j<m;j++)
            if(Free[j]&&fcmp(A[i][j],0.0)!=0)
                cnt++,pos=j;
        if(cnt==1)//一個未知變數,可計算
        {
            Free[pos]=0;
            X[pos]=A[i][m]/A[i][pos];
        }
    }
    return 1;
}
int main()
{
    int n;scanf("%d",&n);
    for(int i=1;i<=n;i++) scanf("%lf",&aa[i]);
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            scanf("%lf",&bb[j]);
            A[i][j]=2*(bb[j]-aa[j]);
            A[i][n+1]+=(bb[j]*bb[j]-aa[j]*aa[j]);
        }
    }
    Gauss(n,n+1);
    check(n,n+1);
    for(int i=1;i<=n;i++)
        printf("%.3lf ",X[i]);
    return 0;
}