ラベル algorithm の投稿を表示しています。 すべての投稿を表示
ラベル algorithm の投稿を表示しています。 すべての投稿を表示

2019年6月2日日曜日

[C++][Google CodeJam] C++20 rangeでGCJ2019 R2-Aを解く

5年空いた更新でのっけから言い訳で始まりますが、C++20どころかC++14から周回遅れ状態で突っ込みどころが多々あるかと思いますがそれ相応にお願いします。

さて2008年以降毎年参加しているGoogle CodeJamですが、最新のC++で頑張る(特に時間に余裕のあるQualification Roundでは)という裏目的は2018年にコンテストシステムがローカル実行からリモート実行へ変更されて以降難しくなってしまいました(2019年ではgcc-6.3.0で-std=c++14になっています)。ということでC++欲が高まっているところにGCJ R2-Aが割ときれいに書けた&cmcstl2を見つけた(遅い)ということでC++20 rangeで書いてみたものです。

該当の問題はNew Elements: Part 1。実装した解法はanalysisに書いてあるのと同じだと思います。競技時間中に書いていたものをC++20 rangeを使って書き直したものが↓でwandboxによる実行結果がこれです。

#include "cmcstl2/include/meta/meta_fwd.hpp" // force to use meta in cmcstl2
#include "cmcstl2/include/meta/meta.hpp" // same as above
#include <experimental/ranges/ranges>
#include <experimental/ranges/iterator>
#include <experimental/ranges/algorithm>
#include <vector>
#include <set>
#include <iostream>
#include <utility>
#include <cstdint>

namespace ranges = std::experimental::ranges;

// for ADL for operator>>
struct vec2 : std::pair<int,int> { using std::pair<int,int>::pair; };
inline std::istream& operator>>(std::istream& is, vec2& p) { return is >> p.first >> p.second; }
inline constexpr vec2 operator-(const vec2 &v1, const vec2 &v2) { return { v1.first - v2.first, v1.second - v2.second }; }
inline constexpr vec2 operator-(const vec2 &v1) { return { -v1.first, -v1.second }; }
inline constexpr auto ratio_less = [](const std::pair<int,int> &v1, const std::pair<int,int> &v2){
    return static_cast<std::int64_t>(v1.second) * v2.first > static_cast<std::int64_t>(v1.first) * v2.second;
};

// from http://ericniebler.com/2018/12/05/standard-ranges/
inline constexpr auto for_each = [] <ranges::Range R, ranges::Iterator I = ranges::iterator_t<R>, ranges::IndirectUnaryInvocable<I> Fun>
                                    (R&& r, Fun fun) /* requires ranges::Range<ranges::indirect_result_t<Fun, I>> */ {
    return std::forward<R>(r) | ranges::view::transform(std::move(fun)) | ranges::view::join;
};
inline constexpr auto all_pairs = [](int N) {
    return for_each(ranges::iota_view(0, N), [=](int n){ return ranges::view::transform(ranges::iota_view(n+1, N), [=](int m){ return std::make_pair(n, m); }); });
};
inline constexpr auto all_pairs_r = [](auto &&r){
    return for_each(
        ranges::iota_view(ranges::begin(r), ranges::end(r)) | ranges::view::transform([e=ranges::end(r)](auto it){return ranges::subrange(it, e);}),
        [](auto &&rr){
            return ranges::view::transform(rr.next(), [pp1=ranges::begin(rr)](auto p2){ return std::make_pair(*pp1, p2); });
        }
    );
};

