憂鬱な午後のひととき


11月29日_
VCでclock()はCLOCKS_PER_SEC = 1000なのでmsec単位だけど, 実際の分解能は十数msec程度(環境によってはmsecレベルのときもある). それは_ftime_s()やGetSystemTimeAsFileTime()を使っても変わらない. それより細かい精度の時間を計りたいならrdtsc系のQueryPerformanceCounter()を使いましょうという話. ところで, rdtscの拡張でプロセッサごとの識別情報を同時に取得できるrdtscpという命令(AMD系はわりと昔からあった)がi7でも使えるのを初めて知った. VC9でも__rdtscpを持っていた.
11月24日_
昼休み, 同僚に「高速化のためには条件分岐を減らした方がいいのか」と聞かれて「確率が見積もれるかそうじゃないかでケースバイケース」と答えた. 自分の予想通りか実際に測ってみる. 環境はどれも32bit exeで(jmp-cmov.cpp).
case 1. 配列の最大値を求める

この場合, データがランダムだったとしても, その時点での最大値との比較になるので後半ほど最大値を更新しない確率が高い. つまり分岐予測があたりやすくなると思われる.

    /* cmov版のコア */
    cmp(eax, ptr [edx + ecx * 4]);
    cmovl(eax, ptr [edx + ecx * 4]);

    /* jmp版のコア */
    cmp(eax, ptr [edx + ecx * 4]);
    jge(".skip");
    mov(eax, ptr [edx + ecx * 4]);
L(".skip");
データの形を変えたときの最大値取得処理の時間の変化 [clk]
乱数 最初が最大値 単調増加 だいたい単調増加
Core2Duo + Xp jmp 2.04 2.01 2.02 4.02
cmov 5.10 5.10 5.11 5.10
Core i7 + Win7 jmp 1.77 1.84 1.85 3.54
cmov 2.68 2.78 2.66 2.58
Opteron 2378 + Win2008 jmp 5.06 5.05 3.03 4.05
cmov 2.02 2.02 2.02 2.02

cmovを使うとデータによらず実行時間は一定. jmpは分岐予測のあたり外れにより時間が変動する. 最初に最大値があったり(=殆ど分岐する), 単調増加(=殆ど分岐しない)だったりするとIntel系では速い.

Opteronは分岐すると常にコストがかかる模様. ということは分岐の確率が偏っている場合, このCPUでは確率が高い方を分岐しないように配置したほうがいい(昔のIntel CPUもそうだった). cmovはかなり速い. 最近のAMD系CPUだとどうだろう. 知りたいけど手近なマシンが無かった.

データが"だいたい単調増加"のときは, 最大値の更新がときどき起こり, しかも最後まであるので分岐予測が難しい. したがって分岐予測が外れる率も高くなり, jmp版のコストが高い.
case 2. 二つの配列を比較して左が大きいときの個数を数える

int countMax(const int *a, const int *b, size_t n)
{
    int ret = 0;
    for (size_t i = 0; i < n; i++) {
        if (a[i] > b[i]) ret++;
    }
    return ret;
}

0. 条件分岐を使う, 1. setg + movzxを使う, 2. carryを使ったややトリッキーなコード, 3. setgで余計にレジスタを使う, で試してみた.

    /* 0. 条件分岐 */
    mov(ebx, ptr [esi + ecx * 4]); // a
    cmp(ebx, ptr [edx + ecx * 4]); // b
    jle(".skip");
    add(eax, 1);
L(".skip");

    /* 1. setg + movzx */
    mov(ebx, ptr [esi + ecx * 4]);
    cmp(ebx, ptr [edx + ecx * 4]);
    setg(bl);
    movzx(ebx, bl);
    add(eax, ebx);

    /* 2. carryを使う */
    mov(ebx, ptr [esi + ecx * 4]);
    cmp(ebx, ptr [edx + ecx * 4]);
    sbb(eax, -1);

    /* 3. setg + レジスタ */
    mov(edi, ptr [esi + ecx * 4]);
    cmp(edi, ptr [edx + ecx * 4]);
    setg(bl);
    add(eax, ebx);
右が大きい個数の割合を変えたときの処理時間の変化 [clk]
割合[%] 0 25 50 75 100
Core2Duo + Xp jmp 2.28 9.65 12.82 9.44 2.27
setg + movzx 3.02 3.02 3.02 3.02 3.02
carry 2.40 2.41 2.41 2.40 2.40
setg 2.35 2.35 2.35 2.35 2.35
Core i7 + Win7 jmp 2.62 9.33 12.71 9.15 1.85
setg + movzx 2.61 2.89 2.65 2.75 2.60
carry 2.78 2.88 2.70 2.61 2.64
setg 2.64 2.65 2.72 2.72 2.65
Opteron 2378 + Win2008 jmp 3.04 7.91 10.90 8.28 4.04
setg + movzx 2.88 2.94 2.95 3.15 2.90
carry 3.03 3.03 3.02 3.03 3.03
setg 2.25 2.04 2.04 2.04 2.04

