多重背包问题及其优化

        本文探讨多重背包问题。要认识这个问题需要先熟悉01背包和完全背包的状态转移方程。

• 01背包

$$dp(a,p)=\begin{cases} \max\{dp(a-1，p)，dp(a-1，p-v_a)+w_a\}\ \ \ & p\geq w_a\\ dp(a-1,p) \end{cases}$$

        数组降维时需要从大到小遍历来更新状态。

• 完全背包

$$dp(a,p)=\begin{cases} \max\{dp(a-1，p)，dp(a，p-v_a)+w_a\}\ \ \ & p\geq w_a\\ dp(a-1,p) \end{cases}$$

        数组降维时需要从小到大遍历来更新状态。

        多重背包是指每一个物品有着数量限制时的背包问题，这也是一类经典动态规划问题，解法主要有以下几种。

• 朴素解法

        如果第i个物品有c_i个，我们可以将这c_i个物品看做不同的物品，从而化为01背包问题。这种方法十分简单，但时间复杂度为$O(V\displaystyle\sum_{i=1}^n w_ic_i)$，时间消耗很大。

• 二进制拆分

        二进制拆分的做法也是将多重背包转化为01背包，但是二进制拆分更为高效，它基于以下数学原理：
        任何一个自然数x都可以写成如下形式：

$$x=2^0+2^1+\cdots+2^k+p,\ p<2^{k+1}$$

        可以理解成将x首先写成2的连续方幂的和，再加上剩余的部分。这k+1个数（如果p=0则为k个数）有一个性质：它们可以表示出0~x中所有的整数。下面给出简单证明：
        一个数都不加就可以得到0，全部加上可以得到x，p可以由单独的p表示出来。对于比p小但比0大的数，它一定小于2^k+1。根据二进制的性质，这个数一定可以通过2⁰~2^k中某些数的组合来表示出来。对于比p大比x小的数，这个数减去p仍然小于2^k+1，这个差仍可以由若干2的方幂来组合表示，最后再加上p即可。证毕。
        那么我们可以将c_i拆成上面的形式，再令得到的每一个数乘以物品价值，问题便转化为等价的01背包问题。这种方法时间复杂度为$O(V\displaystyle\sum_{i=1}^n w_ilog(c_i))$。二进制拆分将物品数量对数化，相比朴素解法是一个很大的改进。这种解法可以过HDU 2844 Coins，但是过不了POJ 1742的相同题目（数据加强）。
        这里附上Coins题目的大致题意：

        给出一些硬币的价格和它们的数量，判断1~m中有多少价格可以用这些硬币组合表示出来。

        题目可以转化为多重背包问题，这里的m为背包容量，硬币为物品（价值与体积相同）。

• 单调队列优化

        单调队列优化可以说是较难理解的一种方式，这种优化思路与前两个不同，它并不是将多重背包问题转化为等价的01背包问题，而是直接从多重背包的状态转移方程入手：

$$dp(a,p)=\max\{dp(a-1，p-kv_a)+kw_a\},\ k\geq 0,p\geq kv_a,k\leq c_i$$

        状转方程是易于理解的，它与01背包其实只有k的差别。常规做法是for一遍a，for一遍p，还要for一遍k，这与朴素解法一样低效。现在拿出v_i的一个剩余系来探讨转移的优化方案。
        所谓剩余系就是在模v_i的基础上，由同余关系建立起的一个个等价类。比如若v_i=5，那么0、5、10在剩余系[0]中，1、6、11在另一个剩余系[1]中。
        对于第二次循环中的p，将其划分为若干个剩余系（显然最多v_i个）。先拿出一个剩余系，假设它们都同余于e（e < v_i）。根据状转方程可以发现：状态的转移都局限于同一个剩余系中。这是因为，转移的状态都是p减去v_i的整数倍，它们显然都是同余的。
        对于剩余系[e]中的某一个元素p，它可以写成$p=uv_i+e$的形式，u为商。将这个式子代入状态转移方程：

$$dp(a,p)=\max\{dp(a-1，uv_a+e-kv_a)+kw_a\},\ k\geq 0,p\geq kv_a,k\leq c_a$$

        再将max里面的每一个元素都减去uw_a，得到：

$$dp(a,p)=\max\{dp(a-1，uv_a+e-kv_a)+kw_a-uw_a\}+uw_a,\ k\geq 0,p\geq kv_a,k\leq c_a$$

        也就是：

