★【动态规划】【斜率优化】【平衡树维护决策序列】【NOI2007】货币兑换

问题描述
小 Y 最近在一家金券交易所工作。该金券交易所只发行交易两种金券:A 纪念券(以下简称 A 券)和 B 纪念券(以下简称 B 券)。每个持有金券的顾客都有一个自己的帐户。金券的数目可以是一个实数。
每天随着市场的起伏波动,两种金券都有自己当时的价值,即每一单位金券当天可以兑换的人民币数目。我们记录第 K 天中 A 券和 B 券的价值分别为 AK 和BK(元/单位金券)。
为了方便顾客,金券交易所提供了一种非常方便的交易方式:比例交易法。比例交易法分为两个方面:
a) 卖出金券:顾客提供一个[0,100]内的实数 OP 作为卖出比例,其意义为:将 OP%的 A 券和 OP%的 B 券以当时的价值兑换为人民币;
b) 买入金券:顾客支付 IP 元人民币,交易所将会兑换给用户总价值为IP 的金券,并且,满足提供给顾客的 A 券和 B 券的比例在第 K 天恰好为 RateK;
例如,假定接下来 3 天内的 Ak、Bk、RateK 的变化分别为:

假定在第一天时,用户手中有 100 元人民币但是没有任何金券。
用户可以执行以下的操作:

注意到,同一天内可以进行多次操作。
小 Y 是一个很有经济头脑的员工,通过较长时间的运作和行情测算,他已经知道了未来 N 天内的 A 券和 B 券的价值以及 Rate。他还希望能够计算出来,如果开始时拥有 S 元钱,那么 N 天后最多能够获得多少元钱。
输入文件
第一行两个正整数 N, S、 分别表示小 Y 能预知的天数以及初始时拥有的钱数。
接下来 N 行,第 K 行三个实数 AK、BK、RateK,意义如题目中所述。
输出文件
只有一个实数 MaxProfit,表示第 N 天的操作结束时能够获得的最大的金钱数目。答案保留 3 位小数。
输入样例
3 100
1 1 1
1 2 2
2 2 3
输出样例
225.000
样例说明

评分方法
本题没有部分分,你的程序的输出只有和标准答案相差不超过0.001时,才能获得该测试点的满分,否则不得分。
数据规模和约定
测试数据设计使得精度误差不会超过 10-7。
对于 40%的测试数据,满足 N ≤ 10;
对于 60%的测试数据,满足 N ≤ 1 000;
对于 100%的测试数据,满足 N ≤ 100 000;
对于 100%的测试数据,满足:
0 < AK ≤ 10;
0 < BK ≤ 10;
0 < RateK ≤ 100
MaxProfit ≤ 109;
提示
输入文件可能很大,请采用快速的读入方式。
必然存在一种最优的买卖方案满足:
每次买进操作使用完所有的人民币;
每次卖出操作卖出所有的金券。

朴素的方程应该很简单（考试的时候没想出来），设f[i]表示第i天最多能得的金钱数量。
fA, fB为第i天最多能得的A、B券数量。
那么有：
f[i] = max{f[i - 1], fA[j] * A[i] + fB[j] * B[i]} (j < i)
代码：

/**************************\
 * @prob: NOI2007 cash    *
 * @auth: Wang Junji      *
 * @stat: TLE: 60         *
 * @date: May. 24th, 2012 *
 * @memo: DP              *
\**************************/
#include <cstdio>
#include <algorithm>
const int maxN = 100010; int n, m;
double f[maxN], fA[maxN], fB[maxN], A[maxN], B[maxN], Rate[maxN];
int main()
{
    freopen("cash.in", "r", stdin);
    freopen("cash.out", "w", stdout);
    scanf("%d%lf", &n, f);
    for (int i = 1; i < n + 1; ++i)
        scanf("%lf%lf%lf", A + i, B + i, Rate + i);
    for (int i = 1; i < n + 1; ++i)
    {
        f[i] = f[i - 1];
        for (int j = 0; j < i; ++j)
            f[i] = std::max(f[i], fA[j] * A[i] + fB[j] * B[i]);
        fB[i] = f[i] / (A[i] * Rate[i] + B[i]);
        fA[i] = fB[i] * Rate[i];
    }
    printf("%.3lf\n", f[n]); return 0;
}

