憂鬱な午後のひととき

11月29日_

・4年振りにboost.Pythonを触る. 以前触ったときはVC2005でエラーになったりしていろいろしないといけなかったけど, 今回VC2008を使った場合はすんなりいった. ただドキュメントなどもろもろあんましメンテされてない感は漂ってますが. とりあえずメモ.

環境はVC2008(32bit) + Python 2.7.1 + boost-1.48.0. まずboostのPythonライブラリをbuildしておく. stage/lib/以下ライブラリができるので適当にコピー. 

bootstrap
b2 --with-python --toolset=msvc-9.0 address-model=32 cxxflags="-D_SECURE_SCL=0 -GS-" \
  threading=multi link=static runtime-link=static,shared variant=release stage

-D_SECURE_SCL=0はVC2010なら不要. 他に必要なライブラリがあれば追加. bjamコマンドがいつのまにかb2になっていた. 古いbjam.exeがパスに入ってると誤動作するかもしれない.

C++でモジュールを書く. 

// python-test.cpp
#define BOOST_PYTHON_STATIC_LIB /* 定義するとboost.pythonのDLLが不要になる */
#include <boost/python/module.hpp>
#include <boost/python/def.hpp>
#include <string>

std::wstring sayHello(const std::wstring& name)
{
  return L"こんにちは，" + name + L"さん";
}

BOOST_PYTHON_MODULE(mie)
{
  boost::python::def("sayHello", sayHello);
}

setup.pyを用意する. パスやオプションは適宜修正すること. 

#!/usr/bin/env python
# setup.py
from distutils.core import setup, Extension

module = Extension(
    'mie',
    sources = ['python-test.cpp'],
    include_dirs = ['c:/p/boost'],
    library_dirs = ['c:/p/boost/lib'],
    extra_compile_args = ['-D_CRT_SECURE_NO_WARNINGS', '/EHsc'], # for vc9
)

setup(
    name = 'mie',
    version = '1.0',
    ext_modules = [module],
)

これでpython setup.py installとしてエラーがでないことを確認する. そしたらPython側で使ってみる. 

# python-test.py
# -*- coding: cp932 -*-
import string
import mie

print mie.sayHello(u"太郎")

>python python-test.py
>こんにちは，太郎さん

さくっとユニコードを扱える. boost::python::long_という多倍長整数の型もある. 

#include <boost/python/long.hpp>
/* あくまでサンプル. 効率の悪いアルゴリズム */
boost::python::long_ power(boost::python::long_ x, int n)
{
    boost::python::long_ y(1);
    for (int i = 0; i < n; i++) {
        y *= x;
    }
    return y;
}
BOOST_PYTHON_MODULE(mie)
{
    boost::python::def("power", power);
}
--
>>>import mie
>>>mie.power(1234L, 56)
12991190255487145194103208439623513775465782010127392384379012704624259433055094
64892567848536247290201061395156473849109449211865238658490562753590662623529116
82504769929216L

よい感じ.

・眠たすぎ. 朝, 飯田橋を乗り過ごした. 顎のあたりに腫れ物ができて痛い. 体調悪いなあ.

