ミケ猫の小部屋

情報学、数学、語学および料理（暫定）について発信します

ICPC2020国内予選のFに挑んでみた！

11/06にICPC国内予選に参加して、三問完で何とかギリギリ予選突破しそうなみけ（チームxjubichanx）です。試合最中にFが昔のJOIのある問題とよく似ていることに気づきながらも解けなかった。

挑戦の第一歩：基本方針

始まる前に、まずは問題を目に通してください。

icpc.iisf.or.jp

一番ナイーブな方法は、災難レベルを0から徐々に上げていく。そして災難レベルごとにグラフのコンポーネントを計算して、サイズが一番大きいコンポーネント以外のコンポーネントをフィルターアウトしていく。災難レベルごとにDFSする必要があるため、計算量は恐らくO(M(N+M))になる。通称TLEってことです。

では、コンポーネントを増量的に探していこう。災難レベルを段々に上げて、そして災難レベル以下のエッジを削除してコンポーネントを段々小さいコンポーネントに分けていくことは難しいだろうと思う。なぜなら、一つのエッジを削除しても他のエッジが二つのサブコンポーネントの間にあるかもしれない。一つのエッジがある二つのコンポーネントを繋ぐ最後の「橋」を確認することは時間がかかる（DFSする必要があるので、計算量はさらにN+Mがかけられる）。

そこで、私は閃いた！災難レベルを初めから高く設定して、各頂点を最初から孤立する。その後、災難レベルを徐々に下げて頂点たちを繋いでいく。なので、推定強度の高いエッジから見に行って、エッジをグラフに加入していく。最初はN個の頂点それぞれ自分のコンポーネントにいるけど、エッジの加入によって小さいコンポーネントは段々大きいコンポーネントになる。 加入によるコンポーネントの融合タイミングはO(1)で判別できる

もし一つのエッジの二つの頂点は違うコンポーネントにいるなら、そのエッジの加入は二つのコンポーネントの融合に繋ぐ。つまり頂点が同じコンポーネントにいるかを判別し、融合をするタスクがメインです。この場合はUnion-FindでO(1)の判別と融合ができる。

挑戦の第二歩：エッジを昇順で見に行く

先までは推定強度の高いやつからコンポーネントを探していましたが、その方法で探したコンポーネントには致命的な冗長がある。

バックから見たコンポーネントは必ずしも存在しない*1

その問題を解決するため、最初は後ろから見て各辺は二つのコンポーネントをつなぐかどうかそして二つのコンポーネントのサイズを記録する。昇順でエッジを削除する時、コンポーネントのマージは分裂になる。その記録により、分裂のタイミングと分裂したコンポーネントのサイズはすぐにわかる。

昇順でエッジを見に行く最中、各災難レベルでの最大コンポーネントサイズはすぐにわかる。その後すぐにその推定強度のエッジをもう一回見て、分裂すべきエッジ*2なら分裂した後のコンポーネントのサイズによる頂点の削除を行う。分裂した後のコンポーネントのサイズがこの災難レベルでの最大コンポーネントのサイズより小さいだったらすべての点にタイムスタンプを付ける。

挑戦の第三歩：チェックメート

残ることはどうやら全ての頂点を見てタイムスタンプのもっとも大きいやつを出力するだけだと思ったら、サンプル出力に合わないところが現れた！もし、答えとなる集合は推定強度のもっとも高いエッジで繋がっているならどうなる？たぶん最後までも「削除」されないだろう。そのエッジを潰したら、残るコンポーネントはどれもサイズ１なので、最後までも削除されない。そのため、タイムスタンプが0（削除されないためタイムスタンプがおかしい）の頂点には一番大きいタイムスタンプを付けておきたい。

その最後の一歩が完成すれば、残りは小さい順から頂点を見てタイムスタンプの一番大きいやつを出力する。

コード

コンテスト中には解けなかったので、正解かどうかはわからないまま。ここでコードを貼り付けて、間違いがあれば指摘お願いします。

#include <cstdio>
#include <queue>
#include <utility>
#include <memory>
#include <map>
#include <algorithm>

using namespace std;

const int MAXN=1e5+7;

int UF[MAXN];
int sz[MAXN];

int parent(int n){
    // printf("%d\n",n);
    if (UF[n]==0) return n;
    else {
        auto ret=parent(UF[n]);
        UF[n]=ret;
        return ret;
    }
}

using E=pair<int,pair<int,int>>;

int N,M;

vector<pair<int,int>> graph[MAXN];
E edges[MAXN];//s a b
pair<int,int> partition_sz[MAXN];//sza szb

int height[MAXN];
int h,w;