仔细观察前面的方程，忽略掉f[i - 1]项，设j为i的最佳决策，k为异于j的一个决策，则有：
fA[j] * A[i] + fB[j] * B[i] >= fA[k] * A[i] + fB[j] * B[i]，即：
 fA[j] - fA[k]        B[i]
--------------- <= - ------  (fB[j] > fB[k])
 fB[j] - fB[k]        A[i]

 fA[j] - fA[k]        B[i]
--------------- >= - ------  (fB[j] < fB[k])
 fB[j] - fB[k]        A[i]

 fA[j] - fA[k] >= 0          (fB[j] = fB[k])

那么设X[i] = fB[i], Y[i] = fA[i]，（也可以设X[i] = fA[i], Y[i] = fB[i]，依据个人喜好而定。）以上就变成了标准的斜率的表达式。
但是可以发现，其中的X并不单调，而我们需要维护的是一个上凸的凸壳！并且每一次就用斜率为- B[i]/A[i]的直线来从右边无限远处向左边逼近，与此凸壳的第一个交点即为当前最优决策。

于是可以用SBT来维护这个决策序列（以X为关键字），每插入一个点之前先判断它是否应该被插入（若它被包含在之前求出来的凸壳内（含边界）则不插入，否则插入），然后如果插入则巩固一下凸壳的性质，即去掉一些被包含在新的凸壳内的一些点，具体细节见程序。

代码：

/****************************************\
 * @prob: NOI2007 cash                  *
 * @auth: Wang Junji                    *
 * @stat: TLE: 80                       *
 * @date: May. 26th, 2012               *
 * @memo: 斜率优化DP，SBT维护决策序列       *
\****************************************/
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cstring>
#include <string>
#include <cmath>
#define K(a, b)  ((Y[a] - Y[b]) / (X[a] - X[b]))

const int maxN = 100010; const double INF = 1e198, zero = 1e-12;

