一类动态规划的四边形优化

前段时间看到csdn论坛的数据结构与算法版上有人问了这么一个问题：
给定n个石子，其重量为a1,a2...,an,要求将其划分为m堆，每一份划分的费用定义为这堆石头中最大重量与最小重量的差的平方。总划分费用等于各堆费用之和。 输入：n m 及a1,a2...,an ，输出：sum

我搜了一下，发现同样的问题问过不止一遍：
http://topic.csdn.net/u/20091216/00/34ca0d8d-9d56-4ed4-821b-5171e1a67679.html?50332
http://topic.csdn.net/u/20091212/01/bcd82df3-ef2a-4ccd-bd1f-af8588e8e114.html

这道题显然用DP可以解决：
对所有石子的重量从小到大进行排序，设排序后的重量为a[1], a[2]....a[n]。

设f[n,m]为n个石子分成m堆后的最小费用，我们可以列出DP转移方程：
f[j,m] = min_i{f[i,m-1] + (a[j]-a[i])^2}  (i<j)
这样我们得到了一个O(m*n^2)的算法。（这里小小的解释一下为什么要排序：因为每个堆的费用只记录最大值和最小值的平方，中间的值完全没有用处，所以应该把所有中间的值都放到这个堆中，使它们不会增大其他堆的费用。当排完序后，每一堆就是序列中的某一段。）

然而，这样的复杂度并不是很好，考虑到m最大可到n，当n大到1000级别的时候，程序运行将会是很慢，我们需要优化它。事实上，这道题中每一堆的费用是满足一些特殊性的。

设从i到j的石子的费用为w[i,j]=(a[j]-a[i])^2。对于这个w[i,j]，可以验证，当i<=i'<=j<=j'时，有w[i,j']+w[i',j]>=w[i,j]+w[i',j']。这就是四边形不等式。有了这个性质，我们该如何用它来对DP进行优化呢？

我们看到每个f[i,m]的计算是根据f[x,m-1]的值进行运算的，我们可以把转移方程看成对下面这个方程进行了m次操作：E[j]=min_i{D[i]+w[i,j]} (i<j)。对于这个方程，我们可以进行下面的优化：

首先，设T[i,j]=D[i]+w[i,j]，对于每个位置，我们定义B[i]为“最小的j，使得对于所有的i'<i有T[i',j]>=T[i,j]”。我们看一下B数组的意义：B[i]里面存的是一个界值，当我们把T[i,j]向右扩展的时候，只要j越过了这个界值，我们可以肯定的说，i一定比之前的所有i'都更优。

其次，有了B[i]，我们发现，有些决策是没用的。哪些呢？当i<j且B[i]>=B[j]的时候，i就是没用的。为什么？假设我现在在求E[k] (i<j<k)，如果k<B[j]<=B[i]，我们可以说，i和j都不会是决策，因为在i和j之前还有其它决策比i和j更优；如果B[j]<=k<B[i]，我们可以说，j是一个可能的决策，因为j之前没有比j更优的了，但i一定不是决策；如果B[j]<=B[i]<=k，那么可以说i和j都是可能的决策，但是j至少不会比i更差，因为i在j之前。我们最终要维护一个决策序列，在这个序列中i_1<i_2<...<i_d，B[i_1]<B[i_2]<...<B[i_d]。那么在求E[k]的时候，序列里面哪个是它的决策呢？设B[i_x]<=k<B[i_(x+1)]，那么i_x即是决策！这是因为：(1)i_x前面的所有的i都没有i_x更优；(2)i_x和i_(x+1)之间的i都是无用的。

再次，我们该如何寻找决策序列呢？我们可以用一个栈，按照B值从小到大存储(i_x, B[i_x])，对于当前求的k，假设我们已经得到了B[k]，我们可以跟栈顶的B值比较，如果栈为空或者B[k]比它大，则(k, B[k])进栈，否则栈顶元素退栈，继续拿B[k]跟当前的栈顶比较。那么对于k，我们如何找i_x使得B[i_x]<=k<B[i_(x+1)]呢？这里无需二分查找，只需要从k-1时的决策开始找即可。

最后，我们说说该如何求B[k]。这里就要用上四边形不等式了。根据四边形不等式，设i<i'<j<j+1，有w[i,j+1]+w[i',j]>=w[i,j]+w[i',j+1]。
根据此式，我们可以得出，T值也符合四边形不等式，即
        T[i,j+1]+T[i',j]>=T[i,j]+T[i',j+1]
把这个式子换个写法：
        T[i',j]-T[i,j]>=T[i',j+1]-T[i,j+1]
也就是说，当i<i'时，随着j的增长，T[i',j]和T[i,j]之间的差距是不增的。当它们的差值减为0的时候，i'就不会比i差了。那么我们要保证i'之前的i都不会更优该怎么办呢？我们只需要让之前的最优值不会更优即可，就是前面B[i_x]<=k<B[i_(x+1)]里面的i_x。现在有了T[i',j]和T[i,j]差值的单调性，我们就可以通过二分法来找到第一个小于等于0的位置，从而求出B[k]。当然，也有可能找不到差值小于等于0。对这种情况，令B[k]=n+1即可。这样，我们就可以通过O(log(n))的时间求出每个k的B值来。

再看看总的时间复杂度，我们发现对于E[j]=min_i{D[i]+w[i,j]} (i<j)这样一个DP转移方程，用O(n*log(n))的时间即可得到。而对于原问题，我们就得到了一个时间复杂度为O(m*n*log(n))的算法。