11月28日_
・ここ数日, 明らかなスパムメールがスパムでないと判定されてるなあ. なぜだろう.
・眠たすぎて全然プログラムできない. Grundy数ってなんかよう分からんな.
11月25日_
・恐る恐るジョギングしてみる. これぐらいは大丈夫なのか. しかし朝6時, 流石に半袖はさぶい.
・社内PRML読書会が終わったので今度は私はしばらく『入門 自然言語処理』を読んでみる. Python上で対話的にNLTKを触る. このいたれりつくせり感は, 古き良き時代のBASICを思い出す. 楽しい.
11月22日_
・国際宇宙ステーションから帰還の古川さん, 着地したところはカザフスタンの零下20度の雪原だそうな. えらいところに降りるもんやな.
・こつこつと働いて稼いでたお金を社長にカジノでぱっと使われたらやりきれない. 何で誰も止められなかったのかな.
11月18日_
無事(?)二人分の論文チェックが終わった. 多分ちゃんと提出してくれるだろうと期待して寝る.
11月17日_
・ここ何日か学生の方の論文推敲を手伝ってる. でも, それより自分自身のもしないとやばい.
・何カ月も前にNASMのバグ報告をメールしたはずなのだけどなかなか直らない. パッチ(ものすごい簡単)送ったらすぐ直してもらえた. やっぱり些細なことでもちゃんとアクションしないと駄目だなあ.
11月11日_
AVX2のmulx
以前ちょっとだけマニュアルを見たときは気が付かなかったけど使いやすい整数乗算命令も増えてた. 
mulx dest1, dest2, src
で64bit x 64bit → 128bitを[dest1:dest2] = src * rdxの形で求められるようだ. しかもフラグ変化を起こさない. これは便利だ. 今までのmulはdestinationが[rdx:rax]に決まっていて更にフラグが変化するので多倍長演算のときにレジスタの使い回しがややこしくてしかたなかった. 早く試してみたいな. ってこのCPUはいつ出るの?
11月10日_
『舟を編む』(三浦しをん)
こないだまとめて買った本の一冊. 25歳ぐらいの男性が, 突如辞書を編纂する部署に抜擢(?)される. ある言葉をどう簡潔に説明するかとか, 辞書のためのよい紙を作ろうとする執念や情熱がすごい. 非常におもしろい. 今年読んだ中で一番か.
話の中で出てきた「左」をどう説明するかは昔よく話題にしてたなあ. 方角を使うと, 今度は方角をどう説明するかで困るんだよな. 大学生のときは冗談で通常観測されるニュートリノのスピンの向きにすればいいなんて言ってたこともある. ところで「一番最初」ってベテラン辞書編集部員が使うものかな. 私はつい使ってしまうけど.
11月9日_
↓嘘の嘘書いた. 頭痛が酷くて寝られなかったのでバファリン飲みながら夜中3時に書くと間違える. 
lea rax, qword [32bit定数]
は32bit定数を正しく符号拡張する. 
#include <stdio.h>
#include <stdint.h>
#include "xbyak/xbyak.h"

const size_t codeSize = 4096;

struct Code : public Xbyak::CodeGenerator {
    Code(uint32_t x, void *userPtr)
        : Xbyak::CodeGenerator(codeSize, userPtr)
    {
        printf("addr=%p\n", getCurr());
        lea(rax, ptr [(void*)(uint64_t)x]);
        ret();
    }
};

int main()
{
    printf("main=%p\n", main);
    using namespace Xbyak;
    uint8 buf1[codeSize + 16];
    uint8 *stackBuf = CodeArray::getAlignedAddress(buf1);
    CodeArray::protect(stackBuf, codeSize, true);

    static uint8_t buf2[codeSize + 16];
    uint8 *staticBuf = CodeArray::getAlignedAddress(buf2);
    CodeArray::protect(staticBuf, codeSize, true);

    const struct {
        const char *name;
        void *ptr;
    } tbl[] = {
        { "heap", 0 },
        { "stack", stackBuf },
        { "static", staticBuf },
    };
    for (size_t i = 0; i < sizeof(tbl) / sizeof(*tbl); i++) {
        printf("%s:\n", tbl[i].name);
        const uint32_t valTbl[2] = { 0x12345678, 0x87654321 };
        for (int j = 0; j < 2; j++) {
            Code c(valTbl[j], tbl[i].ptr);
            void* x = ((void*(*)())c.getCode())();
            printf("%x -> %p\n", valTbl[j], x);
        }
    }
}

// 64bit Linux
main=0x400c98
heap:
addr=0x2299000
12345678 -> 0x12345678
addr=0x2299000
87654321 -> 0xffffffff87654321
stack:
addr=0x7fff1eb57790
12345678 -> 0x12345678
addr=0x7fff1eb57790
87654321 -> 0xffffffff87654321
static:
addr=0x606220
12345678 -> 0x12345678
addr=0x606220
87654321 -> 0xffffffff87654321

// 64bit Windows
main=000000013F901267
heap:
addr=00000000004EB000
12345678 -> 0000000012345678
addr=00000000004EB000
87654321 -> FFFFFFFF87654321
stack:
addr=000000000025EC60
12345678 -> 0000000012345678
addr=000000000025EC60
87654321 -> FFFFFFFF87654321
static:
addr=000000013F6AC2E0
12345678 -> 0000000012345678
addr=000000013F6AC2E0
87654321 -> FFFFFFFF87654321
Linuxのstack, Windowsのコードとstatic変数は4GB以降の領域に置かれる.
11月8日_

32bit定数cを64bitレジスタに代入するいくつかのやりかた. ちょっと面白かった.

cの値によらず0拡張した場合 

mov rax, 0x12345678 ; 7byte
mov rax, 0x87655431 ; 10byte