double f[maxN], A[maxN], B[maxN], Rate[maxN], X[maxN], Y[maxN]; int n;
class SBT
{
private:
    int key[maxN], sz[maxN], lc[maxN], rc[maxN], T, tot;
    void Zig(int &T)
    {
        int tmp = lc[T]; lc[T] = rc[tmp]; rc[tmp] = T; sz[tmp] = sz[T];
        sz[T] = sz[lc[T]] + sz[rc[T]] + 1; T = tmp; return;
    }
    void Zag(int &T)
    {
        int tmp = rc[T]; rc[T] = lc[tmp]; lc[tmp] = T; sz[tmp] = sz[T];
        sz[T] = sz[lc[T]] + sz[rc[T]] + 1; T = tmp; return;
    }
    void maintain(int &T, bool flag)
    {
        if (!T || (!lc[T] && !rc[T])) return;
        if (!flag)
        {
            if (sz[lc[lc[T]]] > sz[rc[T]]) Zig(T);
            else if (sz[rc[lc[T]]] > sz[rc[T]]) Zag(lc[T]), Zig(T);
            else return;
        }
        else
        {
            if (sz[rc[rc[T]]] > sz[lc[T]]) Zag(T);
            else if (sz[lc[rc[T]]] > sz[lc[T]]) Zig(rc[T]), Zag(T);
            else return;
        }
        maintain(lc[T], 0); maintain(rc[T], 1);
        maintain(T, 0); maintain(T, 1); return;
    }
    void Ins(int &T, int v)
    {
        if (!T) {sz[T = ++tot] = 1, key[T] = v; return;}
        ++sz[T]; Ins(((X[v] - X[key[T]] < -zero) ? lc[T] : rc[T]), v);
        maintain(T, X[v] - X[key[T]] > -zero); return;
    }
    int Del(int &T, int v)
    {
        --sz[T]; double delta = X[v] - X[key[T]];
        if (fabs(delta) < zero ||
            (delta < -zero && !lc[T]) ||
            (delta > zero && !rc[T]))
        {
            int tmp = key[T];
            if (!lc[T] || !rc[T]) T = lc[T] + rc[T];
            else key[T] = Del(lc[T], key[T]);
            return tmp;
        }
        return Del(((delta < -zero) ? lc[T] : rc[T]), v);
    }
    int pred(int v)
    {
        X[0] = -INF; 
        int ans = 0;
        for (int t = T; t;)
        {
            if (X[key[t]] - X[v] < -zero && X[key[t]] - X[ans] > zero) ans = key[t];
            t = (X[v] - X[key[t]] < zero) ? lc[t] : rc[t];
        }
        X[0] = 0;
        return ans;
    }
    int succ(int v)
    {
        X[0] = INF; 
        int ans = 0;
        for (int t = T; t;)
        {
            if (X[key[t]] - X[v] > -zero && key[t] - v //避免找后继时找到本身。
                && X[key[t]] - X[ans] < -zero) ans = key[t];
            t = (X[v] - X[key[t]] < -zero) ? lc[t] : rc[t];
        }
        X[0] = 0;
        return ans;
    }
public:
    SBT(): tot(0), T(0)
    {
        memset(key, 0, sizeof key); memset(sz, 0, sizeof sz);
        memset(lc, 0, sizeof lc); memset(rc, 0, sizeof rc);
    }
    void Ins(int v)
    {
        int pre = pred(v), suc = succ(v);
        //在插入之前先找到离这个横坐标最近的两个点
        do
        {
            if (!pre || !suc) break; //若前驱不存在或后继不存在则直接插入。
            if (fabs(X[v] - X[suc]) < zero)
            {
                if (Y[v] - Y[suc] > -zero) Del(T, suc), suc = succ(v);
                else return;
            }
            //若要插入的点v与之前凸壳上某一点的横坐标相同，则舍弃纵坐标较小的点。
            if (K(pre, v) - K(v, suc) < zero) return;
            //若在v处向内凹陷，则放弃插入。
        } while (0);
        Ins(T, v);
        for (int ppre = pred(pre); ; pre = ppre, ppre = pred(pre))
        {
            if (!ppre) break;
            if (K(ppre, pre) - K(pre, v) < zero) Del(T, pre); else break;
        } //将v左侧向内凹陷的点一一去掉。
        for (int ssuc = succ(suc); ; suc = ssuc, ssuc = succ(suc))
        {
            if (!ssuc) break;
            if (K(v, suc) - K(suc, ssuc) < zero) Del(T, suc); else break;
        } //将v右侧向内凹陷的点一一去掉。
        return;
    }
    int Find(double k) 
    {
        for (int t = T; t;)
        {
            int pre = pred(key[t]), suc = succ(key[t]);
            double k1 = pre ? K(pre, key[t]) : k + 1,
                   k2 = suc ? K(key[t], suc) : k - 1;
            if (k1 - k > -zero && k - k2 > -zero) return key[t];
            t = ((k1 - k < -zero) ? lc[t] : rc[t]);
        }
        return 0;
    } //在平衡树中找到与斜率为k的直线相切的切点。
} mp;

int main()
{
    freopen("cash.in", "r", stdin);
    freopen("cash.out", "w", stdout);
    scanf("%d%lf", &n, f);
    for (int i = 1; i < n + 1; ++i)
    {
        scanf("%lf%lf%lf", A + i, B + i, Rate + i);
        int pos = mp.Find(-B[i] / A[i]);
        f[i] = std::max(f[i - 1], Y[pos] * A[i] + X[pos] * B[i]);
        X[i] = f[i] / (A[i] * Rate[i] + B[i]);
        Y[i] = X[i] * Rate[i];
        mp.Ins(i);
    }
    printf("%.3lf\n", f[n]); return 0;
}