carryを使ったコードはsetgよりも常によいだろうと思ったが, Opteronではそうとも言い切れなかった. Opteronではcmovと同様setgは速い.

setg + movzxとsetgを比べると, Core2DuoとOpteronで気持ち速くなった. レジスタが余ってるならcarryよりいいかも.

jmpでは"割合"が大きくなるほど個数をincrementする部分のコードはスキップされる. 割合が50の時のjmpの遅さは注意すべき. 本当に外れるときはこのぐらいのコストがかかっている.
#1/2のならし平均なので恐らく20clk以上.

これを見るとcase1の最大値のときの"だいたい単調増加"では実は結構, 分岐予測があたっていたと推測される.

ところで0以上n未満の二つの自然数a, bの和がnを超える確率はほぼ1/2なので, (a + b) mod nの計算で分岐予測はできない. よって上記の"割合=50"に相当する. というわけでateペアリングの実装ではcmovを使っていた. で, そもそもmodulo演算をできるだけしないようにロジックを変えたのであった.

11月20日_
昨日の夜中ぐらいから, 胃のあたりが痛くて何度ももどす. 朝病院に行ったらえらく混んでて3時間近く待った. 水を飲んでも吐く. どうしようもなくしんどい. せっかく連休に実家に戻って遊ぼうと思ってたのに.
11月19日_
暗号理論研究の現場からを途中から聴講. まさに「企業の現場」のお話だった. 古川さんの選挙のための暗号や, グループ署名がどう「うまくいかなかったか」という話はとても興味深かった. 高島さんの企業としてのアウトプットをどうしていくかという話もよかった. ゆっくりお話を聞きたかったけどちょっとしんどかったので先に帰る. 残念.
11月17日_
昔書いてたboost::spirit(::classic)を使っていたプログラムをspirit::qiで書き直してみる. コンパイル速度が尋常じゃなく遅くなった. 300行足らずのコードに20秒近くかかる. classicも遅かったけどせいぜい数秒. ヘッダのプリコンパイルを使っても重たい. これじゃあtemplateネタプログラミング(さむ)としても使う気になれない.
11月16日_
amachangたちと夜ご飯.
11月12日_
『死の影の谷間』(ロバート・C. オブライエン 越智 道雄訳)
1976年に書かれた児童向け小説の再翻訳版. 核戦争後, "谷間"に残された少女が一人. そこへ別のところで生き残っていた科学者がやってきて共同生活が始まる. 穏やかな日々になるのかと思いきや, あるものが壊れていく様をたんたんと描写する. ラストは"敵"がそんな甘い台詞をいうかなあとも思って釈然としない気がしたけど, 全編を漂う恐怖感は凄いものがあった.
11月11日_
サーバ上の暗号化されたデータを利用するにはを聞きに行く. 当たり前だけど聴講者が社会人の方が多くて実際にどの程度使えるのか, すぐ使えるのか的な質問が多かった. まだまだ研究段階で難しい技術が多そうだ.
11月10日_
Project Euler 250
最初見たときは部分集合って2250250個だから無理だろっと思ったけど, やってみたら意外と簡単だった. C++で0.66秒.
const uint64_t n16 = 10000000000000000ULL;
const int N = 250250;
const int M = 250;
static int vec[N + 1];
static uint64_t num[M];

uint32_t ipowMod(uint32_t x, uint32_t a, uint32_t m)
{
    if (a == 0) return 1;
    if (a == 1) return x % m;
    uint64_t t = x;
    uint64_t out = 1;
    while (a > 0) {
        if (a & 1) out = (uint64_t(out) * t) % m;
        t = (uint64_t(t) * t) % m;
        a >>= 1;
    }
    return out;
}

template<class T>
T addm(T a, T b, T c) { T s = a + b; return (s >= c) ? s - c : s; }

int main()
{
    for (int i = 1; i <= N; i++) {
        vec[i] = ipowMod(i,i,M);
    }
    num[0] = 1;
    for (int i = 1; i <= N; i++) {
        int v = vec[i];
        uint64_t ad[M];
        std::copy(num, num + M, ad);
        for (int j = 0; j < M; j++) {
            int p = addm(v, j, M);
            num[p] = addm(ad[p], ad[j], n16);
        }
    }
    printf("%llu\n", num[0] - 1);
}
でHaskellでも書いてみた.
import Data.List
import Data.Array
import Data.Word

n=250250
m=250
n16::Word64
n16=10^16

pow mul x n = f x (n-1) x
 where f _ 0 y = y
       f x n y = g x n where
                 g x n | even n = g (mul x x) (n`quot`2)
                       | otherwise = f  x (n-1) (mul x y)

{- x^n mod m -}
ipowMod x n m = pow mul x n
 where mul x y = rem (x * y) m

vec::Array Int Int
vec = listArray (1,n).map f $ [1..n]
 where m'::Word64
       m' = fromIntegral m
       f i = fromIntegral $ ipowMod i' i' m'
          where i'::Word64
                i' = fromIntegral i

