Didn't know will it be useful for someone and is it on the right path to solve the problem, but I wish to publish this code for SSP solver. To invent good heuristics for MSSP is utmost harder.

using I = short;

template<typename T = I, typename S = I>
struct DDP // dynamic dynamic programming solver for subset sum problem O(n^3 * 2^sqrt(n)) (n is number of input bits)
{
    T terms[2][10001];
    S sums[2][10001];

    template<typename Term, typename Sum>
    std::pair<Term, Sum> solve(const Term termsBeg, const Term termsEnd, const Sum sumsBeg, const Sum sumsEnd)
    {
        static_assert(std::is_same<typename std::iterator_traits<Term>::value_type, T>::value, "!");
        static_assert(std::is_same<typename std::iterator_traits<Sum>::value_type, S>::value, "!");

        auto sumMax = std::max_element(sumsBeg, sumsEnd);

        auto termsSrc = terms[0];
        auto termsDst = terms[1];

        auto sumsSrc = sums[0];
        auto sumsDst = sums[1];

        *termsSrc = {};
        *sumsSrc = {};

        auto sumsSrcEnd = std::next(sumsSrc);
        for (auto t = termsBeg; t != termsEnd; ++t) {
            auto termsSrcBeg = termsSrc;
            auto termsDstBeg = termsDst;

            auto sumsSrcBeg = sumsSrc;
            auto sumsShiftedBeg = sumsSrc;
            auto sumsDstBeg = sumsDst;

            // merge unique
            assert(std::size_t(std::distance(sumsSrcBeg, sumsSrcEnd)) <= std::extent<std::remove_reference_t<decltype(sums[0])>>::value);
            while (sumsSrcBeg != sumsSrcEnd) {
                if (sumsShiftedBeg == sumsSrcEnd) {
                    break;
                }
                I shiftedSum = *sumsShiftedBeg + *t;
                if (shiftedSum < *sumsSrcBeg) {
                    *sumsDstBeg++ = shiftedSum;
                    ++sumsShiftedBeg;
                    *termsDstBeg++ = *t;
                } else {
                    if (shiftedSum == *sumsSrcBeg) {
                        ++sumsShiftedBeg;
                        assert(*sumsShiftedBeg + *t != *sumsSrcBeg);
                    }
                    *sumsDstBeg++ = *sumsSrcBeg++;
                    *termsDstBeg++ = *termsSrcBeg++;
                }
                assert(std::size_t(std::distance(sumsDst, sumsDstBeg)) <= std::extent<std::remove_reference_t<decltype(sums[0])>>::value);
            }
            if (sumsShiftedBeg == sumsSrcEnd) {
                sumsSrcEnd = std::copy(sumsSrcBeg, sumsSrcEnd, sumsDstBeg);
                termsDstBeg = std::copy_n(termsSrcBeg, std::distance(sumsDstBeg, sumsSrcEnd), termsDstBeg);
            } else {
                assert(sumsSrcBeg == sumsSrcEnd);
                sumsSrcEnd = std::transform(sumsShiftedBeg, sumsSrcEnd, sumsDstBeg, [t] (S s) { return s + *t; });
                termsDstBeg = std::fill_n(termsDstBeg, std::distance(sumsDstBeg, sumsSrcEnd), *t);
            }

            std::swap(termsSrc, termsDst);
            std::swap(sumsSrc, sumsDst);

            for (auto sum = sumsBeg; sum != sumsEnd; ++sum) {
                sumsSrcBeg = std::lower_bound(sumsSrc, sumsSrcEnd, *sum);
                if ((sumsSrcBeg != sumsSrcEnd) && (*sumsSrcBeg == *sum)) {
                    sumsSrcEnd = sumsSrcBeg;
                    assert((*std::next(termsSrc, std::distance(sumsSrc, sumsSrcEnd)) == *t));
                    auto termsDstEnd = termsDst;
                    do {
                        auto term = *std::next(termsSrc, std::distance(sumsSrc, sumsSrcEnd));
                        *termsDstEnd++ = term;
                        sumsSrcEnd = std::lower_bound(sumsSrc, sumsSrcEnd, *sumsSrcEnd - term);
                    } while (sumsSrcEnd != sumsSrc);
                    assert((std::accumulate(termsDst, termsDstEnd, S{}) == *sum));
                    termsDstBeg = termsDst;
                    auto termsSrcEnd = t;
                    assert(*termsSrcEnd == *termsDstBeg);
                    auto termsTail = termsSrcEnd;
                    while (++termsDstBeg != termsDstEnd) {
                        while (*--termsSrcEnd != *termsDstBeg) {
                            *termsTail-- = *termsSrcEnd;
                        }
                    }
                    while (termsBeg != termsSrcEnd) {
                        *termsTail-- = *--termsSrcEnd;
                    }
                    assert(std::distance(termsBeg, termsTail) + 1 == std::distance(termsDst, termsDstEnd));
                    auto term = termsTail;
                    termsDstBeg = termsDst;
                    assert(termsDstBeg != termsDstEnd);
                    do {
                        *termsTail-- = *termsDstBeg++;
                    } while (termsDstBeg != termsDstEnd);
                    assert(termsTail == std::prev(termsBeg));
                    return {term, sum}; // not past the end element for found range but exactly last element
                }
            }
            sumsSrcEnd = std::lower_bound(sumsSrc, sumsSrcEnd, *sumMax); // crop
        }
        return {termsEnd, sumsEnd};
    }
};

auto x = solve(std::begin(ships), std::end(ships), std::cbegin(rows), std::cend(rows)) if succeeded, it gives rearranged ships range such, that range [std::begin(ships), x.first] is a solution and sum of all ships is *x.second. [std::begin(ships), std::end(ships)) still consists of exactly the same ships. If not succeeded then it return a pair of std::end(ships) and std::cend(rows)

Edited by author 11.06.2020 12:53