可以发现，用SBT维护决策序列的方法中，每插入一个点需要逐一进行检查左右的凹陷处，效率偏低，所以考虑使用伸展树进行优化。

自然就是将刚刚插入的点作为根，在左右两边分别找到凹陷的转折点然后一并删掉。

还有一个改进就是没插入一个点的时候就事先把它与前驱结点的斜率和与后继结点的斜率分别计算出来，这样即使树的结构改变前驱后继也不会变。

代码：

/*************************************************\
 * @prob: NOI2007 cash * @auth: Wang Junji       *
 * @stat: Accepted.    * @date: May. 27th, 2012  *
 * @memo: 动态规划、斜率优化、Splay维护决策序列        *
\*************************************************/
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cstring>
#include <string>
#include <cmath>
#define K(a, b)  ((Y[a] - Y[b]) / (X[a] - X[b]))

const int maxN = 100010; const double INF = 1e198, zero = 1e-12;

double f[maxN], A[maxN], B[maxN], Rate[maxN];
double X[maxN], Y[maxN], lk[maxN], rk[maxN]; int n;
class SplayTree
{
private:
    int key[maxN], F[maxN], lc[maxN], rc[maxN], T, tot;
    void Zig(int T)
    {
        int P = F[T], tmp = rc[T];
        if (lc[F[P]] == P) lc[F[P]] = T; else rc[F[P]] = T;
        F[T] = F[P]; F[P] = T; lc[P] = tmp; F[tmp] = rc[T] = P;
        return;
    }
    void Zag(int T)
    {
        int P = F[T], tmp = lc[T];
        if (lc[F[P]] == P) lc[F[P]] = T; else rc[F[P]] = T;
        F[T] = F[P]; F[P] = T; rc[P] = tmp; F[tmp] = lc[T] = P;
        return;
    }
    void Splay(int T, int t)
    {
        if (!t) return; int tmp = F[T];
        while (F[t] - tmp)
        {
            int P = F[t];
            if (F[P] == tmp) (lc[P] == t) ? Zig(t) : Zag(t);
            else
            {
                if (lc[F[P]] == P) (lc[P] == t) ? Zig(P) : Zag(t), Zig(t);
                else (lc[P] == t) ? Zig(t) : Zag(P), Zag(t);
            }
        }
        if (!tmp) this -> T = t; return;
    }
    int FixL()
    {
        int tmp;
        for (int t = lc[T]; t;)
            if (lk[key[t]] - K(key[t], key[T]) < zero) t = lc[t];
            else tmp = t, t = rc[t];
        return tmp;
    }
    int FixR()
    {
        int tmp;
        for (int t = rc[T]; t;)
            if (K(key[t], key[T]) - rk[key[t]] < zero) t = rc[t];
            else tmp = t, t = lc[t];
        return tmp;
    }
    int Union(int T1, int T2)
    {
        int tmp = T1; while (rc[tmp]) tmp = rc[tmp];
        rc[tmp] = T2; F[T2] = tmp;
        return T1;
    }
    void update()
    {
        int t = T;
        lk[key[t]] = lc[t] ? (rk[key[lc[t]]] = K(key[t], key[lc[t]])) : INF;
        //注意计算根的前驱的同时顺便把左子结点的后继求出来。
        rk[key[t]] = rc[t] ? (lk[key[rc[t]]] = K(key[t], key[rc[t]])) : -INF;
        //注意计算根的后继的同时顺便把右子结点的前驱求出来。
        if (lk[key[t]] - rk[key[t]] < zero)
        {
            F[T = Union(lc[t], rc[t])] = 0;
            rk[key[T]] = lk[key[rc[T]]] = K(key[T], key[rc[T]]);
        } //若根处不满足凸性，则删掉根并重新计算后继。
        return;
    }
public:
    SplayTree(): tot(0), T(0)
    {
        memset(key, 0, sizeof key); memset(F, 0, sizeof F);
        memset(lc, 0, sizeof lc); memset(rc, 0, sizeof rc);
    }
    void Ins(int v)
    {
        int tmp = T, t = T;
        while (t)
        {
            if (X[v] - X[key[t]] < -zero) tmp = t, t = lc[t];
            else if (X[v] - X[key[t]] > zero) tmp = t, t = rc[t];
            else if (Y[v] - Y[key[t]] < zero) return;
            else {Y[key[t]] = Y[v], v = key[t]; break;}
        }
        //找到插入点，若v在某个点的正上方
        //则直接用v将其取代并进行下一步的维护。
        if (!t)
        {
            key[t = ++tot] = v;
            (X[v] - X[key[tmp]] < -zero) ? (lc[tmp] = t) : (rc[tmp] = t);
            F[t] = tmp;
        } //插入操作。
        Splay(T, t); //将t旋转至根。
        if (lc[t]) {int tmp = FixL(); Splay(lc[T], tmp); rc[tmp] = 0;}
        //维护左边的凸壳，即找到第一个满足凸性的点
        //并把它旋转至左子树的根，再把它的右子树删掉，
        //于是这时根的前驱一定就是左子结点。
        if (rc[t]) {int tmp = FixR(); Splay(rc[T], tmp); lc[tmp] = 0;}
        //维护右边的凸壳，即找到第一个满足凸性的点
        //并把它旋转至右子树的根，再把它的左子树删掉，
        //于是这时根的后继一定就是右子结点。
        update(); return; //计算出与前驱和后继的斜率并维护v处的凸性。
    }
    int Find(double k)
    {
        for (int t = T; ;)
        {
            if (lk[key[t]] - k < -zero) t = lc[t];
            else if (rk[key[t]] - k > zero) t = rc[t];
            else return key[t];
        }
    }
} mp;

