9/14更新：コメント欄にてご教示いただいた実装例を本記事末尾に掲載しました。（C++, Python）

記事の概要

　慣れないC++で重複を許す組み合わせを作りたいと思ったので色々調べたのですが、なかなか見つからないですね。重複を許す「順列」とか、重複組み合わせの総数を計算する公式（中学とかで習った$\small_n H_r$）はあるのですが、重複組み合わせを列挙できるコードはなかなか転がっていませんでした。
　調べてみると、Pythonとかだとitertoolsに普通にcombination_with_replacement()という関数として入ってるようなので、そのコード（ https://docs.python.org/ja/3/library/itertools.html ）を参考にC++で実装し直してみました。
　何かミスや改善点などございましたらご教示くださると幸いです。

※ちなみに重複を許す「順列」の列挙というのが、自分の実装の中でいうデカルト積に相当するのですが、汎用的なものが見当たらなかったので実装し直しています。

アルゴリズムの流れ

　まずは例として、集合$\left\{ 1, 2, 3 \right\}$から重複を許して$2$つの要素を取り出す場合の組み合わせ列挙を考えます。流れ自体は簡単で、以下のようになります。

$\left\{ 1, 2, 3 \right\}$を$2$組用意し、それらのデカルト積（直積集合）を求める。
$$
(1, 1), (1, 2), (1, 3),
(2, 1), (2, 2), (2, 3),
(3, 1), (3, 2), (3, 3)
$$
デカルト積から、ユニークな組み合わせを取り出す。
$$
(1, 1), (1, 2), (1, 3),
(2, 2), (2, 3),
(3, 3)
$$

　つまり一般に、集合$\Omega$から重複を許して$n$個の要素を取り出す場合には、$\Omega$を$n$個に複製し、それらのデカルト積からユニークな組み合わせを取り出せばよいことになります。

デカルト積の生成

　まずはデカルト積を生成するところから実装します。こちらもPythonだとitertoolsのproduct()という関数があるようですが、ソースコードで２重のリスト内包表記を使っていて頭が追い付かなかったので、Golangでデカルト積を実装されている方（ https://qiita.com/tesujiro/items/2e41e1159948dc90c3d9 ）のコードを参考にさせていただきました。こちらの方は一般に$\Omega_1, \Omega_2, \Omega_3, \ldots$という複数集合間でのデカルト積を求めるコードを書かれています。

　今回の自分の実装では、「$\Omega$から重複を許して$n$個取り出す」ということで、一般的なデカルト積でいう$\Omega_1 = \Omega, \Omega_2 = \Omega, \Omega_3 = \Omega, \ldots, \Omega_n = \Omega$に相当します。$n$の値が大きくなってくると、$n$の数だけ$\Omega$を複製して引数に代入するのが面倒になってくるのでrepeatという引数で$n$を数えています。まずはコードを以下に掲載します。　

#include <bits/stdc++.h>
using namespace std;

template<typename T>
vector<vector<T>> productRecurse(vector<T>, int, vector<vector<T>>);

//デカルト積
template<typename T>
vector<vector<T>> product(vector<T> choices, int repeat) {
    vector<vector<T>> emptyResult(1, vector<T>());//組み合わせを格納するための空のリスト

    return productRecurse(choices, repeat, emptyResult);
}

template<typename T>
vector<vector<T>> productRecurse(vector<T> choices, int repeat, vector<vector<T>> result) {
    if (repeat == 0) {return result;}//デカルト積が完成すれば答えを返す
    
    vector<vector<T>> provisionalResult;//作成途中のデカルト積
    
    for (vector<T> re: result) {
        for (T c: choices) {
            vector<T> temp = re;
            temp.push_back(c);
            provisionalResult.push_back(temp);
        }
    }
    
    return productRecurse(choices, repeat - 1, provisionalResult);
}

　上のコードを解説するにあたって、ここでも例として集合$\left\{ 1, 2, 3 \right\}$を$2$つに複製してデカルト積を作成することを考えます。product(choices = {1, 2, 3}, repeat = 2)を呼び出した時の流れは以下のようになっています。

product(choices = {1, 2, 3}, repeat = 2)を呼び出す。
productRecurse(choices = {1, 2, 3}, repeat = 2, result = {{}})が呼び出される。
provisionalResult = {{1}, {2}, {3}}となる。
productRecurse(choices = {1, 2, 3}, repeat = 1, result = {{1}, {2}, {3}})が呼び出される。
provisionalResult = {{1, 1}, {1, 2}, {1, 3}, {2, 1}, {2, 2}, {2, 3}, {3, 1}, {3, 2}, {3, 3}}となる。
productRecurse(choices = {1, 2, 3}, repeat = 0, result = {{1, 1}, {1, 2}, {1, 3}, {2, 1}, {2, 2}, {2, 3}, {3, 1}, {3, 2}, {3, 3}})が呼び出される。
repeat == 0なので引数のresultが最終的なデカルト積となる。

