放牧代码和思想
专注自然语言处理、机器学习算法
    愛しさ 優しさ すべて投げ出してもいい

POJ 3689 Equations 题解 《挑战程序设计竞赛》

目录

POJ 3689 Equations

对偶线性规划:n维参数向量满足两个等式,求另一个参数向量的极大值?

3.6与平面和空间打交道的计算几何 

凸包 

首先形式化描述该线性规划(线性规划与对偶问题的一般化详见附录):

原始问题:

max  Z = C1X1 + C2X2 + C3X3 + … + CnXn

st       A1X1 + A2X2 + A3X3 + … + AnXn = Si

         B1X1 + B2X2 + B3X3 + … + BnXn = Ti

         Xi ≥ 0

等效为:

max  Z = C1X1 + C2X2 + C3X3 + … + CnXn

st       A1X1 + A2X2 + A3X3 + … + AnXn ≤ Si

         -A1X1 – A2X2 – A3X3 – … – AnXn  –Si

         B1X1 + B2X2 + B3X3 + … + BnXn  Ti

         -B1X1 – B2X2 – B3X3 – … – BnXn Ti

         Xi ≥ 0

但愿没人问出“将一个不等式左右两边乘以-1的话,不等符号应该反过来,第2个不等式为什么不是≥”这种问题,这么做是因为,a≥b && -a≥-b等效于a≥b && a≤b等效于a==b。

对偶形式:

min  G = SiW1 – SiW2 + TiW3 – TiW4

st       A1W1 – A1W2 + B1W3 – B1W4 ≥ C1

         A2W1 – A2W2 + B2W3 – B2W4 ≥ C2

         …

         AnW1 – AnW2 + BnW3 – BnW4 ≥ Cn

         Wi ≥ 0

令X=W1 – W2 ,Y=W3 – W4,则对偶形式简化为:

min  G = SiX + TiY

st       A1X + B1Y  ≥ C1

         A2X + B2Y ≥ C2

         …

         AnX + BnY ≥ Cn

于是我们得到n条直线,每一条都是一个约束。每一条直线将平面分割为两部分,只有其中一个半平面才是合乎约束的。这n个半平面的交集构成解空间,剩下的就是初中水平的一元二次线性规划了

我是很想这么说,但是写起来真的好麻烦

首先,一个严峻的问题是,半平面交集算法。

看到一种算法是,先按直线的方向向量角度和截距排序,然后维护一个栈,先压入前两条线,从第3条线开始,判断这条线与栈顶下一条线的交点是否满足栈顶的线,如果满足,弹出栈顶,继续判断,直到栈中元素不足两个,插入当前直线,继续加入下一条直线。画个图示意一下:

如图,直线标号表示加入顺序,红色表示弹出,绿色表示不弹出。

然后半平面交集有了后,对输入直线集合也排个序,在交集第一条直线左边的肯定无解,在最后一条直线右边的也无解。夹在中间的话,找角度最相似(最接近平行,也就是叉积最接近0)的那一条边,这条边方向上最远的顶点就是极值点(这里有点像凸包的对踵点算法)。

#include<iostream>
#include <cmath>
#include <algorithm>
#include<vector>
using namespace std;

#define MAX_N 100000 * 2
typedef long long  LL;

inline int compare_to_zero(long double x)
{
	if (fabs(x) < 1e-12)return 0;
	return x < 0 ? -1 : 1;
}

struct Line
{
	LL a, b, c;
} 
ls[MAX_N], // 输入直线
st[MAX_N]; // 输入st对

long double ans[MAX_N];
bool impossible[MAX_N];

// 优先斜率,截距次之的排序
inline bool cmp(const Line &a, const Line &b)
{
	LL t = a.a*b.b - a.b*b.a;
	if (t != 0)return t > 0;
	t = a.c*b.b - a.b*b.c;
	return t > 0;
}
struct Point
{
	long double x, y;
};
int stack[MAX_N];
Point p[MAX_N];

// 直线交点
inline Point intersection(const Line &la, const Line &lb)
{
	LL a = la.a, b = la.b, c = la.c;
	LL d = lb.a, e = lb.b, f = lb.c;
	LL der = a*e - b*d;
	Point ret;
	ret.x = c*e - b*f;
	ret.y = a*f - c*d;
	ret.x /= der;
	ret.y /= der;
	return ret;
}

// 点是否满足直线的约束
inline bool is_point_fit_line(const Line &l, const Point &p)
{
	return compare_to_zero(l.a * p.x + l.b * p.y - l.c) >= 0;
}