でかい. 実は32bitレジスタをいじると自然に0拡張されるので 

mov eax, c

とすれば5byteですむ(rexプレフィクスがあれば+1byte).

cの最上位ビットが立っていて符号拡張したい場合 

mov rax, (0xffffffff00000000 | c)

は10byte. ならmov eac, cしてからmovsx rax, eaxかと思ったけどそんな命令はない. movsx rax, eax(REX + 0F BE/r)だとalを符号拡張してしまう. ...っていいながらマニュアルを読み返すとmovsxd rax, eaxなんてあったのか. 知らなかった. それを使うと5 + 3 = 8byte. 

lea rax, qword [c]

なら1命令で8byte. cをint32_tとして扱い符号拡張する.
嘘書いた. 実行時のeipに依存する. これはその目的には使えない.

11月7日_
「代数学と計算」研究集会参加.
橋本駅で乗り換えるつもりが, ふと見ると相模原駅で「しまった, 乗り過ごした」と慌てて降りる. 反対ホームに行って戻ろうとしたけど橋本に行かない. 「ん?」とよく見ると相原駅で一つ前で降りてたのだった. 相模原 - 橋本 - 相原 早めに出てたので助かった.
11月6日_
加算数回分減らして気持ち速くなった. といってももはや計測誤差の範囲内に入ってる(つまり無駄).
11月5日_
最近ちょいちょいちまちまと以前作った暗号のプロトタイプ実装の改良をしてる. ただかなり時間をかけても0.1%速くなるかどうかという状態が続いていて煮詰まりすぎ. Aranhaさんたちのを抜き返したとは思うがあんまし差はなさそう. うーん, どこか改良点が無いかなあ.
11月3日_
『マスカレード・ホテル』(東野圭吾)
こないだ本屋に行ったときにまとめて単行本を何冊か買ったときの一冊. ちょっと技巧的すぎてがんばったなあという印象を受けるけどまあ悪くない. よくこういう水準でたくさん本を出せるなあ. 著者は一体何人いるのだろうとつい思ってしまう(ほめ言葉). 自分の仕事に誇りを持ってるプロフェッショナルという話では『使命と魂のリミット』が好みかな.
11月2日_

MurMurHash3という128bit値hash関数を調べてた.

いくつかのテストを見る限り, 分布はかなり一様でかつ速い. Jenkinsの96bit値hash関数よりも速い.

Jenkinsのは12byteあたりの処理でadd, sub, xorが21個, rolが6個なのにそれよりも速いのはすごい. MurMurHash3(64bit用)は16byteあたりの処理で, なんと64bit x 64bit→64bit乗算を4回使ってる. その代わりadd, xor, leaは8回ぐらいしかない. 64bit乗算命令が2～3クロックしかかからないこと前提なのか. おもろいな. SHA-1に比べると20倍以上速い.
#32bit版は2倍ちょい遅くなる.

なお, MurMurHash3は暗号用のhash関数ではないので注意. 逆関数は比較的容易に構成できる. なんの役に立つか分からないけどせっかく作ったので書いておく.

最終部分の関数 

uint64_t fmix(uint64_t k)
{
    k ^= k >> 33;
    k *= BIG_CONSTANT(0xff51afd7ed558ccd);
    k ^= k >> 33;
    k *= BIG_CONSTANT(0xc4ceb9fe1a85ec53);
    k ^= k >> 33;
    return k;
}
void final(uint64_t& h1, uint64_t& h2, int len)
{
      h1 ^= len; h2 ^= len;

      h1 += h2;
      h2 += h1;

      h1 = fmix(h1);
      h2 = fmix(h2);

      h1 += h2;
      h2 += h1;
}

の逆関数は 

uint64_t rev_fmix(uint64_t k)
{
    k ^= k >> 33;
    k *= BIG_CONSTANT(0x9cb4b2f8129337db);
    k ^= k >> 33;
    k *= BIG_CONSTANT(5725274745694666757);
    k ^= k >> 33;
    return k;
}
void inv_final(uint64_t& h1, uint64_t& h2, int len)
{
    h2 -= h1;
    h1 -= h2;
    h1 = rev_fmix(h1);
    h2 = rev_fmix(h2);
    h2 -= h1;
    h1 -= h2;
    h1 ^= len;
    h2 ^= len;
}

前回のハッシュ値を(ph1, ph2), 16byteデータを(k1, k2)として(h1, h2)を出力する16byte単位の処理関数をhash_lpとする. 