int main()
{
    int P; std::cin >> P;
    ranges::for_each(ranges::iota_view(1, P+1), [&](int p){
        int N; std::cin >> N;
        std::vector<vec2> v(N);
        // At least, ranges::copy_n can't be used because input iterator is incremented for a result value
        // see also http://www.open-std.org/jtc1/sc22/wg21/docs/lwg-active.html#2471
        std::copy_n(ranges::istream_iterator<vec2>(std::cin), N, ranges::begin(v));
#if 1
        auto r = all_pairs_r(v) | ranges::view::transform([](auto p){ return p.first - p.second; })
                                | ranges::view::filter([](auto p){ return p.first * p.second < 0; })
                                | ranges::view::transform([](auto p){ return p.first > 0 ? p : -p; })
                                | ranges::view::common;
#else
        auto r = all_pairs(N) | ranges::view::transform([&v](auto p){ return v[p.first] - v[p.second]; })
                              | ranges::view::filter([](auto p){ return p.first * p.second < 0; })
                              | ranges::view::transform([](auto p){ return p.first > 0 ? p : -p; })
                              | ranges::view::common;
#endif
        // Use deduction guide (gcc8 required for libstdc++ part)
        std::set s(ranges::begin(r), ranges::end(r), ratio_less);
        std::cout << "Case #" << p << ": " << s.size() + 1 << std::endl;
    });
    return 0;
}

以下簡単な説明。

最初の2行とコンパイルオプション-I/opt/wandboxはwandboxだとcmcstl2内のmetaとrange-v3のmetaが衝突するのでその回避用です。まず必要なヘッダを#includeして適当にnamespace aliasを作成。vec2を定義しているのはiostream用のoperator>>がADLで引っ掛けられるようにするためです(std::pairのままだとstd名前空間内にoperator>>を定義しないとADLで引っ掛けられないが未定義動作になる)。using部分はinheriting constructorsです。このvec2に対して、iostream用operator>>、減算、単項マイナス、有理数として見た場合の比較関数(の符号手抜き版)を定義します。比較関数は32bit環境向けに64ビット変数へのキャストを入れています。

for_eachはRangeの神様 Eric Niebler 氏の記事からぱくってきたものです。Constraint lambdaがとてもえぐいですがシグネチャとしては明確になっていると思います(でもあまり書きたいとは思えない)。コメントアウトしてあるrequires部分はtrailing requires-clauseとしてここにかけるはずなのですがgcc-9.1.0では通りませんでした(HEADだとコンパイラが落ちる)。

このfor_eachを利用して、N未満の2整数の組のrangeを生成するall_pairsと、渡されたrange中の2要素の組のrangeを生成するall_pairs_rを定義しています。all_pairs_rについてはとりあえず動いてるように見えますが書いた当人がまじで動いてるのこれ?状態です。

最初の ranges::for_each にはあまり意味はなく、せっかくだら for なくそうというだけでしかありません。

最初にデータ数を読みだしてきてstd::copy_nを使ってstd::cinからvec2の列を読みだしてきます。このとき少なくともranges::copy_nだと意図通り動作しません。ranges::copy_nはstd::copy_nと違い入力側のiteratorについても処理後の値を返り値で返すのですがいわゆるone-past-endの値になっています。一方でistream_iteratorはインクリメント時に次のデータを読みだしてくるのでN個コピーするときにはN+1個読みだされた状態になってしまいます(実質1個分捨てられる)。実際にはstd::copy_nについても現状は明確な記載がなくlwg-2471がOpenになっていますが既知の実装ではN-1回インクリメントになっているようです。

次がメインの処理で、all_pairs起点の場合、N未満の2整数の組を生成→これを添え字としたときのvec2列の要素間の差を生成→firstとsecondの積が負になる場合(=傾きが負になる場合)のみに限定→first(分母)が正になるように正規化しています。all_pairs_r起点の場合には添え字の組ではなく最初から要素の組が生成されています。最後のranges::view::commonは次にsetのコンストラクタ(Iteratorの組を引数にとるもの)に渡すためにbeginとendが同じ型になるようにしています。これを有理数としてみた場合のsetに突っ込んで独立な傾きの数を求めて+1して終了です。setではクラステンプレートのテンプレート引数推論で型引数を省略しています。ranges::to: A function to convert any range to a containerが入れば

auto s = all_pairs_r(v) | ... | ranges::to<std::set>(ratio_less);

みたいに書けるようになるかもしれません(多分)。Input Range Adaptorsも加えて

auto v = ranges::istream_view<vec2>(std::cin) | ranges::view::take(N) | ranges::to<std::vector>;

とでも書けるともっといいと思いますが、copy_nの時と同じ話になってしまいます。せっかく無限Rangeが表現できるようになったので++時には入力を読まない版ができればいいのかもしれません。

