競プロをはじめた家事手伝いロボットのブログ

競技プログラミングをしている家事手伝いロボットのブログです

凍ってこまったこと #6

凍りました(もう凍ってから 1.3 年くらい経ってそうですね)。みなさんは今は Twitter (現 X) にいるのでしょうか。他の SNS への流出もありますか?まとめます。

実質 ARC175-E の解法記事なので、タイトル詐欺かもしれません

Twitter (現 X) 言及?

どしどし

コンテスト感想

楽しかったです(楽しかったので) 今回 A と C に五千年かけていたので E が通っても負けだと思っていたんですけど(C までを通したとき (85:53 (+1)) に D と E をちらっとみて E のほうが解けそうだったので E に特攻)、E があんまり通されていなかったのであたたまりでした

C はひとしきり迷走したあとにすっきりした解法になったので楽しかったです 実装もお気に入り

C の実装

A は迷走したまま最後まで走り抜きました 実装もお気に入らず(ぽか) 問題(と想定解)はすきです

D と F は読んでません(ぽか)

E の解法

あんまり他の人が書いてなさそうな解法で解いていても Twitter にいないと気軽に書けなくてこまります(こまるので)

次のような図形を構築することを考えます(ここで、\(N\) は与えられる \(N\) と関係ありません。制約を満たす任意の \((N,K)\) で構築できるので、なるべく詰めておけば入力の \(N\) は見なかったことにしてよいためです)。

N×N の正方形と、1×S の長方形と、S×1の長方形を 1 つずつと、1×1 の正方形をたかだか 1 つ作る
構築の方針

どのような \(K\) も \(N ^ 2+2S\) もしくは \(N ^ 2+2S+1\) として書けるので、上の形で構築ができればよいです。

ぐっと睨むと、次の \(3\) 集合で \(N ^ 2+2S\) が実現できることがわかります。

\[\lbrace (x,y,z)\mid x+y+z=S-2\wedge0\leq x,y,z\lt N\rbrace\] \[\lbrace (x,y,z)\mid x+y+z=S-1+N\wedge0\leq x,y,z\leq N\rbrace\] \[\lbrace (x,y,z)\mid x+y+z=S-1+2N\wedge0\leq x,y,z\lt N\rbrace\]

お気持ちとしては、\(\lbrace (x,y,z)\mid x+y+z\equiv S-1\pmod{N}\wedge0\leq x,y,z\lt N\rbrace\) に対して、\(\lbrace (x,y,z)\mid x+y+z=S-1+N\wedge N\in\lbrace x,y,z\rbrace\rbrace\) を付け加える代わりに \(\lbrace (x,y,z)\mid x+y+z=S-1\wedge0\leq x,y,z\lt N\rbrace\) の部分を \(1\) つずらす感じです(発想としては、\(S=1,2\) の構築を考えることで出ました。最初から見通しよくこの形だったわけではなくて、降順の列を rotate してがちゃがちゃしていました)。

\((N,N,N)\) は使われないことがわかる(\(S\leq N\) より、\((N,N,z)\) は入れることができない)ので、\(1\) として \((N,N,N)\) をいつでも入れることができます。

実装は以下のようになります。

#include <iostream>
#include <cmath>
#include <ranges>

int main() {
    using namespace std;
    // 和が sum で各座標が 0 以上 limit 未満である (x, y, z) をすべて出力
    constexpr auto sum_triplet{[](const unsigned sum, const unsigned limit){
        for (const auto x: views::iota(0U, limit))
            for (const auto y: views::iota(0U, limit))
                if (const auto z{sum - x - y}; z < limit)
                    cout << x << " " << y << " " << z << endl;
    }};
    const unsigned K{*next(istream_iterator<unsigned>{cin})}, N(sqrt(K)), S{(K - N * N) / 2};
    sum_triplet(S - 2, N);
    sum_triplet(S - 1 + N, N + 1);
    sum_triplet(S - 1 + N + N, N);
    if (1 & (K + N))
        cout << N << " " << N << " " << N << endl;
    return 0;
}