void DFS(int cur,int c=0){
    if (c>=3) return;
    auto &vis=height[cur];
    if (vis) return;

    vis=h;
    // printf("height[%d]=%d h=%d\n",cur,height[cur],h);
    for(auto adj:graph[cur]) {
        if (adj.second<=w) continue;
        DFS(adj.first,c+1);
    }
}

bool judge(const map<int,int>& m){
    return m.crend()->second==1;
}

int main(){

    while(scanf("%d%d",&N,&M)&&N){
        
        map<int,int> count;
        count[N]=1;

        h=0;
        for(int i=1;i<=N;i++){
            sz[i]=1;
            UF[i]=0;
            graph[i].clear();
            height[i]=0;
        }

        for(int i=1;i<=M;i++){
            int a,b,s;
            scanf("%d%d%d",&a,&b,&s);
            edges[i]=E{s,{a,b}};
            graph[a].push_back(pair<int,int>(b,s));
            graph[b].push_back(pair<int,int>(a,s));

            partition_sz[i]=pair<int,int>(0,0);
        }

        sort(edges+1,edges+M+1);

        int i=M;
        while(i>0){
            auto weight=edges[i].first;

            int j=i;
            while(j>0&&edges[j].first==weight){
                auto edge=edges[j].second;
                auto pa=parent(edge.first),pb=parent(edge.second);
                auto &part=partition_sz[j];
                j--;
                if (pa==pb&&(pa!=0)) continue;


                auto sza=sz[pa],szb=sz[pb];
                
                // printf("sza=%d szb=%d\n",sza,szb);
                part=pair<int,int>(sza,szb);
                
                UF[pb]=pa;
                sz[pa]+=sz[pb];
            }
            i=j;
        }


        //昇順
        i=1;
        while(i<=M){
            auto weight=edges[i].first;
            w=weight;
            h++;

            int j=i;
            // for(auto kv:count){
            //     printf("k=%d v=%d\n",kv.first,kv.second);
            // }
            while(j<=M&&edges[j].first==weight){
                auto edge=edges[j].second;
                auto a=edge.first,b=edge.second;
                auto &part=partition_sz[j];
                j++;
                //このエッジがコンポーネントを二つに分けていないなら次
                if (!part.first||!part.second) continue;
                if (height[a]||height[b]) continue;

                auto sza=part.first,szb=part.second,sum=sza+szb;
                
                // printf("a=%d b=%d sza=%d szb=%d\n",a,b,sza,szb);
                if (--count[sum]==0) count.erase(sum);
                count[sza]++,count[szb]++;

            }

            auto max_size=count.rbegin()->first;
            
            // printf("max_size=%d\n",max_size);
            // puts("loop");
            for(int k=i;k<j;k++){
                auto edge=edges[k].second;
                auto a=edge.first,b=edge.second;
                auto const &part=partition_sz[k]; 

                if (!part.first||!part.second) continue;
                if (height[a]&&height[b]) continue;
                auto sza=part.first,szb=part.second,sum=sza+szb;
                
                // puts("dfs");
                if (sza<max_size&&!height[a]) 
                    DFS(a);
                if (szb<max_size&&!height[b]) 
                    DFS(b);
            }

            auto occurences=count.rbegin()->second;
            count.clear();
            count[max_size]=occurences;
            i=j;
        }


        for(int i=1;i<=N;i++) if (height[i]==0) height[i]=h;
        auto max_height=(*max_element(height+1,height+N+1));
        
        
        for(int i=1;i<=N;i++) if (height[i]==max_height) printf("%d ",i);
        putchar('\n');

    }

    return 0;
}

*1:

例として、1級災難AとBの二つのコンポーネントがいるとする。Aは7個の頂点からなる。Bは6個の頂点の頂点からなる。なのでBはAより不安全。でも、災難レベルがさらに上げると、Aは縮小して3個の頂点になる。Bはまた6個の頂点からなる。

バックから見れば、Bはどのコンポーネントの一部だったとかの情報はわからない。なのでBが存在しないコンポーネントということは知らない。

*2: バックからのトラバーサルで分裂するとマークされる以外も、削除された二つの頂点の間にいないことも重要。前言った存在しないコンポーネントを排除するためです

【自称】世界一分かりやすいBITの解説

こんにちは！一ヶ月ぶりのブログ更新です～今度はBITを解説しに行きます。実は一年前ほど研究室でBITとセグツリーを解説しましたが、あの頃解説はちょっと不明なところもあったので、このブログのきっかけになります。

