[CEOI2004]锯木厂选址 斜率优化DP

斜率优化DP

先考虑朴素DP方程，

f[i][k]代表第k个厂建在i棵树那里的最小代价，最后答案为f[n+1][3];

f[i][k]=min(f[j][k-1] + 把j+1~i的树都运到i的代价）

首先注意到“把j+1~i的树都运到i的代价”不太方便表达，每次都暴力计算显然是无法承受的，

于是考虑前缀和优化，观察到先运到下一棵树那里，等一会再运下去，和直接运下去是等效的。

设sum[i]代表1 ~ i的树都运到i的代价，

于是根据前缀和思想，猜想我们可以用1 ~ r 的代价与 1 ~ l-1的代价获取l ~ r的代价,

所以要做的就是吧1 ~ l-1 对 1 ~ r产生的贡献给算出来，然后减掉，

考虑先把1 ~ l-1的树都运到l-1,所以这部分的代价是sum[l-1],

然后再把树一次性运到r，那么代价是sum_weight[l-1] * (sum_len[r] - sum_len[l-1]);

                                    总的重量 * 现在要再次运的路程

这里为了表示方便，用$sw$代表sum_weight,用$sl$代表sum_len;

于是用$sum[r]$ 减去这两部分代价就可以得到$l ~ r$ 的代价(把$l ~ r$的树都运到$r$)

代价(l ~ r)$ =  sum[r] - sum[l-1] - sw[l-1] * (sl[r] - sl[l-1]);$

那么如何计算sum ?

也是一样的思想，用前面的推后面的，先得到前面的代价，再加上新增的代价即可

$sum[i]=sum[i-1] + swt[i-1] * len[i-1];$//len[i-1]代表i-1到i的距离

于是我们就得到了DP方程：

当$k==1$时,$f[i][k]=sum[i]$;

else

        $f[i][k]=min(f[j][k-1] + sum[i] - sum[j] - sw[j] * (sl[i] - sl[j]));$

但是可以发现，由于k最大就是3，而且3必须是n+1才可以取，

而且当$k==1$时，$f[i][k]$就等于$sum[i]$，

所以考虑优化维数：

当$k==1$时，不用求，因为有$sum$了

当$k==2$时，调用的$f[j][k-1]$替换为$sum[j]$,并且还可以发现由于后面有一个$-sum[j]$,所以可以直接消掉

当$k==3$时，由于只有$n+1$可以取，所以直接在外面多写一个循环，相当于最后统计答案即可

转移方式同朴素方程

但是这样是$n^2$的DP，而$n$有20000，那怎么办呢？

考虑斜率优化。

首先我们用暴力打表可以发现，决策是单调的，

打表代码（朴素DP）：

#include<bits/stdc++.h>
using namespace std;
#define R register int
#define AC 20100
int n, ans;
int sum[AC], sum_weight[AC], sum_len[AC], f[AC];
int weight[AC], len[AC];
inline int read()
{
    int x = 0; char c = getchar();
    while(c < '0' || c > '9') c = getchar();
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x; 
}

void pre()
{
    n = read();
    for(R i = 1; i <= n; i ++)
        weight[i] = read(), len[i] = read();
}

void getsum()
{
    for(R i = 1; i <= n + 1; i ++)//山脚的也要求
    {
        sum_len[i] = sum_len[i - 1] + len[i - 1];
        sum_weight[i] = sum_weight[i - 1] + weight[i];
        sum[i] = sum[i - 1] + sum_weight[i - 1] * len[i - 1];
    //    printf("%d : %d
",i,sum[i]);
    }
}

void work()
{
    for(R i = 1; i <= n; i ++)
    {
        int tmp = 0;
        f[i] = INT_MAX;
        for(R j = 1;j < i;j ++)
        {
            if(sum[i] - sum_weight[j] * (sum_len[i] - sum_len[j]) < f[i])
            {
                f[i] = sum[i] - sum_weight[j] * (sum_len[i] - sum_len[j]);
                tmp = j;
            } 
        }
        printf("%d --- > %d
", tmp, i);//打表验证决策单调性
    }
    ans = INT_MAX;
    for(R i = 2; i <= n; i ++)//注意应该是n+1,因为山脚是在下面
        ans = min(ans, f[i] + sum[n + 1] - sum[i] - sum_weight[i] * (sum_len[n + 1] - sum_len[i]));
    for(R i = 2; i <= n; i ++) printf("%d : %d
", i, f[i]);
    printf("%d
", ans);
}

int main()
{
    freopen("in.in", "r", stdin);
    freopen("out.out", "w", stdout);
    pre();
    getsum();
    work();
    fclose(stdin);
    fclose(stdout);
    return 0;
}

于是我们推斜率优化方程:

设有 $k < j < i$,且$j$优于$k$（相当于$j$是后面来的）,则有：

$sum[i] - sw[j] * (sl[i] - sl[j]) < sum[i] - sw[k] * (sl[i] - sl[k])$

$sw[j] * (sl[i] - sl[j]) > sw[k] * (sl[i] - sl[k])$

$sw[j] * sl[i] - sw[j] * sl[j] >  sw[k] * sl[i] - sw[k] * sl[k]$

$sw[k] * sl[k] - sw[j] * sl[j] > sw[k] * sl[i] - sw[j] * sl[i]$

$sw[k] * sl[k] - sw[j] * sl[j] > sl[i] * (sw[k] - sw[j])$

$frac{(sw[k] * sl[k] - sw[j] * sl[j])} {(sw[k] - sw[j])} < sl[i]$ //注意sw[k] - sw[j]小于0，要变号

所以令$K = frac{(sw[k] * sl[k] - sw[j] * sl[j])}{(sw[k] - sw[j])}$;

则    while(head < tail && k(q[head],q[head+1]) < sum_len[i])  ++head;

while(head < tail && k(q[tail-1],q[tail]) > k(q[tail],i)) --tail;

最后上代码：

#include<bits/stdc++.h>
using namespace std;
#define R register int
#define AC 20100
int n, ans;
int sum[AC], sum_weight[AC], sum_len[AC], f[AC];
int weight[AC], len[AC];
int q[AC], head, tail;
inline int read()
{
    int x = 0; char c = getchar();
    while(c < '0' || c > '9') c = getchar();
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x; 
}

inline double k(int x, int y)
{
    double a = sum_weight[x] * sum_len[x] - sum_weight[y] * sum_len[y];
    double b = sum_weight[x] - sum_weight[y];
    return a / b;
}

void pre()
{
    n = read();
    for(R i = 1; i <= n; i ++) weight[i] = read(), len[i] = read();
}

void getsum()
{
    for(R i = 1; i <= n + 1; i ++)//山脚的也要求
    {
        sum_len[i] = sum_len[i - 1] + len[i - 1];
        sum_weight[i] = sum_weight[i - 1] + weight[i];
        sum[i] = sum[i - 1] + sum_weight[i - 1] * len[i - 1];
    //    printf("%d : %d
",i,sum[i]);
    }
}

void work()
{
    head=1;
    for(R i = 1; i <= n; i ++)
    {
        f[i] = INT_MAX;
        
        /*int tmp = 0;
        f[i] = INT_MAX;
        for(R j = 1; j < i; j ++)
        {
            if(sum[i] - sum_weight[j] * (sum_len[i] - sum_len[j]) < f[i])
            {
                f[i] = sum[i] - sum_weight[j] * (sum_len[i] - sum_len[j]);
                tmp = j;
            } 
        }
        printf("%d --- > %d
", tmp, i);//打表验证决策单调性*/
        
        while(head < tail && k(q[head], q[head + 1]) < sum_len[i]) ++ head;
        int now = q[head];
    //    printf("%d --- > %d
",now,i);
        f[i] = sum[i] - sum_weight[now] * (sum_len[i] - sum_len[now]);
        while(head < tail && k(q[tail - 1], q[tail]) > k(q[tail], i)) -- tail;
        q[++tail] = i;
    }
    ans = INT_MAX;
    for(R i = 2; i <= n; i ++)//注意应该是n+1,因为山脚是在下面,注意要从2开始，因为这是在枚举第2个厂在哪
        ans = min(ans, f[i] + sum[n + 1] - sum[i] - sum_weight[i] * (sum_len[n + 1] - sum_len[i]));
    printf("%d
", ans);
}

int main()
{
//    freopen("in.in", "r", stdin);
    pre();
    getsum();
    work();
//    fclose(stdin);
    return 0;
}

 ---------------2018.10.12--------------优化了代码格式