SnackDown Elimination 2017 PLUSMUL
問題要約(相方担当)
- N個の数列が与えられる
- それぞれの数の間に"*" か "+"を入れた式を作る(
通りの式ができる)
個の式の総和 MOD
を出力せよ．
概要
なんとなくDPで解けそうだったので適当にやったら遷移できた
実験
を，数列を左から
個読み込んだ時の式の総和とする．
とする
1つの要素の数列が与えられ，その間に"*"か"+"を挿入すると，
またこの総和をとする．
2つの要素の数列が与えられ，その間に"*"か"+"を挿入すると，
またこの総和をとする．
3つの要素の数列が与えられ，その間に"*"か"+"を挿入すると，
またこの総和をとする．
0->1の遷移に注目すると
1->2の遷移に注目すると
2->3の遷移に注目すると
との部分とそれ以外の部分(
)に分けることができた．
それ以外の部分もDPで求めることが出来る 具体的な証明や導出は省略するが，今回の場合は
と新たに増えた要素を使って簡単に表記できる．に係数2が付いているが，これは注目している数列の個数を
とすると，
となる．
大体法則性がわかったので一般化する．N個の数列をではなく,
と表記する．
に関するDP遷移は以下のようになる．
では
に関するDP遷移は以下のようになる．
では
ここで，に
かかってしまうので，計算しながらメモすることによって
にする．
総計算オーダーはとなり，これで間に合う．
ソースコード
実は実装ミスって上のやつと添字がちょっと違う
void solve() { int N; cin >> N; vector<int> v; REP(i, N) v.push_back(in()); vector<int> twoPow(N); vector<int> m(N); twoPow[0] = 1; FOR(i, 1, N) twoPow[i] = twoPow[i - 1] * 2 % MOD; vector<int> dp(N), dpsum(N); dp[0] = v[0]; dpsum[0] = v[0]; m[0] = v[0]; FOR(i,0,N-1){ m[i + 1] = ((m[i] * v[i + 1])%MOD + (twoPow[i] * v[i + 1])%MOD)%MOD; dp[i + 1] = (dpsum[i] + m[i+1])%MOD; dpsum[i + 1] = (dpsum[i] + dp[i + 1])%MOD; } cout << dp[N - 1] << endl; }
いろんな言語でSI単位系の計算をやってみよう(Cpp, python, rust, Julia)(WIP)
C++
C++なのでコンパイル時になんやかんやできるやつを探してたらBoostであった
Boost.Units を使えば良い．(http://www.boost.org/doc/libs/1_64_0/doc/html/boost_units.html)
sample
使用例として，簡単な問題を解いてみる．
h=634m地点から初速v0=2.4m/sでm=5kgのボールを落とす． ボールにかかる力は重力のみとすると，何秒で地面に到達するか．
#include <iostream> #include <boost/units/systems/si.hpp> #include <boost/units/io.hpp> #include <boost/units/pow.hpp> #include <boost/units/cmath.hpp> #include <boost/mpl/arithmetic.hpp> using namespace boost::units; // create new unit!! typedef boost::mpl::times<length_dimension,mass_dimension>::type LM_type; typedef unit<LM_type, si::system> kurenaif; int main() { const quantity<si::length> height = 634.0 * (si::meter); //x0 = 634m const quantity<si::mass> m = 5 * (si::kilogramme); const quantity<si::velocity> v0 = 2.4 * (si::meter / si::second); // const quantity<si::force> F = m*v0; compile error!!! because of invalid unit type const quantity<si::acceleration> g = 9.8 * (si::meter/pow<2>(si::second)); //create new unit!!! const quantity<kurenaif> test = height*m; const quantity<si::force> F = m*g; const quantity<si::time> t = -(v0 - sqrt(v0*v0+2.0*g*height))/g; std::cout << "height(x0) = " << height << std::endl; std::cout << "weight = " << m << std::endl; std::cout << "init. velocity = " << v0 << std::endl; std::cout << "g: " << g << std::endl; std::cout << "kurenaif:" << test << std::endl; std::cout << "force:" << F << std::endl; std::cout << "time:" << t << std::endl; return 0; } /* output height(x0) = 634 m weight = 5 kg init. velocity = 2.4 m s^-1 g: 9.8 m s^-2 kurenaif:3170 m kg force:49 m kg s^-2 time:11.1326 s */
memo
- 次元間違えはコンパイル時エラー
- auto型使えます
Juliaで競技プログラミング(もうしない, ABC061)
Motivation
研究でちょっと触っているのでABCで使ってみた
解説は色々みたらわかるので割愛
A問題
http://abc061.contest.atcoder.jp/tasks/abc061_a
A,B,C = parse.(split(readline(STDIN))) if A <= C && C <= B println("Yes") else println("No") end
http://abc061.contest.atcoder.jp/submissions/1283186
B問題
http://abc061.contest.atcoder.jp/tasks/abc061_b
N,M = parse.(split(readline(STDIN))) # ブロードキャスト 参考: http://bicycle1885.hatenablog.com/entry/2016/12/13/205646
count = zeros(Int64, N)
for i in 1:M
a,b = parse.(split(readline(STDIN)))
count[a] += 1
count[b] += 1
end
for value in count
println("$value")
end
http://abc061.contest.atcoder.jp/submissions/1284568
C問題 (TLE)
mapのかわりにdictを使ってみたがTLEしてしまった．
http://abc061.contest.atcoder.jp/tasks/abc061_c
N,M = split(readline(STDIN))
N = parse(N)
M = parse(M)
dict = Dict{Int64,Int64}()
for i in 1:N
a,b = split(readline(STDIN))
a = parse(a)
b = parse(b)
dict[a] = get(dict, a, 0)
dict[a] += b
end
for key in sort(collect(keys(dict)))
M -= dict[key]
if M <= 0
println("$key")
break
end
end
http://abc061.contest.atcoder.jp/submissions/1285715
所感
A問題すらすごい時間がかかって謎だった C問題は解法はO(NlogN)だしWAだとしてもTLEはしないはずなんだけど… 原因が分かるまで競技プログラミングでjuliaは封印します．
rustをインストールから行列で連立方程式を解くところまで
Motivation
なぜrust?
- 最近C++がつらい
- rustが流行っているらしい
- https://imoz.jp/note/rust-functions.html
環境
- ubuntu 16.04 64bit
rustのインストール
このコマンドで簡単にインストールする
curl https://sh.rustup.rs -sSf | sh
以下のようなメッセージができる． 初心者なのでdefaultのstableをinsatllする．
default host triple: x86_64-unknown-linux-gnu
default toolchain: stable
modify PATH variable: yes
1) Proceed with installation (default)
2) Customize installation
3) Cancel installation
インストールが終わると以下のメッセージが出てくるので
To configure your current shell run source $HOME/.cargo/env
.bashrcかなんかに
source $HOME/.cargo/env
を追加する
rustc --version
を入力して，それらしいものがでてきたらとりあえずOK．
References
crateの作成
Rust では create (木箱) 単位でプログラムを開発するらしい
今回は簡単な行列計算をやってみるので，mat_calcという名前でcrateを作る
cargo new mat_calc --bin
なんかできた
output: Created binary (application) `mat_calc` project
ディレクトリを覗いてみると
$ ls mat_calc . .. Cargo.toml .git .gitignore src
どうやら.gitも生成されているらしいぱっと見た感じ.gitignoreもいい感じに定義されてた
実行テスト
src/main.rsにはすでにHello Worldを出力するプログラムが入っているのでとりあえず実行してみる
とりあえずディレクトリに移動して
$ cd mat_calc $ cargo run Compiling mat_calc v0.1.0 (file:///home/kurenaif/Desktop/mat_calc) Finished dev [unoptimized + debuginfo] target(s) in 0.53 secs Running `target/debug/mat_calc` Hello, world!
無事出力に成功した
References
dependenciesの設定
行列を扱いたいけどFortranでもないしrustにはそんなライブラリは標準搭載されていない. rustで行列が扱えるライブラリを探したら以下の3つが見つかった
- https://github.com/sebcrozet/nalgebra
- https://github.com/termoshtt/ndarray-linalg
- https://github.com/masonium/linxal
一番使いやすそうで，commitが盛んなnalgebraを使うことにした．
速度面はよくわからなかった．
ここで, Cargo.tomlを編集し，nalgebraを使用できるようにする(pip installに近いイメージ)
[dependencies] nalgebra = "0.11.0"
versionは2017-05-05の時点での最新バージョンを使用した.
githubなりなんなりをみてうまいこと指定する．
githubのURLでの指定も可能だけど https://crates.io に登録されているのでこんな感じでURLなしでもcrateを入れることが出来る.
References
連立方程式を解いてみる
行列の表示
以下のソースコードをsrc/main.rsに実装し，実行してみる
今後の流体計算を考慮して普通のMatrixよりも実行時に長さを決められるDMatrixのほうが都合がよさそうなのでDMatrixを使用する.
使い方その他もろもろはソースコードみて察していただきたい．
個々で気をつけなければならないのが行列の数字はすべて実数で表記しないとちゃんとprintlnしてくれない.
つまり 1,2,3… のような整数を書く場合でも, 1., 2., 3. と書いたほうが後々楽になる．
rustは整数と小数，また同じ整数どうしや小数どうしてもbit幅が違えばコンパイルエラーを出すほど厳密．
extern crate nalgebra as na; fn main(){ let mat = na::DMatrix::from_iterator(3,3,[ 1.,2.,3., 4.,5.,6., 7.,8.,9. ].iter().cloned()); println!("{}",mat); let b = na::DVector::from_iterator(3,[ 1.,2.,3. ].iter().cloned()); println!("{}",b); }
$ cargo run ┌ ┐ │ 1.000 4.000 7.000 │ │ 2.000 5.000 8.000 │ │ 3.000 6.000 9.000 │ └ ┘ ┌ ┐ │ 1.000 │ │ 2.000 │ │ 3.000 │ └ ┘
- めっちゃいい感じの出力をしてくれる.(DMatrixじゃなくて普通のMatrixはこうはいかない)
- どうやらソースコードに書いたものに比べて転置したものが出てくるらしい 縦ベクトルが3つ横に並んでるイメージ.
逆行列
方程式を解くためには逆行列が必要不可欠 さっき用意したmatは実は逆行列が存在しない．
面白そうなのでさっきのmatの逆行列を計算してみる．
nalgebraでは，.try_inverse()を用いて逆行列を計算する.
無理だったら無理なりの挙動をしてくれる.
try_inverseを出力するプログラム．
println!("{}") と println!("{:?}")では呼びされるメソッドが微妙に異なっており，後者の方がちょっと詳しいメッセージが出てくることが多い(デバッグ用出力?)．
今回，前者を使うとエラーを吐くので後者を使っている．
extern crate nalgebra as na; fn main(){ let mat = na::DMatrix::from_iterator(3,3,[ 1.,2.,3., 4.,5.,6., 7.,8.,9. ].iter().cloned()); println!("{:?}",mat.try_inverse()); }
$ cargo run None
Noneが出力された．
リファレンスに書いているが，try_inverse()をしてなかった時はNoneが返ってくるらしい
では逆行列が存在する時はどんな出力をするのだろうか?
今回はtry_inverseの結果を"{}"と"{:?}"の二つを出力している．
最初のいもすさんの記事でも書いてある通り，rustでは基本所有権がmoveしちゃうので最初のprintln!でmatの所有権が切れる．
そこで，clone()することで所有権を継続させている．
extern crate nalgebra as na; fn main(){ let mat = na::DMatrix::from_iterator(3,3,[ 1.,3.,5., 2.,3.,8., 2.,1.,5. ].iter().cloned()); println!("{:?}",mat.clone().try_inverse()); println!("{}",mat.try_inverse().unwrap()); }
Some(Matrix { data: MatrixVec { data: [1.4, -2, 1.8, 1.2, -1, 0.4, -0.8, 1, -0.6], nrows: Dynamic { value: 3 }, ncols: Dynamic { value: 3 } }, _phantoms: PhantomData })
┌ ┐
│ 1.400 1.200 -0.800 │
│ -2.000 -1.000 1.000 │
│ 1.800 0.400 -0.600 │
└ ┘
一つ目の出力がtry_inverse()の中身を詳細に出力したものである．
Some()がついていて，実はこのままでは{}の出力ができない．
ここでSome()は値がNullになるかもしれない的なやつらしい．
Some()を取るためには，unwrap()をすると良い．
詳しくはここを参照する:
これの結果があってるかわからないので元のmatにかけて単位行列が出ること確認する
rustでは演算子のオーバーロードができる．
このライブラリでは行列の掛け算もちゃんと定義されているので，以下のソースコードみたいな感じで表記できる.
extern crate nalgebra as na; fn main(){ let mat = na::DMatrix::from_iterator(3,3,[ 1.,3.,5., 2.,3.,8., 2.,1.,5. ].iter().cloned()); println!("mat\n{}",mat.clone()); println!("inv(mat)\n{}",mat.clone().try_inverse().unwrap()); println!("mat*inv(mat)\n{}",mat.clone().try_inverse().unwrap()*mat); }
$ cargo run mat ┌ ┐ │ 1.000 2.000 2.000 │ │ 3.000 3.000 1.000 │ │ 5.000 8.000 5.000 │ └ ┘ inv(mat) ┌ ┐ │ 1.400 1.200 -0.800 │ │ -2.000 -1.000 1.000 │ │ 1.800 0.400 -0.600 │ └ ┘ mat*inv(mat) ┌ ┐ │ 1.000 -0.000 0.000 │ │ 0.000 1.000 0.000 │ │ 0.000 0.000 1.000 │ └ ┘
ちゃんと逆行列が計算されているっぽい
連立方程式を解いてみる
いよいよ本題
とりあつかう問題はこれの例1: 連立方程式の解き方
連立方程式を $$ x = A^{-1} b$$ の形にしてinvして解くだけ
extern crate nalgebra as na; fn main(){ let a = na::DMatrix::from_iterator(3,3,[ 1.,1.,0., 2.,-1.,1., 3.,2.,1. ].iter().cloned()); println!("A\n{}",a); let b = na::DVector::from_iterator(3,[ 0.,3.,-1. ].iter().cloned()); println!("B\n{}",b); let invA = a.try_inverse().unwrap(); println!("invA\n{}",invA); println!("x\n{}", invA*b) }
A ┌ ┐ │ 1.000 2.000 3.000 │ │ 1.000 -1.000 2.000 │ │ 0.000 1.000 1.000 │ └ ┘ B ┌ ┐ │ 0.000 │ │ 3.000 │ │ -1.000 │ └ ┘ invA ┌ ┐ │ 1.500 -0.500 -3.500 │ │ 0.500 -0.500 -0.500 │ │ -0.500 0.500 1.500 │ └ ┘ x ┌ ┐ │ 2.000 │ │ -1.000 │ │ 0.000 │ └ ┘
ページの答えと一致してるしちゃんと解けてるっぽい．
おまけ
今回色々やったリポジトリをここにおいておきます github.com
次回は簡単な熱拡散問題とか解いてみるかも
OpenFOAM-2.4.0をmacのgccでmakeする
Motivation
Macはclangをつかってコンパイルをしている。 どうやらコンパイラによって微妙に挙動が違うらしくMacでgccが使いたいと先生に言われて色々やった。
OpenFOAM-2.4.0向け
概要
流れ的には
といった順番になる。
一つづつやっていこう
gccを入れる
自分の環境と合わせたかったので最新のものではなくgcc-5.4.0をぶっこんだ
ソースコードからbuildする
まずはソースコードを入手し、解凍
$ wget http://ftp.tsukuba.wide.ad.jp/software/gcc/releases/gcc-5.4.0/gcc-5.4.0.tar.gz $ tar zxvf gcc-5.4.0.tar.gz
解凍したディレクトリに移動して依存ファイル取得
$ cd gcc-5.4.0 $ contrib/download_prerequisites
buildディレクトリを作成し、configure, build
$ mkdir build $ cd build $ ../configure --enable-languages=c,c++,fortran --prefix=/usr/local/gcc-5.4.0 --disable-bootstrap --disableultilib
$ make -j4 $ make install
/usr/local/binにgccをln -s
$ ln -s /usr/local/gcc-5.4.0/bin/gcc /usr/local/bin/gcc $ ln -s /usr/local/gcc-5.4.0/bin/g++ /usr/local/bin/g++
以下をbashrcに追記
export PATH=/usr/local/bin:$PATH
確認
$ gcc --version gcc (GCC) 5.4.0 Copyright (C) 2015 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE
参考: http://qiita.com/liveralmask/items/6ed4a98ebb3bf6b7f707
openmpiを入れる
新しすぎるopenmpiはなんかエラー吐いてめんどくさいのでちょっと古いの使う
$ wget https://www.open-mpi.org/software/ompi/v1.10/downloads/openmpi-1.10.2.tar.gz $ tar zxvf openmpi-1.10.2.tar.gz
参考: https://www.open-mpi.org/software/ompi/v1.10/
buildディレクトリ作ってconfigureしてmake
$ mkdir build $ ../configure --prefix=/usr/local/openmpi $ make $ make install
以下を.bashrcに追加
export MANPATH=/usr/local/openmpi/share/man:$MANPATH export LD_LIBRARY_PATH=/usr/local/openmpi/lib:$LD_LIBRARY_PATH export PATH=/usr/local/openmpi/bin:$PATH
チェック
$ mpicc --version
参考: http://qiita.com/himo/items/c30d83d0f7642fb3af57
flexを入れる
これはmake – installでもbrewでもなんか無理だった macportを使う
sudo port install flex
OpenFOAMを入れる
なんかいろいろする
$ hdiutil create -size 8.3g -type SPARSEBUNDLE -fs HFSX -volname OpenFOAM -fsargs -s OpenFOAM.sparsebundle $ mkdir -p OpenFOAM $ hdiutil attach -mountpoint $HOME/OpenFOAM OpenFOAM.sparsebundle $ cd OpenFOAM
OpenFOAM2.4.0をダウンロード、解凍する
$ wget http://downloads.sourceforge.net/foam/OpenFOAM-2.4.0.tgz?use_mirror=mesh $ tar zxvf OpenFOAM-2.4.0.tar.gz
参考: https://openfoam.org/download/2-4-0-source/ パッチを当てる
$ cd OpenFOAM-2.4.0 これをダウンロード -> https://sourceforge.net/p/openfoam-extend/svn/HEAD/tree/trunk/Breeder_2.4/distroPatches/MacPatch/OpenFOAM-2.4.x-Mac.patch $ patch -p1<OpenFOAM-2.4.x_Mac.patch
make
./Allwmake
エラーが出るのでsrc/OpenFOAM/lnIncludeにFlexLexer.hを無理やりリンクさせて通す
$ ln -s /opt/local/include/FlexLexer.h $FOAM_SRC/OpenFOAM/lnInclude
もう一度make
./Allwmake
チェック
which icoFoam
decomposeParでmanual切りをする
注意: バージョンによってファイル形式が異なる可能性があります.(OpenFOAM-2.4.0)
Motivation
需要は少ないですが，manualでdecomposeParをする必要がありました． その手法について資料がないので解説します
ケースファイル
tutorials/incompressible/icoFoam/cavity
メッシュ作成
blockMesh
decomposePar(自動領域分割)
tutorials/incompressible/pisoFoam/les/motorBike/motorBike/systemを少しいじったものです
以下を使用します
system以下に入れておくと良いです．
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.4.0 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object decomposeParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
numberOfSubdomains 4;
// numberOfSubdomains 3;
method simple;
simpleCoeffs
{
n ( 2 2 1 );
delta 0.001;
}
manualCoeffs
{
dataFile "cellDecomposition";
}
// ************************************************************************* //
領域分割を実行します cellDistを指定すると後でどういうふうに領域分割したのかを見れます．
decomposePar -cellDist
OpenFOAMで領域分割結果を確認します．(cellDistにチェックを入れて見ると良いです)
使用するソースコード
ここに書いてるソースコードをwmakeしてcavityのディレクトリで実行したら動きます．
https://github.com/kurenaif/manualDecomposeTutorial
これを実行すると，constant/cellDecompositionが生まれます．
その中にどのセルがどのprocessorに該当するのかが一セルずつ記述されています．
また，manualでdecomposeをするためには，system/decomopseParを弄る必要があります．
具体的には，
- 今回は3分割のため，numberOfSubdomainsを3に
- methodをmanual
に変更します．変更したら以下のようになります
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.4.0 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object decomposeParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
numberOfSubdomains 3;
method manual;
simpleCoeffs
{
n ( 2 2 1 );
delta 0.001;
}
manualCoeffs
{
dataFile "cellDecomposition";
}
// ************************************************************************* //
ソースコードの解説
argsの準備
コマンドライン引数をOpenFOAMのものに変換する．
Foam::argList args(argc,argv); if(not args.checkRootCase()) Foam::FatalError.exit();
runTimeの準備
meshを用意するために必要． argsはここで必要
Foam::Time runTime(Foam::Time::controlDictName, args);
meshの準備
Foam::fvMesh mesh
(
Foam::IOobject
(
Foam::fvMesh::defaultRegion,
runTime.timeName(),
runTime,
Foam::IOobject::MUST_READ
)
);
各セルがどのprocessorに該当するのかを記述するための変数の用意
OpenFOAMのListはconstructorで渡してやるとそのサイズを用意してくれる．
mesh.cells().size()はメッシュのセルのサイズを返してくれるので，この容量があればそれぞれのセルがどのprocessorに該当するのかを収容できる．
labelList procIds(mesh.cells().size());
各セルをループで回しそれぞれのセルにprocessorIDを振り分ける
回すだけ 斜めと円形に切ってます
forAll(mesh.cells(), cid){
const vector &C = mesh.C()[cid];
if(C[0]*C[0]+C[1]*C[1] < 0.01){
procIds[cid] = 1;
}
else{
procIds[cid] = 0;
}
if(C[0] + C[1] < 0.1){
++procIds[cid];
}
}
出力
IOobjectを使います． ちょっと詳細はよくわかりませんがprocIdsを渡してやるとどうやら出力してくれるらしいです．
labelIOList cellDecomposition
(
IOobject
(
"cellDecomposition",
mesh.facesInstance(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
procIds
);
cellDecomposition.write();
References
decomopseParのソースコード
