集合分割(斜率优化 和 四边形不等式优化)

题目描述
如果 T 是一个整数集合，另 MIN 为 T 中最小的元素，MAX 为 T 中最大的元素，定 义 ????(?) = (??? − ???)2。  现在有一个集合 S，我们要找到 M 个 S的子集??, ??, ?? ⋯ ??，使得?? ∪ ?? ∪ ?? ⋯ ?? = ?，同时，使得 ∑????(??)最小。
 输入格式
 输入第一行两个整数 n(1≤n≤2000)，m(1≤m≤min(n,1000))。
 接下里一行输入 n 个 106以内的整数。
 输出格式
 输出答案。
 样例输入
 4 2
 4 7 10 1
 样例输出
 18

初步分析
首先需要求每个集合的(MAX-MIN)^2的最小值，如果这n个数是乱序的话我们是很难写出状态方程的，因为它没有任何规律。如果一个集合S里面含有x个数。那么这个集合的cost怎样才能更小呢？这n个数肯定满足一个必要条件就是这x个数在(排好序的n个数)里是连续的。倘若任意一个S不满足这个条件。那么肯定有一个数在另外的一个集合里。那么对于这两个集合他们的cost肯定都增大了。（比如s1:1 2 3 s2:7 8 9 10如果两个中任意交换了对他们的cost肯定都增大了）所以这也是他们的充分条件。所以我们先将n个数排好序。再在这个序列上进行dp。

过程分析
设
f[i][j]为前i个数分为j个集合的cost之和对最小值。
a[i]为排好序后的第i个数

对于这个状态方程，我想大家都是能看懂的f[i][j]=min{f[k][j-1]+(a[i]-a[k+1])^2}(j-1=<k<i).从中我们看得出来k的遍历是在f数组的第一维的。所以对于每一个j值我们都得先遍历i再遍历k。由于k的范围在遍历i的时候有许多重复的地方。所以我们考虑是否可以用斜率优化dp

我们设k2>=k1且k2优于k1.那么
f[k2][j-1]+(a[i]-a[k2+1])^2 >= f[k1][j-1]+(a[i]-a[k1+1])^2;
f[k2][j-1]+a[k2+1]^2-2a[i]a[k2+1] >= f[k1][j-1]+a[k1+1]^2-2a[i][k1+1]
设g[i]=f[i][j-1]=a[i+1]^2
a[i]>=(g[k2]-g[k1])/(a[k2+1]-a[k1+1])*1/2

所以我们将a[k+1]当做横坐标，1/2*g[k]当成纵坐标。那么当k1,k2的斜率小于等于a[i]的时候k2就优于k1.并且随着i的增大a[i],k都是递增的。所以我们就可以用一个模拟队列来创建一个关于最优k点的凹图。(关于斜率优化dp这里的原理如果不懂的话可以看看我的另一篇博客，自己总结的一点证明和理解)
传送门：斜率优化DP的(转移，出队)的可行性证明

这样子的话我们就可以来写代码了。需要注意的地方有几点
一 ：每次在遍历一次j的时候都需要将模拟队列重新清空，因为j的值不同前面的队列成员没有任何参考的价值
二 ：因为k的范围是j-1=<k<i，所以我们最开始的时候就将j-1放入队列中，并且k是取不到i的所以我们在遍历凹包更新队首最优点的时候始终是q[head]与q[head+1]比较之后，从尾端维护凹性之后才将此时的i值放进队列中的，这样才能确保k取不到i
三 ：但是在从末尾进行维护的时候我们是把当前的i的这个点放进去一起维护的，实则这个i值还没有进入队列中，并且取最优点是在维护尾端之前，所以与第二点并不冲突

应该差不多了叭，大家结合代码再看看

#include <iostream>
#include <math.h>
#include <algorithm>
using namespace std;
const int MAX_N=2010;
const int MAX_M=1010;
const int inf=0x3f3f3f3f;
int n,m;

int f[MAX_N][MAX_M];
int g[MAX_N];
int a[MAX_N];
int q[MAX_N];
int head=1,tail=1;
int getan(int k,int j,int i){
    return f[k][j-1]+pow(a[i]-a[k+1],2);
}

int compare(int third,int second ,int first,int j){
    g[first]=f[first][j-1]+pow(a[first+1],2);
    g[second]=f[second][j-1]+pow(a[second+1],2);
    g[third]=f[third][j-1]+pow(a[third+1],2);

    return (g[first]-g[second])*(a[second+1]-a[third+1])<=(g[second]-g[third])*(a[first+1]-a[second+1]);
}

