草体にぼ日記

だらだらと

宣誓 私が大統領になったら bitDPを使えないプログラマーを処刑します.

2021年8月12日

宣伝
にぼしの解いた問題たち,こっちもちゃんと更新してるよ.

今日のいっきょく

やまないかな
あめ,やまないと
いや,だろ
梅雨の一句
そういえば「常体(タメ口?)」と「敬体(ですます調?)」を混ぜた日本語は読みづらいらしいです.かめさんが言ってました.
明日から意識します.
僕はブログを上から下まで書いた後推敲しません.めんどくさいので.
だから読みづらいかもしれません.

目次

  1. Coprime Present
  2. Joisino's traveling
  3. Traveling Salesman among Aerial Cities

今日は…

昨日・すぬけそだての記事に引き続き,すぬけそだてーーごはんーーと,その類題,あとは他に最短経路を求める問題をbitDPで解いてみようと言ったことをします.

ABC195 F- Coprime Present

[嬉しい,素数の個数,bitDP,UnionFindではない]

昨日やったすぬけそだてーーごはんーーの制約が厳しい版です(情報提供・あさか,こいつは本番で通した.犯罪だね).

bitDPは 「何個目まで見て」っていうのが意外といらなかったりする この問題では関係ないけどpopcount(s)で今まで何個追加してきたがわかったりするから この問題もそうで 素因数の集合がsでそこにxを足せるかどうかはsに何個足してきたかが大事ではなく互いに素かどうかが大事だから (引用元:あさか)

やばい,解いてから1日経ったから少し記憶が抜けてる.
とりあえず,前回のすぬけそだてーーごはんーーで使ったコードはTLEになります.

前回のコードだと,-Aが72で,72以下の素数は20個あります.計算回数は最大で,220 * 72 * 20 かな?まぁ,109超えるから無理って感じか.


あさかのヒントを元に,とりあえず「前から何個みたか重要じゃない」らしいから,dp[i] [bit] みたいな感じだったやつを, dp[bit] っていう形にします.
DP[bit] := 選んだカードの集合にbitが示す素数で割り切れるもの『しか』ないときの場合の数
『しか』ないとき,みたいな感じでめんどくさいことになってる理由は,僕がそうしないと実装出来なかったからです.

ll dp[1<<20];
int main()
{
    cin.tie(0);
    ios_base::sync_with_stdio(false);
    ll a, b;
    cin >> a >> b;
    ll n = b - a + 1; 
    vector<int> prime(0);
    for(int i = 2; i < n; i++){
        bool p = true;
        for(int j = 2; j * j <= i; j++){
            if(i%j == 0) p = false;
        }
        if(p) prime.push_back(i);
    }
    int psz = prime.size();

    vector<ll> A(n);
    rep(i,n){
        ll num = i + a;
        ll sum = 0;
        rep(j,psz){
            if(num % prime[j] == 0){
                sum += (1<<j);
            }
        }
        A[i] = sum;
    }

    dp[0] = 1;
    vector<ll> mul((1<<psz),1);
    rep(j,n){
        rep(i,(1<<psz)){
            ll tmp = dp[i];
            //debug(i);
            if((i & A[j]) == 0){
                if(( i | A[j] )== i){
                    mul[i] *= 2;
                } else dp[i|A[j]] += tmp;
            }
        }
    }
    ll ans = 0;
    //cout << endl;
    rep(i,(1<<psz)){
        ans += dp[i] * mul[i];
    }
    cout << ans << endl;
}

最後のans に+= dp[i] * mul[i]みたいなことをしてますが,こいつのせいでDPの定義の日本語がわけわからなくなりました.
とりあえず,管理する素数が{2, 3, 5} みたいなとき,bit が2進数で011 のときは,DP[110] := 選んだカードの集合には,3で割り切れる数と5で割り切れる数しかないときの場合の数って感じになってます.
それを満たすようなカードの集合は, {3,5} {15} {5, 9}みたいなものがあります.A,Bは今は考えるのがめんどくさいので考えていませんが,まぁ,上で上げたような3つの集合は,3で割り切れる数と5で割り切れる数しか入ってないですね.
このようにDPを定義すると,各bitについて,bitの変更をせずに加えられるものというものが何個あるかで,各bitのときに何通りあるかが求められます.
つまり,dp[110]に対しては, 7や11,13と言ったような3つは,いずれも2,3,5,で割り切ることが出来ないため,{3,5},{15}{5,9}といったカードの選び方に7,11,13は自由に加えることが出来ます.(すぬけそだて風にいうと,7,11,13を食べてもすぬけくんは悲しまない)
3で割り切れる数と5で割り切れる数だけカードを選ぶときのカードの選び方は,DP[110] * 23になります.これについて説明をすると,{3,5}, {15}, {5,9}のどの集合を使うかで3通り.次に,7,11,13はそれぞれこの集合に加えるか加えないかで2通りあります.だから,3 * 23 で求まる.多分.


やばい,ブログ書いてるけどかくきが起きない

