No Such Blog or Diary

ちゃんと資料を追わなかったので分からんのだけど，元の式の対称性を使って変数を減らしたときに正定値条件は必要十分でカバーできるけど目的関数の最適性って保存されるのだろうか？　変数を減らすプロセスに目的関数の係数行列ってでてきてたっけ？　気になってきたのでそのうち確認しよう．

もう少し頭使えよと言いたくなる結果．まあいいだろうとか考えて放っておいた部分が響いて正解なしに…

250点問題： price を最大までインクリメントしながらチェックしてたら時間オーバーだった．五千万位のループは大丈夫かと思ったけどだめみたい．馬鹿なことせずに与えられた配列内の値だけチェックすれば通っていたものを（誰かの max までは単調増加なので）．

500点問題：図形っぽいからパスした．直方体を２の冪のサイズの立方体で埋める問題．正解を見てわかったけど，２の冪しかないのでな 1x1 で埋めなければならない部分を埋めて再帰的に1/2サイズにした問題を解けばよかったらしい．

1000点問題：結局のところ行列式に (-1)^(N-1) を掛けたものが答えになる．理由は以下のとおり．まず，行列に吐き出し操作をしても答えは変わらない（自明）．そして，吐き出し後には，積が非ゼロになるドミノの選択は全てが対角要素である場合しかなく，その積は行列式に等しい．同時に，ドミノの連結グループ数は行列の大きさである N に等くなっているので，答えは行列式に(-1)^(N-1) を掛けたものに等しい（偶数なら符号反転）．ということで，行列式を求めればおわる．幸いなことに modulo 121547 の答えでよく 121547 が素数で有限体上の計算になるので普通のLU分解とかが使える．

以上を踏まえ，有限体の逆元を求めるための拡張ユークリッドの互除法を実装し，演算を有限体の演算に置換したピボットつき LU 分解を実装して提出したのだが…　正規化のためのルーチンを次のように書いていたため TLE 喰らった．

    int norm(long a) {
        while(a < 0) a+=P;
        return (int)(a % P);
    }

掛け算してんだから数万回ループすることもあるよなぁ．もう少し進化して下のように書いたら普通に通ったし．

    int norm(long a) {
        if(a<0) a -= (a/P-1)*P;
        return (int)(a % P);
    }

あー，もったいねぇ．解けてたのに…

別に痛くないけど気になる部分があるのでいってみた．とりあえず気になった部分は虫歯＋前の詰め物が少しとがっていたと判明．そして余計なものもいくつか見つかった．一番面倒そうなのが左上親知らずの外側の生え際の虫歯．なんとなく前からなんか違和感あると思っていたら虫歯だったと．とりあえずどうなるかわからんけど全部片付けるのに一ヶ月～一ヵ月半くらいらしい．

EuroPar2007で発表のあった The Multi-Core Standard Template Library (MCSTL) が GCC に統合されている（されつつある？）らしいので，セミナーのネタにと試用してみた．GCCのオフィシャルページを見ると libstdc++ parallel mode という名前で統合中らしい．こいつは既存の STL を OpenMP による並列実装の STL に置換しましょうというもので，既存の STL で書かれたコードがそのままで並列プログラムになってくれる．

んで，試そうと思ったら最新のGCCリリースであるところの 4.2.2 にも入っていないということに GCC4.2.2 のビルドが終わった時点で気づき…　しょうがないので svn のリポジトリから最新のソースを落としてきて GCC をビルドした（８コアくらいあるマシンで make -j8 とかやるとビルドがとても速い）．とりあえずソースの状態で libstdc++-v3/include/parallel というディレクトリがあれば parallel mode が使える．

準備できたので次の std::sort を使ったプログラムを試してみた．

#include <iostream>
#include <algorithm>
#include <vector>
const int count = 1000000, num_tests = 10;
int main(int argc, char **argv)
{
    std::vector<int> v(count);
    for (int i = 0; i < num_tests; i++)
    {
        std::generate(v.begin(), v.end(), rand);
        std::sort(v.begin(), v.end());
        std::cout << "." << std::flush;
    }
}