凍ってこまったこと #5

Twitter 言及

こちらこそです お礼に

これを書いておきました(奪いぽかいこ)

これ、計算量の規定が見当たらなかったんですよね... 直前で \(O(\log w)\) と言ってしまっているのでそういう感じの実装にしたんですけど、関数があることには言及しておいたほうがよかったかもしれません

おわり

おわりです。最近は bluesky っていうところにみなさんがいたりするんでしょうか?

凍ってこまったこと #4

タイトル詐欺じゃないですか?新鮮なぽかいこをお届け

ぽか

誤差です(誤差なので)

今週の ABC は私の担当じゃなかったので見ていなかったんですけど(ふまじめぽかいこ)、誤差トピックが出ているらしいと話を聞いて飛び起きて問題を解きました

$\dfrac1D\sum x _ i ^ 2-\left(\dfrac1D\sum x _ i\right) ^ 2$ をそのまま計算すると困って、$D ^ 2$ 倍した整数を求めてから $D ^ 2$ で割ると困らない話をします。 特に、割り算のあたりから double とか long double で計算していると思って話をします。

みなさんももう浮動小数点数の演算と誤差に慣れてきたので困る・困らない理由は説明できちゃいますか?そうかも

そのまま計算すると困るのはどうして

$\dfrac1D\sum x _ i ^ 2$ には少なくとも $\dfrac12\operatorname{ulp}$ くらい相当の相対誤差が乗っています(正規化数が動く範囲を超えることはないですが、ちょうどの値は表せないかもしれないので)。 double だと $10 ^ {-15}$ くらい、long double だと $10 ^ {-18}$ くらいだと思います。 $\dfrac1D\sum x _ i ^ 2$ の最大値は小さく見積もっても $10 ^ {16}$ くらいの値になるので、最悪ケースでは $10$ とか $10 ^ {-2}$ くらいの絶対誤差が生じうることになります。

$\left(\dfrac1D\sum x _ i\right) ^ 2$ についても同様です。

すると、求める答え $\dfrac1D\sum x _ i ^ 2-\left(\dfrac1D\sum x _ i\right) ^ 2$ のスケールがある程度小さいとうまくいかない気がしてきますね。 実際、$10 ^ {-2}$ くらいの絶対誤差のある値どうしの演算結果が $10 ^ 4$ 以下の値になった場合、相対誤差が $10 ^ {-6}$ 以下であることを保証するのは難しいです。

逆に、演算結果が $4.5\times 10 ^ 6+\varepsilon$ 以上なら OK とかも言えそうな感じがしますよね($\dfrac1D\sum x _ i ^ 2$ や $\left(\dfrac1D\sum x _ i\right) ^ 2$ は $2.25\times10 ^ {18}$ とかでおさえられるはずなので) 実際、答えが $4.5\times 10 ^ 6+10$ 以下っぽいときは abort するようにすると、それ以外のケースでは AC することを確認できます(このことはこの評価がバグっていないことを意味しません。誤差が偶然いいかんじになっているかもしれないので)。

atcoder.jp

整数にすると困らないのはどうして

$D\sum x _ i ^ 2-\left(\sum x _ i\right) ^ 2$ の正しい値が求まっているとき、これを浮動小数点数にキャストする際には誤差が $\dfrac12\operatorname{ulp}$ くらい乗ります。

これを $D ^ 2$(double で正しい値が表現できる)で割るので、誤差が $\operatorname{ulp}$ くらいになって、double や long double を使っている場合はもちろん十分小さいので、困りません。

おわり

おわりなので。 最近浮動小数点数トピックがおおめな気がしてうれしいです かなしいのは凍っているので浮動小数点数トピック TL を見逃してしまいがちなこと

かなしいこと

なかまいり(最悪)

凍ってこまったこと #3

凍りました。最近はあんまりこまっていません? まとめます。

こまったこと

Twitter (現 X) 言及