BITはRSQ問題を解決するデータ構造で、O(logn)の計算量で区間の要素の総和を求められる。そして単独な要素をO(logn)計算量で更新できる。一回だけある区間の総和を求めるなら、愚直に足し合わせるだけで済みますが、もし数万回異なる区間の総和を問われったら、BITは凄く早い。

BITの基本概念は分割統治

一番簡単な考え方は区間を二つの部分区間に分けて、それぞれ総和を管理する。ある区間の総和を求めるとき、部分集合の総和を加算して素早く求まる。でも、区間総和の特徴を考えれば、本当は二つの部分区間の総和を記録する必要はない。

例えば、区間[1,8]上の総和と[1,4]上の総和が知ったら、[5,8]の総和は[1,8]の総和から[1,4]の総和を引けば求まる。だから、BITが管理する各区間は下図に示したようになる。

f:id:cre_chan:20200808225945p:plain
BITは半分の区間しか記録しない

配列のインデクスで区間を決める

BITは処理を単純化するため、長さが２のべきの区間だけを扱える。つまり、区間の長さが1、２、４のような区間ばっかです。インデクス1からこうスペース空けて区間を配置していく。区間の右側はそれぞれ違うので、一次元配列で表現できる。配列のi番目のはiを右側とする区間の総和を格納する。

f:id:cre_chan:20200808230003p:plain
配列要素の管理区間ビジュアル

そして、i番目の要素が管理する区間の長さは必ずiの因数の中で一番大きい2のべきだ。これからそのことを証明して行く。簡単な例で説明する。12は二進数表現で00001100になる。もしこれを参照としてBITをアクセスしたら、長さ4の区間[9,12]の総和は得られる。それはなぜでしょうか？2のべきは二進数で表現したら、一ビットだけが1になる。2のべきである最大因数はiを二進数に変換したら、最下位の1まで切り取って得られる。12の場合区間の長さは4で、二進数表現は00000100になる。4は12の中に三個含まれる。１や２は12に偶数個含まれる。

前説明したが、BITはそれぞれの長さの区間をスペース空きで配置する。偶数番目の区間は飛ばされるので、iは必ず奇数番目の区間を管理する。こうして、証明できた。これから実用篇に入ります！

区間の長さを求める

i番目要素は長さがi&-iの区間を管理している。証明は単純なのでこのことは読者に証明してもらう。i&-iをlen(i)と呼ぶ。

BITの構築

BITを構築する時はまず配列で累積和を求める。i番目に1からiまでの要素の総和を入れる。でもそれじゃ余計な部分も入ってしまうので、その余計な部分を引く必要がある。i番目の要素は[i-len(i)+1,i]を管理する。それて、[1,i-len(i)]の部分は余計だ。その部分の総和は今i-n番目の要素だ。それて、配列の後ろから頭へスキャンしてn個前の要素を引けばBITを構築できる。

f:id:cre_chan:20200808230021p:plain
BITの構築プロセス

BITで1からiまでの総和を求める

[1,i]までの総和を求めるには、まず[1,i-len(i)]の総和を再帰的に求める。その結果をBIT配列のi番目の要素に足すと結果になる。例え[1,13]の総和を求める時に、順次に[1,12]、[1,8]の総和を求めていく。従ってO(logn)で求まることはわかる。

f:id:cre_chan:20200808230052p:plain
1から5までの総和を求める

[1,i]の総和を求めれば、[l,r]の総和を求めるのも簡単だ。[1,r]から[1,l-1]を引けばわかる。

ある要素を更新

i番目の要素を更新したとき、i番目の要素を含む区間も更新されるべき。BITのi番目要素を更新した後、再帰的にi+len(i)番目を更新する。一番目の要素を更新するビジュアルを作った。

f:id:cre_chan:20200808230036p:plain
一番目の要素の更新

待ちくたびれたコードセッション

このセッションは作成したコードを上がる。コンテストなどで利用することはできる。

template<class T>
void build(T arr[],T bit[],unsigned int n){
    bit[0]=0;
    for(auto i=2u;i<n;i++) bit[i]=bit[i-1]+arr[i];
    for(auto i=n-1;i>0;i--) bit[i]-=bit[i-(i&-i)];
}

template<class T>
void update(T bit[],unsigned int n,unsigned int index,T diff){
    for(auto i=index;i<n;i+=(i&-i)) bit[i]+=diff;
}

template<class T>
T prefix(T bit[],unsigned int bound){
    T ret=0;
    while(bound>0){
        ret+=bit[bound];
        bound-=(bound&-bound);
    }

    return ret;
}

template<class T>
T query(T bit[],unsigned int l,unsigned int r){
    return prefix(bit,r-1)-prefix(bit,l-1);
}

BITが解決できる問題

