GYM 100827 D.Function（高斯消元）

Description

求z(t)
Input
首先输入一整数T表示用例组数，每组用例七个整数n,m,a,b,c,d,r
(1<=T<=100,1<=n<=m<=100,1<=b,d<=r,1<=a,c,r<=1000)
Output
输出z(t)的值，保证有解，要求相对误差和绝对误差均不超过1e-6
Sample Input
3
1 1 1 1 1 1 1
2 3 1 2 3 4 10
2 3 5 6 7 8 9
Sample Output
-1
0.4225806452
4.111111111
Solution

问题转化为求g(n)，那么只要把所有跟求g(n)有关的x找到，然后把这些x值作为变量解方程即可
Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
#define INF 0x3f3f3f3f
#define maxn 222
int T,N,M,A,B,C,D,R;
#define eps 1e-9
double a[maxn][maxn],x[maxn];//方程的左边的矩阵和等式右边的值，求解之后x存的就是结果
int equ,var;//方程数和未知数个数
/*
*返回0表示无解，1表示有解
*/
int Gauss()
{
    int i,j,k,col,max_r;
    for(k=0,col=0;k<equ&&col<var;k++,col++)
    {
        max_r=k;
        for(i=k+1;i<equ;i++)
            if(fabs(a[i][col])>fabs(a[max_r][col]))
                max_r=i;
        if(fabs(a[max_r][col])<eps)return 0;
        if(k!=max_r)
        {
            for(j=col;j<var;j++)
                swap(a[k][j],a[max_r][j]);
            swap(x[k],x[max_r]);
        }
        x[k]/=a[k][col];
        for(j=col+1;j<var;j++)a[k][j]/=a[k][col];
        a[k][col]=1;
        for(i=0;i<equ;i++)
            if(i!=k)
            {
                x[i]-=x[k]*a[i][k];
                for(j=col+1;j<var;j++)a[i][j]-=a[k][j]*a[i][col];
                a[i][col]=0;
            }
    }
    return 1;
}
int flag[maxn];
int main()
{
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d%d%d%d%d%d",&N,&M,&A,&B,&C,&D,&R);
        memset(flag,0,sizeof(flag));
        memset(a,0,sizeof(a));
        memset(x,0,sizeof(x));
        int k=N,res=0;
        while(1)
        {
            flag[M+k]=1;
            if(k<=0)
            {
                a[res][M+k]+=1.0*R*M*M;
                a[res][M-M-2*k]-=1.0*(B*M*M+R*k*k);
                x[res]=1.0*A*M*M;
                k=-M-2*k;
            }
            else
            {
                a[res][M+k]+=1.0*R*M*M;
                a[res][M+M-2*k]-=1.0*(D*M*M+R*k*k);
                x[res]=1.0*C*M*M;
                k=M-2*k;
            }
            res++;
            if(flag[M+k])break;
        }
        for(int i=0;i<=2*M;i++)
            if(!flag[i])
                a[res][i]=1,x[res++]=0;
        equ=var=2*M+1;
        Gauss();
        printf("%.10f\n",x[N+M]);
    }
    return 0;
}