void init(){
    for(int i=0;i<=m;i++)
        f[0][i]=inf;

    for(int i=0;i<=n;++i)
        f[i][0]=inf;

    f[0][0]=0;

}
int main() {

    scanf("%d %d",&n,&m);
    for(int i=1;i<=n;++i){
        scanf("%d",&a[i]);
    }
    sort(a+1,a+n+1);

    init();

    for(int j=1;j<=m;++j){
        head=tail=1;
        q[tail++]=j-1;

        for(int i=j;i<=n;++i){
            while(head+1<tail && getan(q[head+1],j,i)<=getan(q[head],j,i)) head++;

            f[i][j]=getan(q[head],j,i);

            cout<<"i j "<<i<<" "<<j<<" "<<f[i][j]<<endl;
            while (head+1<tail && compare(q[tail-2],q[tail-1],i,j)) tail--;
            q[tail++]=i;

        }
    }

    cout<<f[n][m]<<endl;
    return 0;
}

四边形不等式优化
我们知道类似于dp[i][j]={dp[i][k]+dp[k+1][j]+w[i][j]}这种类型的状态方程可以用四边形不等式优化。其实这道题的dp[i][j]=min{dp[i-1][k]+w[k+1][j]}也是可以的。方法是一样的，也是先证明w满足四边形不等式和区间包含不等式然后推出dp也满足四边形不等式，最后的到决策数组s是单调的。说实话，做题的时候很少去证明，都是通过O(n^3)的复杂度去打表观察dp和s数组是否满足。不过再这之前可以自己推一下w数组是否满足。比如这道题
我们知道要证明w数组满足四边形不等式等价于证明f(j)=w[i+1][j]-w[i][j]是关于j的一个单调递减的函数(或者关于i的也行)因为
f(j)=w[i+1][j]-w[i][j]
=(a[j]-a[i+1])^2 - (a[j]-a[i])^2
=(a[i]-a[i+1])(2*a[j]-a[i]-a[i+1])
所以f(j)是一个关于j单调递减的函数。又因为w[i][j’]<w[i’]j这个很容易想到想想w函数的定义。所以w函数也满足区间包含不等式。所以我们可以猜到多半这道题可以用四边形不等式去优化，有耐心和兴趣的朋友也可以自己去推导一下dp和s是否满足。方法和这个博客的方法类似
动态规划加速原理之四边形不等式

注意
其实这些题的思路很值，几乎想到了四边形不等式优化dp就能做出来了。虽然每个算法的考题很多，但都有突破的地方，和那些一定要注意的小地方。希望自己能多做并多汲取教训和总结。
先上代码我在代码中讲哪些地方需要注意

#include <iostream>
#include <algorithm>
#include <math.h>
using namespace std;

const int MAX_N=1010;
const int MAX_M=2020;
const int inf=0x3f3f3f3f;

int n,m;
int dp[MAX_M][MAX_N];
int a[MAX_N];
int s[MAX_M][MAX_N];

void init(){
    for(int i=1;i<=n;++i) s[1][i]=1;
    for(int i=1;i<=n;++i) dp[1][i]=(a[i]-a[1])*(a[i]-a[1]);
}//初始化很重要，特别是开头一般都需要自己初始化

int main() {
    cin>>n>>m;
    for(int i=1;i<=n;++i) scanf("%d",&a[i]);
    sort(a+1,a+1+n);
    init();

    //s[i-1][j]=<s[i][j]<=s[i][j+1]
    for(int i=2;i<=m;++i){//注意这里的i是从第2项开始的，因为dp[1][j]和s[1][j]都是初始化了的。
        s[i][n+1]=n;
        for(int j=n;j>=i;--j){
            int temp;
            dp[i][j]=inf;
     
            for(int k=s[i-1][j];k<=s[i][j+1];++k){
                if(dp[i][j]>dp[i-1][k]+pow((a[j]-a[k+1]),2)){
                    dp[i][j]=dp[i-1][k]+pow((a[j]-a[k+1]),2);
                    temp=k;
                }
            }

            cout<<i<<" "<<j<<" "<<dp[i][j]<<endl;
            s[i][j]=temp;
        }
    }

    cout<<dp[m][n];
    return 0;
}