現代の競技プログラミングは ranges::max で常勝!(素振り)

Twitter 以外言及

rsk0315.hatenablog.com

浮動小数点数のはなしがされつつあり、うれしいです うれしいので

こまってないこと

インターネットをはじめてしまう

mstdn.maud.io

はじめました。はじめてから

misskey.kyoupro.com

の存在を知りました(ぽか)。

サブこまってないこと

Twitter (現 X) と同じような広い海には漕ぎ出せない感じがするので、見ないほうがいいインターネットを見ることはすくないです。

サブこまったこと

まだみなさんが Twitter (現 X) にいるので、むにゃむにゃむにゃ…

凍ってこまったこと #2 (WIP)

凍りました。こまりました。まとめます。(2) これから追記するかもしれません?

こまったこと

ひとのライブラリにバグを見つけたとき

https://dr0gsk0l.github.io/library/superstd/Multiset.cppdr0gsk0l.github.io

これの erase_k

  void erase_k(const T&a,int k){
    if(map<T,int>::count(a))return;
    at(a)-=k;
    if(at(a)<=0)erase(a);
  }

if(count(a))return したらよくない気がします

こういうとき、Twitter がないと言いにくいですね。 github でプルリクエストを出したらいいんじゃないですか?ぽか。

rate limit にひっかかったとき

なんだかすぐに制限がかかる気がします 気のせいですか?そうかも?

ぽか

書きません(??) 編集前のコードはこんな感じでした

#include <iostream>
#include <vector>

int main() {
    using namespace std;
    unsigned N, T, M;
    cin >> N >> T >> M;

    // hate[i] の j ビットめが 1 ⟹ i 番目の選手と j 番目の選手の相性が悪い (0-indexed)
    vector<unsigned> hate(N);
    for(unsigned i{}, a, b; i < M; ++i){
        cin >> a >> b;
        hate[--b] |= 1U << --a;
    }

    // 再帰関数で計算した値を出力
    cout << [dfs{[N, T, &hate](auto&& f, vector<unsigned>& teams, unsigned now) -> unsigned {
        // 最後の選手まで見て T チームに分かれていれば OK
        if(now == N)
            return size(teams) == T;

        unsigned ans{};

        // すでにあるチームに now 番目の選手を追加する
        for(auto&& team : teams)
            // チームに now 番目の選手と相性の悪い選手がいなければ
            if(!(team & hate[now])){
                team ^= 1U << now;
                ans += f(f, teams, now + 1);
                team ^= 1U << now;
            }

        // チーム数に余裕があるとき、新しいチームを作る
        if(size(teams) < T){
            teams.emplace_back(1U << now);
            ans += f(f, teams, now + 1);
            teams.pop_back();
        }

        return ans;
    }}, T]{
        vector<unsigned> team;
        team.reserve(T);
        return dfs(dfs, team, 0);
    }() << endl;

    return 0;
}

「無名関数再帰をキャプチャで受けてどうこうするのはイディオムすぎる」「unsigned じゃなくて int でもよくない?」「整数の初期化に一様初期化を使う意味はない」などのご意見があります。

  • 無名関数で再帰をする
    • 実際かなり知らないと読みにくそうな気がします。ごめんなさいです。グローバルに変数とか関数とかを置いて再帰を書くことをしないまま時を過ごしているのでラムダ式再帰を書きがちなんですけど、やさしくないな〜と言われると返す言葉もございません(ぽか)。
  • unsigned じゃなくて int でもよい
    • まあそうなんですけど、そうなんですけど、むん。符号なしでいいなら符号なしでよくないですか?オーバーフローも未定義にならないですし。
  • 整数の初期化に一様初期化を使うな
    • これはわたしもなんで使ってるのかわからないんですけど、vector も set も modint も {} で初期化するのに慣れちゃいました。統一感。= 0 と打つより shift を押しながら中指を滑らせて {} と打鍵するほうが速いです(?)。

