くれなゐの雑記

例を上げて 自分で手を動かして学習できる入門記事を多めに書いています

Educational Codeforces Round 53 E. Segment Sum

問題

codeforces.com

問題概要

以下の条件を満たす[l, r]の範囲の数字の和を求める問題。(場合の数ではなく、数字そのものを足す点に注意する。)

条件: 10進数表記で、使用されている数字の種類がk 以下

方針

f(A): ある数字A以下の、条件を満たす数字の和 とすると、f(r)-f(l-1)でこの問題の答えは求められる。

こういう条件を数字の数え上げは桁DPがいい。そう古事記にも書かれている。というわけで桁DPで方針を立てる。

そもそも、数え上げDPは条件を満たすDAGのパスの総数を求める問題とも言いかえることができる(参考)。 桁DPの解説は個人的にこのページが分かりやすかった。

pekempey.hatenablog.com

ということで、上記のページの要領でグラフを構築し、各辺に対し、その数字の重み(右から3桁目の4なら400の重みがつく)を付け、条件に沿うDAGのedgeの合計値を求めれば良い。

状態としては、桁数、使用した数字(10桁のbitで管理)、A未満であることが確定しているかどうか の3次元。 10桁の数字のうち、どれを使用したかをメモするためには1024必要なので、合計で18*1024*2で間に合う。

最後の桁数に関して、使用した数字の数がk以下のもののみを数え上げ、出力する。

数え上げならば、これだけで良いのだが、今回は数字の合計を求めなければならないのでもうひと工夫が必要。 ちゃんと各数字ごとに使用される回数をかけてやらなければならない。 上流側は、数字を流せば自動的に分裂して数え上げしてくれるが、下流側は分裂する回数が足りないのでキチンと数え上げをしてくれない。 例えば、一番右のケタの1は10回しか呼ばれないので雑に計算すると桁数にかかわらず1*10で10になってしまう。 でも実際には、1000までの数字の和だったら100回 1 が呼ばれる。 下流側では、上流側の数え上げた分を反映しなければならない。

なので、DPはシンプルに数え上げるだけのものと、数字そのものの合計を求めるものの2つを用意して、下記ソースコードのように遷移した。

コーナーケースとして、最上位の桁に対して0を選択し続けた場合に関して、それは0を選択していないことになる。 (000023は23扱いなので、0は含まれていない)ので、注意する。

int dp[22][(1<<10)][2];
int sums[22][(1<<10)][2];
int solve(const string str, const int K) {
    memset(dp, 0, sizeof(dp));
    memset(sums, 0, sizeof(dp));
    vector<int> nums(str.size());
    REP(i, str.size()) {
        nums[i] = str[i] - '0';
    }

    vector<int> weights(nums.size());
    weights[0] = 1;
    FOR(i, 1, nums.size()) {
        weights[i] = (weights[i - 1] * 10LL)%MOD;
    }
    reverse(ALL(weights));

    dp[0][0][0] = 1;
    sums[0][0][0] = 0;
    REP(i, nums.size()){
        int weight = weights[i];
        REP(j, (1<<10)) {
            REP(k, 2) {
                int lim = k ? 9 : nums[i];
                REP(d, lim + 1) {
                    int jj = j;
                    if (j == 1) {
                        jj = 0;
                    }
                    (dp[i + 1][jj | (1 << d)][k || d < lim] += (dp[i][j][k])) %= MOD;
                    (sums[i + 1][jj | (1<<d)][k || d < lim] += (sums[i][j][k] + (d * (weight * dp[i][j][k])%MOD)%MOD)%MOD) %= MOD;
                }
            }
        }
    }

    int res = 0;
    REP(j, (1 << 10)) REP(k, 2) {
        bitset<10> bits = j;
        if (bits.count() <= K) {
            (res += sums[nums.size()][j][k])%=MOD;
        }
    }
    return res;
}

void solve() {
    int l, r, k; cin >> l >> r >> k;
    cout << ((solve(to_string(r), k) - solve(to_string(l-1), k)) + MOD)%MOD << endl;
}

OpenFOAMをairflowでジョブスケジューリングしてみる

動機

普段はCentOSにtorqueを使ってジョブスケジューリングして計算していた。 ある日普段使用しているPCも空いてる時間は計算を回そうと思い、torqueをインストールしようとしたが、Ubuntuのaptで入らなかった。 入れようと思えば入れれるが、せっかくなので最近のツールも使ってみようと思い、試してみた。

ジョブスケジューリングするツールの選択

