melancholic afternoon

Finite formula found for partition numbers
3ヶ月ほど前にProject Eulerでもよくでてくる整数の分割関数p(n)の新しい関係式が見つかったそうな. 漸化式や無限和じゃなくて, 代数的数の有限和の形でかけるとか(an algebraic formula for the partition function). 残念ながら, ぱっとみ定理をp(1) = 1で納得するのも難しい….

・東京に来たときに自分の部屋があまりに殺風景なので買ってみた観葉植物(幸福の木とか言ったっけ)がもう枯れかけで, 先週末駄目元で土を替えてみた. そこから更に枯れてきたのだけど3枚の葉っぱだけぎりぎりのところで残ってる. このまま復活してくれるのだろうか.
・査読一つ出した. さて, もう一つ. GWにすればいいだろうと思ってたら明日からGWだということに気が付いた.

社内勉強会でPRMLの3章が終わったのでまとめてみた. prml3.pdf

『タロットの迷宮』(小笠原慧)
『DZ』がおもしろかったので読んでみた. 読み始めてしばらくしてシリーズものの2冊目だと分かる. どこかに書いといて欲しい. 内容は期待外れ. 一部SF的なアイテムが出てくる割りにはなんとなく古くさい. 私はタロットといった占い的なものはどうも好きじゃないようだ.

自然言語処理勉強会参加
n_shuyoさんの話は昼休みにちょいちょいと聞いてた. このあたりの話は, 検証できなくて仮説が多いので, みんな自己主張するする(苦笑).

20日の配列の最大値のコード, 64bitのVCだと酷いコードが生成されていることに気が付いた. えっ, なんで. もちろん32bitのと同じオプションの/Oy /Ox /Ob2 /GS-.

  mov    ebx, 1
  mov     edi, ecx
  mov     dword a$3841[rsp+4], 5
  mov     dword a$3841[rsp], ebx
  mov     edx, ebx
  mov     r8d, ebx
  mov     dword a$3841[rsp+8], 9
$LL6@main:
  mov     eax, dword PTR a$3841[rsp+r8*4]
  cmp     eax, edx
  cmovg   edx, eax
  inc     r8
  cmp     r8, 3
  jb  SHORT $LL6@main
  lea     rcx, "%d\n"
  call    qword __imp_printf

数論暗号勉強会に参加. グレブナー基底入門. ブッフバーガーさんは3回しか指導教授のグレブナーさんに会ってないというエピソード. 3回しか会ってないのによく教授の名前をつけるつもりになったなあ.

initializer_listが内部的には一端配列になるから昨日のようなコードになるのかと思って別のコードを試してみた.

#include <stdio.h>

template<class T, size_t N>
T my_max(const T (&x)[N])
{
    T m = x[0];
    for (size_t i = 1; i < N; i++) {
        if (x[i] > m) m = x[i];
    }
    return m;
}

int main(int argc, char *argv[])
{
    {
        int a[] = { 1, 5, 9 };
        printf("%d\n", my_max(a));
    }
    {
        int a[] = { 3, 4, argc, 5, 3 };
        printf("%d\n", my_max(a));
    }
}

このコードをgcc 4.6.0でコンパイルすると

    movl    $9, 4(%esp)
    movl    $.LC0, (%esp)
    call    printf
    movl    $4, %edx
    movl    $5, %eax
    cmpl    $4, %ebx
    cmovge  %ebx, %edx   // edx = (argc >= 4) ? argc : 4
    cmpl    $5, %edx
    cmovge  %edx, %eax   // eax = (edx >= 5) ? edx : 5
    movl    $3, %edx
    cmpl    $3, %eax
    cmovle  %edx, %eax   // eax = (eax >= 3) ? eax : 3
    movl    $4, 32(%esp) // (A)
    movl    $5, 40(%esp) // (A)
    movl    $3, 44(%esp) // (A)
    movl    %eax, 4(%esp)
    movl    $.LC0, (%esp)
    call    printf

