斜率优化DP

引用一片经典的文章：http://blog.sina.com.cn/s/blog_508dd1e60100tvk0.html

  最近发现HDU上的题目ms比POJ上的好些似的， 因为每次都是一道题搞一天。 也说明我真的还是很菜啊。 昨天搞了一道题，一道最短路的题目，弄了一天。其实最主要的原因是没有搞清楚Dijkstra、SPFA、Bellman_ford这几个最短路算法的复杂度。这里作个总结，数组实现的Dijkstra，复杂度为O(N^2)；用优先队列优化的Dijkstra复杂度为O(ElogE)，但是图需要用邻接矩阵实现；Bellman_ford，O(VE)；SPFA,O(KE),K<<V; HDU 3349，给出一张图，题意是可以选择一条边进行减半，其他的不能变，求从S到T的最短路。直接枚举每条边，然后进行Dijkstra,由于题意中的N给了1e5，N^2的算法显然会超时，分析一下可知，超时的原因是每次只枚举了一条减半的边，但是重复求了其他点的最短路，一个优化是将S到每个点的最短路保存在一个数组里，然后求出每个点到T点的最短路。依次枚举每条边，min = MIN(sDis[i] + tDis[j] + cost[i][j]/2 ),<i,j> ;求两次最短路的时候，可以用Heap+Dijkstra, 当然更容易实现的自然是SPFA。 一开始没有想到用SPFA，一直过不了，调了半天，一直WA。 然后就是看到了网上有人用SPFA，换了SPFA，但是没有考虑起点和终点是同一个点的情况，各种WA。 最后msWA了20+次，才AC， 汗。。。。。

       今天做了一道DP的题目，题意如下：

【题目大意】

有N个数，现要将它们分成连续的若干段，每段的代价为(∑Ci)^2+M，求最小的代价。

【题目分析】

容易得到这样的一个动态规划算法：

令dp[i]表示前i个数分成若干段的最小代价，能得到一个经典的动态转移方程： 

dp[i]=  Min  (dp[j] + Cost(j+1,i))+M，0<=j<i。 

其中Cost(a,b)表示把a..b这些数分成一段的费用。

这个算法的时间复杂度是O(N^2)。 这题N有500000，必须优化。

我们看到Cost(j+1,i)。如果预处理出1..i的和，记为s[i]，那么Cost可以表示为(s[i]-s[j])^2。

将它展开并且代入转移方程中可以得到： dp[i]= Min (dp[j] + s[i]^2 -2s[i]*s[j] +s[j]^2)+M。

因为s[i]^2是只和i有关的，可以移到Min的外面，得到

dp[i]= Min (dp[j] -2*s[i]*s[j] +s[j]^2) + M + s[i]^2。

现在我们单独考虑Min()里面的内容。

如果我们设 k=2*s[i],x=s[j],y=dp[j]+s[j]^2的话（注意将和i有关的量以及和j有关的量结合到一起）

再令G=dp[i]，那么将得到：

G=-k*x+y。

移项得：

y=k*x+G。

看到这个式子，知道斜率优化的朋友显然可以做出来了。

这里我讲讲何谓斜率优化。

得到这个式子之后，我们可以看到，k是一个常数（由当前枚举的i在O(1)时间内计算得出）。将这个式子看成是一个直线的函数表达式的话，k就是斜率，也就是说这是一个斜率固定的直线。

y和x则是和j有关的常量。而j的这些值应该都已经在之前计算过了。（因为j<i）

这个式子中，G是未知的，G和y以及x有关。

如果我们把每个j对应的x和y值看成一个坐标系中的点的话。

那么当我们枚举到i时，坐标系中就有一系列的点。

对于每个点，做一条斜率为k的直线，就能得到一个G的值。G的值为这条直线与y轴交点的纵坐标。

我们可以看到，实际上，如果我们让G的值由负无穷变化到正无穷，相当于一条直线，它满足斜率为k，然后从坐标系的下方慢慢地向上平移到坐标系的上方。

那么，我们要找到，G的最小值，就是在这个过程中，这条直线所碰到的第一个点！

而这个点，必然是这个点集的凸包上的点。（不知道凸包概念的去baidu一下好了）

再加上很关键的一点。随着i的增加，点的坐标是单调不下降的（s[i]），直线的斜率也是单调不降的(2s[i])！

满足这两个单调性，我们就可以利用单调队列，来维护一个凸壳。因为我们要找G的最小值，所以要维护一个下凸壳。方法和之前那篇数形结合题一样， 类似Garham求凸包的算法 。

之所以可以用单调队列是因为坐标和斜率都满足单调性的话，可以证明每个点如果不是某个i的最优决策，也不能会是之后的i的最优决策，可以被抛弃。

因为每个i只会进出队列一次，所以时间复杂度降为O(N)，只需要保存单列就可以，空间复杂度为O(N)，比较完美地解决了这个问题。 

附图：