どうでしょうか。もちろんこの処理なら普通に書いてもよっぽどシンプルですがまるでアルゴリズムの教科書に書いてあるかのごとく処理内容・意図が簡潔・明確に示されている感じがしないでしょうか。C++11以降とそれより前でC++は大きく変わっているわけですが、Conceptが入ることもありC++20はそれ以上に大きな変革になりそうです。

2014年5月29日木曜日

[C++] next_combination() 書いてみた

C++ 標準ライブラリの <algorithm> には next_permutation() という関数があります。辞書順で順列を列挙してくれる大変便利な関数で競プロ関連でお世話になったことのある人も結構いるのではないかと思います。順列があるなら組み合わせはないの?ということで next_combination() で検索をかけると nyama++ さんの 全ての組み合わせを作る【C++】 - Programming Magic という記事が引っかかります(他の検索結果からも結構参照されているようです)。この記事のリンクで指し先が消失しているものがありますが、内容自体は ISO C++ 標準化委員会(JTC1/SC22/WG21) に提案されているペーパー n2639 Algorithms for permutations and combinations, with and without repetitions です(実装は Reference implementation のところ)。next_pertial_permutation() (next_permutation() が n! の列挙なら、これは nPr を列挙する) の実装が

template <class BidirectionalIterator>
bool next_partial_permutation(BidirectionalIterator first, BidirectionalIterator middle, BidirectionalIterator last)
{
  reverse(middle, last);
  return next_permutation(first , last);
}

というのも私的に感動ポイントですが、next_combination() は確かに訳が分かりません。ってことで一度自分で書いてみたところ最終的にほぼ同様のアルゴリズムとなり理解が進みましたので解説記事っぽいものを書いてみようと思います。なお図としては全ての要素が1ずつ異なるようなイメージで書いていますが、値が跳んでいてもいいですし重複があっても問題ないはずです。

まず最初に意味不明と書かれていた prev_combination() について。処理上何かと便利なので [first, middle), [middle, last) がソートされている、という前提を置いておきます。このとき、[first, middle) に対して next_combination() した場合の例は以下のようになります。

見ての通り、[middle, last) に対しては prev_combination() になっています(両方の領域がソート済み、かつ、辞書式順序で列挙なのでこうなる)。ということで先の前提の下で next_combination() が正しく実装できれば他方の領域に対して next_combination() をかけてやれば元々の領域に対する prev_combination() になります。

さて、では本論のコードについて。さすがに真っ白の状態だと書けないので n2639 で参考文献として挙げられている Knuth 先生の The Art of Computer Programming. Volume 4, Fascicle 3 GeneratingAll Combinations and Partitions を見ると(※立ち読み。分冊が一冊になったら買います、多分)確かに列挙アルゴリズムが記述されてはいますが添字の列挙としてのアルゴリズムであって実際に要素の並び替えを行う next_permutation() 風の next_combination() の場合だとそのままでは使えなさそうです。そうすると指針としては、汎用的な方法として挙げられている「最も右端の増加可能な要素を見つけて増加、以降の要素は設定可能な最小値を設定していく」になります。これを実装してみましょう。

まずは「最も右端の増加可能な要素」を見つける必要があります。右端(middle - 1)が最大値でない場合はそこが増加可能です。では右端ではない部分が「最も右端の増加可能な要素」となる場合を考えてみましょう。この時その要素(*targetとしておきます)の右側には、右端(*(middle-1))が最大値となる状態で連続して昇順に要素が並んでいることになります。さもなくば、より右側の要素が増加可能になるからです。[middle, last) もソート済であることと合わせると、*target < *(last-1) となる最も右端の target を見つければいいことになります。これが存在しない場合には列挙終了です。

次に、*target をどの値に変更すれば良いかを考えましょう。target の右側が最大値から連続で詰まっているので、[target, middle) の間の値は考慮する必要がありません(target から middle まで埋めるには要素が足りない)。ということで、[middle, last) のうちで *target より大きくなる最小値に設定すれば良いことになります(next とします)。*target < *(last-1) なので [middle, last) で必ず見つかります。

