rickytheta's blog

競プロ

AlienDPを理解するだけのログ

理解するだけのログシリーズ

ざっくり

ABC218 HがAlienDPの説明として非常に詳しいです:

https://atcoder.jp/contests/abc218/editorial/2621

問題の整理

  • 個数制約付きの最大化問題がある
  • 個数 → スコア のグラフが上に凸になっている (徐々に伸びが鈍くなる)
    • 個数を変数 x と置いたときのスコアを f(x) と書く
    • ぱっと見で貪欲っぽい雰囲気の問題がこうなりがち
  • 求めたいものはある定数 K=O(N) に対して f(K) だが、普通にDPをすると O(N^2) などかかる
    • 「個数がぴったり K になる」という条件を考えようとすると、DPのキーに個数を持つことになり、更新が O(N) 回なのでトータルで O(NK)=O(N^2) となる

緩和問題と幾何学的解釈

  • 個数を気にする代わりに「1個使うと定数 C だけペナルティがかかる」という緩和問題を考える
    • 元の問題が \max f(K) だったのに対し、 \max f(x)-Cx となっている
    • 最終的な解に使われた個数制約に気をつける必要がなくなり、スコアさえ管理できれば良いので、 o(N^2) で解ける(ことがある)
  • この時の最適解が x=x' に対して f(x')-Cx'=B と求まったとする
    • 移項すると f(x')=Cx'+B となり、これは凸グラフ y=f(x) と直線 y=Cx+B の交点 (x',y':=Cx'+B) を導く
  • B は最大化問題の解なので、この交点 (x',y') は凸グラフ f(x) と交点を持つような傾き C の直線の交点として一番切片の大きい直線上にある
    • 凸グラフと直線の交点であることを考えて、より具体的に書くと、
    • 解周りの3点を P_{-1}(x'):=(x'-1,f(x'-1)), P_0(x'):=(x',f(x')), P_{+1}(x'):=(x'+1,f(x'+1)) として、
    • 線分 P_{-1}(x')P_{0}(x') の傾きを \Delta_-f(x'):=f(x')-f(x'-1)、線分 P_0(x')P_{+1}(x') の傾きを \Delta_+f(x'):=f(x'+1)-f(x') と定義すると、
    • \Delta_-f(x')\ge C \ge \Delta_+f(x') が成り立つ
    • (最大化のイメージ図)

二分探索の構成

  • C をうまく選んで緩和問題を解くことで、 f(K) を求めることを考える
    • \Delta_-f(K)\ge C \ge \Delta_+f(K) となる傾き C を見つけられれば良い
    • この [\Delta_-f(x),\Delta_+f(x)] という範囲は、上に凸なので x に対して広義単調減少 (等号がありえる点に注意)
      • \Delta_+f(x-1)=\ \Delta_-f(x)\ge\Delta_+f(x)\ =\Delta_-f(x+1)
    • ある C に対して \Delta_-f(x)\ge C\ge\Delta_+f(x) を満たす値 x の範囲 [L(C),R(C)] は緩和問題を解くときについでに計算可能(スコアが同一になった際に最小個数・最大個数の2つを保持する)
      • \Delta_-f(L(C))>C>\Delta_+f(R(C)) が成り立つ
      • (以下図であれば L(C)=7,R(C)=10 となる)
    • 元の問題の答えは緩和問題の答えが B=f(x)-Cx であることを考えると f(x)=B+Cx と復元できる
    • つまり、最小個数・最大個数を管理したDPを行うことで、凸包上である傾きの直線が触れるような領域の x の範囲、及びその時の評価値 f(x) の範囲を計算することができる
      • これをクエリとして用いて二分探索を構成する
  • 各 C \ge 0 に対して二分探索することで、 L(C')\le K\le R(C') なる C' を探す
    • 適切な C_{low}\lt C_{high} からスタートし、 L(C_{low})\le R(C_{low})\ \lt K\lt\ L(C_{high})\le R(C_{high}) を維持するように更新する
    • (以下図であれば [L(C_{low}),R(C_{low})]=[0,2],\ [L(C_{high}),R(C_{high})]=[10,\infty]; N=19 とすれば右端は19で切れる)

応用

  • 上記は f(K) を求めるために x の範囲を狭めているが、同時に f(x) の範囲も求められるため、 f(x)=V となるような x を求める、 f(x)\le V となる最大の x を求める、ということも可能
    • 上図であれば [L(C_{low}),R(C_{low})] = [0,2] の時 [f(L(C_{low})),f(R(C_{low}))]=[0,16] のような形
    • f:\mathrm{cost}\mapsto \mathrm{score} だと凸性(連続性)が成り立たない場合、 f の代わりに f':\mathrm{score}\mapsto \mathrm{cost} と言い換えて立式する、としても通せることがある
    • ちなみに f(L\le x\le R) の最大値、という問題は凸性から f(L) または f(R) に帰着されるので、この形は考えなくて良い

提出例

その他

  • 制約の代わりにペナルティ項を添加するという考え方はラグランジュ緩和/ラグランジュ双対が有名
  • LPにおけるラグランジュ緩和・ラグランジュ双対はその先の「任意のペナルティ係数における最大値の最小化問題」を考え、それが元の問題の最適解と一致する、というところが本質
    • 強双対性
  • LPの実行可能解は凸包を成すので、緩和問題の解を考えるのはある傾きの半平面で凸包の表面をなぞる……というところまではイメージは同じ
  • AlienDPは「この傾きの範囲をコントロールすることで、凸包のある条件を満たす領域を二分探索で求める」というのが本質
    • 基本的に凸包は二次元空間で、凸包全体に興味がなく(そもそも凸包でなく凸関数)、一部の定義域または値域にのみ興味がある
  • LPのラグランジュ双対は「この傾きをすべて試して凸包の表面を全部なぞると、元の問題の大域解と一致する」というのが本質
    • 一般に幾何学的解釈はあまり意味がなく、緩和(双対)された結果の問題が既存のグラフ問題で解けるような形になることが重要
  • 制約がスコアに移動して解ける形になる、緩和に使うペナルティ係数の範囲を二分探索/三分探索等でコントロールする、という点では似ている
    • AlienDPはある条件を満たす範囲を求める目的のため
    • LPのラグランジュ双対は \min\max f の \max f がペナルティ係数に対して凸性を持つため(確か)
  • LPの話は以下を参照
  • 本来のAlien Techniqueはグラフの辺の重みにペナルティを添加して、グラフ問題の解に使われる辺の数をコントロールするテクニックを指していたらしい
    • noshi91.github.io
    • グラフ上でこれを行う場合、上記ならmonge性、ICPC2015国内予選Fならマトロイド性、といった「辺の数に対してスコアが凸」という良い性質が必要となる
    • 上記ページではラグランジュ双対として定式化されているが、本ブログやABC解説で書かれている形式は上記ページの「その他」にあたるものと思われる (相互に解釈を変換ができるかは不明……？)
    • (ICPC2015、懐かしすぎますね)

Define-By-Run風味の全探索ライブラリを作る

これ何？

Define-By-Run とは

ChainerやOptunaで使われているAPIの作られ方……という認識

例えばChainerであれば「入力画像/ベクトル にいくつか層を適用して出力を計算する」という(ようにみえる)コードを書くことで、内部的にネットワーク構造が作られ、その学習を行うことができる(確か)

例えばOptunaであれば「乱数(っぽいもの)を振って計算し、それを利用して評価値を出力する」という(ようにみえる)コードを書くことで、内部的にその評価値の分布をTPEやCMA-ESにかけて最適化できる(確か)

特にOptunaは例えば「最初に乱数でステップ数を決定して、そのステップ数に応じて処理のクオリティや複雑度を変える」といったことも可能

要は「手続き的に書いているだけで、そのコードの構造や性質を具体的に定義せずとも、ライブラリ側で抽出できる」という手法……と考えられる

今回はこの考えをもとに「Define-By-Run風味の全探索ライブラリ」を作ります。

ただこれをDefine-By-Runと名付けて良いものか正直怪しいので、以降はすべてDBRと略記し区別することにします……(すでに名のあるテクニックだったら教えてください)

作りたいもの

例えば以下の問題を考えましょう(つい最近の問題なのでupsolveしたい人は気をつけてください、このページでは解法には触れません):

https://atcoder.jp/contests/arc168/tasks/arc168_c

小さいケースで検算するために、操作後にありえる文字列を全探索することを考えるとDFSを書きたくなります:

vector<string> bf;
void dfs(string s, int k) {
    bf.push_back(s);
    if (k == 0) return;
    int n = (int)s.size();
    for (int a = 0; a < n; ++a) {
        for (int b = 0; b < a; ++b) {
            swap(s[a], s[b]);
            dfs(s, k-1);
            swap(s[a], s[b]);
        }
    }
}
void main() {
    int n,k;
    string s;
    cin>>n>>k>>s;

    dfs(s,k);
    sort(ALL(bf));
    bf.erase(ALL(bf), bf.end());
}

これ結構冗長だな……って思いますよね、それでこうです:

void main() {
    int n,k;
    string s;
    cin>>n>>k>>s;

    auto f = [&](DBRSource& source) {
        string t = s;
        int swap_num = source.choice(0,k);
        for (int times = 0; times < swap_num; ++times) {
            int a = source.choice(1,n-1);
            int b = source.choice(0,a-1);
            swap(t[a], t[b]);
        }
        return t;
    };
    vector<string> bf = DBR<string>(f);
    sort(ALL(bf));
    bf.erase(ALL(bf), bf.end());
}

割と素直に書ける気がしてきませんか……？

作るもの

変数を生成できるオブジェクト を受け取り、値を返す関数 f(source) を定義します。

関数 f では source を使って「呼び出しごとに変わる変数」を生成します。

この変数を使って値を作り、返します。これが全探索の内の1回の探索に相当します。

この関数 f を受け取って、DBRSource を渡しながら繰り返し呼び出して値を記録する関数 DBR を設計します。

作る

source から出力できるのはひとまず (l,r) を受け取って l以上r以下の数 を返すものだけにします。

最終的な実装は実は結構簡単で、こんな感じになります:

class DBRSource {
public:
    int choice(int l, int r) {
        if (query_times >= (int)history.size()) {
            // new query
            Query q;
            q.type = QueryType_Choice;
            q.query_value[0] = l;
            q.query_value[1] = r;
            history.push_back(q);
        }
        auto& q = history[query_times++];
        assert(q.type == QueryType_Choice && q.query_value[0] == l && q.query_value[1] == r);
        return q.query_value[0] + q.internal_value;
    }
//private:
    enum QueryType {
        QueryType_Choice,
        QueryTypeNum
    };
    struct Query {
        QueryType type = QueryType_Choice;
        int query_value[4];
        int internal_value = 0;
    };
    std::vector<Query> history;
    int query_times = 0;
};

template<class X, typename F>
std::vector<X> DBR(F f) {
    std::vector<X> ret;
    DBRSource source;
    source.history.clear();
    source.query_times = 0;
    do {
        source.query_times = 0;
        ret.push_back(f(source));
        while (source.history.size() > 0u) {
            auto& q = source.history.back();
            q.internal_value++;
            bool is_end = false;
            switch (q.type) {
            case DBRSource::QueryType_Choice:
                if (q.internal_value > q.query_value[1] - q.query_value[0]) is_end = true;
                break;
            default:
                // expected error
                is_end = true;
            }
            if (!is_end) break;
            source.history.pop_back();
        }
    } while (source.history.size() > 0);
    return ret;
}

もうちょっとインターフェースとか整理したほうが良いんですが、それは追々……

基本的には DBRSource に問い合わせされた履歴を持って、前の問い合わせより1多い値を返す、というだけでOKです。

同じ問い合わせかどうかという判定は「何回目に呼ばれたか」で管理すれば良いです。

が、ARC168 Cの例のように、スワップ回数自体が可変の場合や、nC2選択のために値の範囲自体が可変の場合があるため、実際には返した値に応じて分岐する木構造を考える必要があります。

(木構造という観点では、これはDecision TreeをDefine-By-Runで構築している、とも考えられそうです)

ただ実際には木構造のあるパスにしか興味が無いので、結局「値を払い出し終わったクエリをpop_backする」という形でクエリの履歴を管理すれば十分です。

実践

ARC168 Cで検算に使いました (以下提出の95行目のnaive関数):

https://atcoder.jp/contests/arc168/submissions/47835839

展望

今回は一様乱数のようなインターフェースしか無いですが、以下のような拡張を作ると便利かもしれません:

  • permutation
  • combination

また、DBRSource を受け取って値を返す(DBRSource以外に副作用を産まない)関数を作ることで、その場で拡張することも容易そうです(Reactのカスタムフックみたいな考え方)。

まとめ

全探索してる時点で早解き的には負けなので、使わないに越したことは無いけど、たまにあれば便利そうですね

元ツイート

(現ポスト)

平日の深夜2時に何してるんでしょうね

TTPC 2023 をupsolveするだけのログ

atcoder.jp

TTPC2023、かなり面白い問題が多い(と思った)ので、コンテスト終了後も結構粘って解くなどしていました(もはやコンテスト5時間とか何も考えずに……　ほぼProject Euler)

一部問題は解説・解法を見てしまいましたが、基本的に解説を見ずに解いています

ネタバレ・解法バレを多分に含むので解きたい人は読まないでください

submissions: atcoder.jp

A - Numerous Elimination

  • コンテスト中に解いた
  • 頭が回ってなくてOEISで解いたそうです

B - Almost Large

  • コンテスト中に解けず、解説を読んだ
  • 完全に頭が回ってなくて、解説を読んで悲しい気持ちになりました
    • 僕はこの手の経路探索を改札ダイクストラと呼んでいるのですが、完全にそれでした
    • 由来はARC061です (当時すごく感動した記憶があります)
    • atcoder.jp

C - Yet Another Simple Math Problem

  • コンテスト中に解いた
  • 頭が回ってなくてパターンを列挙したそうです

D - Spacecraft

  • コンテスト中に解けず、後で解いた
  • ハイパー3D幾何
  • 実は無限遠だけ考えれば良いことが分かれば、あとは球面幾何でホイ
    • 3D幾何の回転とかが行けるなら割と球面幾何方針もきれいに書けると思う
    • 詳細はhackmdに置いた
    • hackmd.io

E - R-Connected Components

  • コンテスト中に解いた
  • x2 + y2 = R なるベクトル (x,y) を集めて、その線形和が張る格子の格子定数が答え
  • 格子の縮約にLLL……という気持ちをグッとこらえ、格子のサイズを減らすという操作は
    • 2ベクトルのなす平行四辺形の面積は外積
    • 外積には線形性が成り立つ
    • 使えるのは整数倍と足し算引き算
  • となり、
  • なので、任意の2ベクトルの外積を集めて、そのGCDを取れば答え
  • あまり厳密な証明はしていないし、そもそも計算量も(任意の2ベクトルの外積が必要なので)結構悪い
  • 公式解説が賢いが、まだ理解していない……

F - Na (log N)b

  • コンテスト中に解いた
  • 構文解析だ～～～
  • 構文解析をバグらせたり、log多重周りでバグらせたりした
    • 数式の意味単位のBNF表記が与えられているとBNFどおりに書けば良いので楽 (なんか再帰下降パーサが流布され始めた頃は雰囲気で定義されていた……気がする……)
      • 掛け算と足し算をどっち先に処理しなきゃいけない……みたいなのが良いBNFだとそのまま書くだけで行ける (今回のは良いBNF)
      • <expr> ::= <factor> | <expr> “+” <factor> | <expr> “*” <factor> みたいな定義をされると考えないといけなくなる
    • 例: https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=0109&lang=jp
  • 多WA枠として結構面白い問題でした

G - Cola

  • コンテスト中に解けず、後で解いた
    • これは解説を読んでないが、TLで五角数定理が流れてきてしまったのを見ています
  • 前から一個ずつ特定して、当てたら次のを1回だけ0コストでクエリできる、と言い換えると問題がシンプルになる
    • Pr[rand(1) + rand(2) + … + rand(N) ≤ M-1] みたいな形になった
  • あとは(五角数定理カンニングもあり)多項式表現に落として数式をこねると答えが出てくる
    • rand(k) を (1/k + 1/k x + … + 1/k xk-1) = 1/k (1-xk)/(1-x) と表現して計算
    • 分子側で五角数定理の形が出るので変換、分母側は負の二項係数で分子に持っていって、必要な係数の部分だけ畳み込んで終わり
    • 五角数定理、今更ながら初めて触ったかもしれない
    • あとFPSをこねるのもそこまで慣れてない

H - 404 Chotto Found

  • コンテスト中に解けず、後で解いた
    • 解説を読んでおらず、TLでSAが流れてきてしまったのを見ていますが、そもそもどうせSAだろうな……と思ってコンテスト中に問題をそっ閉じした記憶があります
  • SAとLCPを丁寧にやってセグ木に乗せると解ける
  • ACLにSAがあるから書かなくて良いの、いい時代だな……と思います
    • 自前ライブラリをちゃんと整備しろという話ではある

I - T Tile Placement Counting

  • コンテスト中に解けず、後で解いた
  • これね、すごい
  • 結論から言えば公式解説にもあるように、Tタイル敷き詰めの論文を見つけて、それを実装すればOK
  • Wがめちゃ長いので、遷移行列を求めて行列累乗する方針になるが……
    • 遷移行列の計算に論文知識を入れる必要がある
    • 公式解説はchain graphの方でやってるみたいだけど、僕はもっと数式が直接書かれている "Tilings of rectangles with T-tetrominoes" の9章のTutte多項式を使いました
    • で、これを数え上げようと思うと連結成分を管理する必要があり……
    • 完全にフロンティア法(フカシギの数え方)じゃん！！！！
    • となりました
      • (興奮するソース)
  • あとは行列累乗をすれば終わり……ではなく、logWが異常に大きく間に合わないので、Berlekamp-Masseyするか定数倍鬼速行列乗算を持ってくるかしましょう、で終わり
    • ……ちなみにBerlekamp-Masseyすると何やらカタラン数より小さい次数が出てくるんですが、そもそも上下対称な構造であることを思い出すとステートを上下ひっくり返して同一視して良くて、それを使って状態圧縮した行列を考えると実は小さくできる (BMの次数と一致する)
    • この圧縮版行列を使えば普通の行列累乗でも間に合う
  • (まだ完全に理解してないけど、行列累乗の高速化といえば Black box linear algebra だと思っていてそれを持ち出したが、ただのBMで良かった)

J - Set Construction

  • コンテスト中に解けず、後で解いた
  • パズルか？と思い色々実験してそのまま通した
  • 0、2N-1 が必ず入ることに気をつけると、ビット排他なパターンを組み合わせて M-1 か M-2 を作る問題になる
    • ビット排他なものを組み合わせる場合、それまで使ったビット位置をすべて1で塗りつぶす、をすることでAND/ORに影響が出なくなるので足し算になる
  • 使ったパターンは以下
    • 2冪 (2k - 1 通り)
    • n=3,4,5 の全探索
      • n=5 は 22n-2 = 230 がさすがに間に合わないので最低限の表現力になる 2min(2n-2,17) とかでやった
    • 2冪+α
      • 0~k-1 ビットで2冪パターンを作ったところに、k-b~k+c ビットを1にした値を添加して推移的に得られる集合
      • これは 2k + 2k-b - 1 みたいな値が作れる
        • これを繰り返したら純粋な2冪を総和できるかも……？と思ったけど試してないしあまりうまく行く気もしない
      • cは余ったビットを1に使い切るためだけの調整用
    • 乱択埋め込み
      • 上記で n=7,m=28 と n=8,m=30 以外はカバーしたので、あとはこのパターンだけ乱択で探して埋め込みました
  • ソース atcoder.jp

K - Dense Planting

  • コンテスト中に解けず、後で解いた
    • これは解説を読んでいないが、TLに流れてきた天才グラフ形状を見てしまった
  • 天才グラフ形状を頭に入れた状態で探索
    • ただこの探索も結構簡単にはできない
    • ……と見せかけて、実は普通に上限を丁寧にやる全探索がかなり小さな探索空間になるので、余裕で間に合う
    • ソース atcoder.jp

L - Next TTPC 3

  • コンテスト中に解いた
  • BもCもわからん……をしていた時にふと Next TTPC 3 が問題リストにいるのを見て、あわよくばFA取れそうと思い着手
    • 残念ながら2番手……
  • 問題自体は 周期で剰余を取る・二分探索で数え上げに落とす・境界を丁寧にやる と比較的素直で教育的な問題という印象、良い
    • あのNext TTPCシリーズがこんなつよつよ問題になってしまい、感動しています

M - Sum is Integer

  • コンテスト中に解いた
  • ユーザー解説の方の解き方をしていた
  • この手のハッシュ系乱択問題、すき

N - Bracket Sequestion

  • コンテスト中に解けず、後で解いた
  • まず部分点がきつい
    • はてなを括弧に一対一対応させる数え上げ典型をすると、O(N2)の経路DPっぽいのが手に入る
    • 素直に書くと2べきのかかったbinomial convolutionで、よくわからない
      • NxNグリッドで、左下領域を通らず、ある縦遷移を固定してそこで初めて閉じ括弧確定のはてなを使う、とすると、その縦遷移までのdyck pathと、縦遷移後のdyck path (にはてな使う使わないの2べきが掛かった形)として計算できる
      • 最終的に欲しいのは(combの中身が±1ぐらいする) 2i+n-j comb(i+j,i) comb(2n-i-j-1,n-j-1) みたいな感じ
    • 2べきの項が一致する箇所をまとめた総和の値を眺めると…………
      • k=i-j として k ごとに見てみる
        • この時カタラン数みたいに斜め線ができて反対に折り返すと comb(2n-1,…) みたいな式になりそうだなぁ……という気持ちになっておく
      • n=5, i-j=5: comb(5+0,5)comb(9-5-0,4-0) = 1
      • n=5, i-j=4: comb(4+0,4)comb(9-4-0,4-0) + comb(5+1,5)comb(9-5-1,4-1) = 11
      • n=5, i-j=3: comb(3+0,3)comb(9-3-0,3-0) + comb(4+1,4)comb(9-4-1,4-1) + comb(5+2,5)comb(9-5-2,4-2) = 56
      • n=5, i-j=2: … = 176
    • ……パスカルの三角形のある三角形の部分和っぽくないか……？(エスパー)
      • comb(2n-1,…) = comb(9,…) に注目していたので、そこから上を眺めると
        • 11 = (1) + (1+9)
        • 56 = (1) + (1+8) + (1+9+36)
        • 176 = (1) + (1+7) + (1+8+28) + (1+9+36+84)
    • → hockey stick identity が成り立つのであるラインのbinomialの総和！
      • つまりcomb(2n,…) = comb(10,…) を見れば
        • 11 = 1 + 10
        • 56 = 1 + 10 + 45
        • 176 = 1 + 10 + 45 + 120
      • (ここのエスパー、結局なんで成り立つのかまだ理解できてないが……)
    • binomial を累積すればいいので O(N\log p) とかそれぐらいになった
  • で、満点を探す
    • 上記は本当は dyck path で (binomial - binomial) * (binomial - binomial) みたいな形をしていたので、諸々展開して各二項係数を整理すると \displaystyle\mathrm{Catalan}(n) + \sum_{i=0}^{n-1} 2^{\max(2n-2-i,n)} \binom{2n}{i} になる
    • とりあえずカタラン数を求めたい
    • そもそも N! \bmod p で O(N) で計算できるのかな？ → できない……(コードテストで確認)
      • 階乗だけなら埋め込みで高速化できるらしい……というのは後で知った
    • 高速に計算する手段があるはず…… → P-recursive sequenceに出会う
    • カタラン数が高速に計算できるのは分かったけど、2冪×二項係数の総和も P-recursive ってやつなのか……？ → そうっぽい
    • で、
    • r-th.hatenablog.com
    • こうです。
  • この問題が一番時間かかった
    • 5時間コンテストとは……？
    • ちなみに「これP-recursiveだな、あとはどうやって実装するか……」を調べる際に “P-recursive AtCoder” で検索したらTTPC2023の解説がトップヒットしてしまい、半分ネタバレをくらいました(結局P-recursiveなんだな……という情報しか得られてないのでセーフ)

O - 2D Parentheses

  • コンテスト中に解けず、後で解いた
  • 包含関係を考えるとタイブレークできて、同じ座標の点を無くすことができる
  • 端の方から閉じ括弧を見ていくことを考えると、一番近い開き括弧を選ぶとしてよい
    • 近いのを選ばず遠いのを選ぶと、後々その近い開き括弧を他の閉じ括弧が選ぼうとした時に邪魔になる可能性がある
  • あとはこれが判定できればOK
    • 二次元包含クエリで、長方形の四隅にxorが0になるような4つの乱数を振っておけば、長方形領域のxor総和が0か？で半端な包含を判定できる
    • ……が、どうやらこれはオーバーキルらしい

P - Bridge Elimination

  • コンテスト中に解けず、後で解いた
  • 対称性がすごい
    • i 個の値を乗算した項の総和は (1+A_ix) の乗算結果(ほぼDP)でOK
  • あとは「 i 個の二重辺連結成分になるようなグラフの数(厳密には対称式 [x^i]\prod_{i=1}^N(1+A_ix) の寄与)」を数え上げればOK
    • これは寄与典型
  • 二重辺連結成分はOEISに転がってる
    • A095983
    • なんかMathematicaでFPSいじれば行けるって書いてあるのでやると、なんか求まる
  • あとは k 個の二重辺連結成分に分割されたグラフを一成分ずつ決めていくDPで求めていく
    • ただ最終的な k 成分の木の繋ぎ方がCayleyの公式でうまくいかない(各成分の次数に応じてどの頂点と繋ぐかが累乗される)のでなんとか数えたい
    • “count tree fixed degree” とか調べると
    • math.stackexchange.com
    • \displaystyle\frac{(k-2)!}{\prod_i(d_i-1)!} らしいです
  • 各成分の木としての次数の総和を管理しながらやれば行けそうですね。
    • dp[k][v][d] := k成分をv頂点で木次数総和dで作る数え上げ
    • 状態数O(N3) 遷移O(N2) でTLE……
    • 2次元FFTすれば O(N3 logN) も見えそうだけど N=400 はきつい……
  • 各成分の構成頂点数が m_1,\dots,m_k で固定の場合に、次数制約とか気にせず数え上げられれば良さそう
    • もう一回考える
    • e_i := d_i-1 として、dp[i][e] := i 成分まで見て e_i の総和が e
      • 最終的に欲しいのは dp[k][k-2]
    • 使う次数を f として dp[i+1][e+(f-1)] ← dp[i][e] * m_i^f\cdot\frac{1}{(f-1)!}
    • 多項式にすれば \displaystyle f_{i+1}(x)=f_i(x)\cdot\sum_{f=1}m_i^f\frac{x^{f-1}}{(f-1)!}=f_i(x)\cdot m_i\sum_{e=0}m_i^e\frac{x^e}{e!}
      • \displaystyle g_i(x):=\sum_{e=0}m_i^e\frac{x^e}{e!} としましょう
      • 最終的に欲しいのは \left(\prod_{i=1}^k m_i\right)[x^{k-2}]\prod_{i=1}^k g_i(x) ですね
  • これ指数型母関数ってやつですよね？入門しちゃいましょう
    • 37zigen.com
    • 集合と関連付けられるとすると m_i^e なので e 個要素を選んで m_i 色に塗り分ける……と考えられそう
    • つまりこれは「 k-2 個のボールがあって、 e_i 個ずつ k 個のグループに分けて、グループごとに決まった m_i 色の内の色で各ボールを塗り分ける」という通り数……？
    • もはやグループごとに色種類を分ける必要はなくてどのグループであっても \sum_{i=1}^k m_i 色から1色選ぶ……として重複なく合算できるので、結局これは N^{k-2}
    • ほしかった値が \left(\prod_{i=1}^k m_i\right)N^{k-2} になったので、次数制約をDPから省いてトータル O(N^3)
      • ちなみに O(\log p) がつくとTLE(1敗)

おまけ: 解いた順

(コンテスト開始)

10/14: A → L → C → F → E → M

(コンテスト終了)

10/17: G

10/19: H

10/21: D

10/23: I

11/03: O → B → K

11/12: N

11/18(11/19): P → J

時間かかりすぎている

P-recursive sequenceの第N項の平方分割計算アルゴリズムの理解ログ

これ↓ を理解するまでのログです。

https://nyaannyaan.github.io/library/fps/find-p-recursive.hpp.html

問題定義

扱う数は素数 p を法とした剰余環 Z/pZ とし、以降は単に整数・数と呼びます。

求めたい数列を \{a_i\}_i とし、最終的には a_N を求めます。

数列 \{a_i\}_i は、整数を入力とする高々 d 次の多項式 n+1 個の列 \{p_j(x)\}_{j=0,\dots,n} を用いて、以下の関係が成り立ちます。

 \displaystyle \sum_{j=0}^{n} a_{i+j}p_j(i)=a_{i+n}p_{n}(i)+\dots+a_{i+1}p_1(i)+a_ip_0(i)=0

行列積としての表現

線形漸化式における行列累乗のような考えで、この関係式を行列積として表現し直します。

まず関係式を少し整理します:

 \displaystyle \begin{align*}a_{i+1}p_{n}(i+1-n)&=-\sum_{j=0}^{n-1}a_{i+1-n+j}p_j(i+1-n) \\ a_{i+1}&=\frac{\sum_{j=0}^{n-1}a_{i+1-n+j}p_j(i+1-n)}{-p_n(i+1-n)}\end{align*}

よって、 v_{i} = (a_{i},a_{i-1}, \dots, a_{i-n+1})^T から v_{i+1} = (a_{i+1},a_{i},\dots,a_{i-n})^T に変換するには、 n\times n行列

 \displaystyle \begin{align*}M(x)&=\begin{pmatrix}\{p_{n-1-j}(x+1-n)\}_{j=0,\dots,n-2}&p_0(x+1-n)\\-p_n(x+1-n)\ I_{n-1}&\mathbf0\end{pmatrix}\\&=\begin{pmatrix}p_{n-1}(x+1-n)&p_{n-2}(x+1-n)&\cdots&p_{1}(x+1-n)&p_0(x+1-n)\\-p_n(x+1-n)&0&\cdots&0&0\\0&-p_n(x+1-n)&\cdots&0&0\\\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&\cdots&-p_n(x+1-n)&0\end{pmatrix}\end{align*}

を用いて、 v_{i+1}=\displaystyle\frac{M(i)}{-p_n(i+1-n)}v_{i} と書けます。

求めたいものは a_N、つまり v_N の最初の要素なので、 v_N を求めようと思うと、( p'(i):=-p_n(i+1-n) と定義して)

 \displaystyle v_N=\left(\prod_{i=N-1}^0\frac{M(i)}{-p_n(i+1-n)}\right)v_0=\frac{\left(\prod_{i=N-1}^0M(i)\right)v_0}{\prod_{i=N-1}^0p'(i)}

が計算できれば良いです。(行列積が非可換であるため、総乗は逆順である点に留意します。逆順総乗は記法あってるか分かりませんが雰囲気で……)

なお、ここでは a_i=0 (i\lt 0) としても関係式が成り立つとしていますが、そうでなければ b_i:=a_{i-n+1} と適切にシフトした数列について考え、 v_0 = (a_{n-1},a_{n-2},\dots,a_0) として v_{N-n+1} を求めることを考えれば良いです。

さて、分母については 1\times1 行列 \begin{pmatrix}p'(i)\end{pmatrix} の総乗だと思うと、これは M(i) で n=1 とした場合と構造としては同じになります。よって \prod_{i=N-1}^0 \mathit{Matrix}(i) を計算するアルゴリズムがあれば、これを分母・分子それぞれに1回ずつ適用することで v_N が求まります。

(天下り的な)平方分割の構成

n\times n 行列同士の積の計算量を O(MM(n)) とします。

P([a,b)) := \prod_{i=b-1}^a M(i) = M(b-1)M(b-2)\cdots M(a+1)M(a) と定義します。これは b-a 個の行列の積です。見通しの良さのため半開区間を用います。

ブロックの分割数 B を B:=\lfloor\sqrt{dN}-(d-1)\rfloor とします。この数字は調整可能です。

ブロックの幅 W を W:=\lfloor N/B\rfloor とすると BW\le N なので、 B 個の「 W 個の行列積のブロック」に分けます。

 \displaystyle \begin{align*}&\prod_{i=N-1}^0 M(i)\\=&P([0,N))\\=&P([BW,N))P([(B-1)W,BW))\cdots P([W,2W))P([0,W))\\=&P([BW,N))\prod_{i=B-1}^0P([iW,(i+1)W))\end{align*}

前半にかかっている P([BW,N)) は W>N/B-1 より N-BW\lt B なので愚直に計算すれば O(B\cdot MM(n)) です。よって後半の B 個の総乗 \{ P([iW,(i+1)W)) \}_{i=0,\dots,B-1} を高速に求められれば良いです。

ここまでの図はこんな感じです(説明用の例: N=22, B=5):

この各 P([iW,(i+1)W)) を、以下の図のようなイメージのダブリングによってまとめて計算します(説明用の例: N=22, B=5, d=1):

このダブリングをうまく動かすには、これらの値(行列積)がある多項式によって表現されていることを意識する必要があります。

行列積の多項式表現

乗算される行列の多項式表現を再掲します:

 \displaystyle \begin{align*}M(x)&=\begin{pmatrix}\{p_{n-1-j}(x+1-n)\}_{j=0,\dots,n-2}&p_0(x+1-n)\\-p_n(x+1-n)\ I_{n-1}&\mathbf0\end{pmatrix}\\&=\begin{pmatrix}p_{n-1}(x+1-n)&p_{n-2}(x+1-n)&\cdots&p_{1}(x+1-n)&p_0(x+1-n)\\-p_n(x+1-n)&0&\cdots&0&0\\0&-p_n(x+1-n)&\cdots&0&0\\\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&\cdots&-p_n(x+1-n)&0\end{pmatrix}\end{align*}

p_j(x) は高々 d 次の多項式なので、 p_j(x+1-n) も高々 d 次の多項式です。

例えば M(0) は M(x) に x=0 を代入した時の行列です。

では M(4) について考えると、これも x=4 とすれば良いですが、ここであえて x_4:=x-4 とシフト(変数変換)した行列 M_4(x_4)=M(x+4) を考えます:

 \displaystyle \begin{align*}M_4(x_4)&=\begin{pmatrix}\{p_{n-1-j}(x_4+5-n)\}_{j=0,\dots,n-2}&p_0(x_4+5-n)\\-p_n(x_4+5-n)\ I_{n-1}&\mathbf0\end{pmatrix}\\&=\begin{pmatrix}p_{n-1}(x_4+5-n)&p_{n-2}(x_4+5-n)&\cdots&p_{1}(x_4+5-n)&p_0(x_4+5-n)\\-p_n(x_4+5-n)&0&\cdots&0&0\\0&-p_n(x_4+5-n)&\cdots&0&0\\\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&\cdots&-p_n(x_4+5-n)&0\end{pmatrix}\end{align*}

すると、 M_4(x_4) の各要素も高々 d 次の多項式で、また M(4)=M_4(0) です。

このようにシフトされた行列表現 M_i(x):=M(x+i),M_i(0)=M(i) を考えると、「ある多項式行列の積の x=0 の標本点」と表現し直すことが可能です:

 \displaystyle P([a,b))=\prod_{i=b-1}^a M(i)=\left(\prod_{i=b-1}^a M_i(x)\right)(0)

多項式としての行列積にも表記を与えておきます:

 \displaystyle \begin{align*}P([a,b);x)&:=\prod_{i=b-1}^a M_i(x)\\P([a,b);0)&=P([a,b))\end{align*}

シフトされた行列 M_i(x) の各要素は高々 d 次の多項式だったので、その積である行列 P([a,b);x) はの各要素は高々 d(b-a) 次の多項式になります。

(ここまで次数を執拗に確認していますが、(上界の)次数のみが重要で、実際にどのようなシフトをしてどう展開され係数がいくつになるか、ということには興味が無い点に留意してください。)

また M(x) から M_i(x)=M(x+i) とシフト出来るように、 M_s(x) から M_t(x)=M_s(x+(t-s)) とシフトすることも可能です。これを利用すると、 P([a,b);x) についてもシフトが可能です:

 \displaystyle P([a,b);x+s)=\prod_{i=b-1}^a M_i(x+s)=\prod_{i=b+s-1}^{a+s} M_i(x) = P([a+s,b+s);x)

これを考えると、例えば P([4,8))=P([4,8);0)=P([0,4);4) となります。

ダブリングの構成

さて、行列積の多項式表現と、そのシフトを考えながら図の例に戻してみると、要はこういった構造になっています:

つまり、各ステップで P([0,w); x) という多項式を考え、その各点の標本点(行列)を管理する、という形になっています。 P([0,w); x) の各要素は高々 dw 次の多項式なので、 dw+1 個の標本点を管理すれば、多項式としての情報量を損ないません(ref: ラグランジュ補間)。

よって作られるダブリングアルゴリズムは以下の通りで、最終的に P([0,W); x) の x=0,W,\dots,(B-1)W での標本点が分かれば良いです。

  • 幅 w の行列積 P([0,w); x) の x=0,W,\dots,dwW での標本点 P([0,w);iW) が求まっていると仮定する。
    • 標本点の数は dw+1 個。
    • w=1 のケースは愚直に d+1 個の行列を評価をすることで計算可能。
  • 幅を1増やすには:
    • i\le dw の場合: 各標本点 P([0,w);iW) に単純に行列を掛ける。
      • P([0,w+1);iW)=M(iW+w)P([0,w);iW)
      • dw+1 回の行列評価・行列乗算
    • dw+1\le i\le dw+d の場合: 単純に行列を掛ける。
      • P([0,w+1);iW)=P([iW,iW+w+1))
      • d(w+1) 回の行列評価、 dw 回の行列乗算
    • (動作例: 次数 d=2。w=1 (標本点3個) から w=2 (標本点5個) を w+=1 で作る)

  • 幅を2倍に増やすには:
    • i\le dw の場合: 各標本点を w だけシフトした値を考える。
      • P([0,2w);iW)=P([0,w);iW+w)P([0,w);iW)
      • \{P([0,w);iW+w)\}_{i=0,\dots,dw} は \{P([0,w);iW)\}_{i=0,\dots,dw} に関して、多項式 P([0,w);x) から多項式 P([0,w);x+w) へとシフトした値を考えれば良い。
      • これは多項式の標本点シフトとして知られており、 O(dw\log dw) 回の演算で可能。
      • 標本点シフトに O(n^2 dw\log dw) (行列の各要素の多項式に注目しているので、行列としての演算は発生しない)
      • シフト後の行列を乗算するのに dw+1 回の行列乗算
    • dw+1\le i\le2dw の場合: +(dw+1)W と +(dw+1)W+w だけシフトした値を考える。
      • P([0,2w);iW)=P([0,w);iW+w) P([0,w);iW)
      • j=i-(dw+1) とすると 0\le j\le dw-1
        • P([0,w); iW)=P([0,w);jW+(dw+1)W)
        • P([0,w);iW+w)=P([0,w);jW+(dW+1)W+w)
      • i\le dw と同様に多項式の標本点シフトを2回行うことで計算可能。
      • 標本点シフトに O(n^2 dw\log dw)
      • シフト後の行列を乗算するのに dw 回の行列乗算
    • (動作例: 次数 d=2。w=1 (標本点3個) から w=2 (標本点5個) を w×=2 で作る)

あとは実装するだけです。

多項式の標本点シフト

標本点シフトは以下の問題を解くアルゴリズムです:

  • 高々 d 次の多項式 f(x) の標本点 f(0),f(1),\dots,f(d) が与えられる
  • 整数 c が与えられるので、 f(c),f(c+1),\dots,f(c+d) を求めよ

行間埋めは以下が非常に参考になります:

https://mathlog.info/articles/3351

最終的に3回のFFT畳み込みによって O(d\log d) で高速に計算が可能なことが分かります。

本問題は library-checker-problems (yosupo judge) に収録されているので、具体的な実装はその想定解ファイルを読むのが一番簡便かと思われます。

https://github.com/yosupo06/library-checker-problems/blob/master/math/shift_of_sampling_points_of_polynomial/sol/correct.cpp#L196

今回の問題では f(0),f(W),\dots,f(iW),\dots,f(dwW) という形の標本点になっているので、そのまま使うことは出来ません。

g(x):=f(xW^{-1}) とすると、 f(iW)=g(i) となります。

g(i) として管理をした場合、標本点シフトは f(iW+c)=g(i+cW^{-1}) と考えられ、結局 g(i) に関して cW^{-1} のシフトを実行すれば良いことになります。

残り

なお、「最終的に必要なブロックの数は B」「幅 W の行列積の多項式の次数は dW で、必要な標本点は dW+1」ということを考えると、ダブリングの結果が余ったり足りなかったりします。

  • 余るケース dW+1> B
    • 何も考えずに、最終的に捨てれば良いです。
  • 足りないケース dW+1\lt B
    • k 回の標本点シフトによって標本点の数を k(dW+1) にすることが出来ます。

W:=\lfloor N/B\rfloor だったことを思い出すと、この条件は以下と等価です:

  • 余るケース
    • d\lfloor N/B\rfloor +1\ge d(N/B-1)+1>B
    • B^2+(d-1)B-dN\lt 0
    • B\lt \frac{-(d-1)+\sqrt{(d-1)^2+4dN}}{2} (が成り立つならば余る)
  • 足りないケース
    • d\lfloor N/B\rfloor +1\le d(N/B)+1 \lt B
    • B^2-B-dN \gt 0
    • B>\frac{1+\sqrt{1+4dN}}{2} (が成り立つならば足りない)

特に B\le\sqrt{dN} - (d-1) の時、

 \displaystyle B \le \sqrt{dN}-(d-1) = \frac{-(d-1)+2\sqrt{dN}-(d-1)}{2} = \frac{-(d-1)+\sqrt{(d-1)^2-4(d-1)\sqrt{dN}+4dN}}{2} \le \frac{-(d-1)+\sqrt{(d-1)^2+4dN}}{2}

が成り立つので、 B=\lfloor\sqrt{dN}-(d-1)\rfloor と設定することで、足りないケースは回避することが出来ます。

(もちろん実装時には不等式 d\lfloor N/B\rfloor + 1 \ge B を満たす最大の B を全探索/二分探索で求めた方が効率的です。ここではオーダー把握のため簡易な形式を採用します。)

計算量

最後に、計算量の見積もりを行います。

ある行列 M(i) を評価するのに必要な計算量は、多項式が一番上の行と対角線一個下にしか無いことに注意すると、 O(n^2+nd) です。

行列の乗算に必要な計算量を O(MM(n)) と置いていました。ここではシンプルに乗算するとして O(n^3) ということにします。

\{ P([0,w);iW) \}_{i=0,\dots,dw} を求める計算量を T(w) と置きます。

T(1)=O(d(n^2+nd))=O(dn^2+d^2n) です。

w が奇数の場合は w-1 に再帰します。

T(w)=T(w-1)+O(dw(n^3+nd))=T(w-1)+O(n^3dw+nd^2w) です。

w が偶数の場合は w/2 に再帰します。

T(w)=T(w/2)+O(n^2dw\log dw+ dwn^3)=T(w/2)+O(n^3dw+n^2dw\log dw)) です。

w が高々2ステップで w/2 になることを考えると T(w)=O(n^3dw+n^2dw\log dw+nd^2w) です。

よって T(W)=O(dW(n^3+nd+n^2\log dW)) です。

B=\lfloor\sqrt{dN}-(d-1)\rfloor とすればこのあとに追加で評価点シフトする必要も無いです。

あとは B 個の P([iW,iW+W)) と高々 W 個の M(BW+i) を乗算することで、 \prod_{i=N-1}^0 M(i) が求まります。

この行列積を求める総計算量は

O(dW(n^3+nd+n^2\log dW)+(B+W)n^3)=O(n^2dW\log dW+n^3(B+dW)+nd^2W)

となります。 n,d \ll W,N とすれば O(n^2dW \log dW) で、 W=\lfloor N/B\rfloor=O(\sqrt{N/d}) を思い出すと O(n^2\sqrt{dN}\log dN) です。

(なお、NTT-friendly でない素数の場合、畳み込みに garner's algorithm 等を用いた任意mod畳み込みを使用する必要があります。garner's algorithm を用いる場合 \bmod\ p での逆数計算が必要になるため、 O(n^2\sqrt{dN}(\log p+\log dN)) になります。)

その他

  • 本アルゴリズムでは「たかだか dw 次の多項式は dw+1 個の標本点を保存すれば良い」という性質を用いて、標本点の数自体も幅 w と一緒にダブリングしていますが、最初から最終的な多項式の次数 dW に対応した dW+1 個の標本点を持ち続けるアルゴリズムも設計可能です。
    • この場合、標本点の数が最初から最後まで O(dW) なので、段数がかかって管理する標本点の総数が O(dW\log\cdot) になります。
    • 標本点の数自体もダブリングする方式では級数を考えることで、標本点の総数が O(dW) になるので、この分の \log が落ちることになります。
  • また、多項式 P([0,W);x) を陽に求め、多点評価によって P([0,W);iW) を求めるアルゴリズムも設計可能です。

後記

  • お久しぶりです
    • また O(\sqrt N) の記事書いてる……
    • 関係ないですがTTPC良かったですね。僕はTTPC upsolveに睡眠時間を削られています。来年はTTPC(という名前のコンテスト)が開催されるのでしょうか……
  • P-recursive がちょうど流行り始めてそうな2020年～2021年頃にぬるっと引退したので存在を知らず、すげ～となった
    • でも N! \bmod p 自体はもっと昔から存在はしていたらしいので、単純に精進が足りなかっただけでもありそう
    • 雑に Berlekamp-Massey に投げる感覚で P-recursive 判定に投げる世界、ヤバい
  • 階乗を o(N) で解けるのだいぶ驚きだったんだけど、 o(\sqrt N) はやばそうだな……という感覚がある
    • 素因数分解が解けてしまうので……
    • と思ったらそれが Strassen’s Factoring algorithm だった
  • FPSの文脈で出がちなイメージだけど、おそらくコレ自体は”FPS”とは関係ない
    • x に値をバリバリ代入するし、やってることはただの多項式演算
    • 係数に意味があるわけでもない(多分、おそらく; 研究的価値はありそう)
    • とはいえ結局数列 \{a_i\}_i を錬成するのにFPS使いがちなのでジャンルとしてはFPSの延長っぽい

参考文献

yukicoder No.2116 Making Forest Hard を平方分割+rerootingで解く

お久しぶりです。最近はTTPC2022のtesterやったりICPC観戦したりして急に競プロ熱が再燃しています。

ところで今回はyukicoder No.2116 “Making Forest Hard” を平方分割+rerootingで解こうと思います。

https://yukicoder.me/problems/no/2116

問題概要

  • N 頂点の木 T がある
  • 各頂点には整数 A_i が書かれている
  • 木 T の任意の連結成分 U について、 U 内の頂点数を S(U) 、 U 内の頂点の最大の A_i の値を M(U) と定義し、連結成分 U のスコアを g(U) = S(U)M(U) と定義する
  • 木 T の辺集合 E の任意の部分集合 E' を考え、 T から E' の辺をすべて削除した木(連結成分)の集合を T'(E') とする
  • T'(E') のスコアを、任意の T'(E') 内の連結成分 U のスコアの総和 f(E') = \sum_{U\in T'(E')} g(U) と定義する
  • 任意の部分集合 E'\subseteq E に対する T'(E') のスコアの総和 \sum_{E'\subseteq E}f(E') を MOD で剰余を取って出力せよ

基本考察: 寄与

最大値(×サイズ)の任意の操作における総和を考える際、典型的な方法の1つが寄与の計算です。 今回はある頂点 i に書かれている値 A_i がどの程度スコアに寄与するかを考えます。

なお、 M(U) の計算のために使う最大値計算を、頂点IDによってタイブレークしてどの頂点の値かを一意に定めることにします。つまり M(U) = A_{\arg\max_{i \in U} (A_iN + i)} のようにします。

以降はこのタイブレークされた値を考え、 A_i は全て相違なるものと仮定し、単純に A_i の大小で記述します。

木DP O(N2)

さて、頂点 i に着目するため、一度 i を根として木を見、これを R_i とします(元々の木 T は根の決まっていない無向木です)。

例として適当な木 R_i を挙げます。 i=9 に着目しています。

なお簡単のため、この図では頂点ID = A_i としています。実際の問題ではそのような制約は無いので留意してください。

またグレーに塗られた頂点は A_v>A_i となる頂点です。これらの頂点は連結成分に含まれないようにする必要があります。

この根付き木 R_i において、各頂点 v について

  • v を根とした R_i の部分木(連結成分) U で、 M(U) \le A_i であるようなものの通り数 X_v
  • そのような U についての S(U) の総和 Y_v
  • v を根とした R_i の部分木に含まれる辺の数 Z_v

を求める木DPを考えます。

まず辺の数 Z_v はそのまま辺の数を足し上げるだけで良いです。

次に X_v, Y_v について、頂点 v に対して子 w の部分木の結果を順番にマージした際の値の変化を考えます。 子 w の手前までマージしていった現在の頂点 v の値を (X'_v,Y'_v) とし、子 w をマージした後の値を (X''_v,Y''_v) とします。子 w の値は (X_w,Y_w,Z_w) とします。

  • 辺 (v,w) を加え入れる場合
    • 通り数は v 側と w 側の値を乗算することになり X''_v \leftarrow X'_vX_w となる
    • サイズの総和は v 側の総和 Y'_v が X_w 通り寄与し、 w 側の総和 Y_w が X'_v 通り寄与するため、 Y''_v \leftarrow Y'_vX_w+X'_vY_w となる
  • 辺 (v,w) を取り除く場合
    • 子 w の部分木内に含まれる辺はどのようにしても良くなる
    • そのため X''_v \leftarrow X'_v\mathrm{pow}(2,Z_w) 及び Y''_v \leftarrow Y'_v\mathrm{pow}(2,Z_w) となる

以降、簡単のため W_v := \mathrm{pow}(2,Z_v) と表記します。また毎回 pow を計算するのも面倒なので実装上も基本的にはこの形で持つことにします。

この2通りを足し合わせれば良いですが、 A_w>A_i となる場合は辺 (v,w) を取り除かなければなりません。よって

  • A_w \le A_i の時
    • X''_v=X'_vW_w+X'_vX_w
    • Y''_v=Y'_vW_w+Y'_vX_w+X'_vY_w
  • A_w>A_i の時
    • X''_v=X'_vW_w
    • Y''_v=Y'_vW_w

となりますが、 A_w>A_i の時に X_w=Y_w=0 とすれば上側の計算だけで済みます。 またこれらの計算は全て X'_v,Y'_v の線形和になるので、初期化のタイミングで X_v=Y_v=0 とすれば場合分けが除去できます。 ちなみに A_w\le A_i の場合は X_v=Y_v=1 と初期化します。

木DPの結果、 Y_i に A_i の寄与する乗数が求まるため、答えに A_iY_i を足し合わせます。

以上で、各頂点 i に対して O(N) のDFSで寄与が求まるため、全ての頂点に対して行えばトータル O(N^2) の計算量で求まりました。

発展考察: 平方分割+rerooting

さて、木DPを高速化するにあたり、典型の一つである全方位木DPを考えてみます。が、上述の木DPは根の値 A_i に応じてグレーに塗られた切らなければいけないノード A_v>A_i が変わるという問題点があります。

逆に A_v>A_i という制約に対処しようと考えると A_i の昇順もしくは降順で木DPすることで、変化する頂点は根に来るため O(1) の変化でほとんど同じ木DPを解くことが出来るようになりますが、逆にrerootingのコストが1回につき最大で直径分かかります。

そこで、平方分割によって木の直径を O(\sqrt N) に圧縮し、rerootingのコストを1回あたり O(\sqrt N) とし、トータルコストを O(N\sqrt N) にすることを目指します。

平方分割

平方分割した状態でrerootingとクエリ対象の頂点編集を行いたいので、次のような性質が必要です:

  • クエリする頂点集合を D={i_1,i_2,\dots,i_B} とする
  • 頂点 i\in D は圧縮しない
  • 平方分割後の木の直径が O(\sqrt N) である
  • 木DPの値の整合性を保ったまま回転が出来る

このように平方分割するには次のようにするのが良いです:

  • i_1 を根とするような木を考える
  • 任意の i_x, i_y \in D について、 \mathrm{lca}(i_x,i_j) も圧縮しない
    • このとき、圧縮しない頂点を分割点と呼び、その集合を D' とする
    • このとき i_1 が根であることに注意すると \lvert D'\rvert \le 2B-2
  • D' によって分割された連結成分について、それらは1個か2個の分割点と繋がっている
    • 0個: 元々その連結成分が木全体である場合のみ起こり得るが B\ge 1 のため存在しない
    • 3個以上: そのような連結成分があると仮定すると、連結成分内のある1つの頂点から相違なる3つの分割点 i_x,i_y,i_z\in D' への交わらないパスが取れるが、このときその頂点は3つの分割点のどれかのLCAであり、 D' に属さない頂点から構成されるという仮定に矛盾するため、存在しない
    • (これはつまるところtop treeにおけるclusterとboundary vertexである)
  • 2個の分割点と繋がっている連結成分はたかだか \lvert D'\rvert -1 個
    • つまり \le 2B-3 個

下図はクエリする頂点集合を 9~12 としたときの木のサンプルです。

青い 9~12 の頂点が D、赤い 7 の頂点が D' \setminus D、太めの黒い曲線に囲まれた連結成分が2点の分割点に隣接する(圧縮対象の)連結成分です。

分割点では通常の木DPを計算することとし、連結成分に対する木DPを効率化します。

  • 全ての分割された連結成分について、隣接する分割点を根方向としたときの木DPの値を計算できるようにする
    • 隣接1個: rerootingによって根方向が変わることは無いため、特に何もしなくて良い(極論圧縮もしなくて良い)
    • 隣接2個: rerootingによって根方向が変わることがあり、たかだか1個の分割点を部分木に持つ。この分割点以下の木DPの値を入力とし、もう一方の分割点を根方向とした連結成分の木DPの値を出力する作用素を構成する必要がある

この作用素について考えます。

木DPの圧縮

元の木DPを再掲します:

  • W''_v=2W'_vW_w
  • X''_v=X'_vW_w+X'_vX_w
  • Y''_v=Y'_vW_w+Y'_vX_w+X'_vY_w

葉方向の分割点を u とすると、木の回転などによって (W_u,X_u,Y_u) の値は変化するため、これらの値を自由変数と考えると、

  • 分割点 u : (X_u,Y_u,W_u) という値は (W_u,X_u,Y_u) の線形結合(恒等関数)で書ける
  • 分割点 u を子孫に含まない頂点 : (W_u,X_u,Y_u) に依存しない定数値で書ける
  • 分割点 u を子孫に含む頂点 : 分割点 u を含む部分木が (W_u,X_u,Y_u) の線形結合で書けるなら、定数とのマージなので計算結果も (W_u,X_u,Y_u) の線形結合で書ける
    • u の深さに関する帰納法で証明できる

となるため、この線形結合の係数行列さえ保存すれば木の回転によって根方向が変わった場合でも木DPの値を高速に更新することが出来ます。

木DPの計算式から上三角行列となるため、係数保持に必要な値は6個で十分です(実際には5個でも大丈夫です)。

木DPのrerooting

次に木の回転(rerooting)を考えます。

rerootingの典型として「逆演算を定義し、回転する部分木を打ち消す」という手法が考えられます。一方で今回は剰余乗算を含んでいるため、0が乗算される際の扱いが困難です。また剰余除算に O(\log\mathrm{MOD}) がかかってしまうのもネックです。

今回は平方分割によって主要な木の頂点数が O(B) となっているため、そもそも任意のパス上の次数の総和も O(B) になりそうです。

ただし分割点との隣接が1点のみの連結成分の数は抑えられていません。これらの部分木については値の更新が行われないため、木の回転のたびに計算しないようにすることで、1回の rerooting が O(B) で収まります。

実装方針

  • 根 i_1 を決める
    • 最初のクエリは A_{i_1} について行うことに留意
  • D の頂点を塗っておく
  • 圧縮のDFSを行う: 各頂点 v について
    • v の子孫に D' の頂点を含むか、 D' の頂点を含む子の数を求める
      • D' の頂点を含む子の数が1の場合は圧縮対象、2以上の場合は分割点となる
    • 通常の木DPの値 (W_v,X_v,Y_v) 及び D' の頂点を含まないstaticな木DPの値 (W^s_v,X^s_v,Y^s_v) を求める
    • 圧縮対象の分割点へのリンクと係数行列を計算する
      • 圧縮対象のパスの長さが1の場合は圧縮しないことにするとちょっと楽
  • i\in D について
    • i を根にする: i の親を根にする(再帰) → i が浅くなるように回転する
      • 回転は圧縮ノードを間に挟むか否かで若干変わる
    • i についてクエリ結果を求める
    • i を変化させる (今回は A_i > A_{\mathrm{next}\ i} となるように変化させる等)

解答ソース

まとめ

ちなみにオフライン限定ですが、直径が O(\sqrt N) のためパス系クエリや、圧縮をうまく取り扱えば遅延伝播系クエリ、更に森に対して圧縮すればLink Cutにも対応出来ると思います。

まぁ多分top-treeで良いんですが、しいて言えば以下のメリットを見いだせると思います:

  • 「木を木で取り扱う」のに対し、このアルゴリズムは連結成分の圧縮しかしない
  • 定数倍は軽い(はず)
    • オーダーはめちゃめちゃ悪い、まだ最適化の余地がありそうですが N=100000 が限度だと思います
  • 実装も軽い(はず)
    • 多分……オフラインだからその取り扱いが煩雑かも

補遺: 経緯

今回のアルゴリズムはTTPC2022のtesterをしている際に思いついたものになります。TTPC2022 testerに呼んでくださった運営の方ありがとうございました。今更ですが参加してくださった皆様もありがとうございました。

まぁTTPC2022の木の問題はD: XOR Tree Pathだけだったんですけどね。

参考文献

ARC068 E Snuke Line O((N+M)√M + NlogN)解

みんなFenwick Tree使ってて泣いたので。

問題概要

arc068.contest.atcoder.jp

N個の区間[l,r]でお土産が売っている。全ての1<=d<=Mに対して、dの倍数の地点を全て訪れることで入手できるお土産の種類を答えよ。

1 \le N \le 3\times10^5, 1 \le M \le 10^5, 1 \le l_i \le r_i \le M

解法

dの値に関して平方分割をします。

d \le \sqrt M の時、全ての区間に対してdの倍数が含まれているかをチェックします。具体的には \lceil \frac{l_i} d \rceil\times d \le r ならば含まれています。よって計算量は O( N \sqrt M )です。

d > \sqrt M の時、訪れる地点はたかだか \sqrt M です。この訪れる地点を各dに対して更新するという方針を取ります。訪れる地点を指し示す構造体を便宜上ポインタと呼びます。

0番地側から1-basedでx番目のポインタは d \times xの地点を指していることになります。dを1増やすとポインタの位置はxだけ動きます。

ポインタが動くと区間に入ったり出たりするイベントが発生するので、ポインタには「区間左側ソート順でどこまでの区間に入ったか」「区間右側ソート順でどこまでの区間から出たか」という情報を持たせておきます。

区間の出入りで情報が更新される度に、区間を指すポインタの数を配列で管理し、インクリメント・デクリメントします。また、区間を指すポインタの数が0から1になったら答えを1増やし、1から0になったら答えを1減らします。

このようにポインタを動かすと、1つのポインタについては全体でたかだかMしか動かず、区間イベントもたかだかNで、そのポインタが \sqrt M個なので、合計して計算量は O( (N+M) \sqrt M )で、区間のソートがあるのでここに O( N \log N ) が加わります。

以上を合算して O( (N+M) \sqrt M + N \log N ) の計算量を達成します。

提出

arc068.contest.atcoder.jp

1276ms。TLEが2secなので厳しい。

感想

はてなでmarkdown+TeXってありえんつらい

CODE FESTIVAL 2016 本戦 参加記 #CODE_FESTIVAL

去年に引き続き2回目の参加でした。楽しかった。

朝

去年同様の若干の遅刻でした。

着いてすぐにお昼ごはんを食べました。今年はTシャツの説明をちゃんと聞けたためTシャツが恵まれました。

本戦

今年のパーカーボーダーは1600点。

結果はABCDF+E部分+H部分で3700点(全体33位、国内15位)でした。去年に比べたら成長を感じられる結果。以下解答方針。

A

書きました。cin>>sしながらチェックするだけでとても楽。

B

N<=107なのでO(NlogN)が怪しい気がしたんですけど大丈夫でした。多分「メモリアクセスは√Nと考えるべき」みたいなのが効いてるんだと思う(つまりメモリアクセスが無い分若干大きい計算量見積もりでも大丈夫という意味)。

使える数の上限を決め打ちすると貪欲に取れるので、それで二分探索しました。

C

結局「同じ言語話せる人たちをuniteしたUFで全部連結になるか」ってことなのでUFを書いて提出。誤読が多かったらしい(この日は珍しく誤読しなかった)。

D

サンプルでも分かる通り、mod MがiとM-iのものをできるだけ組み合わせるほうが良い。

まずmodMで分けて、組み合わせて、余ったやつの中で同じ数のやつを組み合わせる、ということを貪欲っぽくやる。4WA。

E

部分点のみ。dp[i] := 速度iで始めてN枚クッキーを焼けるまでの最短時間 というDPを考える。

dp[i] = max( (N+i-1)/i, max{ j+A+dp[j*i] } ) という更新式を考えるとO(NlogN)。対数は自然数の逆数和。

F

DP。「A:頂点0から行き来可能な頂点」「B:未訪問の頂点」「C:訪問済みだけど行き来可能にはなってない頂点」の状態で考える。

この3つをキーにすると空間O(N3)時間O(MN3)っぽく見えるけど、1つのキーはNから他の2つのキーの値を引けば求まるのでNが一個落ちる。

Aに移動する場合はCの頂点がAに追加される。Bに移動する場合はBが1つCに遷移する。Cに移動する場合は変化なし。という感じでDPの漸化式を求めるとOK。

H

この時点でEの満点を考えつつ順位表を見るとHの部分点があることに気付いたので着手。

部分点は数列A上でのゲームなのでDP。

dp[i] := 自分がi、相手がi+1の場所にいる場合の最適なスコア と定義して漸化式を考える。

j>=i+2にジャンプすることが出来るので、大きい方から求めていってmaxを取っていけば良い。segment treeは不要。

残り時間

E満点、H満点、Gを考えていました。

E満点については「数を掛け算してN以上を作る」という発想はあったのですが、DPの考えが抜けず割り算をするDPだと思いこんでいたのが失策。

H満点はDPの式をもうちょっといじると満点に繋がったらしいので考察力が足りないんだなぁと実感しました。

Gは考察はだいぶいい線いってたと思うんですけど考察と実装どっちがミスっていたのか今でもはっきりわかってないです(Gを通せていたら20位or21位だったのでかなりへこんでいる)。

無事パーカーをもらえました。今回は事故が多発したらしくパーカーボーダーが高かった。

夜

7歳で競プロ始めたガチプロやばい。

講演の後は東工大15勢で🍕🍕🍣🍣🍕🍣🍣したり書道したりしていました(エンディングムービーに書道の光景が若干写っていた)(phi氏は賞ももらっていた)。

エキシビションパラレルも一応出ていたのですが、脳が停止していたので考察漏れ多数でダメでした。

その後はリレーで英語話せない辛いを多発しつつホテルに行きました。

今年は短縮王も無くやることがTwitterとTypingWarぐらいしか無かったのでポットの水を沸かしてお茶をいれることにしたのですが、

WAをたくさん出しました。

2日目朝

最近はリポビタンDにはまっています。

トーナメント

今年はあさプロではなくトーナメントが行われました。本戦順位が近い人同士で刺し合う血祭り。

30分2問なので早解き有利でした。

Round1及びRound2で0完でした。部分点すらない……

Round1 A

どう見てもMST+ダブリングLCAクエリです本当にありがとうございました。とD: たこ焼き屋とQ人の高橋君 - AtCoder Regular Contest 048 | AtCoderのソースを持ってきます。

めっちゃ古かったので使い方を完全に間違えてバグらせました。終了。

その後はパラレルで並走してましたが、Round2Aを眺めてはバグらせ、Round3Aを眺めてセグ木を書いてはTLE、スライド最小値を書いてはDPの初期化をミスってWAと散々でした。Round3Bは終わってから読んだのですがパット見D: Median Pyramid Hard - AtCoder Grand Contest 006 | AtCoder感がある。

昼

いくらと聞いてテンションが上がっていましたが速攻で完売。競プロ勢は🍣といくらが好物。

ご飯(とんかつ)を食べた後はchokudaiさんの生ハンドスプリングを見てrngさんの謎クイズを見て昼の部を終了しました。

リレー

お　　　ま　　　た　　　せ

今回は33位で海外勢抜いてチーム内最高順位なので一体どんな難しい問題が来るんだとめっちゃ緊張してました。

私はJを担当しました。リレー中の流れとしては、Jを考えた後、Iをみんなで考えていました。

J

N*Nのチェスボード( (i,j)はi+jが偶数なら黒、奇数なら白に塗られている)があるので、適切に白いマスを黒く塗って黒いマスを連結にしろという問題。N<=1000、黒く塗れるマスの数の制限が170000個。

これは割と簡単で、一番上の行と一番右の列をまず黒く塗って(N個ぐらい)、それから3列に1列を黒く塗っていけばOK(N2/6個ぐらい)。合計すると1000+166666=167666で問題なし。(行と列の内半分はすでに黒く塗られている点に注意)

I

N<=100で、条件を満たす写像をN-1個作れという問題。f_i(x)を写像とすると(1<=i<=N-1, 1<=x<=N)、

  • 全てのxについて、 f_i(x) != x
  • 全てのxについて、 i != j ならば f_i(x) != f_j(x)
  • 全てのiについて、 f_i(f_i(x)) != x

というのが条件。

解けてないのと疲れたので詳細には書かないのですが、素因数分解して分割統治すると解ける形にはなりましたが、実装量がめちゃくちゃ重い……

後で聞いたところもっと簡単な方法があるらしいです。

Kは海外勢の方におまかせしていたのですがなんかバグってしまったらしい。もうちょい議論しておくべきだった……

コンテンツ表彰

書家さんが今年はめちゃくちゃおもしろいタイプの人でした。今後も腕立てしてから書いて欲しい。

と思ったら芸人のおさるって人だったんですね……気付かなかった……(小さい頃に見た記憶がある)

解散

今年は懇親会的なやつないんですね……

感想

去年と比べるとかなり成長できたと感じられた点と、まだまだ詰めが甘いと感じた点があって、微妙な気持ちになりました。特に国内だけなら賞金圏内だったことを考えるとかなり悲しい。

来年こそは20位に入れる実力を持って挑みたいです。

運営の方々お世話になりました。ありがとうございます。

ちなみに来週はDDCC2016の方に行く予定です。よろしくお願いします。