特に A-B 問題ではなるべく癖を抜くようにしてるんですけど、C 問題以降だとなかなか抜けていないことも多いかんじがします。 気をつけます。 ところで今日の B 問題の解説は入れ子になった any_of があるんですけど、ちゃんと癖抜けてますか??気をつけます(でもこっちも意味論的に読みやすい気がするんですけどね(?)。std::ranges::any_of が使えるようになったらもっと読みやすくなる気がするんですけど(??)。for 文とにゃんにゃんするプログラムも併記するかもしれません)

ちゃんと文を読むだけでわかってもらえる解説が書けるようになりたいですね。〜完〜

制限を満たしている状態で writer になったあと制限を満たさなくなると、こうなります(ぽか)。

mod 998244353 で有理数を答える問題で出た答えが有理数で何か知りたいですね

ですね。

クイズです。$p/q\equiv 709787742\pmod{998244353}$ を満たす互いに素な正整数のペアはなんでしょうか。

答えは $(p,q)=(25,30233088)$ です。 $(p,q)=(23583,26948)$ だと思った人もいるかもしれませんが、サイコロを $10$ 回だけ振ってちょうど $9$ 回だけ $1$ の目が出る確率 $\pmod{998244353}$ を想定していたので、ちがいます(ちがいはしませんけど)。1

こんなふうに、競技プログラミングの問題を解いていて出てくるおおきめ有理数 $\pmod p$ の値からもとのおおきめ有理数がわかると、サンプルがわかりやすくなったりどこで壊れているのか確認しやすくてうれしいですね。 知る方法について書きます。

元ネタ?はこちらです。いっしょにどうぞ?

rsk0315.hatenablog.com

主に C++競技プログラミングをしている人に向けて書いています。 別の言語でプレイしている人がどういう感じになるのかはわかりません。ぽか

まとめ

最初から有理数で計算してしまえばいいですね。

日記要素ですか?

ちがいます。

まじめなはなし

modint などを使って dp をしたりするとき、boost::rationallong double と併用したり切り替えたりできると、サンプルを試してみたときに「なんだか逆数を求めていた」とか、「なんだか確率じゃなくて期待値を計算していた」みたいなミスに気づきやすそうですね。

わたしはこんなときには main 関数の頭に