$$dp(a,p)=\max\{dp(a-1，e+v_a(u-k))-w_a(u-k)\}+uw_a,\ k\geq 0,p\geq kv_a,k\leq c_a$$

        用k代替u-k，立即得到：

$$dp(a,p)=\max\{dp(a-1，e+kv_a)-kw_a\}+uw_a,\ l_a\leq k\leq r_a\ \cdots (*)$$

        l和r是k取值的限制条件，易知：

$$l_i=\max(0,u-c_i)=u-\min(c_i,u),\ r_i=u$$

        可以发现在同一个剩余系中，状态转移方程(*)max中的元素其实是不变的，可以看做一段区间。每一次转移相当于在某一子区间中内取最值，而区间约束变量k又受p的影响。优化思路就很明显了，这便是如何快速求区间最值的问题。可以发现l和r随着p的增大单调递增，于是就可以考虑用单调队列优化。在实现时，需要将p分成若干个剩余系依次解决。
        以HDU 2844 Coins作示例（会TLE，进一步优化见后文）：

#include<iostream>
#include<cstdio>

using namespace std;
int n, m, A[105], C[105], dp[2][100005], temp[100005], que[100005], key = 0, l, r, ans;//dp是滚动数组

inline int read() {//读入优化
    char e = getchar();
    while ((e < '0' || e > '9') && (e != '-'))e = getchar();
    bool k = false;
    if (e == '-')k = true, e = getchar();
    int s = 0;
    while (e >= '0' && e <= '9')s = s * 10 + e - '0', e = getchar();
    return k ? -s : s;
}

inline void add(int x) {//入队
    while (r >= l && temp[que[r]] <= temp[x])r--;
    que[++r] = x;
}

inline int find(int x) {//找到最大值
    while (que[l] < x)l++;
    return temp[que[l]];
}

int main() {
    while (true) {
        n = read(), m = read(), key = 0, ans = 0;
        if (!n && !m)return 0;
        for (int i = 1; i <= n; i++)A[i] = read();//价格+体积
        for (int i = 1; i <= n; i++)C[i] = read();//数量
        for (int i = 0; i <= m; i++)dp[0][i] = min(i / A[1], C[1]) * A[1];//递推初始化
        for (int i = 2; i <= n; i++) {
            key = !key;
            for (int j = 0; j < A[i]; j++) {//剩余系枚举
                l = 0, r = -1;
                for (int z = 0; j + z * A[i] <= m; z++) {//同一个剩余系内元素枚举
                    temp[z] = dp[!key][j + z * A[i]] - z * A[i], add(z);
                    dp[key][j + z * A[i]] = find(z - min(C[i], z)) + z * A[i];
                }
            }
        }
        for (int i = 1; i <= m; i++)if (dp[key][i] == i)ans++;
        cout << ans << endl;
    }
}

        这种做法会T，这是因为即使用单调队列优化，也不能达到01背包和完全背包的O(VN)复杂度。但是可以发现，对于不同的物品，如果它的数量为1，则可以按01背包处理；如果v_ic_i≥V，则可以按完全背包处理。只有剩余的情况才按照多重背包加单调队列处理，下面的代码还在上面代码基础上做了一些小的常数优化。

#include<iostream>
#include<cstdio>

using namespace std;
int n, m, A[105], C[105], dp[2][100005], temp[100005], que[100005], key = 0, l, r, ans;

inline int read() {
    char e = getchar();
    while ((e < '0' || e > '9') && (e != '-'))e = getchar();
    bool k = false;
    if (e == '-')k = true, e = getchar();
    int s = 0;
    while (e >= '0' && e <= '9')s = s * 10 + e - '0', e = getchar();
    return k ? -s : s;
}

