hdu4870 高斯消元

#include <iostream>
#include <algorithm>
#include <string.h>
#include <cstdio>
#include <cmath>
using namespace std;
const int maxn=25;
const double eps=0.000000001;
int sgn(double f){
  if(fabs(f)<eps)return 0;
  return f>0?1:-1;
}
int id[maxn][maxn],cnt;
int equ,var;//方程数和未知数个数
double P;
double a[250][250],x[250];
void perId()
{
    cnt=0;
    for(int i=0; i<20; i++)
        for(int j=0; j<=i; j++)
            id[i][j]=cnt++;
    id[20][19]=cnt++;
    equ=var=cnt;
}
void init()
{
    memset(a,0,sizeof(a));
    for(int i=0; i<20; i++)
        for(int j=0; j<=i ;j ++)
         {
              int u=id[i][j];
             x[ u ]=a[ u ][ u ] = 1.0;
               int nx,ny;
               nx=max( i , j+1 ),ny=min( i , j+1 );
               a[ u ][ id[nx][ny] ]-=P;
               nx=i; ny=max(j-2,0);
               a[ u ][ id[nx][ny] ]-=(1-P);
         }
    x[cnt-1]=0.0;
    a[cnt-1][cnt-1]=1.0;
}

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(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][col];
              for(j=col+1;j<var;j++)a[i][j]-=a[k][j]*a[i][col];
              a[i][col]=0;
          }
    }
    return 1;
}
int main()
{
    perId();

    while(scanf("%lf",&P)==1)
    {
        init();
        Gauss();
        printf("%.6lf\n",x[0]);
    }
    return 0;
}