となった. { 1, 5, 9 }の即値の配列に対しては9という即値を得ている. 素晴らしい. 二つ目に関しては, (A)の4, 5, 3をスタックに積んでる部分が配列を作りかけて途中でやめてしまったゴミのようだ. 使わないからコード自体削除してもよさそうなので残念. ただ簡単な配列だったら作らないという最適化はある程度するようだ.
VC9だと

  push     9         // 即値9
  push     offset "%d\n"
  call     esi       // printf

  mov     eax, dword _argc$[esp+8]
  add     esp, 8
  cmp     eax, 4
  jle     .L1 // argcが4より小さいならL1へ
  cmp     eax, 5
  jge     .L2 // 5以上ならL2へ
.L1:
  mov     eax, 5    // 5に確定
  push     eax
  push     offset "%d\n"
  call     esi       // printfして終了
  add     esp, 8
  xor     eax, eax
  pop     esi
  ret
.L2:
  cmp     eax, 3
  jge     .L3       // 3より大きいならL3を (B)
  mov     eax, 3    // そうじゃないなら3を (B)
.L3:
  push     eax
  push     offset "%d\n"
  call     esi

(B)の部分はもう一声欲しかったが, 配列は生成していない.

C++0xで追加されたvariadic templateとinitializer_list, どちらがよいのかなと試してみた.

#include <stdio.h>
#include <algorithm>

template<class T>
T my_max(const T& a)
{
    return a;
}
template <class T, class... Args>
int my_max(const T& x, const Args&... xs)
{
    return std::max(my_max(xs...), x);
}

int main(int argc, char *argv[])
{
    printf("%d\n", std::max({1, 5, 9}));
    printf("%d\n", std::max({3, 4, argc, 5, 3}));
    printf("%d\n", my_max(1, 5, 9));
    printf("%d\n", my_max(3, 4, argc, 5, 3));
}

試したのはgcc 4.6.0(32bit)
まずはinitializer_listを使った方の出力.

   movl    $1, %eax
   movl    %ebx, %edx
   movl    $1, 40(%esp) // 1をスタックに積む
   movl    $5, 44(%esp) // 5をスタックに積む
   movl    $9, 48(%esp) // 9をスタックに積む
   jmp .L7
  .p2align 4,,7
  .p2align 3
.L8:
   movl    4(%edx), %edx // 一つずつ取り出して
   cmpl    %edx, %eax
   cmovl   %edx, %eax    // 比較して
   movl    %ecx, %edx
.L7:
   leal    4(%edx), %ecx
   cmpl    %esi, %ecx
   jne .L8               // ループ
   movl    %eax, 4(%esp)
   movl    $.LC0, (%esp)
   call    printf
    (略)

ださい, ださすぎる. いちいちデータをスタックにコピーしてから比較してる. 一度配列に変換されている.
で, variadic templateの場合.

    movl    $9, 4(%esp)   // max(1, 5, 9)
    movl    $.LC0, (%esp)
    call    printf
    cmpl    $5, %edi
    movl    $5, %eax      // eax = 5;
    cmovge  %edi, %eax    // eax = (argc <= 5) ? 5 : argc
    movl    %eax, 4(%esp)
    movl    $.LC0, (%esp)
    call    printf

うって変わってこちらの生成コードのなんと見事なことか. 定数は全てコンパイル時に確定. max(3, 4, argc, 5, 3)も

 max(3, max(4, max(argc, max(5, 3))));

に展開された後, 一番大きな5とだけ比較すればよいと判断したコード生成になっている.

int func1(int a, int b, int c, int d, int e)
{
    return std::max({a, b, c, d, e});
}
int func2(int a, int b, int c, int d, int e)
{
    return my_max(a, b, c, d, e);
}

についても同様にinitializer_listの場合は一端引数をスタックにコピーしていたが, variadic templateの場合は直接見るコードが生成された. 素晴らしい.
VCも早くvariadic templateに対応して欲しいものだ.

「マイケル・サンデル　究極の選択」大震災特別講義

もう少し踏み込んだ議論を期待したけど, 題材が生々しすぎて腫れ物に触るような印象. 日本の参加者もちょっと優等生気味な嫌いがあったように思った.

アメリカの災害のときに起きた暴動について, 個人主義だから自分の身は自分で守るという意見があった. 言葉の定義だけども個人主義をそういうところで使うのは好きじゃない. 個人主義は各自の自律(もちろん他人のも)を尊重することであって, 自分の好きなようにするというのは単なる利己主義だ. ああいう議論をするときは明確に区別すべきじゃなかろうか.

