2014-team3/code/lp

从 Trac 迁移的文章

这是从旧校内 Wiki 迁移的文章，可能存在一些样式问题，您可以向 memset0 反馈。

原文章内容如下：

{{{
#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstring>
using namespace std;

typedef double dbl;
const dbl inf = 1e99,eps = 1e-8;
struct LP{
    static const size_t N=20,M=20,L=N+M*2,faux=~0; 
    dbl A[M][L],b[M],c[N+M],c1[L],x[N],d[L],zm;
    size_t bs[M],n,m,l;
#define go(i,n) for(i=0;i<n;++i) 
    dbl simplex(dbl *c){
        for(size_t j,k,i,o;;bs[o]=i){
            dbl dmax=0,u=inf,t;
            i=o=-1;
            go(k,l){
                d[k]=c[k];
                go(j,m)d[k]-=c[bs[j]]*A[j][k];
                if(d[k]>dmax+eps){dmax=d[k];i=k;}
            }
            if(!~i){
                zm=0;fill(x,x+n,0);
                go(j,m)if(bs[j]<n)
                    zm+=c[bs[j]]*(x[bs[j]]=b[j]);
                return zm;
            }
            go(j,m)if(A[j][i]>eps&&b[j]/A[j][i]+eps<u)
                u=b[o=j]/A[j][i];
            if(!~o)throw "unbounded!";
            b[o]*=t=1/A[o][i];
            go(k,l)A[o][k]*=t;
            go(j,m)if(j!=o&&fabs(A[j][i])>eps){
                t=-A[j][i];b[j]+=t*b[o];
                go(k,l)A[j][k]+=t*A[o][k];
            }
        }
    }
    void init(int _n){
        memset(A,0,sizeof(A));
        memset(b,0,sizeof(b));
        memset(c,0,sizeof(c));
        l=n=_n;m=0;
    }
    void add_contraint(dbl* a,dbl _b,bool eq=false){
        copy(a,a+n,A[m]);b[m]=_b;
        if(eq||_b<0)bs[m]=faux;
        else A[m][bs[m]=l++]=1;
        m++;
    }
    dbl solve(){
        size_t ol=l,j,k;
        go(j,m)if(bs[j]==faux)A[j][bs[j]=l++]=1;
        if(ol!=l){
            go(k,l)c1[k]=k<ol?0:-1;
            if(simplex(c1)<-eps)throw "infeasible";
            l=ol;
        }
        return simplex(c);
    }
};
LP a;
int main(){
        int n,m,i,j;
        while(cin >> n >> m){
                a.init(n);
                dbl t[100];
                go(i,n)cin >> a.c[i];
                go(j,m){
                        go(i,n)cin >> t[i];
                        dbl b;cin >> b;
                        a.add_contraint(t,b);
                }
                long long ans = a.solve()*m+0.9999999;
                cout << "Nasa can spend "<< ans<<" taka."<< endl;
        }
}
}}}

#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstring>
using namespace std;
typedef double dbl;
const dbl inf = 1e99,eps = 1e-8;
struct LP{
    static const size_t N=20,M=20,L=N+M*2,faux=~0; 
    dbl A[M][L],b[M],c[N+M],c1[L],x[N],d[L],zm;
    size_t bs[M],n,m,l;
#define go(i,n) for(i=0;i<n;++i) 
    dbl simplex(dbl *c){
        for(size_t j,k,i,o;;bs[o]=i){
            dbl dmax=0,u=inf,t;
            i=o=-1;
            go(k,l){
                d[k]=c[k];
                go(j,m)d[k]-=c[bs[j]]*A[j][k];
                if(d[k]>dmax+eps){dmax=d[k];i=k;}
            }
            if(!~i){
                zm=0;fill(x,x+n,0);
                go(j,m)if(bs[j]<n)
                    zm+=c[bs[j]]*(x[bs[j]]=b[j]);
                return zm;
            }
            go(j,m)if(A[j][i]>eps&&b[j]/A[j][i]+eps<u)
                u=b[o=j]/A[j][i];
            if(!~o)throw "unbounded!";
            b[o]*=t=1/A[o][i];
            go(k,l)A[o][k]*=t;
            go(j,m)if(j!=o&&fabs(A[j][i])>eps){
                t=-A[j][i];b[j]+=t*b[o];
                go(k,l)A[j][k]+=t*A[o][k];
            }
        }
    }
    void init(int _n){
        memset(A,0,sizeof(A));
        memset(b,0,sizeof(b));
        memset(c,0,sizeof(c));
        l=n=_n;m=0;
    }
    void add_contraint(dbl* a,dbl _b,bool eq=false){
        copy(a,a+n,A[m]);b[m]=_b;
        if(eq||_b<0)bs[m]=faux;
        else A[m][bs[m]=l++]=1;
        m++;
    }
    dbl solve(){
        size_t ol=l,j,k;
        go(j,m)if(bs[j]==faux)A[j][bs[j]=l++]=1;
        if(ol!=l){
            go(k,l)c1[k]=k<ol?0:-1;
            if(simplex(c1)<-eps)throw "infeasible";
            l=ol;
        }
        return simplex(c);
    }
};
LP a;
int main(){
        int n,m,i,j;
        while(cin >> n >> m){
                a.init(n);
                dbl t[100];
                go(i,n)cin >> a.c[i];
                go(j,m){
                        go(i,n)cin >> t[i];
                        dbl b;cin >> b;
                        a.add_contraint(t,b);
                }
                long long ans = a.solve()*m+0.9999999;
                cout << "Nasa can spend "<< ans<<" taka."<< endl;
        }
}