#include <iostream>
#include <cstdio>
#include <cstring>
#include <queue>
#include <cmath>
#include <algorithm>
using namespace std;
#define N 1002
#define M 1020
#define inf 0x3f3f3f3f
#define eps 1e-8
struct SIMPLEX {
// AX <= B
// minimize C^T * X
double A[N][N], b[N], c[N];
int x[N], y[N];
double X[N];
int n, m;
double ans;
int vis[N];
int sim_type;
void init(int n, vector<double> &t) {
this -> n = n;
for(int i = 1; i <= n; ++i) {
c[i] = t[i];
}
m = 0;
}
void pivot(int l, int e) {
b[l] /= A[l][e];
for(int i = 1; i <= n; ++i) {
if(i != e) {
A[l][i] /= A[l][e];
}
}
A[l][e] = 1 / A[l][e];
for(int i = 1; i <= m; ++i) {
if(i != l && fabs(A[i][e]) > eps) {
b[i] -= A[i][e] * b[l];
for(int j = 1; j <= n; ++j) {
if(j != e) {
A[i][j] -= A[i][e] * A[l][j];
}
}
A[i][e] = -A[i][e] * A[l][e];
}
}
ans += c[e] * b[l];
for(int i = 1; i <= n; ++i) {
if(i != e) {
c[i] -= c[e] * A[l][i];
}
}
c[e] = -c[e] * A[l][e];
}
void simplex() {
while(1) {
int e = -1;
for(int i = 1; i <= n; ++i) {
if(c[i] > eps) {
e = i;
break;
}
}
if(e == -1) return;
double tmp = 1e10;
int l = -1;
for(int i = 1; i <= m; ++i) {
if(A[i][e] > eps && b[i] / A[i][e] < tmp) {
tmp = b[i] / A[i][e], l = i;
}
}
if(l == -1) return;
pivot(l, e);
swap(y[l], x[e]);
}
}
void dual() {
int t = max(n, m);
for(int i = 1; i <= t; ++i) {
for(int j = i + 1; j <= t; ++j) {
swap(A[i][j], A[j][i]);
}
}
for(int i = 1; i <= t; ++i) swap(b[i], c[i]);
swap(n, m);
}
void solve() {
ans = 0;
dual();
for(int i = 1; i <= n; ++i) x[i] = i;
for(int i = 1; i <= m; ++i) y[i] = i + n;
simplex();
memset(vis, 0, sizeof vis);
for(int i = 1; i <= n; ++i) vis[x[i]] = i;
for(int i = 1; i <= m; ++i) {
if(vis[i+n]) X[i] = -c[vis[i+n]];
else X[i] = 0;
printf("ans %d %.10lf\n", i, X[i]);
}
}
// tpye = 0, <=
// type = 1, >=
// type = 2, ==
void addLimit(vector<double> &t, double lm, int type) {
if(type != 2) {
if(type == 0) {
for(int i = 1; i <= n; ++i) t[i] = -t[i];
lm = -lm;
}
++m;
for(int i = 1; i <= n; ++i) A[m][i] = t[i];
b[m] = lm;
}
else {
addLimit(t, lm, 0);
addLimit(t, lm, 1);
}
}
}go;
int n, m;
int main() {
printf("input n, m\n");
scanf("%d%d", &n, &m);
vector<double> t;
t.push_back(0);
for(int i = 1; i <= n; ++i) {
double c;
scanf("%lf", &c);
t.push_back(c);
}
go.init(n, t);
for(int i = 1; i <= m; ++i) {
int type;
scanf("%d", &type);
t.clear();
t.push_back(0);
for(int j = 1; j <= n; ++j) {
double c;
scanf("%lf", &c);
t.push_back(c);
}
double v;
scanf("%lf", &v);
go.addLimit(t, v, type);
}
go.solve();
printf("%.10lf\n", go.ans);
return 0;
}
单纯形模板
最新推荐文章于 2018-11-27 16:58:12 发布