/*
    input : (h1, h2) := (ph1, ph2), (k1, k2)
    output: (h1, h2)
*/
void hash_lp(uint64_t k1, uint64_t k2, uint64_t& h1, uint64_t& h2)
{
    uint64_t c1 = BIG_CONSTANT(0x87c37b91114253d5);
    uint64_t c2 = BIG_CONSTANT(0x4cf5ad432745937f);
    k1 *= c1; k1  = ROTL64(k1,31); k1 *= c2; h1 ^= k1;

    h1 = ROTL64(h1,27); h1 += h2; h1 = h1*5+0x52dce729;

    k2 *= c2; k2  = ROTL64(k2,33); k2 *= c1; h2 ^= k2;

    h2 = ROTL64(h2,31); h2 += h1; h2 = h2*5+0x38495ab5;
}

(k1, k2), (h1, h2)が与えられたときに(ph1, ph2)を求める逆関数をinv_h_hash_lp, (ph1, ph2), (h1, h2)が与えられたときに(k1, k2)を求める逆関数をinv_k_hash_lpとすると 

void inv_k_hash_lp(uint64_t& k1, uint64_t& k2, uint64_t ph1, uint64_t ph2, uint64_t h1, uint64_t h2)
{
    uint64_t rc1 = BIG_CONSTANT(0xa98409e882ce4d7d);
    uint64_t rc2 = BIG_CONSTANT(0xa81e14edd9de2c7f);

    h2 -= 0x38495ab5;
    h2 *= 0xcccccccccccccccd;
    h2 -= h1;
    h2 = ROTL64(h2,33);
    k2 = h2 ^ ph2;
    k2 *= rc1;
    k2 = ROTL64(k2, 31);
    k2 *= rc2;

    h1 -= 0x52dce729;
    h1 *= 0xcccccccccccccccd;
    h1 -= ph2;
    h1 = ROTL64(h1, 37);
    k1 = h1 ^ ph1;
    k1 *= rc2;
    k1 = ROTL64(k1, 33);
    k1 *= rc1;
}

void inv_h_hash_lp(uint64_t& ph1, uint64_t& ph2, uint64_t k1, uint64_t k2, uint64_t h1, uint64_t h2)
{
    uint64_t c1 = BIG_CONSTANT(0x87c37b91114253d5);
    uint64_t c2 = BIG_CONSTANT(0x4cf5ad432745937f);
    uint64_t rc1 = BIG_CONSTANT(0xa98409e882ce4d7d);
    uint64_t rc2 = BIG_CONSTANT(0xa81e14edd9de2c7f);

    k2 *= c2; k2 = ROTL64(k2, 33); k2 *= c1;

    h2 -= 0x38495ab5;
    h2 *= 0xcccccccccccccccd;
    h2 -= h1;
    h2 = ROTL64(h2,33);
    ph2 = h2 ^ k2;

    k1 *= c1; k1 = ROTL64(k1, 31); k1 *= c2;

    h1 -= 0x52dce729;
    h1 *= 0xcccccccccccccccd;
    h1 -= ph2;
    h1 = ROTL64(h1, 37);
    ph1 = h1 ^ k1;
}

これらの補助関数を使えばseed = 0としたときの16byteのMurmurHash3_x64_128(key, 16, 0, h);の逆関数は 

void inv_murmurHash3_16byte(const uint64_t *h, uint64_t *key)
{
    const int len = 16;
    uint64_t h1 = h[0];
    uint64_t h2 = h[1];

    inv_final(h1, h2, len);
    inv_k_hash_lp(key[0], key[1], 0, 0, h1, h2);
}

とかける. 32byte以上はinv_k_hash_lpとinv_h_hash_lpを組み合わせればよいだろう.

11月1日_
『飲めば都』(北村薫)
彼の作品で読後感がそうだったのは殆どなかったのだけど, 今回はちょっといまいち. 大事な原稿を扱う編集者が記憶を失うまでお酒を飲むかなあ. しかも1回や2回じゃないし. 一人だけじゃないし. そんなことしてたら気づいてないだけで絶対無くしてると思う. 王妃の髪飾りの勘違いも強引すぎると思った. もし, そんなの渡されても他人に見せへんだろう. 駄洒落の応酬やうんちく部分は好きなんやけどね.
昔
[Back]
御意見は光成滋生<herumi@nifty.com>までお願いします