单纯形模板

#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;
}
		

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值