{- a + b mod c -}
addm a b c = let s = a + b in if s >= c then s - c else s

up num i = array (0,m-1) $ map f [0..m-1]
 where f j = (p, addm (num ! p) (num ! j) n16)
        where p = addm (vec ! i) j m

num::Array Int Word64
num=listArray (0,m-1).take m $ 1:[0,0..]
calc = foldl' up num [1..n]

ans = (calc ! 0) - 1
main = print $ ans
n=25025と一桁落とした場合は2秒ぐらいだったのでこのままでもいけるかと思ったらn=250250ではメモリ不足で落ちた. foldl'でstrict評価してるのになんでメモリが足りなくなるのだろう. GC回収してくれ. 悩むのも面倒なのでSTArrayを使ってup以下を書き直す.
import Control.Monad.ST
import Data.Array.ST

{- 上と同じ -}

up num i = do
 let v = vec ! i
 ad <- newArray_ (0,m-1) :: ST s (STArray s Int Word64)
 mapM_ (cp ad) [0..m-1]
 mapM_ (f v num ad) [0..m-1]
 where cp a j = do
          x<- readArray num j
          writeArray a j $! x
       f v ao ai j = do
          let p = addm v j m
          x <- readArray ai p
          y <- readArray ai j
          writeArray ao p $! (addm x y n16)

calc = do
 num <- newArray (0,m-1) 0 :: ST s (STArray s Int Word64)
 writeArray num 0 1
 mapM_ (up num) [1..n]
 getElems num

ans = runST calc
main = print $ (ans !! 0) - 1
これなら7.36秒でC++の10倍遅い程度ですんだ. でもどうみてもC++のループの方が簡単. 速くてわかりやすく書く方法が知りたいなあ. [2011/4/13 追記高速化遍]
11月9日_
Project Euler 249, 250をひよってC++で解いてしまった. またshuyoさんに怒られる. Haskellで書き直さなきゃ.
11月7日_
自然言語処理勉強会参加. 統計的係り受け解析の話がとても面白かった. 係り受けのパターンが何通りあるか?なんてProject Eulerに出そうだなとかどうでもいいこと思ってしまったのだけど. 終わったらジュンク堂行ってみるかと考えていたら, がっつり8時過ぎまであって疲れて帰る. お腹減りすぎで渋谷でパンかって電車の中で食べてしまった. 頭も働かないし6時半過ぎたらさくっと帰ったほうがいいかな.
11月6日_
風邪引いてくしゃみと鼻水が酷い. 今日は早く寝る.
11月5日_
Project Euler 193
250未満のsquarefreeなものの個数. 最初見たときは, 素数を数えるだけでもでかすぎてお手上げ. しばらくして逆から数えてがんばってみても235で1分超えて諦める.
今日, もうちょっとがんばってようやく解けた. Haskellで28秒. C++なら2.6秒. 解いたあとでフォーラムを見たらメビウス関数がうんたらかんたら. そんなものもあったなあ(遠い目).
import Data.List
import Data.Array

limit=2^50

{- primesは素数の無限列, isqrtは整数の平方根(略) -}

ps = map (^2).takeWhile(<isqrt limit) $ primes
pa = listArray (0,(length ps)-1) ps

f::Integer->Int->Integer
f 0 _ = 0
f n j | n == (fromIntegral j-1) = n
      | otherwise = n - sum [f n' (i+1)|(n',i)<-sub]
 where sub = takeWhile((/=0).fst) [(quot n (pa ! i), i)|i<-[j..(length ps)-1]]

main = print $ f limit 0
typedef std::vector<uint64_t> IntVec;

/* primesは2^25までの素数の列(略) */

IntVec ps(primes.size());

uint64_t f(uint64_t n, size_t j)
{
    if (n == 0) return 0;
    if (n == j-1) return n;
    uint64_t sum = 0;
    for (size_t i = j; i < ps.size(); i++) {
        uint64_t x = n / ps[i];
        if (x == 0) break;
        sum += f(x, i + 1);
    }
    return n - sum;
}

int main()
{
    for (size_t i = 0; i < primes.size(); i++) {
        uint64_t p = primes.get(i);
        ps[i] = uint64_t(p) * p;
    }
    printf("%llu\n", f(limit, 0));
}
11月4日_
・季節の変わり目で眠たいのにProject Eulerやってるので全然時間が無い. ようこんな問題, さくさく解いてるなあ. 電車の中で考えてると本も読めない. 192番は単なる連分数の問題かと思ったら微妙に微調整が必要だった. 補正方法を適当に書いたらうまくいった. ただし理屈はちゃんと考えてない.
クラウド時代の暗号解読を聴講. 攻撃側の話とか, EC2を使って解くとどれぐらいお金がかかるとかの話が面白かった. 終わってから講師の方とお知り合いと食事.

御意見は光成滋生<herumi@nifty.com>までお願いします