俺の中で110っていうのは,5で割り切れる,3で割り切れる,2で割り切れるものは入ってない.っていうのを示すものなんだけど,見てる人には暗号になっているだろうね.

まぁDP[110]は,で割り切れる,3で割り切れる,2で割り切れるものは入ってない.また,余計なもの(7や11,13)も入れないで今は考える.みたいな感じ・・?

とりあえず,dp[110]のときにやった動きを実装するために動いていきましょう.

まずは,A以上B以下を割り切る数として,どのような素数を考えればいいかをprimeで持ちます.

 ll a, b;
    cin >> a >> b;
    ll n = b - a + 1; // a ~ bの数字はb-a+1個
    //2~n-1の素数を管理すれば良い.
    vector<int> prime(0);
    for(int i = 2; i < n; i++){
        bool p = true;
        for(int j = 2; j * j <= i; j++){
            if(i%j == 0) p = false;
        }
        if(p) prime.push_back(i);
    }
    int psz = prime.size();

で,そう,A以上B以下の数について,それぞれ,何で割り切れて何で割り切れないのかを簡単に見る方法がほしいので,A以上B以下の数の変換を行います.

例えば,見る素数が{2,3,5}のときは,次のように変換されます.

5 -> 100 (4)
6 -> 011 (3)
7 -> 000 (0)

8 -> 001 (1)

9 -> 010 (略) 10 ->101
11 -> 000
12 -> 011
13 -> 000
14 -> 001
15 -> 110

こんな感じです. X -> x2, x1, x0 みたいな感じに変換されますが,xiが1のとき,prime[i]の数でXが割り切れることを示しています.

 vector<ll> A(n);
    rep(i,n){
        ll num = i + a;
        //debug(num);
        ll sum = 0;
        rep(j,psz){
            if(num % prime[j] == 0){
                sum += (1<<j);
            }
        }
        A[i] = sum;
    }

A[0]はAを変換した姿.A[n-1]はBを変換した姿が入ってます.

DPします.
まず,dp[000] は 1 です.(何も選ばないという1通り).

 dp[0] = 1;
    vector<ll> mul((1<<psz),1);
    rep(j,n){ // A ~ Bを回す
        rep(i,(1<<psz)){
            ll tmp = dp[i];
            //debug(i);
            if((i & A[j]) == 0){
                if(( i | A[j] )== i){
                    mul[i] *= 2;
                } else dp[i|A[j]] += tmp;
            }
        }
    }

まぁ,わけわからんっすよね.
例を出します.001の状態から,何かしら状態が変わって遷移できる先とし考えられるものは,111や,101,011があります
意味は,2で割り切れるものしか入っていない状態から遷移できる状態が「2,3,5で割り切れるものしか入っていない状態(ただし,すぬけくんは悲しまないようになってる)」「2,5で割り切れるものしか入ってない状態」「2,3で割り切れ…」

って感じです.2で割り切れるものしか入ってない状態から(010)のように,2で割り切れず,3で割り切れるものが入っている状態にはなりませんよね.

例えば,選んだカードの集合が{2}なのに,新たにカード増やして{2,X}になったら,{2,X}の中身に2で割り切れるものがない状態になった!なんてならんよな.(2があるんだから)
で,ある状態bitから,使える数字についてだけど,まずprimeに含まれていない素数(7,11,13)と言ったものは使えなくて,コッチはmulで管理しないといけません.というのも,7,11,13と言ったものは状態bitに変更を加えません.これをDPで管理できるようには今出来ていないので,仕方ないのでmulで管理します.

次に,(001)という状態から,2で割り切れる数字も使えません.この判定は変換をしたおかげでO(1)で出来ます.

if( (i&A[j]) == 0 ) {
    // true なら使える.
}

例えば,3は変換後の姿が010で,bitが001のとき,010 AND 001 は 0になる.

6は変換後の姿が011で,bitが001のとき, 001 AND 011 は 001(not equal 0)で,使えません.
ここで,bit AND primeに含まれてない素数 のときに問題が発生する(俺の解き方では,使っちゃいけないはずなのに使われてしまう)ので,どうにかします.

  rep(j,n){ // A ~ Bを回す
        rep(i,(1<<psz)){
            ll tmp = dp[i];
            //debug(i);
            if((i & A[j]) == 0){
                if(( i | A[j] )== i){ // このif文がどうにかした結果
                    mul[i] *= 2; // 使っちゃいけないときはこっち
                } else dp[i|A[j]] += tmp; // 使うときはこっち
            }
        }
    }

primeに含まれていない素数は,変換後の姿が000です.ANDを取ると必ず0になってしまう.変換後の姿が000ってことは i と ORを取れば iになるね,だからそれが判定だね.っていうか

if( A[j] == 0) 

でいいじゃん.草.

これをすると何が起こるかというと,j=7,11,13のときには,mul[bit] *= 2;が行われるようになります.
これによって,全てのbitについてmul[bit] = 8; (2倍にする.を3回行ったから)

 else dp[i|A[j]] += tmp; // 使うときはこっち tmp = dp[i] 

