Pisinger 快速解决子集和算法

Question

这是我之前的question 的后续。我仍然觉得这是一个非常有趣的问题，因为有一种算法值得更多关注，所以我将其发布在这里。

来自Wikipedia：对于每个 xi 为正且以相同常数为界的情况，Pisinger 找到了线性时间算法。

有一篇不同的论文似乎描述了相同的算法，但对我来说有点难以阅读，所以请 - 有人知道如何将第 4 页 (balsub) 中的伪代码翻译成工作实现吗？

以下是我目前收集到的几点建议：

http://www.diku.dk/~pisinger/95-6.ps（论文）
https://stackoverflow.com/a/9952759/1037407
http://www.diku.dk/hjemmesider/ansatte/pisinger/codes.html

PS：我并不真正坚持这个算法，所以如果你知道任何其他类似性能的算法，请随时在下面提出建议。

编辑

这是 oldboy 在下面发布的代码的 Python 版本：

class view(object):
    def __init__(self, sequence, start):
        self.sequence, self.start = sequence, start
    def __getitem__(self, index):
        return self.sequence[index + self.start]
    def __setitem__(self, index, value):
        self.sequence[index + self.start] = value

def balsub(w, c):
    '''A balanced algorithm for Subset-sum problem by David Pisinger
    w = weights, c = capacity of the knapsack'''
    n = len(w)
    assert n > 0
    sum_w = 0
    r = 0
    for wj in w:
        assert wj > 0
        sum_w += wj
        assert wj <= c
        r = max(r, wj)
    assert sum_w > c
    b = 0
    w_bar = 0
    while w_bar + w[b] <= c:
        w_bar += w[b]
        b += 1
    s = [[0] * 2 * r for i in range(n - b + 1)]
    s_b_1 = view(s[0], r - 1)
    for mu in range(-r + 1, 1):
        s_b_1[mu] = -1
    for mu in range(1, r + 1):
        s_b_1[mu] = 0
    s_b_1[w_bar - c] = b
    for t in range(b, n):
        s_t_1 = view(s[t - b], r - 1)
        s_t = view(s[t - b + 1], r - 1)
        for mu in range(-r + 1, r + 1):
            s_t[mu] = s_t_1[mu]
        for mu in range(-r + 1, 1):
            mu_prime = mu + w[t]
            s_t[mu_prime] = max(s_t[mu_prime], s_t_1[mu])
        for mu in range(w[t], 0, -1):
            for j in range(s_t[mu] - 1, s_t_1[mu] - 1, -1):
                mu_prime = mu - w[j]
                s_t[mu_prime] = max(s_t[mu_prime], j)
    solved = False
    z = 0
    s_n_1 = view(s[n - b], r - 1)
    while z >= -r + 1:
        if s_n_1[z] >= 0:
            solved = True
            break
        z -= 1
    if solved:
        print c + z
        print n
        x = [False] * n
        for j in range(0, b):
            x[j] = True
        for t in range(n - 1, b - 1, -1):
            s_t = view(s[t - b + 1], r - 1)
            s_t_1 = view(s[t - b], r - 1)
            while True:
                j = s_t[z]
                assert j >= 0
                z_unprime = z + w[j]
                if z_unprime > r or j >= s_t[z_unprime]:
                    break
                z = z_unprime
                x[j] = False
            z_unprime = z - w[t]
            if z_unprime >= -r + 1 and s_t_1[z_unprime] >= s_t[z]:
                z = z_unprime
                x[t] = True
        for j in range(n):
            print x[j], w[j]

Answer 1

// Input:
// c (capacity of the knapsack)
// n (number of items)
// w_1 (weight of item 1)
// ...
// w_n (weight of item n)
//
// Output:
// z (optimal solution)
// n
// x_1 (indicator for item 1)
// ...
// x_n (indicator for item n)

#include <algorithm>
#include <cassert>
#include <iostream>
#include <vector>

using namespace std;

int main() {
  int c = 0;
  cin >> c;
  int n = 0;
  cin >> n;
  assert(n > 0);
  vector<int> w(n);
  int sum_w = 0;
  int r = 0;
  for (int j = 0; j < n; ++j) {
    cin >> w[j];
    assert(w[j] > 0);
    sum_w += w[j];
    assert(w[j] <= c);
    r = max(r, w[j]);
  }
  assert(sum_w > c);
  int b;
  int w_bar = 0;
  for (b = 0; w_bar + w[b] <= c; ++b) {
    w_bar += w[b];
  }
  vector<vector<int> > s(n - b + 1, vector<int>(2 * r));
  vector<int>::iterator s_b_1 = s[0].begin() + (r - 1);
  for (int mu = -r + 1; mu <= 0; ++mu) {
    s_b_1[mu] = -1;
  }
  for (int mu = 1; mu <= r; ++mu) {
    s_b_1[mu] = 0;
  }
  s_b_1[w_bar - c] = b;
  for (int t = b; t < n; ++t) {
    vector<int>::const_iterator s_t_1 = s[t - b].begin() + (r - 1);
    vector<int>::iterator s_t = s[t - b + 1].begin() + (r - 1);
    for (int mu = -r + 1; mu <= r; ++mu) {
      s_t[mu] = s_t_1[mu];
    }
    for (int mu = -r + 1; mu <= 0; ++mu) {
      int mu_prime = mu + w[t];
      s_t[mu_prime] = max(s_t[mu_prime], s_t_1[mu]);
    }
    for (int mu = w[t]; mu >= 1; --mu) {
      for (int j = s_t[mu] - 1; j >= s_t_1[mu]; --j) {
        int mu_prime = mu - w[j];
        s_t[mu_prime] = max(s_t[mu_prime], j);
      }
    }
  }
  bool solved = false;
  int z;
  vector<int>::const_iterator s_n_1 = s[n - b].begin() + (r - 1);
  for (z = 0; z >= -r + 1; --z) {
    if (s_n_1[z] >= 0) {
      solved = true;
      break;
    }
  }
  if (solved) {
    cout << c + z << '\n' << n << '\n';
    vector<bool> x(n, false);
    for (int j = 0; j < b; ++j) x[j] = true;
    for (int t = n - 1; t >= b; --t) {
      vector<int>::const_iterator s_t = s[t - b + 1].begin() + (r - 1);
      vector<int>::const_iterator s_t_1 = s[t - b].begin() + (r - 1);
      while (true) {
        int j = s_t[z];
        assert(j >= 0);
        int z_unprime = z + w[j];
        if (z_unprime > r || j >= s_t[z_unprime]) break;
        z = z_unprime;
        x[j] = false;
      }
      int z_unprime = z - w[t];
      if (z_unprime >= -r + 1 && s_t_1[z_unprime] >= s_t[z]) {
        z = z_unprime;
        x[t] = true;
      }
    }
    for (int j = 0; j < n; ++j) {
      cout << x[j] << '\n';
    }
  }
}

Answer 2

伟大的代码人，但它有时会在这个代码块中崩溃

    for (mu = w[t]; mu >= 1; --mu) 
    {
        for (int j = s_t[mu] - 1; j >= s_t_1[mu]; --j) 
        {
            if (j >= w.size())
            { // !!! PROBLEM !!!

            }


            int mu_prime = mu - w[j];

            s_t[mu_prime] = max(s_t[mu_prime], j);
        }
    }

...