では、[target+1, middle) (と [next, last) をどのように埋めれば良いでしょうか? これまでの条件から [target, middle) [next, last) 近傍の順序は下図上段のようになっています。[target, middle) はできるだけ最小になるように設定すること、[middle, last) もソートされている必要があることを考えると最終結果は下図下段のようになっている必要があります。これは target と next を iter_swap した後、[target+1, middle), [next+1, last) が一つの領域だとみなして、next+1 が先頭となるように rotate() することに他なりません。

というのをコードにまとめると以下のようになります。C++er な人には言うまでもないとは思いますが、わざわざ !(a < b) みたいな書き方になってるのは operator< のみ使用にしたかったからです。なお、rotate() は en.cppreference.com の参考実装(Forward Iterator 向けですが)に分割2領域用にちょっとだけ修正したものです。

// possible implementation introduced at http://en.cppreference.com/w/cpp/algorithm/rotate with slight modification to handle parted ranges
template<typename FI>
void parted_rotate(FI first1, FI last1, FI first2, FI last2)
{
 if(first1 == last1 || first2 == last2) return;
 FI next = first2;
 while (first1 != next) {
  std::iter_swap(first1++, next++);
  if(first1 == last1) first1 = first2;
  if (next == last2) {
   next = first2;
  } else if (first1 == first2) {
   first2 = next;
  }
 }
}

template<typename BI>
bool next_combination_imp(BI first1, BI last1, BI first2, BI last2)
{
 if(first1 == last1 || first2 == last2) return false;
 auto target = last1; --target;
 auto last_elem = last2; --last_elem;
 // find right-most incrementable element: target
 while(target != first1 && !(*target < *last_elem)) --target;
 if(target == first1 && !(*target < *last_elem)) {
  parted_rotate(first1, last1, first2, last2);
  return false;
 }
 // find the next value to be incremented: *next
 auto next = first2;
 while(!(*target < *next)) ++next;
 std::iter_swap(target++, next++);
 parted_rotate(target, last1, next, last2);
 return true;
}

// INVARIANT: is_sorted(first, mid) && is_sorted(mid, last)
template<typename BI>
inline bool next_combination(BI first, BI mid, BI last)
{
 return next_combination_imp(first, mid, mid, last);
}

// INVARIANT: is_sorted(first, mid) && is_sorted(mid, last)
template<typename BI>
inline bool prev_combination(BI first, BI mid, BI last)
{
 return next_combination_imp(mid, last, first, mid);
}

基本的に n2639 とほぼ同じことをするコードです。nyama++ さんの記事で②③となっているところが実は合わせて rotate() だったわけです。

reverse(first, mid);
reverse(mid, last);
reverse(first, last);

という reverse() 3連打の rotate() 実装を 2 領域対応にした形が②③になります。後の違いは、false を返す時に先に return する(rotate() が 2 ヶ所になる)か、1ヶ所で rotate() するかくらいですね。

ついでなので、計算量(iter_swap の呼び出し回数)について実際に計数させてみたところ、2nCn 列挙時での 1 next_combination() 辺りの平均 iter_swap() 回数は 1.7 くらいに収束しそうな感じです。 もちろん組み合わせの数自体が指数関数的に増えるので全体的にはすぐに苦しくなるのですが、かなり悪くない数字に感じます。

#n,r,the number of iter_swap,the number of combinations,average iter_swap per one next_combination() call
2,1,2,2,1
4,2,8,6,1.33333
6,3,30,20,1.5
8,4,112,70,1.6
10,5,414,252,1.64286
12,6,1540,924,1.66667
14,7,5754,3432,1.67657
16,8,21656,12870,1.68267
18,9,81994,48620,1.68643
20,10,312068,184756,1.68908
22,11,1192954,705432,1.6911
24,12,4577356,2704156,1.69271
26,13,17619034,10400600,1.69404
28,14,68003992,40116600,1.69516
30,15,263097002,155117520,1.69611
32,16,1019997844,601080390,1.69694
34,17,3961678122,2333606220,1.69766
36,18,15412309480,9075135300,1.6983
38,19,60046904394,35345263800,1.69887
40,20,234252753696,137846528820,1.69937