首先把板子放在这儿......不过这个板子好像还是有点不太对，后面有时间再改改吧（
简化版：
namespace Simplex //单纯形法，求得目标函数的max值
{
const int N = 150, M = 150; //N是变量个数，M是条件限制个数
const double INF = 0x3f3f3f3f3f3f3f3f;
const double eps = 1e-7; //精度误差
int n, m;
double a[M][N], b[M], c[N], v;
//max cX
//s.t aX<=b
void pivot(int l, int e)
{
b[l] /= a[l][e];
for (int j = 1; j <= n; j++)
if (j != e)
a[l][j] /= a[l][e];
a[l][e] = 1 / a[l][e];
for (int i = 1; i <= m; i++)
if (i != l && fabs(a[i][e]) > 0)
{
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];
}
v += c[e] * b[l];
for (int j = 1; j <= n; j++)
if (j != e)
c[j] -= c[e] * a[l][j];
c[e] = -c[e] * a[l][e];
}
double simplex()
{
while (true)
{
int e = 0, l = 0;
for (e = 1; e <= n; e++)
if (c[e] > eps)
break;
if (e == n + 1)
return v;
double mn = INF;
for (int i = 1; i <= m; i++)
if (a[i][e] > eps && mn > b[i] / a[i][e])
mn = b[i] / a[i][e], l = i;
if (mn == INF)
return INF; // unbounded
pivot(l, e);
}
}
} // namespace Simplex
几乎完全版
namespace Simplex
{
const int N = 150, M = 150; //N是变量个数，M是条件限制个数
const double INF = 3e9;
const double eps = 1e-7; //精度误差
int n, m;
double a[M][N], b[M], c[N], v, A[M][N];
int idx[N], idy[N];
void pivot(const int &e, const int &l)
{
double tmp = -A[l][e];
for (int i = 0; i <= n; i++)
if (i != e)
A[l][i] /= tmp;
A[l][e] = -1 / tmp;
for (int i = 0; i <= m; i++)
{
if (i == l)
continue;
for (int j = 0; j <= n; j++)
{
if (j == e)
continue;
A[i][j] += A[i][e] * A[l][j];
}
A[i][e] *= A[l][e];
}
swap(idx[e], idy[l]);
}
double ans[N];
double simplex()
{
for (int i = 1; i <= n; ++i)
A[0][i] = c[i];
for (int i = 1; i <= m; ++i)
{
for (int j = 1; j <= n; ++j)
A[i][j] = a[i][j];
A[i][0] = b[i];
}
for (int i = 1; i <= n; i++)
idx[i] = i;
for (int i = 1; i <= m; i++)
idy[i] = n + i;
int x, y;
while (1)
{
x = y = 0;
for (int i = 1; i <= m; i++)
if (A[i][0] < -eps)
{
y = i;
if (rand()
break;
}
if (!y)
break;
for (int i = 1; i <= n; i++)
if (A[y][i] > eps)
{
x = i;
if (rand()
break;
}
if (!x)
return -INF;
pivot(x, y);
}
while (1)
{
x = y = 0;
for (int i = 1; i <= n; i++)
if (A[0][i] > eps)
{
x = i;
break;
}
if (!x)
break;
double temp = INF;
for (int i = 1; i <= m; i++)
if (A[i][x] < -eps)
{
double nt = A[i][0] / -A[i][x];
if (nt < temp)
temp = nt, y = i;
}
if (!y)
return INF;
pivot(x, y);
}
// printf(
return A[0][0];
// if (!t)
// return;
// for (int i = 1; i <= m; i++)
// if (idy[i] <= n)
// ans[idy[i]] = A[i][0];
// for (int i = 1; i <= n; i++)
// printf(
}
} // namespace Simplex
怎么用呢？这个板子实际上求解的是：
$$
\max\space\space c^TX \\
s.t\space\space AX\le b
$$
如果要求 $\min$ 的话，转一下对偶即可：
$$
\max\space\space b^TY \\
s.t\space\space A^TY\ge c
$$
总之要注意所有的条件一个都不能漏！
然后接下来开始线性规划的学习吧~
线性规划：将不等式约束通过引入松弛变量和剩余变量转化为等式约束，后者称为标准型
最优解一定在多面体端点上，检查所有端点的复杂度为指数级别。
松弛变量：$\sum_{j=1}^na_{pj}x_j+x_{n+p}=b_p(x_{n+p\ge 0})$
剩余变量：$\sum_{j=1}^na_{pj}x_j-x_{n+p}=b_p(x_{n+p\ge 0})$
自由变量：$x_i=x_i^{+}-x_i^{-}(x_i^{+},x_i^{-}\ge 0, x_i\in R)$
将目标函数求 max 转为求 min
$$
AX=b(b\ge 0) \
A=[B,N],B可逆(A的基向量组)\
基本解：x=(B^{-1}b,0)^T\ge0\space\space(令N=0)
$$
基可行解：既是基本解，又是可行解（可行解：满足方程组变量条件，$X\ge 0$）
不同的基最多有 $\binom{n}{m}$ 个，而一个基最多对应一个基可行解。
顶点（极点）
设 $X^0\in D\space R^n,D是凸集$，$X^0$ 不能表示为 $D$ 中其他任意两不同点的凸组合，则 $X^0$ 为 $D$ 的顶点（极点）。
在 $D^2$ 中，多边形的顶点就是那些端点。半圆的话，圆弧上的点都是顶点。
或者说这个点不在 $D$ 内的任意一条线段上。
设 $X$ 是线性规划的可行解，则 $X$ 为基可行解充要条件为：$X$ 的非零向量在 $A$ 中对应的列向量组线性无关。
一组基最多对应一个顶点。
线性规划、凸优化的两个最优解的凸组合也是最优解（凹优化也是）
单纯形法
从多面体某一个顶点出发，沿着目标函数值下降的方向寻找下一个顶点（沿着端点挨着跳，一边跳一边检查）。
一般来说 $AX=B$ 好解，$AX\le B$ 很困难。
准备工作
最优解判别准则：何时迭代终止，找到最优解
$$
\min\space f(X)=C^TX \
s.t.\space AX=b \
X\ge 0
$$
设 $A=[I\space\space N]$ 转化为
$$
\min\space f(X)=C_I^TX_I+C_N^TX_N \
s.t.\space [I\space\space N][X_I\space\space X_N]^T=b \
X_I,X_N\ge 0 \
I\space 是单位矩阵
$$
$\bar{x}=\binom{b}{0}\ge 0\space\space(基本可行解，即顶点),\space 则\space f(\bar{x})=C_{I}^Tb\space 为最小值的充要条件为\space C_I^TN-C_N^T\le 0$
证明：$f(x)=C_I^TX_I+C_N^TX_N=C_I^T(b-NX_N)+C_N^TX_N=f(\bar{x})-[C_I^TN-C_N^T]X_N\ge f(\bar{x}) \space (C_I^TN-C_N^T\le 0)$
通常叫 $\space C_I^TN-C_N^T$ 为判别式。
若判别式中某个判别数 $\sigma_j>0$ 但相应的列向量 $P_j\le 0$，则该线性规划无最优解（即没有下界）。
换基运算：从一个顶点转到另一个顶点
利用不交换行的行初等变换，选定其中一个位置（该位置不等于 $0$），该位置变为 $1$，这一列其他位置变为 $0$.
选择 $\sigma_j$ 最大的一列来换，这样函数值下降最多。$b_k'=\frac{b_k}{a_{kl}},\space b_i'=b_i-b_k'a_{il}$，故应选择这一列的正数 $a_{kl}$ 且使得 $\frac{b_k}{a_{kl}}$ 最小。$\frac{b_k}{a_{kl}}$ 称为转轴元素。这一列称为进基列。
单纯形表
用分块矩阵的运算来解释，进行行初等变换即可。
$$
\begin{bmatrix}
A&b \\
-C^T&0
\end{bmatrix}
\to
\begin{bmatrix}
I&N &b \\
-C_I^T&-C_N^T &0
\end{bmatrix}
\to
\begin{bmatrix}
I&N &b \\
0&C_I^TN-C_N^T &C_I^Tb
\end{bmatrix}
$$
最后一个称为初始单纯形表。如下：
$$
\begin{bmatrix}
P_1& P_2&...&P_m&P_{m+1}&...&P_l&...&P_n&b \\
1 & & & &a_{1,m+1}&...&a_{1l}&...&a_{1n}&b_1 \\
...\\
0&0&...&0&\sigma_{m+1}&...&0&...&\sigma_n&f
\end{bmatrix}
$$
两阶段方法
如果 $A=[B\space\space N]$，即不显含单位矩阵 $I$，则构造人工变量 $Y=(y_1,y_2,...,y_m)^T$，有
$$
\min\space\space \sum y_i \\
s.t\space\space Y+AX=b
$$
注意到含 $y_i$ 的系数都是 $1$，此时有一个显含的单位矩阵 $I$，进行 $m$ 次行初等变换之后就变成一个初始单纯形表。
上述问题为辅助线性规划，求得最优解。
若 $\sum Y>0$，则原问题无可行解。
证明：若有可行解，则 $\binom{0}{\bar{x}}$ 一定是辅助线性规划的一个特解（因为这是原问题的一个特解）。此时一定有 $\sum Y=0$。 证毕。
若 $\sum Y=0$，则 $y_1=y_2=...=0$。这是我们希望的。此时约束条件就是原问题的约束条件。
转轴之后单位矩阵的列会转到所有 $x_i$ 列上，此时就找到一个极大无关组，即 $B$。
这种做法比直接行初等变换找极大无关组速度更快。
对偶问题
$$
\max\space\space c^TX \\
s.t\space\space AX\le b
$$
这类问题可以用对偶原则转化为：
$$
\min\space\space b^TY \\
s.t.\space\space A^TY\ge c
$$
证明：写出拉格朗日函数：$L(x,\lambda,y)=c^TX+\lambda^T(AX-b)-\mu^TX $
$q(\lambda,\mu)=\min\limits_{X\in \mathcal{R^{n}}} L(x,\lambda,\mu)=\min\limits_{X\in \mathcal{R^n}}{(c^T-\lambda^TA-\mu^T)X}+\lambda^Tb=\lambda^Tb+ \min\limits_{X\in \mathcal{R^n}}(c-A^T\lambda-\mu)^TX=\lambda^Tb\space\space(c-A^T\lambda-\mu=0) $
故对偶问题可以表示为：
$$
\max\space\space \lambda^Tb \\
s.t\space\space c-A^T\lambda \ge0
$$
至于线性规划如何转化成网络流问题，这个坑先留在这里吧（
下面是一些题目
方格取数问题
这个题用费用流解蛮经典的了，但是实际上用线性规划来表示的话也很直观：
$$
\max\space\space \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} c_{ij}x_{ij} \\
s.t\space\space x_{ij}+x_{i(j+1)}\le 1 \\
x_{ij}+x_{(i+1)j}\le 1
$$
一共 $n\cdot m$ 个变量，$(n-1)(m-1)$ 个约束，实际上直接求解即可......可能是因为单纯形法的复杂度一直都很玄学吧（
#include <bits/stdc++.h>
#define ll long long
#define ls id << 1
#define rs id << 1 | 1
#define mem(array, value, size, type) memset(array, value, ((size) + 5) * sizeof(type))
#define memarray(array, value) memset(array, value, sizeof(array))
#define fillarray(array, value, begin, end) fill((array) + (begin), (array) + (end) + 1, value)
#define fillvector(v, value) fill((v).begin(), (v).end(), value)
#define pb(x) push_back(x)
#define st(x) (1LL << (x))
#define pii pair<int, int>
#define mp(a, b) make_pair((a), (b))
#define Flush fflush(stdout)
#define vecfirst (*vec.begin())
#define veclast (*vec.rbegin())
#define vecall(v) (v).begin(), (v).end()
#define vecupsort(v) (sort((v).begin(), (v).end()))
#define vecdownsort(v, type) (sort(vecall(v), greater<type>()))
#define veccmpsort(v, cmp) (sort(vecall(v), cmp))
using namespace std;
const int N = 150;
const int inf = 0x3f3f3f3f;
const ll llinf = 0x3f3f3f3f3f3f3f3f;
const int mod = 998244353;
const int MOD = 1e9 + 7;
const double PI = acos(-1.0);
clock_t TIME__START, TIME__END;
void program_end()
{
#ifdef ONLINE
printf("\n\nTime used:
system("pause");
#endif
}
namespace Simplex //单纯形法，求得目标函数的max值
{
const int N = 10001, M = 9900; //N是变量个数，M是条件限制个数
const double INF = 2e9;
const double eps = 1e-9; //精度误差
int n, m;
double a[M][N], b[M], c[N], v;
//max cX
//s.t aX<=b
void pivot(int l, int e)
{
b[l] /= a[l][e];
for (int j = 1; j <= n; j++)
if (j != e)
a[l][j] /= a[l][e];
a[l][e] = 1 / a[l][e];
for (int i = 1; i <= m; i++)
if (i != l && fabs(a[i][e]) > 0)
{
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];
}
v += c[e] * b[l];
for (int j = 1; j <= n; j++)
if (j != e)
c[j] -= c[e] * a[l][j];
c[e] = -c[e] * a[l][e];
}
double simplex()
{
while (true)
{
int e = 0, l = 0;
for (e = 1; e <= n; e++)
if (c[e] > eps)
break;
if (e == n + 1)
return v;
double mn = INF;
for (int i = 1; i <= m; i++)
if (a[i][e] > eps && mn * a[i][e] > b[i])
mn = b[i] / a[i][e], l = i;
if (mn == INF)
return INF; // unbounded
pivot(l, e);
}
}
} // namespace Simplex
int n, m, a[N][N];
inline void solve()
{
scanf(
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= m; ++j)
scanf(
Simplex ::n = n * m;
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= m; ++j)
Simplex::c[(i - 1) * m + j] = a[i][j];
int k = 0;
for (int i = 1; i <= n; ++i)
{
for (int j = 1; j < m; ++j)
{
k++;
Simplex::b[k] = 1;
Simplex::a[k][(i - 1) * m + j] = 1;
Simplex::a[k][(i - 1) * m + j + 1] = 1;
}
}
for (int j = 1; j <= m; ++j)
{
for (int i = 1; i < n; ++i)
{
k++;
Simplex::b[k] = 1;
Simplex::a[k][(i - 1) * m + j] = 1;
Simplex::a[k][i * m + j] = 1;
}
}
Simplex::m = k;
int ans = Simplex::simplex() + 0.5;
printf(
}
int main()
{
TIME__START = clock();
int Test = 1;
// scanf(
while (Test--)
{
solve();
// if (Test)
// putchar('\n');
}
TIME__END = clock();
program_end();
return 0;
}
Codeforces 1061E - Politics
这题时限 4s，我愣是 3650ms 卡过去了呃......
反正线性规划式子也很好推，注意要限制每一个点最多被选一次，还有判断无解的情况。
然后直接按照题意来设限制即可，复杂度玄学（
#include <bits/stdc++.h>
#define ll long long
#define ls id << 1
#define rs id << 1 | 1
#define mem(array, value, size, type) memset(array, value, ((size) + 5) * sizeof(type))
#define memarray(array, value) memset(array, value, sizeof(array))
#define fillarray(array, value, begin, end) fill((array) + (begin), (array) + (end) + 1, value)
#define fillvector(v, value) fill((v).begin(), (v).end(), value)
#define pb(x) push_back(x)
#define st(x) (1LL << (x))
#define pii pair<int, int>
#define mp(a, b) make_pair((a), (b))
#define Flush fflush(stdout)
#define vecfirst (*vec.begin())
#define veclast (*vec.rbegin())
#define vecall(v) (v).begin(), (v).end()
#define vecupsort(v) (sort((v).begin(), (v).end()))
#define vecdownsort(v, type) (sort(vecall(v), greater<type>()))
#define veccmpsort(v, cmp) (sort(vecall(v), cmp))
using namespace std;
const int N = 505;
const int inf = 0x3f3f3f3f;
const ll llinf = 0x3f3f3f3f3f3f3f3f;
const int mod = 998244353;
const int MOD = 1e9 + 7;
const double PI = acos(-1.0);
clock_t TIME__START, TIME__END;
void program_end()
{
#ifdef ONLINE
printf("\n\nTime used:
system("pause");
#endif
}
namespace Simplex
{
const int N = 505, M = 2550; //N是变量个数，M是条件限制个数
const double INF = 3e9;
const double eps = 1e-6; //精度误差
int n, m;
double a[M][N], b[M], c[N], v, A[M][N];
int idx[N], idy[M];
void pivot(const int &e, const int &l)
{
double tmp = -A[l][e];
for (int i = 0; i <= n; i++)
if (i != e)
A[l][i] /= tmp;
A[l][e] = -1 / tmp;
for (int i = 0; i <= m; i++)
{
if (i == l)
continue;
for (int j = 0; j <= n; j++)
{
if (j == e)
continue;
A[i][j] += A[i][e] * A[l][j];
}
A[i][e] *= A[l][e];
}
swap(idx[e], idy[l]);
}
double ans[N];
double simplex()
{
srand(1);
for (int i = 1; i <= n; ++i)
A[0][i] = c[i];
for (int i = 1; i <= m; ++i)
{
for (int j = 1; j <= n; ++j)
A[i][j] = -a[i][j];
A[i][0] = b[i];
}
for (int i = 1; i <= n; i++)
idx[i] = i;
for (int i = 1; i <= m; i++)
idy[i] = n + i;
int x, y;
while (1)
{
x = y = 0;
for (int i = 1; i <= m; i++)
if (A[i][0] < -eps)
{
y = i;
if (rand()
break;
}
if (!y)
break;
for (int i = 1; i <= n; i++)
if (A[y][i] > eps)
{
x = i;
if (rand()
break;
}
if (!x)
return -INF;
pivot(x, y);
}
while (1)
{
x = y = 0;
for (int i = 1; i <= n; i++)
if (A[0][i] > eps)
{
x = i;
break;
}
if (!x)
break;
double temp = INF;
for (int i = 1; i <= m; i++)
if (A[i][x] < -eps)
{
double nt = A[i][0] / -A[i][x];
if (nt < temp)
temp = nt, y = i;
}
if (!y)
return INF;
pivot(x, y);
}
return A[0][0];
}
} // namespace Simplex
int n, rx, ry;
vector<int> e1[N], e2[N];
int a[N];
int depx[N], depy[N];
void dfsdepx(int u, int fa)
{
depx[u] = depx[fa] + 1;
for (auto v : e1[u])
{
if (v == fa)
continue;
dfsdepx(v, u);
}
}
void dfsdepy(int u, int fa)
{
depy[u] = depy[fa] + 1;
for (auto v : e2[u])
{
if (v == fa)
continue;
dfsdepy(v, u);
}
}
void dfsx(int u, int fa, int k, int f)
{
Simplex::a[k][u] = f;
for (auto v : e1[u])
{
if (v == fa || depx[v] < depx[u])
continue;
dfsx(v, u, k, f);
}
}
void dfsy(int u, int fa, int k, int f)
{
Simplex::a[k][u] = f;
for (auto v : e2[u])
{
if (v == fa || depy[v] < depy[u])
continue;
dfsy(v, u, k, f);
}
}
inline void solve()
{
scanf(
Simplex::n = n;
for (int i = 1; i <= n; ++i)
{
scanf(
Simplex::c[i] = a[i];
}
for (int i = 1; i < n; ++i)
{
int x, y;
scanf(
e1[x].push_back(y);
e1[y].push_back(x);
}
for (int i = 1; i < n; ++i)
{
int x, y;
scanf(
e2[x].push_back(y);
e2[y].push_back(x);
}
dfsdepx(rx, 0);
dfsdepy(ry, 0);
int K = 0;
for (int i = 1; i <= n; ++i)
{
K++;
Simplex::b[K] = 1;
Simplex::a[K][i] = 1;
}
int qx, qy;
scanf(
while (qx--)
{
int k, x;
scanf(
K++;
Simplex::b[K] = x;
dfsx(k, 0, K, 1);
K++;
Simplex::b[K] = -x;
dfsx(k, 0, K, -1);
}
scanf(
while (qy--)
{
int k, x;
scanf(
K++;
Simplex::b[K] = x;
dfsy(k, 0, K, 1);
K++;
Simplex::b[K] = -x;
dfsy(k, 0, K, -1);
}
Simplex::m = K;
double ans = Simplex::simplex();
if (ans == -Simplex::INF)
return puts("-1"), void();
// printf(
int ans2 = ans + 0.5;
cout << ans2 << endl;
}
int main()
{
TIME__START = clock();
int Test = 1;
// scanf(
while (Test--)
{
solve();
// if (Test)
// putchar('\n');
}
TIME__END = clock();
program_end();
return 0;
}
Comments NOTHING