其中绿色部分：类似Garham求凸包的算法：一个一个按逆时针来，加入j1, j2. j3. j4，再加入j5时，向量差乘，发现j5在s射线j3j4的右边，结果为负值，就把j4去掉，继续j6j7j8....

其中背景绿色的部分的简单证明：G=y-kx，对于两个点G1=y1-kx1, G2=y2-kx2，G2在G1点后面出现，也就是说y1<y2, x1<x2，对于一个较小的k尚且有y1-kx1>y2-kx2，即k(x2-x1)>y2-y1，当k增大时，该式当然成立。

而下面代码中之所以能够通过比较相邻的j方式来实现，是因为对于每一个k（或者说没一个i），关于j的函数图像都是 U 字形的。而随着k（或者说i）的增大，这个U字形的极值点也是在不断增大的。（可以参考http://wenku.baidu.com/view/d3d979dcd15abe23482f4d58.html 里面的第一段证明）

附 HDU 4258：Covered Walkway的题目和代码：

Covered Walkway

Time Limit: 30000/10000 MS (Java/Others)    Memory Limit: 131072/131072 K (Java/Others)
Total Submission(s): 350    Accepted Submission(s): 140

Problem Description

Your university wants to build a new walkway, and they want at least part of it to be covered. There are certain points which must be covered. It doesn’t matter if other points along the walkway are covered or not. 
The building contractor has an interesting pricing scheme. To cover the walkway from a point at  x to a point at  y, they will charge  c+( x- y) ², where  c is a constant. Note that it is possible for  x=y. If so, then the contractor would simply charge  c. 
Given the points along the walkway and the constant  c, what is the minimum cost to cover the walkway?

 

Input

There will be several test cases in the input. Each test case will begin with a line with two integers,  n (1≤ n≤1,000,000) and  c (1≤ c≤10 ⁹), where  n is the number of points which must be covered, and  c is the contractor’s constant. Each of the following  n lines will contain a single integer, representing a point along the walkway that must be covered. The points will be in order, from smallest to largest. All of the points will be in the range from 1 to 10 ⁹, inclusive. The input will end with a line with two 0s.

 

Output

For each test case, output a single integer, representing the minimum cost to cover all of the specified points. Output each integer on its own line, with no spaces, and do not print any blank lines between answers. All possible inputs yield answers which will fit in a signed 64-bit integer.

 

Sample Input

  
  

   
   10 5000
1
23
45
67
101
124
560
789
990
1019
0 0
  
  

 

Sample Output

  
  

   
   30726
  
  

 

Source

The University of Chicago Invitational Programming Contest 2012

#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
typedef long long LL;
template <class _Tp> _Tp inline sqr(_Tp a) { return a*a;}

LL GetLL() {
	char ch = ' ';
	while(ch < '0' || ch > '9')
		ch = getchar();
	LL x = 0;
	while (ch >= '0' && ch <= '9')
		x = x * 10 + ch - '0',ch = getchar();
	return x;
}
#define MAXN 1000005

struct point{
	LL x,  y;
	point() {x=y=0;}
	point(LL a, LL b) {x=a; y=b;}
};
LL multi(point o, point a, point b){
	return (a.x-o.x)*(b.y-o.y) - (a.y-o.y)*(b.x-o.x);
}

LL dp[MAXN];
point q[MAXN]; int head, tail;

int n, c;
LL x[MAXN];

int main()
{
	while( scanf("%d%d", &n, &c), n||c ){
		memset(x, 0, sizeof(x));
		memset(dp, 0, sizeof(dp));

		for(int i=0; i<n; i++)
			//scanf("%I64d", &x[i]);
			x[i] = GetLL();      //这两句话意思一样，解释见文章末尾

		head = tail = 0;
		q[tail ++] = point(x[0], sqr(x[0]));
		dp[0] = c;
		for(int i=1; i<n; i++){
			point pp(x[i], dp[i-1] + sqr(x[i]));
			while( head+1<tail && multi( q[tail-2], q[tail-1], pp ) < 0 )   tail --;   //凸包 加点
			q[tail ++] = pp;

			while( head+1<tail && q[head].y - 2*x[i]*q[head].x >= q[head+1].y - 2*x[i]*q[head+1].x )  // y-kx 去点
				head ++;
			dp[i] = q[head].y - 2*x[i]*q[head].x + x[i]*x[i] + c;
		}
		printf("%I64d\n", dp[n-1]);
	}
	return 0;
}

===============================================================================================================

顺便还有个输入输出的优化：

LL GetLL() {
	char ch = ' ';
	while(ch < '0' || ch > '9')
		ch = getchar();
	LL x = 0;
	while (ch >= '0' && ch <= '9')
		x = x * 10 + ch - '0',ch = getchar();
	return x;
}

用GetLL代替scanf("%lld", &)，将2859MS的代码降低到了1968MS，(读写次数与时间复杂度比例1:1)，而后面玩扑克牌的那道题（读写次数与时间复杂度比例1:800），也能用类似的方法将时间从4281MS降低到4031MS。

利用的就是getchar() 的神速度....