高斯消元.cpp

一些比较难找的算法代码 中国剩余定理.cpp 高斯消元.cpp 红黑树.cpp 排序后第k位置数.cpp 修正单纯形.c

//小心精度
# include <iostream.h>
# include <string.h>
# include <stdlib.h>
# include <math.h>
# define MAX 150
const double eps=1e-8;
int solve(double a[][MAX],double b[],int n,double x[])
{
        int i,j,k,l,top;
        double temp,lamder;
                int Pi[MAX];
        top=0;
        for(j=0;j<n;j++){
                for(i=top;i<n;i++)if(fabs(a[i][j])>eps) break;
                if(i<n){
                        for(k=0;k<n;k++) {
                                temp=a[i][k];
                                a[i][k]=a[top][k];
                                a[top][k]=temp;
                        }
                        temp=b[top];b[top]=b[i];b[i]=temp;
                        for(k=top+1;k<n;k++)if(fabs(a[k][j])>eps){
                                                                lamder=a[k][j]/a[top][j];
                                for(l=0;l<n;l++) a[k][l]-=lamder*a[top][l];
                                b[k]-=lamder*b[top];
                        }
                        Pi[j]=top++;
                }
                else Pi[j]=-1;
        }
        for(j=n-1;j>=0;j--)if(fabs(b[j])>eps)  break;
        if(j+1>top) return 0;
        for(j=n-1;j>=0;j--)if(Pi[j]!=-1){
                for(i=Pi[j]-1;i>=0;i--)if(fabs(a[i][j])>eps) {
                                                lamder=a[i][j]/a[Pi[j]][j];
                        for(k=0;k<n;k++) a[i][k]-=lamder*a[Pi[j]][k];
                        b[i]-=lamder*b[Pi[j]];
                }
        }
        for(j=0;j<n;j++)
                if(Pi[j]!=-1) x[j]=b[Pi[j]]/a[Pi[j]][j];
                else x[j]=0;
        return 1;
}

int main()
{
        double a[MAX][MAX],b[MAX],x[MAX],p;
        int i,j,n;
        cout.setf(ios::fixed);
        cout.precision(6);
        while(cin>>n>>p){
                for(i=0;i<2*n-1;i++) for(j=0;j<2*n-1;j++) {
                        if(i==j) a[i][j]=-1;
                        else if(i-j==1) a[i][j]=1-p;
                        else if(j-i==1) a[i][j]=p;
                        else a[i][j]=0;
                }
                for(i=0;i<2*n-2;i++) b[i]=0;
                b[i]=-p;
                solve(a,b,2*n-1,x);
                cout<<(x[n-1]+eps)<<endl;
        }
    return 0;
}