第2回さくさくテキストマイニング勉強会参加

よい勉強会だった. 進行役の人が熱かった. それぞれの発表が始まる前に, 発表者がどういう人で, この発表に何を期待してるかを話す. やや大袈裟なところもあるのだけど, それが駄目というのでは全然なくて, 「ほほう, 聞いてみようか」という気にさせられる. 終わったら, どういう点がよかったという感想もきちんと述べる. あれは素晴らしい座長だ.

それぞれの発表も, やってみたらおもろいよと思わせる遊び精神が多かった. 最近, お勉強が続いてたので私も簡単なものから手を動かしてみよう.

会場も快適だった. 会場提供のOracleの方, ありがとうございました.

PRML3章やってて派生した問題.
「A, Bを実n次対称行列とするとき, A > B > 0ならばB^-1 > A^-1か」
ここでA > 0は任意のx(≠0)に対して^txAx > 0を意味する.
なんか盛大に勘違いしていたのでやり直す.

示せたかな.
C = A - Bとおくと C > 0. またB > 0からB^-1 > 0.
更に B^-1 - A^-1 = B^-1(A - B)A^-1 = B^-1CA^-1.

すると
A(B^-1 - A^-1)A = AB^-1C = (B + C)B^-1C = C + CB^-1C.

よって任意のx(≠0)についてy = A^-1xとおくと
^tx(B^-1 - A^-1)x = ^ty (C + CB^-1C) y = ^tyCy + ^t(Cy)B^-1(Cy) > 0.

『サウスポー・キラー』(水原秀策)
野球賭博の疑いを掛けられた投手が濡れ衣を張らそうとする. なんか久しぶりにオーソドックスなミステリーを読んだ. 派手さは無いのについつい引き込まれるのはたいしたものだと思う. やな性格のオーナーや周りの選手の中で, 熱くなりすぎることも無く, でもきっちり主張する主人公のキャラが気に入った. これはいい.

@tanakhさんに高速版コードをいただいた. どうもありがとうございました. 彼によるとC++の2倍遅いぐらい(C++ : Haskell = 0.209s : 0.401s)にまではなったそうな.

しかし, 残念ながら私の環境(32bit Windows XP(Core2Duo 2.6GHz) + GHC 6.12.3(Haskell Platform 2010.2.0.0) + vector-0.7.0.1)では実行すると

Segmentation fault/access violation in generated code

と表示されて動かなかった. どうやらIntを64bitを仮定してるからのようだ. HaskellのIntは2²⁹ぐらいまでを持っていることという規約なので私の32bit環境では駄目みたい. コンパイルできてるのに実行時にさくっと落ちるのはunboxedだからか? 実は250を解くときに使ってみようとしたのだけどすぐにはできなくて諦めていたのであった. 今回実例を見せていただいたので再挑戦してみるかも.

というわけで自分の環境でも動くように最小限の修正をしてみた. なんかすぐsegvで死ぬし, 動いたと思っても答えが違ってたりで難しい. n16はWord64に, ipowModは演算途中をWord64に変更. addmはIntとWord64の両方必要.