inline double getreal()
{
    static long long base[18] = {1}; static bool flag = 1;
    if (flag) 
    {
        flag = 0; 
        for (int i = 1; i < 18; ++i) base[i] = base[i - 1] * 10;
    }
    int Int = 0, cnt = 0; long long Dec = 0; char tmp; bool sgn = 1;
    do tmp = getchar();
    while (!isdigit(tmp) && tmp - '-');
    if (tmp == '-') tmp = getchar(), sgn = 0;
    do Int = (Int << 3) + (Int << 1) + tmp - '0';
    while (isdigit(tmp = getchar()));
    if (tmp == '.') tmp = getchar();
    else return sgn ? Int : -Int;
    do Dec = (Dec << 3) + (Dec << 1) + tmp - '0', ++cnt;
    while (isdigit(tmp = getchar()));
    return sgn ? Int + (double)Dec / base[cnt]
        : -Int - (double)Dec / base[cnt];
}
//手写读入实数……有些麻烦，题目中说了要用
//最快的读入方式。（虽然不用这个也能过。）

int main()
{
    freopen("cash.in", "r", stdin);
    freopen("cash.out", "w", stdout);
    scanf("%d%lf", &n, f);
    for (int i = 1; i < n + 1; ++i)
    {
        A[i] = getreal(); B[i] = getreal(); Rate[i] = getreal();
        int pos = i == 1 ? 0 : mp.Find(-B[i] / A[i]);
        f[i] = std::max(f[i - 1], Y[pos] * A[i] + X[pos] * B[i]);
        X[i] = f[i] / (A[i] * Rate[i] + B[i]);
        Y[i] = X[i] * Rate[i];
        mp.Ins(i);
    }
    printf("%.3lf\n", f[n]); return 0;
}