int main(){
    using modint = atcoder::static_modint<998244353>; // ふつうの modint です
//    using modint = boost::rational<boost::multiprecision::cpp_int>; // 有理数です
//    using modint = long double; // 浮動小数点数です

のように型エイリアスをおいています。 これのコメントをつけはずししながらサンプルを実行してみてデバッグをしているというわけですね。 この型エイリアスの名前が modint なのはどうかと思ってたまに dp_value_type とか ans_type みたいにすることもあるんですけど、結局戻ってきてしまっています。 コメントつけはずしなのもどうかと思ってたまに template<typename modint> みたいにした solve 関数を書くこともあるんですけど、結局戻ってきてしまっています。 もうちょっといい名前とか便利な切り替え方法をご存知の方はおしらせください。

きをつけていること

raw とか pow とか inv とか、上の 3 つのうち atcoder::static_modint にしかないものを使っているとコンパイルエラーになりがちです。 使わないようにするか、もう 2 つについても同じ名前の関数を作ってしまうことでなんとかなるかもしれません。 抽象化した関数を作ってあげて関数の中でこれらの違いを吸収してあげてもいいかもしれないですね。 val は入出力の operator<<(std::ostream&, modint) みたいなのをつくっておくほうがおきにいりです。

また、boost::rational<boost::multiprecision::cpp_int> を使っているときに、適当に計算した結果を auto で受けたりすると型が変になったりすることがあるので、modint で受けていい感じにしておくことがだいじかもしれません。

またまた、boost::rational<boost::multiprecision::cpp_int> で分子や分母の値がおおきくなっていくと、遅くなります。 こういうときはあきらめて long double などを使いましょう。

atcoder.jp

たとえばこういうときですね。$10 ^ 8$ 桁くらいになりそうな感じがします。ならないかも。

おわり

Q. ちょっと文章がてきとうじゃないですか?A. 日記要素なので。

Q. ei1333 の日記さんは文章がてきとうというわけではないんじゃないでしょうか?A. ぽか。


  1. ここで確率を想定していたのに出てくる有理数が $1$ より大きいですとかだともうちょっといい感じでしたよね。最初に思いついた有理数が $5\times10/60466176$ だったのでごようしゃいただけるとぽかです。

浮動小数点数の話シリーズ #n

これはなに

お気持ち表明ではないです。

これをみたので。

比較的ねむいのでこのツイートに関する話ではないかもしれません。

以下では「iff で判定可能」というのを、ある数 $x$ が与えられたときにある数 $\alpha$ からの相対誤差もしくは絶対誤差が $\varepsilon$ 以下であるか否かを誤差なしで判定できることとします。

たぷ

解の真の値が 0 である場合、0.0 を AC としないとワクワクになります。1

まじめな話

浮動小数点数を出力しようとすると、符号?十進数列(.十進数列)?(e符号?十進数列)? みたいな形になることが多いと思います(ランダムにビット列を取ってくると nan のことも多いですか?そうかも)。 この形式の値は有限小数なので有理数として表すことができますね。

有理数として値が得られれば、ある代数的数 $\alpha$ からの絶対誤差もしくは相対誤差がある代数的数 $\varepsilon$ 以下であるかを厳密に判定することができますね($\alpha$ を根にもつ多項式はあるとうれしいですね そもそも一般には $\alpha$ を持っておくために必要な気がしますけど)。

atcoder.jp

たとえばこれとかは iff で判定可能な問題です。

でも現実的に iff で判定可能な問題は少なくて、

atcoder.jp

この問題は解の真の値は代数的数になりますが、最小多項式の次数が $2 ^ {\Theta(N)}$ くらいになりそうな気がします。 証明していないのでふんわりしておきます。

解が代数的数にならない問題は iff で判定可能と言いにくいような気がします。 問題を準備する際に十分時間をかけ、許容範囲の上限と下限の十進展開をOLEぎりぎりまで評価してしまうくらいしか有効な手立てがないような気がします。 この手段を取ったとしても、適切な問題とテストケースを用意するとどれだけの精度で計算を始めればいいのかわからないような例が作れるような気がします(OLEぎりぎりの直下が $000000\ldots$ なのか $999999\ldots$ なのか(どちらも長い有限桁ののち $0$ や $9$ でない桁が出現するものとします)の判定が現実的には終わらないような場合があるかもしれません)。
そうでなくても、$10 ^ 8$ 桁以上正しいことを保証しながら計算を行うのは非常に骨が折れそうですね。 1ケースに1時間くらいかかると50ケースでは2日以上かかるので、計画的に準備をすることでなんとかなる場合があるかもしれません。

$10 ^ 8$ 桁までは OLE にならなさそうな話。

そもそも $10 ^ 8$ 文字の文字列どうしを浮動小数点数として大小比較するのも多少時間がかかりそうな気がしますね。

本題

1e100 とか -3.14e-3 みたいな値から有理数を得る関数を書いてみました。 構文解析には自信がないロボットなので、みなさまの助力をおまちしています(?)。

$\dfrac ab\times10 ^ c\ (10\nmid a,b)$ みたいな形式で表すのもいいかなと思ったんですけど、boost::rational の恩恵にあずかりたかったので。

    using big_integer = boost::multiprecision::cpp_int;
    using big_rational = boost::rational<big_integer>;

    const auto decimal_floating_point_number_literal_to_rational{[](const string& s) -> optional<big_rational> {
        vector<pair<unsigned long, unsigned long>> sub_idx(5, {size(s), size(s)});
        unsigned long state{}, index{};
        bool positive{true}, dup_sign{false};
        const auto change_state{[&sub_idx, &state, &index, &dup_sign](unsigned long next_state){
            dup_sign = false;
            for(; state < next_state; ++state){
                sub_idx[state].second = index;
                sub_idx[state + 1].first = index;
            }
        }};
        sub_idx[0].first = 0;
        for(const auto c : s){
            if(state == 0){
                if(!dup_sign && (c == '-' || c == '+')){
                    if(c == '-')positive = false;
                    dup_sign = true;
                }else if('0' <= c && c <= '9')change_state(1);
                else if(c == '.')change_state(2);
                else return {};
            }else if(state == 1){
                if('0' <= c && c <= '9');
                else if(c == '.')change_state(2);
                else if(c == 'e' || c == 'E')change_state(3);
                else return {};
            }else if(state == 2){
                if('0' <= c && c <= '9');
                else if(c == 'e' || c == 'E')change_state(3);
                else return {};
            }else if(state == 3){
                if(!dup_sign && (c == '-' || c == '+'))dup_sign = true;
                else if('0' <= c && c <= '9')change_state(4);
                else return {};
            }else if(state == 4){
                if('0' <= c && c <= '9');
                else return {};
            }
            ++index;
        }
        if(sub_idx[2].first != sub_idx[2].second)++sub_idx[2].first;
        if(sub_idx[3].first != sub_idx[3].second){
            ++sub_idx[3].first;
            if(s[sub_idx[3].first] == '+')++sub_idx[3].first;
        }
        sub_idx[3].second = sub_idx[4].second;
        sub_idx.pop_back();

        big_integer exponent{sub_idx[2].second - sub_idx[2].first};
        exponent *= -1;

        while(sub_idx[1].first != sub_idx[1].second && s[sub_idx[1].first] == '0')++sub_idx[1].first;
        if(sub_idx[1].first == sub_idx[1].second)while(sub_idx[2].first != sub_idx[2].second && s[sub_idx[2].first] == '0')++sub_idx[2].first;
        while(sub_idx[2].first != sub_idx[2].second && s[sub_idx[2].second - 1] == '0'){
            --sub_idx[2].second;
            ++exponent;
        }
        if(sub_idx[2].first == sub_idx[2].second)while(sub_idx[1].first != sub_idx[1].second && s[sub_idx[1].second - 1] == '0'){
            --sub_idx[1].second;
            ++exponent;
        }

        big_integer numerator{1}, denominator{1};
        if(!positive)numerator *= -1;
        numerator *= big_integer{s.substr(sub_idx[1].first, sub_idx[1].second - sub_idx[1].first) + s.substr(sub_idx[2].first, sub_idx[2].second - sub_idx[2].first)};
        const auto pow_ten{[](big_integer n) -> big_integer {
            big_integer r{1}, a{10};
            while(n){
                if(n & 1)r *= a;
                a *= a;
                n /= 2;
            }
            return r;
        }};
        if(sub_idx[3].first != sub_idx[3].second)exponent += big_integer{string_view{s}.substr(sub_idx[3].first, sub_idx[3].second - sub_idx[3].first)};
        if(exponent < 0)denominator *= pow_ten(-exponent);
        else numerator *= pow_ten(exponent);

        return big_rational{numerator, denominator};

    }};

これを使ってジャッジするの、遅くなるだけでうれしくなる人はそんなに多くないと思ってますけど、お好みでお使いください(?) バグってたらごめんなさいです

ところで

こう書いてもいいと思いますけど、この条件を満たすかの判定が煩雑になるわりに問題やジャッジの振る舞い/コンテスタントが得る情報としてはあんまり増えないような気もします… どうでしょう? こう書くのがダメではないので、こう書いてもいいと思います(想定解からの誤差よりよっぽどよいです なぜならば想定解からの誤差は悪いので)

おわりに

洗濯機型音楽ゲームでぽかいこと勝負!

負けました(なんですか?)

🍑🍝🦃(ぽかパスタターキー(982766353型文字列))


  1. ほんとうに?数え上げで答えが 0 のときに 0.0 を AC にしなくても文句は出ないと思います