线性规划模板

前段时间我参加华为比赛学习的线性规划
虽然现在看来没有卵用
将模板贴一下

#include<cmath>
#include<map>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<set>
#include<vector>
#include<algorithm>
using namespace std;
typedef long long ll;
const int MAXN = 1005;
const double INF = 1e10;
const double eps = 1e-7;

int n,m; // 行， 列
double c[MAXN]; double b[MAXN]; double a[MAXN][MAXN]; double x[MAXN*2]; double v; int N[MAXN], B[MAXN];
int ans;

int cnt = 0;
void debug(int pos , int flor) {
    if(cnt < 10) {
        cnt ++;
    }else return ;
    printf("pos = %d, flor = %d\n",pos,flor);
    printf("%d %d %.2f\n",n,m,v);

    for(int i = 1; i <= m; ++i)
        printf("%.2f ", c[i]); printf("\n");

    for(int i = 1; i <= n; ++i) {
        printf("%.2f ",b[i]);
        for(int j = 1; j <= m; ++j) {
            printf("%.2f ",a[i][j]);
        }
        printf("\n");
    }
    for(int i = 1; i <= n; ++i) printf("%d ",B[i]); printf("\n");
    for(int i = 1; i <= m; ++i) printf("%d ",N[i]); printf("\n");
//  for(int i = 1; i <= n+m; ++i) printf("%.2f ",x[i]); printf("\n\n");
}


int equal(double a, double b) {
    double tt = fabs(b-a);
    if(tt < eps) return 1;
    else return 0;
}
void priov(int l, int e) { // 第l行，第e列
    b[l] /= a[l][e];
    for(int i = 1; i <= m; ++i) {
        if(i != e) {
            a[l][i] /= a[l][e];
        }
    }    
    a[l][e] = 1/a[l][e];
    for(int i = 1; i <= n; ++i) {
        if(i != l && fabs(a[i][e]) > eps) {
            b[i] -= b[l] * a[i][e];
            for(int j = 1; j <= m; ++j) {
                if(j != e) {
                    a[i][j] -= a[i][e] * a[l][j];
                }
            }
            a[i][e] = -a[i][e] * a[l][e];
        }
    }
    v += c[e] * b[l];
    for(int i = 1; i <= m; ++i) {
        if(i != e) {
            c[i] -= c[e] * a[l][i];
        }
    }
    c[e] = -c[e] * a[l][e];
}
double simple(){
    while(1) {
        int l,e = 0;
        for(int i = 1; i <= m; ++i) {
            if(c[i] > eps) {
                e = i; break;
            }   
        }
        if(!e) return v;
        double minn = INF;
        for(int i = 1; i <= n; ++i) {
            if(a[i][e] > eps && minn > b[i]/a[i][e]) {
                minn = b[i]/a[i][e]; l = i; 
            }
        }
        if(minn == INF) return minn;
        priov(l, e);
        swap(B[l], N[e]);       
    }
}
int initsimplex() { // 找到一个初始可行解
    // 备份一些信息
    int _m = m; double _v = v;
    int _N[m+5]; int _c[m+5]; 
    for(int i = 1; i <= m; ++i) _N[i] = N[i];
    for(int i = 1; i <= m; ++i) _c[i] = c[i];

    int id, minn = INF;
    for(int i = 1; i <= n; ++i) 
        if(b[i] < minn) {
            id = i; minn = b[i];
        }

    if(b[id] >= 0) return 1;
    //最初形式不满足条件，就需要寻找一个

    m ++;
    v = 0;
    for(int i = 1; i <= m; ++i)  c[i] = 0; c[m] = -1;
    for(int i = 1; i <= n; ++i) a[i][m] = -1;

    N[m] = 0;
//  debug(2, 1);
    printf("%d %d\n",id,m);
    priov(id, m);
    swap(B[id], N[m]);
    debug(3, 1);
    simple();
    debug(4,1);
    if( equal(v, 0) == 0) return 0;
    //现将生成的新松弛型中的x0项删去
    for(int i = 1; i <= n; ++i) {
        if(B[i] == 0) {
            int pos;
            for(int j = 1; j <= m; ++j) {
                if(a[i][j] != 0) {
                    pos = j;
                }
            }
            priov(i, pos);
            swap(B[i], N[pos]);
            break;
        }
    }
    for(int i = 1; i <= n; ++i) {
        int fl = 0;
        for(int j = 1; j <= m; ++j) {
            if(N[j] == 0) fl = 1;
            if(fl) {
                a[i][j] = a[i][j+1]; 
            }
        }
    }

    int fl = 0;
    for(int i = 1; i <= m; ++i) {
        if(N[i] == 0) fl = 1;
        if(fl) N[i] = N[i+1];
    }
    m --;
    map<int, int> m1; m1.clear();
    map<int, int> m2; m2.clear();

    for(int i = 1; i <= n; ++i) m1[B[i]] = i;
    for(int i = 1; i <= m; ++i) m2[N[i]] = i;
    for(int i = 1; i <= m; ++i) c[i] = 0;
    v = _v;
    printf("%.3f\n",v);
    for(int i = 1; i <= _m; ++i) {
        if(m1.count(_N[i])) {
            int tt = m1[_N[i]];
            for(int j = 1; j <= m; ++j) {
                c[j] -= _c[i] * a[tt][j];
            }
            v += _c[i] * b[tt];
        }else {
            c[m2[_N[i]]] += _c[i];
        } 
    }
    return 1;
}
void solve(){
    while(~scanf("%d %d",&n,&m)) {
        ans = -1;
        scanf("%lf",&v);
        for(int i = 1; i <= m; ++i) scanf("%lf",&c[i]);
        for(int i = 1; i <= n; ++i) {
            scanf("%lf",&b[i]);
            for(int j = 1; j <= m; ++j) {
                scanf("%lf",&a[i][j]);
            }
        }
        for(int i = 1; i <= m; ++i) N[i] = i;
        for(int i = 1; i <= n; ++i) B[i] = i+m;
        set<int> st;
    //  for(int i = 1; i <= m; ++i) st.insert(N[i]);

        int tt = initsimplex();
        for(int i = 1; i <= m; ++i) st.insert(N[i]);
        if(!tt) {
            printf("there is none true answer!\n");
            return;
        }
        debug(1,1);
        simple();
        memset(x, 0, sizeof(x));
        for(int i = 1; i <= n; ++i) {
            if(st.find(B[i]) != st.end()) {
                x[B[i]] = b[i];
            }
        }

        debug(1, 1);
        printf("ans = %d \n", ans);
    }
}
int main(){
    solve();
    return 0;
}

/*
   2 2 
   0
   3 13
   40 2 9
   82 11 -8

   5 2 0.00
   3.00 13.00 
   40.00 2.00 9.00 
   82.00 11.00 -8.00 
   9.00 1.00 0.00 
   2.00 0.00 1.00 
   -9.00 -1.00 0.00

   4 2 0.00
   3.00 13.00 
   40.00 2.00 9.00 
   82.00 11.00 -8.00 
   9.00 1.00 0.00 
   -3.00 0.00 -1.00 
*/