int main() {
    while (true) {
        n = read(), m = read(), key = 0, ans = 0;
        if (!n && !m)return 0;
        for (int i = 1; i <= n; i++)A[i] = read();
        for (int i = 1; i <= n; i++)C[i] = read();
        for (int i = 0; i <= m; i++)dp[0][i] = min(i / A[1], C[1]) * A[1];
        for (int i = 2; i <= n; i++) {
            key = !key;
            if (C[i] == 1) {//01背包
                for (int j = 0; j <= m; j++) {
                    if (j >= A[i])dp[key][j] = max(dp[!key][j], dp[!key][j - A[i]] + A[i]);
                    else dp[key][j] = dp[!key][j];
                }
            } else if (A[i] * C[i] >= m) {//完全背包
                for (int j = 0; j <= m; j++) {
                    if (j >= A[i])dp[key][j] = max(dp[!key][j], dp[key][j - A[i]] + A[i]);
                    else dp[key][j] = dp[!key][j];
                }
            } else {//多重背包+单调队列
                for (int j = 0; j < A[i]; j++) {
                    l = 0, r = -1;
                    for (int z = j, p = 0; z <= m; z += A[i], p++) {
                        temp[p] = dp[!key][z] - p * A[i];
                        while (r >= l && temp[que[r]] <= temp[p])r--;//入队
                        que[++r] = p;
                        while (que[l] < p - min(C[i], p))l++;//出队
                        dp[key][z] = temp[que[l]] + p * A[i];
                    }
                }
            }
        }
        for (int i = 1; i <= m; i++)ans += dp[key][i] == i;
        cout << ans << endl;
    }
}

        但是上面的代码仍然都过不去POJ 1742的强化数据，这是单调队列中两个while循环导致的。其实如果只探讨多重背包的话这里的确可以结束了，下面根据POJ 1742 Coins的特点来进一步优化。
        由于题目只是判断哪些可以被组合表示出来，那么dp数组只需要储存0和1即可，这里0不可以被表示，1为可以。这样0必然是最小值、1必然是最大值。这时的单调队列由于只可能有两种值的情况，可以省去上面代码中的两个while循环。实际做法与普通单调队列类似，差异体现在以下几点：

• 当新元素入队时，不需要将比其小的元素弹出，直接加到队尾即可。这样队列中的元素就是按照序号排列的。
• 保证队列的元素个数不大于c_i，这一步相当于将不在目标区间中的元素剔除。做法是判断长度，若长度达到c_i，队首右移一位。由于队列元素是按照序号排列的，因此右移一位剔除的元素一定不在目标区间中。

        这里会有一个疑问：为什么要保证队列元素数目不大于c_i呢？从前文可知k的取值范围长度为min(u,c_i)+1，可以肯定的是长度一定不会大于u+1，所以维护不大于c_i长度的区间可以保证不越界。当然道理上应该维护c_i+1长度，但是k的最后一个取值很特殊：在数组降维时就是本身，不考虑在内。
        还一个疑问：这两处差异已经违反单调队列的定义了，如何判定最值？其实这种做法可以认为是扩展版的单调队列，故操作方式也有所不同。用变量sum记录队列中值为1的数目，当新加元素为1时sum自加1，队首剔除时sum自减1（队首为1时），这样就可以通过sum来判别：sum为0时，最值为0，否则最值为1。
        根据上面的讨论，两个while循环便可以去除，进一步提高效率解决POJ 1742。

#include<iostream>
#include<cstdio>
#include<cstring>

using namespace std;
int n, m, A[105], C[105], que[100005], ans;
bool dp[100005];//dp数组降维

inline int read() {
    char e = getchar();
    while ((e < '0' || e > '9') && (e != '-'))e = getchar();
    bool k = false;
    if (e == '-')k = true, e = getchar();
    int s = 0;
    while (e >= '0' && e <= '9')s = s * 10 + e - '0', e = getchar();
    return k ? -s : s;
}

int main() {
    while (true) {
        n = read(), m = read(), ans = 0;
        if (!n && !m)return 0;
        for (int i = 1; i <= n; i++)A[i] = read();
        for (int i = 1; i <= n; i++)C[i] = read();
        memset(dp, 0, sizeof(dp));
        for (int i = 0; i <= m && i / A[1] <= C[1]; i += A[1])dp[i] = 1;
        for (int i = 2; i <= n; i++) {
            if (C[i] == 1) {//01背包
                for (int j = m; j >= A[i]; j--) {
                    if (dp[j - A[i]])dp[j] = 1;
                }
            } else if (A[i] * C[i] >= m) {//完全背包
                for (int j = A[i]; j <= m; j++) {
                    if (dp[j - A[i]])dp[j] = 1;
                }
            } else {//多重背包
                for (int j = 0; j < A[i]; j++) {//剩余系枚举
                    int l = 0, r = 0, sum = 0;//[l,r)队列
                    for (int z = j; z <= m; z += A[i]) {//剩余系元素枚举
                        if (r - l > C[i])sum -= que[l++];//长度过大，剔除队首，更新sum
                        sum += que[r++] = dp[z];//sum加上新加入的元素
                        if (sum)dp[z] = 1;//更新状态
                    }
                }
            }
        }
        for (int i = 1; i <= m; i++)ans += dp[i];
        cout << ans << endl;
    }
}