{-# Language BangPatterns #-}
{- ghc -O2 --make -fvia-C vec.hs -}

import Control.Applicative
import Data.Vector.Unboxed.Mutable as V
import Data.Word

n :: Int
n = 250250

m :: Int
m = 250

n16 :: Word64
n16 = 10^16

ipowMod :: Int -> Int -> Int -> Int
ipowMod x n m = f x (n-1) x
 where
   f _ 0 y = y
   f x n y = g x n where
     g x n | even n = g (x `mu` x) (n`quot`2)
           | otherwise = f  x (n-1) (x `mu` y)
   mu::Int->Int->Int
   mu x y = fromIntegral $ to64 x * to64 y `rem` to64 m
   to64::Int->Word64
   to64 x = (fromIntegral x)

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

for :: Monad m => Int -> Int -> (Int -> m ()) -> m ()
for !f !t m = go f where
  go !f | f == t = return ()
        | otherwise = m f >> go (f+1)

main :: IO ()
main = do
  num <- V.replicate m 0
  ad <- V.new m
  V.unsafeWrite num 0 1

  for 1 (n+1) $ \i -> do
    let v = ipowMod i i m
    V.unsafeCopy ad num

    for 0 m $ \j -> do
      let p = addm v j m
      V.unsafeWrite num p =<<
        addm <$> V.unsafeRead ad p <*> V.unsafeRead ad j <*> pure n16

  print =<< (pred <$> V.read num 0)

これだと3.81秒になった. 私が示したHaskell版の解答は7.39秒なので約2倍の高速化. それでも, 同環境でC++版は0.66秒なのでまだHaskell版は5.7倍ぐらい遅い. @tanakhさんとの環境での結果と大分違うが, あちらはおそらく64bit版だからだろう(GHC 7.x系だからというのも多少あるのかも). 32bit環境ではWord64だとunboxの効果が薄れるのかな. もう少し速くできないものか.

@shelarcyさんによるとGHC 6.x系はvectorが少し遅いということらしいが私がProject Euler 250の自分の解答を書いてたときはまだGHC 7系は出てなかったし, 今回は触れない. 機会があれば後日調べる.

別の環境で測ってみた.

各コードの実行時間[sec]
	C++	私の解答	@tanakhさんの	@tanakhさんの32bit用修正版
32bit Win Xp(VC9 / GHC 6.12.3)	0.66	7.39	落ちる	3.81
64bit Mac (gcc 4.2.1 / GHC 6.10.4)	0.251	2.571	19.270	19.664
64bit Linux (gcc 4.4.5 / GHC 6.12.1)	0.13	3.01	0.480	0.480

環境によって全然違うなあ. Macのは@nishioさんによると一応portの最新版らしいけどvector版めちゃ遅い. 私の解答の方が速い. Linux版はHaskellがC++の3.5倍遅い程度と@tanakhさんのに近い.

いくつか分かったこと

-fvia-CはCコード経由で最適化するためのオプションらしい. ただ将来の環境では非推奨らしい.
{-# INLINE addm #-}で速くなることがある. 最適化で勝手にinlineしてくれるもんやないのね. C++の感覚と違う.
powのbinary法による実装は, 今回の整数のmodだけでなく行列などのときにも使えるので別のモジュールに作っている. 今回は埋め込んでいたけど, それを分離してみたら4.23秒と遅くなった. これもC++のtemplateの感覚からすると同じ速度であって欲しい.
forとか手続き型っぽいな. 速いHaskellコードは見た目がどんどんCに近づくの仮説.
unboxedは辛い. スペースリークで返って来ないのはまだパラメータを小さくすれば分かるけど, こちらは何が起こってるのか分からず即死なのが辛い.
GHCはバージョン書かないと速い遅いが意味を成さないほど速度が変わることがある.

この結果だけから見ると私のコードも「こういうコードが転がってて、Haskell遅いとか言われるのは問題だな・・・」と言われるほど酷いとは思えないなあ. この問題に関しては, リストを使わず, 副作用ありにするためmonadにして, 遅延評価を止めるためにBangPatternsを使い, 最新の処理系と非標準のライブラリを使って, unboxedを使って64bitでしか動かなくなって, C++の半分の速度なのかという感じ.

Project Euler 250の自分の解答
@tanakhさんに「こういうコードが転がってて、Haskell遅いとか言われるのは問題だな・・・」と言われた. ぜひぜひどういう風に書くのが良いのか知りたい(それなりに調べたり試したけどよい情報を得られなかったので, そう書けばどなたかつっこんでくださるだろうと, 多少上記日記は挑戦的に書いたところもあるのであった. すいません).
ただ一応標準しばりでやってるのでGHC6.12.3では入ってなかったvectorは使わないという方法がありがたいです.

Haskell初心者としては

大抵の入門書がリストしかでてこなくて, 普段はよいけどこういう問題を解こうとしてその遅さに悩む.
で, 配列は無いものかと探してData.Arrayを使ってみて, その遅さに悩む.
で, Data.Arrayの更新がO(n)だと知って衝撃を受ける.
で, 書き換えが必要なんだ, 世の中副作用無しのpureなんて無理なんだと思ってSTArrayに手を出す.
なんて使いにくいんだ, でも結構速くなったぞ.
あれ, それでもでもC++の1/10ぐらい. まだ駄目か….
こうなったら遅延評価も捨ててunbox化を…しようとしたけどちょっと間違えるとすぐランタイムエラーになってよう分からん. もういいや.

リストが駄目で, 副作用ありにしないと駄目で, 遅延評価も捨てないといけないならHaskellって本当にC++並に速いと言っていいの?という心境でして, その自分の評価が変わるなら嬉しい.

こう書くと, すぐ全称記号的に, 「○○の場合はHaskellでも速いからそれは嘘だ」という意見を見るけどいつでも遅いとは言ってない. HaskellがC++よりもよい分野もある. そうじゃなくて配列的なものを操作する計算主体ではやっぱり遅いんじゃないのかという(初心者の)感想. そしてそれが違うのならどうやればいいのか知りたいというだけです.
#欲を言うとProject Euler 258をC++並の速度(数秒)で解く方法も知りたい.

無性にピザが食べたくなって久しぶりに緑が丘の近くのピザ屋さんに行ってみた. 生地がしっかりしてておいしい. 食べ終わったら緑道沿いに桜が咲いてたのでぶらっと歩いて帰る. もう春なんだなあ.

『月の恋人』(道尾秀介)
どんなミステリーかと期待しつつ読み始めたっけど, えらくキャラクターが類型的で, ストーリーや細かいエピソードもあたりきなものばかり. なんか一昔前のトレンディードラマみたいだと思ったら去年の月9用の話だったらしい. しまった.

サイボウズ・ラボユースで林さんが開発したJavaScriptによる型推論器で遊んでみた.
型推論の途中を表示するVisualized type inferer using slideshowはかなりおもろい. Haskellにあったらいいなあ. まあ, それはともかくProject Eulerの最初の方を解いてみたけどまだ機能に制限が多いので辛い.
PE1: 1から999までの数で3か5で割れるものの和

let f = fun n -> fun r -> r * (n/r) * (n/r + 1) / 2
;;
let pe1 = fun n -> f n 3 + f n 5 - f n 15
;;
let ans = pe1 999
;;

PE2: 400万以下のフィボナッチ数列のうち偶数なものの和
まだリストが無いので作りながら足していく

let limit = 4000000
;;
let rec fib_sum = fun sum -> fun x -> fun y ->
 if x >= limit then
  sum
 else
  if x mod 2 = 0 then
   sum + x + fib_sum sum y (x + y)
  else
   sum + fib_sum sum y (x+y)
;;
fib_sum 0 1 2
;;

PE3: 600851475143の最大素因数
うちのマシンは遅いのでもうちょい小さい数字71 * 839 * 1471 = 87625999で試す.

let rec f = fun n -> fun p ->
 if n < p * p then
  n
 else
  if p = 2 then
   if n mod 2 = 0 then
    f (n / 2) 2
   else
    f n 3
  else
   if n mod p = 0 then
    f (n/p) p
   else
    f n (p+2)
;;
let pe3 = fun n -> f n 2
;;
pe3 87625999

100までの素数を表示する

let rec f = fun i -> fun n ->
 if i * i > n then
  true
 else
  if n mod i = 0 then
   false
  else
   f (i+1) n
;;
let isPrime = fun n -> f 2 n
;;
let rec showPrimeTo = fun n ->
 if n = 1 then
  ""
 else
  if isPrime n then
   showPrimeTo (n-1) ^ " " ^ show n
  else
   showPrimeTo (n-1)
;;
showPrimeTo 100
;;
---
結果
- : (int -> (int -> bool)) = <closure>
- : (int -> bool) = <closure>
- : (int -> string) = <closure>
- : string = " 2 3 5 7 11 13 17 19 23 29 31
    37 41 43 47 53 59 61 67 71 73 79 83 89 97"

今朝の朝日新聞の社説余滴という欄で, 著者が, 彼のボスである論説主幹が書いた「新しい日本へ造り直していく」という意見に真っ向から反論していた.

えらく気合の入った文章だったので(そんなところで上司に楯突いてる様が)不謹慎にもちょっと笑ってしまった. 内容はまっとうなもので, リセットなんてできない, 過去の全てが否定されたわけではない, 性急はいけないというものだ.

何故早急に0/1の二値的な判断を下そうとするのか. もっと中庸に物事を決めようとできないのか. 私の世代はセンター試験の○×的な軽さやバーチャルとかデジタル人間などと揶揄されたものだが, そう揶揄していた人たちは, それよりもっと酷くないか.

震災で科学技術への信頼は崩壊したのか

福島第一原発の深刻な事態で, 震災を機に科学技術への信頼は崩壊したという論調の文章をよく見かける. でも本当にそうなのか. あの凄まじい地震でも, 200km/h以上で走っていた新幹線は1台も脱線することなく止まった. 私はあんなに揺れたら絶対脱線してるだろと思っていたので全然そういう事故の報道が無かったことに驚いた. 福島第二原発や女川原発も自動停止して安全な状態になった.

福島原発についてはたとえば福島原発１０基の耐震安全性の総点検等を求める申し入れに

福島原発はチリ級津波が発生した際には機器冷却海水の取水が出来なくなることが、すでに明らかになっている。これは原子炉が停止されても炉心に蓄積された核分裂生成物質による崩壊熱を除去する必要があり、この機器冷却系が働かなければ、最悪の場合、冷却材喪失による苛酷事故に至る危険がある。

などど言われていた, 他にも同様のことを指摘していた人はいるようだ. 全くの「想定外」では無かったわけだ. ということは, そう予測されていたにも関わらず, 対策をとらなかったのが問題だ.

何故正しい対処ができなかったのか, それはコストを見誤ったからか, 一部の利権か, 政治的な駆け引きか, どんな理由があったからかは知らない.

技術が正しい予測をしていてもそれを判断する人間が駄目だったから悲惨なことになったのだ. その中には東電の人たちだけでなく, 政治的な立場の人たちも当然いるだろう. それをきちんと明らかにすべきだと思う.

もちろん, そういう人を選んでしまった責任の一部は国民の側にもある. もしかしたら自分も入ってるのかもしれない. どうみても不要なところに新幹線の駅を誘致してしまった政治家は税金の無駄遣いをして腹立たしいというだけだが, 原発での対策ミスではとんでもないことになってしまった. そういうことだと思う.

fmathのベンチマークを取り直してる. 相変わらずgccのexpが遅いように見える. Macのgccだとそんなに遅くないという報告もある. ということはLinuxのgccやcygwinだけが遅いのか. それとも私がみたマシンだけが遅いのか.
それからAVXは従来のSSEと混ぜてはいけないのだから, 実行時に命令を切り替える場合は, AVX対応オプションをつけてコンパイルされてるかそうでないか知りたいよなあ. gccは__AVX__のありなしで判別できるけどVCには無いみたい(AVX用オプション). それは困る. 早く対応してMicrosoft.

3日前に書いてたpdfの日本語文字の間に妙なすき間ができる現象はTeXの問題ではなくビューワ(foxit)の問題だった. acrobat readerはもう久しく使ってないのだけどそんな微妙な違いがあったとは.

ぼちぼち桜も咲き始めた. 最近吉野山に登ってない. また行きたいなあ.

『DZ』(小笠原慧)
複数の登場人物からの視点で話が同時に進みつつ少しずつ収斂していく. 先が気になってしかたがない. タイトルの意味が分かるのは物語も佳境に入ってから. これは解説を見ないで読んだ方がいいだろうなあ.
一つ難をつけるなら, ラストである人が別の人に打ち明ける事実はちょっと偶然すぎないかなあと思った. 別にそれなくてもいいんじゃないかと.

・ニッポン人脈記の「成長の陰に女性の叫び」
インドにはダウリーという, 花嫁の家族が花婿に年収の2～5倍の金や物を貢ぐ習慣があるらしい. で, それが払えないと花嫁が焼き殺されるとか. 2009年は年間8400人殺されているってこれはエイプリルフールか何かですか. 内容が酷すぎて理解できないのだが.

・これから暗号関係をやめて機会学習に専念しようと今日なら言ってもいい.
さて, 社内でPRMLの2章を読んだときの一部の数学をまとめてみた(prml2.pdf). 間違い指摘, つっこみ歓迎.

自分の環境ではdviからpdfに変換するとき日本語の文字の間に妙なすき間ができてしまった. どうやったら直るのだろう. 自分の環境はTeXインストーラ 3のabtexinst_0_81r4.zipを使ってさっきクリーンインストールしたばかり. dvipdfmxで変換してある. もとソースはgithub.com/herumi/prml/に.

憂鬱な午後のひととき