AOJ DSL_2_B　 直接解決できる問題はRSQ問題で、BITだけで楽勝だ！

ABC 174F Range Set Query 異なる色が出現するところを1にして、ある区間に出現した色の種類はBITで集計できる。

ベクトルから見るフーリエ級数

「フーリエ級数は周期関数を周期の異なる三角関数に分解する」と学部の一年から教え込まれたが、その公式がなかなか覚えられないダメ人間です（笑）。たぶん私と同感してる人はかなり多いと思います。この記事は公式を頭の中に焼き付けるためでも、一時の発想を記録するためでも存在する。

フーリエ級数の復習

随時参照できるように、まずは公式をここに乗せる。周期が Tの関数 f(x)があるとする。 f(x) = \sum_{k=-K}^{K} c_k e^{ i\frac{2k\pi}{T} x } のように展開できる。そこで、 c_k = \frac{ \int f(x)e^{ -i\frac{2k\pi}{T} x }dx }{T} と置く。ブログの表示問題で、ここは実際0からTまでの定積分となる。

なんだその謎の c_k?

その謎を解ける前に、まず線形代数を復習しましょう！平面上に二つのベクトル \vec{a}と \vec{b}があるとすると、下の図はその状態を示した。

f:id:cre_chan:20200630151821j:plain
ベクトルの内積

まずは \vec{a}と \vec{b}を同一の始点に揃えて、 \vec{a}の終点から \vec{b}に垂線を引く。上の図の点線はその垂線を表す。そして、その垂線の切片を \vec{a}が \vec{b}にの投影という。投影の長さと \vec{b}の長さの積はまさに \vec{a}と \vec{b}の内積だ。公式で表すと \vec{a} \cdot \vec{b} = \cos \theta |\vec{a}| |\vec{b}| になる。その公式見て、 \vec{b}の長さを1に置く時に何が起こる？内積は \vec{a}の投影の長さになった！

普段ベクトルを言えば、常に多くの数が縦に並べるイメージですが、それらの数は実このベクトルと各軸の正の方向に向いてる単位ベクトルの内積だ。下の図のようにベクトルは二方向のベクトルに分解された。x個の (1 , 0)とy個の (0, 1)の和に分解した。

f:id:cre_chan:20200630161137j:plain
ベクトルを分解する
(1,0)と (0,1)のような長さが1かつお互い直交するベクトルの組は分解することにぴったりだ。なぜなら、直交することはお互いが表現する部分に重なりがないことを意味する。x軸の単位ベクトルが表現することはy軸の単位ベクトルで表現できない。そしてある次元の成分を弄っても、他の次元に影響がない。そのような分解は人間にとっても理解しやすい。

フーリエ級数の話に戻ります。周期がKの関数 f(x)と g(x)の内積は \int f(x)\bar{g(x)}と定義される。ここでは表示問題のため、実際は0からKまでの定積分になる。その知識があれば、 c_kは実際元関数と \frac{e^{i\frac{2k\pi}{T} x}}{T}の内積を表してる。そして、 l \neq rの場合に&lt; e^{i2l\pi x},e^{i2r\pi x}>=0がある。 e^{i2k\pi x}は周波数の異なる振動を表すため、それは異なる周波数の振動がお互い直交してることを意味する。つまり、下の図のように、その関数を座標軸に見なすことはできる。

f:id:cre_chan:20200630171110j:plain
三倍振動までフーリエ展開する

この図の場合、元関数を何らかのベクトルとして扱う。積分により、各軸での投影が求まって、それらの振動の線形コンビネーションとして表すことこそフーリエ級数の正体だ。さらに、公式の分母にいる Tはベクトル e^{ i\frac{2k\pi}{T} x } の長さを1に抑えるための定数となっている。

難しそうなこと言って何が嬉しい

もし突然公式を思い出せなくなる場合、ベクトル空間を想像してください。それは { \frac{e^{i\frac{2k\pi}{T} x}}{T} }を基底にする空間で、元関数をを各基底に投影し、基底の線形コンビネーションとして表す。そうしたら、ちょっと覚えやすくなるでしょう。少なくとも私はそう思ってる。

雑談：なんだが最小二乗法と似てない？

フーリエ級数は最初二乗誤差を最小化する発想から導かれた級数なので、そのことを思いついた。級数なので、無限に足し続けると誤差はゼロにちかづいていく。最小二乗法の場合は次元数が限られるため、どうしても誤差は生じる。そして、 Ax=bの解を A^{T}Ax=A^{T}bの解で近似することは実際bをまず Aの基底が張る超平面にまず投影して解くことと解釈される。そのような色んな「証拠」からすると、私はフーリエ級数が最小二乗法と似てると思う。