　Pythonだとリスト内包表記で済む部分を再帰的に繰り返しているということになります。では実際に実行してみます。

int main(){
    vector<int> choices = {1, 2, 3};
    int repeat = 2;
    
    vector<vector<int>> answer = product(choices, repeat);
    
    for (int i = 0; i < answer.size(); i++) {
        cout << answer[i][0] << " " << answer[i][1] << endl;
    }
}

　手計算通り、デカルト積が正しく列挙されているのがわかります。

重複組み合わせの抽出

　さて、ここまでくればあとは簡単です。デカルト積からユニークな組み合わせを取り出します。要は例えば$(1, 2)$を取り出す場合、$(2, 1)$は除けば良いわけですので、ある組み合わせcombについてcomb == sort(comb)となるものだけを抽出します。

template<typename T>
vector<vector<T>> combWithReplace(vector<T> choices, int n) {
    vector<vector<T>> answer;
    
    for (vector<T> comb: product(choices, n)) {
        vector<T> toSorted = comb;
        sort(toSorted.begin(), toSorted.end());
        if (comb == toSorted) {
            answer.push_back(comb);
        }
    }
    return answer;
}

実行結果

　さて、重複組み合わせを列挙するプログラムが完成しました。最後に再び、集合$\left\{ 1, 2, 3 \right\}$から重複を許して$2$つの要素を取り出す場合の組み合わせ列挙を試してみます。

int main(){
    vector<int> choices = {1, 2, 3};
    int n = 2;
    
    vector<vector<int>> answer = combWithReplace(choices, n);
    
    for (int i = 0; i < answer.size(); i++) {
        cout << answer[i][0] << " " << answer[i][1] << endl;
    }
}

結果は以下のようになります。手計算と同じ結果が得られていることが分かります。また各組み合わせ内もソートされているのがわかります。

終わりに

　Pythonitertoolsのソースコードを見た限りでは、c++の実装がここまで面倒になるとは思いもしませんでした。リスト内包表記というものが可読性を落としてまで存在しているのには訳があったんですね$\ldots$。
　何かミス等ございましたらご教示くださると幸いです。最後までお読みくださりありがとうございました。

参考文献

追記（9/14）

　以下は@Nabetani 様にご教示いただいた実装例です。パフォーマンス等についてはコメント欄をご覧ください。

　自分の実装が足並み揃えて
{{}}
⇒{{1}, {2}, {3}}
⇒{{1, 1}, {1, 2}, {1, 3}, {2, 1}, {2, 2}, {2, 3}, {3, 1}, {3, 2}, {3, 3}}
とデカルト積を完成させていった後に該当を抽出したのに対し、こちらの実装ではそれぞれの組み合わせを

{}⇒{1}⇒{1, 1}⇒答えに追加
　　　　⇒{1, 2}⇒答えに追加
　　　　⇒{1, 3}⇒答えに追加
{}⇒{2}⇒{2, 2}⇒答えに追加
　　　　⇒{2, 3}⇒答えに追加$\cdots$

という形で完成させては答えに追加しているようです。（構造体の中のappend()関数）
　手計算で重複組み合わせを列挙する際の操作に近いのでこちら方が理解しやすいかもしれません。

c++14

#include <iterator>
#include <vector>

template <typename container_type>
std::vector<std::vector<typename container_type::value_type>>
combWithReplace(container_type const &choices, size_t n) {
  using value_type = typename container_type::value_type;
  using selected_t = std::vector<value_type>;
  using itor_t = typename container_type::const_iterator;
  struct impl { // lambda で再帰は面倒なので クラスにする
    std::vector<std::vector<value_type>> &r_; // コピーを避けるために参照を持つ
    impl(std::vector<std::vector<value_type>> &r) : r_(r) {}
    void append(selected_t &s, itor_t b, itor_t e, size_t n) {
      if (n == 0) {
        r_.push_back(s);
      } else {
        for (auto it = b; it != e; ++it) {
          s.push_back(*it);
          append(s, it, e, n - 1);
          s.pop_back();
        }
      }
    };
  };
  std::vector<std::vector<value_type>> r;
  impl o{r};
  selected_t e;
  e.reserve(n);
  o.append(e, std::cbegin(choices), std::cend(choices), n);
  return r;
}

python3.8


def comb_with_replace(a, n):
    r = []

    def append(e, b, add_count):
        """重複組み合わせの要素を再帰を使って列挙し、見つかっった組み合わせを r に追加する

        Parameters
        ----------
        e : list
            組み立てている途中の重複組合せ
        b : int
            利用可能な要素の先頭のインデックス
        add_count : int
            e に追加すべき要素の数
        """
        if add_count == 0:
            r.append(e)
        else:
            for ix in range(b, len(a)):
                append(e+[a[ix]], ix, add_count-1)
    append([], 0, n)
    return r


def main():
    choices = [1, 2, 3]
    for e in comb_with_replace(choices, 2):
        print(e)


main()