調べてもよくわからなかった。 なんかモダンっぽいGUIで、並列実行ができて、Pythonでいろいろかけそうだったのでこの記事を見てAirflowを選んだ。 Python3に対応しており、最近もしっかりメンテナンスされているのもいい。

qiita.com

実行環境

kurenaif@ubuntu:~$ screenfetch 
                          ./+o+-       kurenaif@ubuntu
                  yyyyy- -yyyyyy+      OS: Ubuntu 18.10 cosmic
               ://+//////-yyyyyyo      Kernel: x86_64 Linux 4.18.0-10-generic
           .++ .:/++++++/-.+sss/`      Uptime: 29m
         .:++o:  /++++++++/:--:/-      Packages: 1637
        o:+o+:++.`..```.-/oo+++++/     Shell: bash 4.4.19
       .:+o:+o/.          `+sssoo+/    Resolution: 1920x983
  .++/+:+oo+o:`             /sssooo.   DE: GNOME 
 /+++//+:`oo+o               /::--:.   WM: GNOME Shell
 \+/+o+++`o++o               ++////.   WM Theme: Adwaita
  .++.o+++oo+:`             /dddhhh.   GTK Theme: Yaru [GTK2/3]
       .+.o+oo:.          `oddhhhh+    Icon Theme: Yaru
        \+.++o+o``-````.:ohdhhhhh+     Font: Ubuntu 11
         `:o+++ `ohhhhhhhhyo++os:      CPU: Intel Core i7-4790 @ 4x 3.6GHz [100.0°C]
           .o:`.syhhhhhhh/.oo++o`      GPU: svgadrmfb
               /osyyyyyyo++ooo+++/     RAM: 952MiB / 3921MiB
                   ````` +oo+++o\:    
                          `oo++.      
mysql: 5.7.24
airflow: v1.10.0
OpenFOAM: dev

Airflowのインストール

Airflowのインストール

リブセンスのAirflowの記事のページと、airflowの公式ドキュメントを参考にしてインストールした。

$ sudo apt install python3-pip
$ export  SLUGIFY_USES_TEXT_UNIDECODE=yes 
$ pip3 install apache-airflow

ここでairflowを実行しようとすると

$ airflow initdb
airflow: command not found

おや、新しい環境なのでパスが通っていないらしい pip3 show apache-airflowコマンドでパスを調べて… パスを通す

$ airflow initdb
$ airflow version
[2018-11-04 10:46:53,290] {__init__.py:51} INFO - Using executor SequentialExecutor
  ____________       _____________
 ____    |__( )_________  __/__  /________      __
____  /| |_  /__  ___/_  /_ __  /_  __ \_ | /| / /
___  ___ |  / _  /   _  __/ _  / / /_/ /_ |/ |/ /
 _/_/  |_/_/  /_/    /_/    /_/  \____/____/|__/
   v1.10.0

v1.10.0 が入ったらしい

mysqlのインストール、準備

今回は、リブセンスのAirflowの記事と同じ構成にするので、続けてmysqlceleryのパッケージをインストールしてゆく。

$ sudo apt install libmysqlclient-dev
$ sudo apt install mysql-server mysql-client
$ sudo mysql_secure_installation
$ pip3 install "apache-airflow[mysql, celery]"

これでとりあえず必要なものはインストールできたので、mysqlの設定をしていく airflowというユーザー名で作ると設定が楽なので、airflowというユーザー名と、airflowというDBを作る。(ちなみに、airflowのパスワードをairflowにすると一番何も考えずに設定ができるが、セキュリティ的にアレ)

雑なパスワード設定の方法は

qiita.com

を参照する。

mysql> create user airflow@localhost identified by {パスワード}
mysql> create database airflow;
mysql> grant all on airflow.* to airflow@localhost;

ここで、dbを見てみよう

$ mysql -u airflow -p
Enter password: 
(メッセージ略)
mysql> show databases;
+--------------------+
| Database           |
+--------------------+
| information_schema |
| airflow            |
+--------------------+
2 rows in set (0.00 sec)

OK ユーザーとデータベースの作成ができた。

airflowの設定

構成として、

とすべてMySQLの構成で行う。 これが一番設定が楽だった。

今回は、以下に示す設定を変更した。

executor = CeleryExecutor
sql_alchemy_conn = mysql://airflow:{パスワード}@localhost:3306/airflow
load_examples = False
catchup_by_default = False
worker_concurrency = 3
executor

今回はCeleryExecutorを使います。 Localでしか使わないので、LocalExecutorと悩んだのですが、なんか動かなかったので今回はこれで。

sql_alchemy_conn

CeleryExecutorではデフォルトのsqliteは使えません。

mysql://{ユーザ名}:{パスワード}@{アドレス}:{ポート}/{DB名}

というフォーマットで指定します。

load_examples

これは趣味なのですが、examplesをloadするとすごいloadされるので僕は切りました。

catchup_by_default

catchupという機能をデフォルトでonにするかどうかです。 リブセンスのAirflowの記事がわかりやすいです。

worker_concurrency

一つのワーカーの並列実行する最大のプロセス数だと思っています。 4スレッドなので3を与えました。合ってるかわかりません。

設定の反映

満足行く設定ができたら、設定を反映します

airflow initdb

このときに

Exception: Global variable explicit_defaults_for_timestamp needs to be on (1) for mysql

のようなエラーが出るかもしれません。その時は以下の記事を読んでみてください。

MySQL - 5.6 系で TIMESTAMP 型デフォルト値警告! - mk-mode BLOG

設定の確認

$ airflow webserver

コマンドでウェブサーバーを起動してみましょう。 ブラウザを立ち上げて、 ipaddress:8080 で管理画面っぽいやつが出るはずです。

f:id:kurenaif:20181105204523p:plain

計算するケースの準備

今回はOpenFOAMチュートリアルから5ケース使います デスクトップに、5ケース分置きます

$ cd Desktop/
$ cp -r ~/OpenFOAM/OpenFOAM-dev/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/ .
$ cp -r ~/OpenFOAM/OpenFOAM-dev/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct .
$ cp -r ~/OpenFOAM/OpenFOAM-dev/tutorials/multiphase/interFoam/RAS/damBreak/damBreak/ .
$ cp -r ~/OpenFOAM/OpenFOAM-dev/tutorials/combustion/chemFoam/gri/ .
$ cp -r ~/OpenFOAM/OpenFOAM-dev/tutorials/compressible/rhoCentralFoam/shockTube/ .

DAGファイルの作成

airflowでは、DAGファイルというpythonスクリプトを書いて、 airflow.cfg にかかれている dags_folder のフォルダに入れます。

airflowの公式ページのチュートリアルに簡単なスクリプトチュートリアルがあるので、 一度目を通してみてください。ここでは、いきなりOpenFOAMの計算を流します。

下記のソースコードを適当なファイル名で保存し、dags_folderに入れます。(デフォルトでは ~/airflow/dags になっていると思います。)

これでファイルの準備は完了です。

2点、チュートリアルと大きく違っており、

  • start_date は本来はn日前みたいな表記は非推奨なのですが、今回は1回実行でキャッチアップなしなので、これで大丈夫です。
  • schedule_interval は今回は一回実行のみなので、 @once を与えています。

また、今回はタスクの依存関係は特にないので、set_upstreamは不要です。 あとは基本的にこのファイルを見ればおおよそやっていることがわかると思います。 ファイルがちゃんとできているかどうかは

python3 file_name.py

とすると、チェックできます。 BashOperatorの変数への代入は、今回のケースは不要です

from airflow import DAG
from airflow.operators.bash_operator import BashOperator
from datetime import datetime, timedelta

default_args = {
    'owner': 'airflow',
    'depends_on_past': False,
    'start_date': datetime.today() - timedelta(7),
    'email': None,
    'email_on_failure': False,
    'email_on_retry': False,
}

dag = DAG('cfd-test', default_args=default_args, schedule_interval = '@once')

# t1, t2 and t3 are examples of tasks created by instantiating operators

aerofoilNACA0012 = """
source ~/OpenFOAM/OpenFOAM-dev/etc/bashrc
cd ~/Desktop/aerofoilNACA0012/
./Allrun
"""
task_angledDuct = BashOperator(
    task_id='aerofoilNACA0012',
    bash_command=aerofoilNACA0012,
    dag=dag)

angledDuct = """
source ~/OpenFOAM/OpenFOAM-dev/etc/bashrc
cd ~/Desktop/angledDuct/
./Allrun
"""
task_angledDuct = BashOperator(
    task_id='angledDuct',
    bash_command=angledDuct,
    dag=dag)

damBreak = """
source ~/OpenFOAM/OpenFOAM-dev/etc/bashrc
cd ~/Desktop/damBreak/
./Allrun
"""
task_damBreak = BashOperator(
    task_id='damBreak',
    bash_command=damBreak,
    dag=dag)

gri = """
source ~/OpenFOAM/OpenFOAM-dev/etc/bashrc
cd ~/Desktop/gri/
./Allrun
"""
task_gri = BashOperator(
    task_id='gri',
    bash_command=gri,
    dag=dag)

shockTube = """
source ~/OpenFOAM/OpenFOAM-dev/etc/bashrc
cd ~/Desktop/shockTube/
./Allrun
"""
task_shockTube = BashOperator(
    task_id='shockTube',
    bash_command=shockTube,
    dag=dag)

airflowの起動、実行

今回の構成では、以下の3つのairflowのツールを実行し、その後タスクをトリガーします。

3つターミナルを起動し、(systemdで裏で実行してもいいです)

airflow webserver
airflow worker
airflow scheduler

を起動します。それぞれ何を表すかは、リブセンスのairflowの記事がとてもわかりやすかったです。

localhost:8080 につなぐと、

f:id:kurenaif:20181105204541p:plain

のような画面があるので、

一番左のOFFになってるやつを ONにしましょう。 これでトリガーするとジョブが流れるようになります。

f:id:kurenaif:20181105204549p:plain

どんなタスクが中に入っているかは、ON/OFFのスイッチの右の cfd-test をクリックすると、Graph ViewやTree Viewでタスクの一覧がみえます

f:id:kurenaif:20181105204603p:plain

f:id:kurenaif:20181105204614p:plain

実行方法は、CUIとウェブページを使う方法があり、

CUIの方法は

airflow trigger_dag cfd-test

GUIの方法は、 ウェブページのLinksの再生ボタンみたいなやつを押せばトリガーされます。

f:id:kurenaif:20181105204620p:plain

実行が終わると、再生ボタンを押したページのRecent Tasksの数字が増えたり、DAG Runsの数字が増えたり、

Treeviewの緑のランプが増えたり

f:id:kurenaif:20181105204628p:plain

ガントチャートを確認できたり(途中でも確認できます)

f:id:kurenaif:20181105204634p:plain

します。

これで実行完了です。

angledDuctをみてみると…

$ ls
0  10  Allrun  constant  log.blockMesh  log.rhoPimpleFoam  system

できていますね。

所感

並列処理なジョブを回すときに、どうやって制限するのかよくわからなかったタスクのresourcesフィールドに、cpusみたいなのがあって、それに3とか入れても平気で他と並列処理してるしウーンという感じ

GUIでジョブの進行状況が見れたり、Pythonスクリプトで簡単にジョブファイルを作れるのはいいですね。 いろいろ応用できそうです。

今後も少し運用してみます。

OpenFOAMをallWmakeなしで必要なものだけwmakeするツール auto_wmakeを作りました

動機

最近、諸事情で様々な環境でOpenFOAMを使うことがおおくなりました。 OpenFOAMでソルバを使えるようにするといえば、Allwmake コマンドですが、このコマンドありえんくらい時間がかかるので、かなりのボトルネックになる可能性があります。

実際、すべてのユーティリティをmakeしなくてもよく、実は一部のコマンドはそこまで多くwmakeしなくてもいいのです。 (例えばicoFoamに燃焼系のライブラリは不要です)

そこで、必要なソルバを使えるようにするために、必要最低限のmakeで済ませるコマンドがほしいなと言うことで作りました。

auto_wmakeの使い方や紹介

github.com

コチラにおいてあります。 README.mdに書いてあるとおりだと思いますので、何かわからない場合はブログのコメントや知ってる方はtwitterからよろしくお願いします。 cargoで管理してるとかんたんにinstallできていいですね

技術的な話

なぜRust?

特に理由はありません。 書いてみたかった go langでも良かったかもしれないですが、個人的にRustのほうが好きなのでRustにしました。 Pythonはグラフ生成が思ったより遅かったのでやめました

生成されるファイル名、依存関係の抽出

各ライブラリどうしの依存関係は get_make_target

Make/options, 各ライブラリの生成されるファイルの名前は Make/files からうまく取ってきています。 例外的な文字列処理の塊で正直気合です

依存関係から依存関係グラフの作成

各ライブラリ同士の依存関係をグラフにいい感じにまとめます。 get_edgesで先程の関数を使って辺を抽出、 出てきた辺は 以下のコードのようにしてグラフにします。 graph が隣接リスト, in_degree は入次数を表しています。

グラフは、(依存する側)->(依存される側) と定義しており、入次数が大きければ大きいほど、いろんなライブラリに依存していることになります。

let mut graph : HashMap<String, Vec<String>> = HashMap::new();
let mut in_degree : HashMap<String, i32> = HashMap::new();
for edge in edges {
    in_degree.entry(edge.1.clone()).or_insert(0);
    let deg = in_degree.entry(edge.0.clone()).or_insert(0);
    *deg += 1;
    graph.entry(edge.0.clone()).or_insert(Vec::new());
    let from = graph.entry(edge.1.clone()).or_insert(Vec::new());
    from.push(edge.0);
}

wmakeをする始点の探索

入次数が0のところがwmakeの始点です。 get_zero_in_degreeがこれに対応しています。

wmakeをする順番

計算が終了したら、終了したディレクトリから辺が伸びている先のノードの in_degree を1減らし、0になったタイミングでキューに追加します。 これを繰り返し、最終的にキューが尽きたら終了です。

while let Some(target) = queue.pop_front() {
    cnt += 1;
    // output progress
    println!("[{}/{}] {}", cnt, size, target);
    wmake(&target, &jobs_string, is_stdout_detail);
    let nexts = graph.get_mut(&target).unwrap();
    for nxt in nexts {
        *in_degree.get_mut(nxt).unwrap() -= 1;
        if in_degree.get(nxt).unwrap() == &0 {
            queue.push_back(nxt.to_string());
        }
    }
}

例外的なwmake

wmakeそのもの、Pstretam, OSspecific に関しては少し特殊なmakeをしなければならないため、例外的にmakeを挟んでいます。 これのせいでおそらくいくらかのバージョンのOpenFOAMでは、auto_wmakeが動かなくなる可能性があります。

例外的なwmakeはinit_buildに対応しています。

その他機能

auto_wmakeは最初はコマンドライン引数に、makeしたいディレクトリのパスを与えていたのですが、icoFoamへのパスが長すぎてめんどくさくなったので、applicationsフォルダの中にあるものは省略可能にしました。 やり方としては、Make/filesの中身を見て、一致したものを始点として auto_wmake するだけです。

また、この際に打ち間違えのリコメンド機能を実装してみました。applicationフォルダの中にあるbinの名前はすべてコチラで持っているので、レーベンシュタイン距離 が最も近いものをリコメンドしてくれます。 結構かんたんなアルゴリズムではありますが、かなり良い精度でリコメンドしてくれます。

終わりに

今回始めてrustでリポジトリを作りましたが、インストールがかんたんで、所有権の移動等も慣れればとりあえずコンパイルは通ってくれるようになってくれました(美しいかどうかはさておき…) とりあえず困ったらC++でしたが少しずつrustに置き換えていく予定です。

もし、このソースコードのだめなところを見つけたら教えていただけると幸いです。 issueプルリク等待っています。よろしくお願いします。

ARC081 E - Don't Be a Subsequence

問題

https://beta.atcoder.jp/contests/arc081/tasks/arc081_c

問題概要

  • 英小文字のみからなる文字列が与えられる
  • 文字列の部分文字列"でない"文字列のうち、

  • 最小の長さのもの

  • 最小の長さのものが複数ある場合は、辞書順最小のもの

を出力する。

考察

  • 英小文字のみからなるので、結構なゴリ押しはできる
  • 答えは必ず、0文字を含む部分文字列+一文字になる。
  • 答えの部分文字列の最後の文字は、それ以降でa-zのうちどれかが出現していないような文字
  • 答えの部分文字列の最後の文字に付け足す文字は、それ以降で出現していないアルファベットのうち、最も小さいもの
  • それぞれの文字列の位置を始点とし、そこから文字X(a-zに関してすべて行う)が初めて出てくる場所に向けて辺を貼る。
  • 最初にa-zが現れる場所を最初の探索点として、幅優先探索をし、「最後の文字」が現れたら探索終了。「追加する一文字」を追加してそれを出力する。
  • 探索は必ずa-zの順番でqueueに入れていき、すでに訪れている位置に関しては無視して良い。(なぜならば、aから順番に探索しているので後から訪れるような状態は必ず先に訪れたものよりも辞書順が後になるから)

オーダーは文字種をMとすると、O(AM)

実装

void solve() {
    string S; cin >> S;
    vector<int> A(S.size());
    for (int i = 0; i < S.size(); ++i) {
        A[i] = S[i] - 'a';
    }
    vector<bool> check('z' - 'a'+1);
    vector<vector<int> > G(A.size());
    vector<int> memo(check.size(), -1);
 
    RREP(i, A.size()) {
        REP(j, check.size()) {
            if (memo[j] == -1) {
                G[i].push_back(-j-1);
                break;
            }
        }
        if (G[i].size() == 0) {
            G[i].resize(memo.size());
            copy(ALL(memo), G[i].begin());
        }
        memo[A[i]] = i;
    }
 
    vector<bool> visited(A.size());
    queue<int> que;
    for (int i = 0; i < memo.size(); ++i) {
        if (memo[i] == -1) {
            cout << char(i + 'a') << endl;
            return;
        }
        que.push(memo[i]);
        visited[memo[i]] = true;
    }
 
    vector<int> rev(A.size(), -1);
    int last = -1;
    string res;
    while (not que.empty()) {
        int node = que.front(); que.pop();
        if (G[node].back() < 0) {
            res.push_back(-G[node].back()-1 + 'a');
            last = node;
            break;
        }
        for (auto &to : G[node]) if(not visited[to]) {
            que.push(to);
            visited[to] = true;
            rev[to] = node;
        }
    }
    while (true) {
        if (last == -1) break;
        res.push_back(S[last]);
        last = rev[last];
    }
    reverse(ALL(res));
    cout << res << endl;

ARC096 E - Everything on It

問題

E - Everything on It

問題概要 & note

  • トッピングが N 種類ある
  • 一つのラーメンにはそれぞれのトッピングを1つか0つ乗せる
  • つまり全部で 2^{N} 個ラーメンができる
  • そのラーメンをいくつか選ぶ組み合わせは全部で 2^{2^{N}}
  • そのようなラーメンを幾つか選んで、それぞれのトッピングが2個以上ラーメンに乗っているようにしたい。
  • 全部で何通りのラーメンの組み合わせがあるでしょうか

考察(解説見た)

全部で  2^{2^{N}} 通りあって、そこから 一つ以上ののトッピングが1個以下 のようなものを除外する方針にする。どのように除外していくうかと言うと、たとえば、トッピングが2種類の場合は

  1. 「トッピングAが1個以下でほかは気にしない(1個以下でも2個以上でもいい)ラーメンの組み合わせ」除外する。
  2. 「トッピングBが1個以下でほかは気にしないラーメンの組み合わせ」を除外する
  3. 「トッピングAとトッピングBの両方が1個以下のラーメンの組み合わせ」を 追加する (なぜなら、Aを除外してBを除外すると、AB両方共通のものは2回減算されることになるから1回たさなければならない)

じゃあ3つの時はというとそれは3つめはまた引き算になる。一般化すると包除原理になる。

ここでよく考えてみると、「トッピングAが1個以下でほかは気にしないラーメンの組み合わせ」も「トッピングBが1個以下でほかは気にしないラーメンの組み合わせ」も数が一緒である「トッピングAB(ry」も「トッピングAC(ry」も数が一緒になるはずである。なのでより一般化すると「トッピング1~jまではi個以下で、ほかは気にしない組み合わせ」 * nCi とすれば、いちいち全部列挙しなくてもよくなる。

さて、それでは「トッピング1~iまでが一個以下ほかはどうでもいい」ラーメンの組み合わせの数は幾つあるのか? ここでways2(i,j)を, iまでのトッピングをjこのグループに分ける組み合わせの数(但し、グループに選ばないものがあっても良い)とする。これは、iまでのトッピングを、j個のラーメンに、2つ以上乗らないようにする組み合わせ数と同じで、\sum_{j=0}^N ways2(i,j) とすることで、iまでのトッピングが一つ以下しか乗っていないものの個数がO(N)で求められる。

ways2は第二種スターリング数に似たdpで解くことができる。第二種スターリングのままでは、選ばないトッピングがものを数え上げてしまうので、選ばないものを許すために、二項目の\tex[k]を\tex[(k+1)]とすることで、うまくいく(どこにも所属しないという新たなグループを作るイメージ)また、第二種スターリング数と違い、0個の選択を許すので、\tex[k=0]の時は1になることに注意する。

あとは、図を書いて説明しないといけないので、こちらの解説放送を見ていただきたい AtCoder Regular Contest 096/ Beginner Contest 095 解説放送 - YouTube

source code

// binomial coefficients O(n^2)
// res[i][j] = C(i,j) (i , j <= n)
// C(i,j) =  C(i-1, j-1) + C(i-1, j)
vector<vector<int> > BinomialCoefficients(int n) {
    vector<vector<int > > res;
    res.resize(max(n+1,2LL));
    res[0].push_back(1);
    res[1].push_back(1);
    res[1].push_back(1);
 
    FOR(i, 2, res.size()){
        res[i].push_back(1);
        FOR(j,1,i){
            res[i].push_back((res[i - 1][j - 1] + res[i - 1][j])%MOD);
        }
        res[i].push_back(1);
    }
    return res;
}
 
// Stirling NUmber of The Second Kind O(n^2)
// note: NOT symmetry!!!
// S(n, k) = S(n, k-1) + k*S(n-1, k)
vector<vector<int> > StirlingNumber2(int n) {
    vector<vector<int > > res;
    res.resize(max(n+1,2LL));
    res[0].push_back(1);
    res[1].push_back(1);
    res[1].push_back(1);
 
    FOR(i, 2, res.size()) {
        res[i].push_back(1);
        FOR(j,1,i){
            res[i].push_back((res[i-1][j - 1] + (j+1)*res[i-1][j]) % MOD);
        }
        res[i].push_back(1);
    }
    return res;
}
 
void solve() {
    int N; cin >> N >> MOD;
    vector<vector<int> > b = BinomialCoefficients(N);
    vector<vector<int> > s2 = StirlingNumber2(N);
    vector<int> pow2(N*N+1);
    vector<int> powpow2(N+1);
    pow2[0] = 1;
    powpow2[0] = 2;
    FOR(i, 1, pow2.size()) pow2[i] = (pow2[i - 1] * 2)%MOD;
    FOR(i, 1, powpow2.size()) powpow2[i] = (powpow2[i - 1] * powpow2[i - 1]) % MOD;
    vector<int> ways(N + 1);
 
    int res = 0;
    REP(i, ways.size()) {
        int w = 0;
        REP(j, i + 1) {
            int temp = (s2[i][j] * pow2[(N - i)*j]) % MOD;
            temp = (temp * powpow2[N - i]) % MOD;
            w = (w + temp) % MOD;
        }
        int temp = (w * b[N][i]) % MOD;
        temp *= (i % 2 == 0 ? 1 : -1);
        res = (res + temp);
        while (res < 0) res += MOD;
        res %= MOD;
    }
    cout << res << endl;
}

note

a^{b} (\mathrm{mod}\  p) = a^{b\ \mathrm{mod}\ (p-1)} (\mathrm{mod}\ p)

フェルマーの小定理より、a^{p-1} = 1 つまり、 a^{p-1} を何度かけても1を掛け算するので答えは変わらない  a^{b} = a^{b/p-1} * a^{b \% (p-1)} = a^{p-1} * a^{p-1} * \cdots * a^{p-1} *  a^{b \% (p-1)} = a^{b \% (p-1)}

AGC021 D - Reversed LCS

問題概要

  • 文字列 Sが与えられる
  •  S' Sを逆から読んだものである。
  •  Sのうち K文字を変えることが出来る
  • 最大 K文字を変更し,  S S'最長共通部分列を求める

考察

  • これってつまり回文を求めろってことでは…?
  • 回分なら典型的な区間DPがある(あとで解説)
  •  Kも一緒にDPテーブルに突っ込んでやる

 \mathrm{dp}(l)(r)(x) x回文字を変更し, [l, r)の間で最長の回文

漸化式

文字を変更しない場合

通常の回文であれば、以下のDPで更新できる。

l, r-1の文字が等しい場合(A???????A みたいな文字列の時は、最初と最後の文字を除いた真ん中の部分+2してやれば良い.)
 \mathrm{dp}(l)(r) = \mathrm{dp}(l+1)(r-1)+ 2
そうではない場合(A???????B みたいな文字列の時は、Aを含みBを含まない文字列とAを含まずにBを含む文字列で比較してやれば良い)
 \mathrm{dp}(l)(r) = \max(dp(l+1)(r), dp(l)(r-1))

文字を変更する場合

今回は、文字の変更ができるので、もう少しアレンジする。

l, r-1の文字が等しい場合(A???????A みたいな文字列の時は、最初と最後の文字を除いた真ん中の部分+2してやれば良い.)
 \mathrm{dp}(l)(r)(x) = \mathrm{dp}(l+1)(r-1)(x) + 2
そうではない場合(A???????B みたいな文字列の時は、Aを含みBを含まない文字列とAを含まずにBを含む文字列で比較してやれば良いし、さらに、xを一回分消費して、最初と最後の文字を強引に合わせても良い。)
 \mathrm{dp}(l)(r)(x) = \max(dp(l+1)(r)(x), dp(l)(r-1)(x), dp(l+1)(r-1)(x-1)+2)

実装

agc021.contest.atcoder.jp

#include <iostream>
#include <queue>
#include <map>
#include <list>
#include <vector>
#include <string>
#include <stack>
#include <limits>
#include <climits>
#include <cassert>
#include <fstream>
#include <cstring>
#include <cmath>
#include <bitset>
#include <iomanip>
#include <algorithm>
#include <functional>
#include <cstdio>
#include <ciso646>
#include <set>
#include <array>
#include <unordered_map>
#include <complex>

using namespace std;

#define FOR(i,a,b) for (int i=(a);i<(b);i++)
#define RFOR(i,a,b) for (int i=(b)-1;i>=(a);i--)
#define REP(i,n) for (int i=0;i<(n);i++)
#define RREP(i,n) for (int i=(n)-1;i>=0;i--)

#define inf 0x3f3f3f3f
#define PB push_back
#define MP make_pair
#define ALL(a) (a).begin(),(a).end()
#define SET(a,c) memset(a,c,sizeof a)
#define CLR(a) memset(a,0,sizeof a)
#define VS vector<string>
#define VI vector<ll>
#define DEBUG(x) cout<<#x<<": "<<x<<endl
#define MIN(a,b) (a>b?b:a)
#define MAX(a,b) (a>b?a:b)
#define pi 2*acos(0.0)
#define INFILE() freopen("in0.txt","r",stdin)
#define OUTFILE()freopen("out0.txt","w",stdout)
#define ll long long
#define ull unsigned long long
#define pii pair<ll,ll>
#define pcc pair<char,char>
#define pic pair<ll,char>
#define pci pair<char,ll>
#define eps 1e-14
#define FST first
#define SEC second
#define SETUP cin.tie(0), ios::sync_with_stdio(false), cout << setprecision(15)

namespace {
	struct input_returnner {
		ll N; input_returnner(ll N_ = 0) :N(N_) {}
		template<typename T> operator vector<T>() const { vector<T> res(N); for (auto &a : res) cin >> a; return std::move(res); }
		template<typename T> operator T() const { T res; cin >> res; return res; }
		template<typename T> T operator - (T right) { return T(input_returnner()) - right; }
		template<typename T> T operator + (T right) { return T(input_returnner()) + right; }
		template<typename T> T operator * (T right) { return T(input_returnner()) * right; }
		template<typename T> T operator / (T right) { return T(input_returnner()) / right; }
		template<typename T> T operator << (T right) { return T(input_returnner()) << right; }
		template<typename T> T operator >> (T right) { return T(input_returnner()) >> right; }
	};
	template<typename T> input_returnner in() { return in<T>(); }
	input_returnner in() { return input_returnner(); }
	input_returnner in(ll N) { return std::move(input_returnner(N)); }
}

const ll MOD = 1e9 + 7;

void solve();

signed main() {
	SETUP;
	solve();
#ifdef _DEBUG
	system("pause");
#endif
	return 0;
}

int dp[302][302][301] = {}; // dp[l][r][x] : x個文字列を変更する権利を使用した時に, [l,r)間で最長となる回文

void solve() {
	string S; cin >> S;
	int K; cin >> K;
	REP(i, S.size()) REP(x,K+1){
		dp[i][i + 1][x] = 1;
	}
	FOR(w, 2, S.size()+1) {
		FOR(l, 0, S.size() + 1 - w) {
			FOR(x, 0, K+1) {
				int r = l + w;
				// cout << S[l] << " " << S[r-1] << endl;
				if (S[l] == S[r-1]) {
					dp[l][r][x] = dp[l + 1][r - 1][x] + 2;
				}
				else {
					dp[l][r][x] = max(dp[l + 1][r][x], dp[l][r - 1][x]);
					if(x-1 >= 0) dp[l][r][x] = max(dp[l][r][x], dp[l + 1][r - 1][x - 1] + 2);
				}
			}
		}
	}
	int res = 0;
	REP(i, K+1) {
		res = max(res, dp[0][S.size()][i]);
	}
	cout << res << endl;
}

感想

  •  xの更新方法をミスったdpをする時は左辺に+1とか描いちゃダメだ
  •  x-1>0 とか描いてしまった配列外参照対策に、あえて変形せずにそのまま描いて保証する書き方をしているが、こういうミスは良くない 反省
  • 今回始めて半開区間区間DPを書いたが、l+1, r-1なテーブルにアクセスする時に便利だった。

AGC012 B - Holes

問題概要

  • めちゃめちゃでかい Rを持つ円が与えられる。
  • その真ん中 \pm 10^6 の範囲内に穴を置く。
  • めちゃめちゃでかい円 R内に点を起き、その点は最も近くにある穴に落ちる。
  • 円の内部すべてに点を起き、各穴に落ちる面積の割合をそれぞれ求める。
  • 同じ座標に2つ置くようなケースはない
  • ちなみに与えられる座標は整数値(今回はこの条件は使わない ゴリ押し解法だと必要かも?)
続きを読む