///////////////////////////SubMain//////////////////////////////////
int main(int argc, char *argv[])
{
#ifndef ONLINE_JUDGE
	freopen("in.txt", "r", stdin);
	freopen("out.txt", "w", stdout);
#endif
	int n, m;
	scanf("%d%d", &n, &m);
	{
		for (int i = 0; i < n; i++)
		{
			scanf("%d%d%d", &ls[i].a, &ls[i].b, &ls[i].c);
		}
		sort(ls, ls + n, cmp);
		int size = 1;
		for (int i = 1; i < n; i++)
		{
			if (ls[i].a* ls[i - 1].b == ls[i - 1].a* ls[i].b)	continue;	// 排重
			ls[size++] = ls[i];
		}
		n = size;
		LL S, T;
		if (n == 1)
		{
			while (m--)
			{
				scanf("%lld%lld", &S, &T);
				if (S* ls[0].b != T* ls[0].a)	puts("IMPOSSIBLE"); // 平行,没有最小值,也就是最小值是无穷小
				else printf("%.5f\n", ls[0].c*S*1.0 / ls[0].a);
			}
		}
		else
		{
			// 半平面交集
			stack[0] = 0;	// 一个栈,储存的是凸包边所在直线的下标
			p[0].x = 0;		// p是最终的交集对应的点,每个点是每条边的最远点,第一个点当然是由x>=0这个半平面和第一个半平面的交点
			p[0].y = ls[0].c / (double)ls[0].b;
			stack[1] = 1;
			p[1] = intersection(ls[0], ls[1]);	// 第二个点
			int stack_size = 2;
			int ls_index;
			for (ls_index = 2; ls_index < n; ls_index++)
			{
				while (stack_size > 1)
				{
					int prev = stack[stack_size - 2];
					int top = stack[stack_size - 1];
					if (is_point_fit_line(ls[top], intersection(ls[ls_index], ls[prev]))) --stack_size;
					else
					{
						break;
					}
				}
				p[stack_size] = intersection(ls[ls_index], ls[stack[stack_size - 1]]);
				stack[stack_size++] = ls_index;
			}
			// 半平面交集求解完毕

			for (int i = 0; i < m; i++)
			{
				scanf("%d%d", &st[i].a, &st[i].b);
				st[i].c = i;	// 用c表示id
				impossible[i] = false;
			}

			sort(st, st + m, cmp);
			int st_index;
			for (st_index = 0; st_index < m; st_index++)
			{
				if (st[st_index].a * ls[0].b > st[st_index].b * ls[0].a) impossible[st[st_index].c] = true;	// 角度比第一条还小,这条线肯定与半平面交集不相交
				else break;
			}
			int cur = 0;
			while (cur < stack_size && st_index < m)
			{
				ls_index = stack[stack_size - 1];
				if (st[st_index].a * ls[ls_index].b < st[st_index].b * ls[ls_index].a) // 比最后一条还大,还是不相交
				{
					break;
				}
				// 角度在半平面交集的范围内,叉积由小增大
				while (cur < stack_size)
				{
					ls_index = stack[cur];
					LL det = st[st_index].a * ls[ls_index].b - st[st_index].b * ls[ls_index].a;	// 叉积
					if (det >= 0)	// 跟这条约束线的方向最相近,在此方向上取最远点作为极值点
					{
						ans[st[st_index].c] = st[st_index].a * p[cur].x + st[st_index].b * p[cur].y;
						break;
					}
					else
					{
						cur++;
						continue;
					}
				}
				st_index++;
			}
			while (st_index < m)
			{
				impossible[st[st_index].c] = true;
				st_index++;
			}
			for (st_index = 0; st_index < m; st_index++)
			{
				if (impossible[st_index])
				{
					puts("IMPOSSIBLE");
				}
				else
				{
					printf("%.5f\n", ans[st_index]);
				}
			}
		}
	}
#ifndef ONLINE_JUDGE
	fclose(stdin);
	fclose(stdout);
	system("out.txt");
#endif
	return 0;
}
///////////////////////////End Sub//////////////////////////////////
13915434 hankcs 3689 Accepted 2888K 610MS C++ 4057B 2015-02-27 02:51:16

其实半平面算法我还找到了一篇文章讲解得不错:http://www.csie.ntnu.edu.tw/~u91029/Half-planeIntersection.html

Reference

对偶线性规划.pdf

http://hi.baidu.com/chenwenwen0210/item/4133c1cbe26308e4994aa050

知识共享许可协议 知识共享署名-非商业性使用-相同方式共享码农场 » POJ 3689 Equations 题解 《挑战程序设计竞赛》

评论 欢迎留言

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址

我的作品

HanLP自然语言处理包《自然语言处理入门》