これは,001から3は使えるよ.みたいなときの遷移ですね.

bit(変数名はi) = 001 (2の倍数のみある)状態が2通り例えば{2},{4}ある時,bit' = 011 に遷移するときを考えると,{2}に3を加えるか,{4]に3を加えるかで二通りの選択肢があるよね.
だから,dp[i|A[j]] += dp[i]にします.(tmpじゃなくてdp[i]でいい)

状態は bit OR A[j] になるよねぇ…

最終的な答えは,全てのbitについて dp[bit] * mul[bit] で求められます.


提出AC
乗り気しないのに記事書いたからちょっと雑になっちゃったかな…ごめんね。

ABC073 D- Joisino's travel

[ワーシャるフロイド法,next_permutation,bitDPによる巡回セールスマンっぽい]

まずは,各頂点間の最小移動コストを計算する.(ワーシャルフロイド法)
ここまでは共通.その後,2つの方法がある.

  1. ネクストパーミュテーションを使って,都市の並びを全探索
  2. bitDPで解く

ネクストパーミュテーションの方は簡単.実際配列が{0,5,2,3,1,4}みたいなときは,「0 -> 5のコスト + 5 -> 2のコスト + ... + 1 -> 4 のコスト」を求めてやればいい. 各配列について「〜のコスト+〜のコスト+...」を求めて,そのうちの最小値を答えとしてやればいい.

 int n, m, R;
    cin >> n >> m >> R;
    vector<int> r(R);
    rep(i,R) cin >> r[i];
    rep(i,R) r[i]--;
    rep(i,n)rep(j,n) dist[i][j] = INF;
    rep(i,n) dist[i][i] = 0;
    rep(i,m){
        int a, b, c;
        cin >> a >> b >> c;
        --a, --b;
        dist[a][b] = c;
        dist[b][a] = c;
    }
//ワーシャルフロイド法. 各頂点間の移動コストを最小にsルウ.
rep(k,n)rep(i,n)rep(j,n)chmin(dist[i][j], dist[i][k] + dist[k][j]);

    SORT(r);
    ll ans = INF;
    do{
        ll now = 0;
        rep(i,R-1){
            now += dist[r[i]][r[i+1]];
        }
        chmin(ans,now);
    } while (next_permutation(r.begin(),r.end()));
    cout << ans << endl;

行数少ないと,コード全体が俯瞰できて嬉しいので,ワーシャルフロイド法は一行で書けるようにしたほうがいいと思います.

俯瞰って日本語の意味がわかりません.僕は「全体を眺めることができる」という意味で使ってみました.

  • bitDP使う場合

ワーシャルフロイドで各長点間の距離を求めるまでは一緒なので省略します.

ll dp[(1<<10)][10]; // 8頂点しかないからまあ,2^8より多けりゃよし.
int main(){
    // フロイドまで省略    
    // DP表の初期化
    rep(i,(1<<n)){
        rep(j,10){
            dp[i][j] = INF;
        }
    }
    rep(j,n){
        dp[0][j] = 0; //どこも訪れていないときjをスタート地点に選ぶ.こうしとけば,dp[0010][2]みたいな感じに,1つ頂点を訪れたときのDPが正しく求まる.
    }
    
    for(int tmp = 0; tmp < (1 << n); tmp++){
        rep(j,n){
            //dp[tmp][j] からの遷移を考える
            // k へ移動する r[j] -> r[k];
            rep(k,n){
                chmin(dp[tmp|(1<<k)][k], dp[tmp][j] + dist[r[j]][r[k]]);
            }
        }
    }
    ll ans = INF;
    rep(j,n){
        chmin(ans, dp[(1<<n)-1][j]);
    }
    cout << ans << endl;
}

DP配列は,DP[訪問済みの国をbit的に管理] [最後に訪れた頂点] って感じにします.DP[1001] [1]は,頂点0,3を(少なくとも)訪問済みで,今頂点0にいることを意味しています.

少なくとも,というのは,頂点0,3を訪問する過程で1や2も訪問している可能性があるからです.

DP[1001] [2]はありえません.今頂点2にいるならDP[1011] [2] なはずだろふざけるなぶち殺すぞ.
(DP[0] [0~n-1]を0としているので,上の文は正しくなかった.)

ABC180 E- E - Traveling Salesman among Aerial Cities

[ワーシャルフロイド法,bitDP,巡回セールスマン]

Joisino's travelのように,next_permutationを用いていたら実行時間が間に合わないっすね.17の階乗は,10を8回以上かけているのでどう考えても間に合わないです.

JoisinoのbitDPの方を意識してやれば解けます.頑張ってください.


万が一俺のJoisinoでbitDPちょっと分かった…って人にとっては,俺のACコードがあると嬉しいと思うので(書き方がJoisinoと似てるから)一応貼っておきます. あ,俺のACコードって部分がリンクになってるからクリックしてね.