コンパイルは parallel mode 用のフラグ（-D_GLIBCXX_PARALLEL -fopenmp）を立てればよいらしい．

 g++-XXX test.cpp -o testS -march=native  # sequential mode
 g++-XXX test.cpp -o testP -march=native -D_GLIBCXX_PARALLEL -fopenmp  # parallel mode

Quad の Xeon x2 という８コアで動かした結果は以下のとおり．

> OMP_NUM_THREADS=8 time ./testS   # sequential mode
..........        6.26 real         6.24 user         0.01 sys
> OMP_NUM_THREADS=8 time ./testP   # parallel mode
..........        1.39 real         5.66 user         0.75 sys

普通に４～５倍速くなってるのを確認．すげー．並列プログラミング簡単．とりあえず sort は簡単に並列に動いてくれた．

ついでに， parallel mode は他にも並列実装を持っているらしいので試してみた．

まず std::partial_sum．4.51秒が1.65秒になった．８コアで約四倍．

#include <iostream>
#include <algorithm>
#include <vector>
#include <numeric>
struct one {
  int operator()() const {
    return 1;
  }
};
const int count = 10000000, num_tests = 10;
int main(int argc, char **argv)
{
    std::vector<int> v(count);
    std::vector<int> r(count);
    std::generate(v.begin(), v.end(), one());
    for (int i = 0; i < num_tests; i++)
    {
        std::partial_sum(v.begin(), v.end(), r.begin());
        std::cout << "." << std::flush;
    }
}

次は std::accumulate．3.20秒が3.79秒に増えた．全要素の和をとるだけなので普通に考えると８コアで８倍いきそうだけど…　速くならない理由は不明．

#include <iostream>
#include <algorithm>
#include <vector>
#include <numeric>
struct one {
  int operator()() const {
    return 1;
  }
};
const int count = 10000000, num_tests = 10;
int main(int argc, char **argv)
{
    std::vector<int> v(count);
    std::generate(v.begin(), v.end(), one());
    for (int i = 0; i < num_tests; i++)
    {
        int sum = std::accumulate(v.begin(), v.end(), 0);
        std::cout << "." << std::flush;
    }
}

最後に試したのは std::for_each．3.81秒が3.75秒．全要素に独立に関数適用するだけなので８コアで８倍いきそうだけど…　速くならない理由は不明．

#include <iostream>
#include <algorithm>
#include <vector>
#include <numeric>
struct f {
  void operator()(int &x) const {
    x = x + 1;
  }
};
struct one {
  int operator()() const {
    return 1;
  }
};
const int count = 10000000, num_tests = 10;
int main(int argc, char **argv)
{
    std::vector<int> v(count);
    std::generate(v.begin(), v.end(), one());
    for (int i = 0; i < num_tests; i++)
    {
        std::for_each(v.begin(), v.end(), f());
        std::cout << "." << std::flush;
    }
}

ということで，簡単に並列化できそうな部分でいろいろと失敗した模様．使い方が悪いのか実装途中なのか？　よくわからんけどうまく動いた部分に関しては効果が確認できたのでよし．

うまく使えずにごちゃごちゃやってたけど最終的に落ち着いた方法：

(setq tramp-default-method "sshx")

を .emacs に追加して ssh 経由をデフォルトにしつつ， cygwin の ssh-agent でキー認証をかけておいて meadow をそのコマンドラインから起動．

落ちとしては，HOME の設定が meadow 用と cygwin 用で別々になっているため .emacs がうまく読まれなかったり，同じく HOME のせいで known_hosts へのキーの追加に失敗したり，キー追加確認用の YES/NO 選択で応答不能になったりと．

とりあえず HOME の設定を整理しようかなと思う．

管理者権限での実行は runas コマンドで実行アカウントを切り換えればできるらしい．でも，Administrator アカウント自体が有効にされていなかったために何度試しても失敗しまくった．

ということで，Administrator アカウントの有効化：スタートメニューとかのコンピュータの右クリックメニューの管理（コンパネからコンピュータの管理でも可）→ システムツールの下のローカルユーザとグループ → ユーザ → Administrator を選択してプロパティ開いて無効化のチェックをはずして閉じる → Administrator を選択してパスワード設定．これでやりたい放題．