憂鬱な午後のひととき

11月30日_
ペアリングフォーラム参加
毎度秘密鍵を渡すなら対称鍵暗号でいいんじゃない(秘密鍵 = hash(マスター秘密鍵 | ID | from | to))ってのは私の質問が悪かったか意図が伝わらなかったようで, 素のタイムリリース暗号のよさをアピールされてしまった. 5年前に実装したしそれは分かってるんだけどなあ.
そのあと会社に戻って22時過ぎまでお仕事(おお珍しい:初めてじゃ).
11月27日_
夜中4時頃目が冴えてしかたなかったのでXbyakで(ある程度)FPUサポートした. これで32bit用のexp()も作れる.
11月26日_
社内でテスト駆動型ペアプログラミングをしてみた. 但しPHPやけど.
11月24日_
少し前だけどゴードン・ベル賞に浜田・長崎大助教ら[NIKKEI] GPUクラスタによる高性能計算技術の実証[理化学研究所]
似鳥さんすごいな. おめでとうございます.
11月23日_
Make:Tokyo Meeting 04に行ってきた.
今回は暗いところでピカピカ光るデモが多くてちょっとしんどかった. 1年半前に行ったときおもしろかった手作りプラネタリウムは人が時間が合わなくて見られなかった. 残念. 次世代自転車を考える会のに試乗したがめちゃ不安定で怖かった. 座ったまま休憩できるとかあればいいのに.
11月18日_
float xに関するexp(x)の実装を考える(6)
小数部のテーブルを3bit分整数部にもっていくことでテーブルサイズを小さくするためにはもとの 
 int n = (int)floor(x) - 1;
 float d = x - n;
する直前でxの値を8倍すればよい. 対応するd0, d1のビットも相応にずらす. expを求めるところは8分の1にする必要があるので 
 e = exp(n - 1.0) * exp(d0 + 1) * exp(d1 - 1);
から 
 e = exp((n - 1.0) / 8) * exp((d0 + 1) / 8) * exp((d1 - 1) / 8);
となる. 近似計算は上はexp(d1 - 1) = 1 + (d1 - 1) = d1とできていたが, 下はexp((d1 - 1) / 8) = 1 + (d1 - 1) / 8 = (d1 + 7) / 8となる. この計算を更に簡略化してみよう. d1は1にすごく近い値, 厳密には d1 = 1 + α/224. α <= 214である. 代入すると(d1 + 7) / 8 = 1 + α/227. つまりd1をつくるところで 
  fi.i =0x3f800000 | α // 0x3f800000は1.0fのバイナリ表現
  d1 = fi.f;
としていたところを 
  fi.i =0x3f800000 | (α >> 3);
  d1 = fi.f;
とすれば(d1 + 7) / 8を計算できる. FPUの足し算1回とFPUの掛け算1回が整数のbitシフト1回に減るわけだ. トータルコストとしては最初の掛け算1回と今回のシフト1回のコストが増えるがテーブルサイズが半分になる(あと情報落ちに関する誤差が若干増える). 場合によってはこちらを採用してもよいだろう. ちなみに 
    for (float x = 0; x <= 20; x+= 1e-5f) {
        sum += exp(x);
    }
を計算してみると, 今までの高速化ではオリジナルexp()による演算結果と厳密一致(!)してたが, テーブルサイズ半分の場合は相対誤差4e-7ほど出た(まあ問題ないレベルだと思うが).
xを0から2まで1e-6ずつ増やしながらexp()の和を計算するベンチマーク(詳細はexp.cpp参照). テーブル半分じゃないバージョンも前回より多少ブラッシュアップした.
Core2duo 2.6GHz
32ibt@VC10 64bit@VC10
org exp 0.828 0.250
my exp 0.219 0.109
my exp(テーブル半分) 0.359 0.156
Opteron 2376 2.3GHz
32ibt@gcc 4.4.1 64bit@gcc 4.4.1
org exp 1.22 4.65
my exp 0.50 0.43
my exp(テーブル半分) 0.560 0.48
32bitVCだとオリジナルの3.5～4倍, 64bitVCで2～3倍, 32bit gccで2～2.5倍, 64bit gccで10倍(?)速い. やっぱりgccのoriginal expは遅い. expを多用する数値計算してる人はgccをVCにするだけで幸せになれるかもしれない.
11月17日_
float xに関するexp(x)の実装を考える(5)
少し間が空いてしまった. 
// __m128 d;
float e;
*(int*)&e = _mm_cvtsi128_si32(_mm_castps_si128(d));

>最後の2行はx86ではfloatはFPUレジスタに入れて返すのでこうするしかない.
と書いていたけどヘッダを見ていて_mm_cvtss_f32という関数があるのに気がついた. これを使えば 
float e = _mm_cvtss_f32(d);
とかける. 32bit環境ではスタックを経由してFPUレジスタにコピーされる. x64環境ではFPUがXMMレジスタを使うので何の命令も生成せずに(見かけ上)変換できる. すばらしい.
ところで下記バージョンではテーブルサイズが整数用に7bit, 小数用に12bit必要だった. 整数と小数は本来つながっているので3bit分の小数部分を整数に持っていけば 整数部7 + 3 = 10bit, 小数部12 - 3 = 9bitとなり, テーブルサイズが4224 dwordから1536 dwordと半分以下になる. そのためにはxを最初に*8して補正すればよい. 詳細は次回.
11月12日_
float xに関するexp(x)の実装を考える(4)
帰りの電車の中でexp()のC版をそのままSISD(not SIMD)を使って実装してみた(exp2.cpp).
core2duo 1.8GHz
VC9(32bit) VC10(32bit) VC10(64bit) gcc 4.3.4(32bit)
org exp 1.625 1.563 0.516 2.72
my C exp 1.484 1.031 0.500 0.75
my SISD exp 0.734 0.359 0.250 0.20
なかなかいい感じ(gcc 4.1.2だとちゃんと動かなかった). VC9からVC10への最適化性能の向上も凄いものがあるな. 次はもうちょいテーブルサイズを小さくできないか考える. SSE4ならもっと簡単にできるんだけどな. 一応コードの説明. 
float Exp2(float x)
{
static const float one = 1.0f;
static const unsigned int mask2 = (mask(8) << 23) | mask(11);

// float xの値をx1に代入
__m128 x1 = _mm_load_ss(&x);
// floor(x1)を求めて1を引く
int n = _mm_cvttss_si32(x1) - 1;
ここで負の無限大への切り捨てが使えれば負の場合を場合分けせずにできるはず. それはSSE4のroundssだと可能だけど, SSSE3までだとSSEの浮動小数制御レジスタを触る必要がありコストが高い. 大量のfloat配列を一度にexpする場合はまだしも1回ごとに触るのは微妙. しかたないのでとりあえず最初に場合分けしてる. 
__m128 d = _mm_setzero_ps();
d = _mm_sub_ss(x1, _mm_cvtsi32_ss(d, n));
整数にしたnを再度浮動小数に戻してx1から引く. _mm_setzero_ps()は無くてもよいが出力コードを見る限りあった方がよいコードが生成された. 依存関係を断ち切れるからかな. 
__m128i di = _mm_castps_si128(d);
unsigned int d0i = (unsigned int)_mm_cvtsi128_si32(di) >> 11 & mask(12);
__m128と__m128iは異なる型なのでキャストが必要. _mm_castps_si128を使って行う. これはあくまでCで必要なだけでSSEとして対応する命令があるわけではない. 
d = _mm_and_ps(d, _mm_load_ss((float*)&mask2));
d = _mm_mul_ss(d, _mm_load_ss(&expTbl0[n + 1]));
d = _mm_mul_ss(d, _mm_load_ss(&expTbl1[d0i]));
float e;
*(int*)&e = _mm_cvtsi128_si32(_mm_castps_si128(d));
テーブル引きして掛け算する. 最後の2行はx86ではfloatはFPUレジスタに入れて返すのでこうするしかない. x64だとSSEレジスタをそのまま返せるので余計なキャストやメモリアクセスは不要.
11月11日_
・カップヌードルが98円で安いと買って帰ってよくみたらスープヌードルだった. 二度目な気がする.
・三井住友VISAカードのポイント交換ページの景品選択方法が酷い. たとえば2130ポイントあったら2000ポイント以上で絞り込みたいのが普通だと思うのに, 選択ダイアログは2001～3000というものしかない. 2000ポイント丁度の景品を見るためだけに1001～2000を選択しないといけない(しかもそれだと最後のページに置かれるので何度もクリックしないといけない). 設計した人絶対自分で試してないだろ. そもそもたいした景品数が無いのだから視認性のよい一覧ページを作ればよいのに. 昔はそうだった気がする.
11月10日_
float xに関するexp(x)の実装を考える(3)
負なら最後に逆数を計算するのでそのためのフラグを保持する. 
float Exp(float x)
{
    bool isPositive;
    if (x >= 0) {
        isPositive = true;
    } else {
        x = -x;
        isPositive = true;
    }
x を 整数部 n と d ∈ [1, 2) に分ける． 
    int n = (int)floor(x) - 1;
    float d = x - n;
浮動小数演のビットパターンを取り出す. *(int*)&dとかするとstrict-aliasing ruleを破るのでunionを使う. 
    union fi {
        float f;
        unsigned int i;
    };
    fi fi;
    fi.f = d;
    unsigned int di = fi.i; // これで di は d のビットパターンになる
dを二つの小数d0, d1に分ける. それぞれの仮数は12bitと11bitで論理和を取るとdの仮数となるようにする. mask(n)はn bitのマスクを作る関数. d0は上位12bit部分. 
    fi.i = di & ((mask(8) << 23) | (mask(12) << 11));
    float d0 = fi.f;
d1は下位11bit分. 
    fi.i = di & ((mask(8) << 23) | mask(11));
    float d1 = fi.f;
ここまででx = n + d0 + d1 - 1となってる(下の(***)の説明の部分相当). よって 
    float e = exp(n) * exp(d0) * exp(d1) / exp(1.0);
で計算可能となった(この式変形は情報落ちを除いて厳密). 最後にもともと負だったら逆数を求める. 
    if (!isPositive) {
        e = 1 / e;
    }
    return e;
}
さて高速化. d1 - 1が0に近いことを利用としてexp(d1 - 1) = 1 + (d1 - 1) = d1とする. nとd0に関するexpはそれぞれ7bit, 12bitのテーブルを用意する. 
static float expTbl0[128 + 1];
const size_t d0MaskLen = 12;
static float expTbl1[1 << d0MaskLen];
#define NUM_OF_ARRAY(x) (sizeof(x) / sizeof(*x))
void InitExp()
{
    for (int i = 0; i < NUM_OF_ARRAY(expTbl0); i++) {
        expTbl0[i] = exp((float)i - 1);
    }
    for (int i = 0; i < (1U << d0MaskLen); i++) {
        fi fi;
        fi.i = (127 << 23) | (i << 11);
        expTbl1[i] = exp(fi.f);
    }
}
d0に関するテーブル引きのための指数(d0i)取り出しはもともとのdiを使えばよい. 
    int d0i = (di >> 11) & mask(12);
これで 
    float e = expTbl0[n + 1] * expTbl1[d0i] * d1;
という計算に置き換えられた. 実際に試してみる. 0～2までのxに関して最大誤差と平均誤差と処理時間@CoreDuo 1.66GHz(exp.cpp). 
// VC9(32bit)
max=4.768372e-007, ave=8.847215e-008
sum=19910788.588494, org exp time=2.000000
sum=19910788.038545, my  exp time=1.453000

// gcc 4.3.4(32bit)
max=4.768372e-07, ave=9.258090e-08
sum=19910788.588641, org exp time=3.150000
sum=19910788.034149, my  exp time=1.490000
おお, この時点で速くなってる(logもそうだがVCの超越関数速いなー). floor()使うなら負数の場合もまとめて(整数部のテーブルを1bit増やして)テーブル引きしたほうがいいかもしれんな. あと, 下に凸なexpを下から追いかけるのでどうしても値が小さい方にいってしまう. 何かうまい方法ないかなあ.
11月9日_
float xに関するexp(x)の実装を考える(2)
最後の部分を図で説明してみる. dが1以上2未満のとき, d = sign ・ 2指数 ・ (1 + 仮数 / 224) と表すとsign = 0, 指数 = 0と固定され, 動くのは仮数fのみとなる. よって, たとえばfを8bit, 8bit, 7bitの論理和に分けてみると 
f = f1 | f2 | f3
f1 = f & 0b11111111000000000000000
f2 = f & 0b00000000111111110000000
f3 = f & 0b00000000000000011111111
d = (1 + f / 224) = (1 + f1 / 224) + (1 + f2 / 224) + (1 + f3 / 224) - 2 とできる. f1, f2, f3に相当する浮動小数をd1, d2, d3とすると, x = n + d1 + d2 + d3 - 2.(***) よってex = en ed1 ed2 ed3 / e2. だいぶいい感じ. 次にどういうビット幅をとればよいかを考える.
11月8日_
float xに関するexp(x)の実装を考える(1)
11月7日_
昔作ってたlogの計算をintrinsicでやってみた(log.cpp). とりあえずasmを焼き直しただけなので改良の余地は多々あるだろうがCore2Duo 1.8GHzでベンチ.
32bit VC9 32bit VC9 /fp:fast 64bit VC 32bit gcc 4.1.2 64bit gcc 4.1.2
log 0.687 0.375 0.188 0.48 0.330
my_log 0.156 0.125 0.130 0.110 0.100
計測がいいかげんなので有効数字2桁なさそうやけど32bitで概ね3～4倍. 64bitのVCのlogはがんばってるなあ. 昔はdivpsよりもrcpps+ニュートン補完が速かったんやけど今は違うのか?(要検証) まあ, 昔の最適化が現在でもそれなりに意味はあると分かったので次はexpを考えたい.
11月6日_
from_time_tの問題は1.41.0で対応するっぽい. と思ったらまだ駄目だった.
11月5日_
・第二回カーネル/VM探検隊
ほどほどに面白かった. ただ16時～20時半の長丁場は辛かった. 途中でおなか空きすぎて頭が半分回ってなかった. MINIX SMPのところで66h, 67hをnopに置き換えてる理由がわからないと言ってはったので32bitモードで16bitコード生成したいからとちゃうって適当に答えてしまったが, ソースのコメントに書いてあった.
・副都心線初めて乗った. 渋谷で上に出ようとして階段使ったら100段あった. なんて深いところにあるんや.
11月4日_
・CANS2009で発表するF21223, F3509上のTateペアリングのソースコードが公開された. 私は途中参加だったのでC++ならもちっと楽できるんだけどなあと思いつつ郷に従う. なにはともあれ, こういう最速狙うコードで公開されてるのは殆ど無いのでめでたい.
・下記のposix_time::from_time_tの件, 2007年に既に指摘されてるのに直ってないということは直す気ないんやろか. trackが見えんでよう分からん. こっちの人もfrom_time_tについて何か言ってるなあ. これも2年前だ.
11月2日_
boost::posix_time::from_time_t(time_t)は64bit用time_tに対応してないぽい(boost 1.40). VC9@32bit, gcc 4.1@64bitで確認. 
#include <stdio.h>
#include <time.h>
#include <iostream>
#include <boost/date_time/posix_time/posix_time.hpp>
#include <boost/date_time/gregorian/gregorian.hpp>

void Put(time_t t)
{
    struct tm tm;
    char buf[64];
#ifdef _WIN32
    _gmtime64_s(&tm, &t);
    asctime_s(buf, &tm);
#else
    gmtime_r(&t, &tm);
    asctime_r(&tm, buf);
#endif
    printf("time=%s", buf);
}

int main()
{
    time_t t;
    boost::posix_time::ptime bt;

    time(&t);
    printf("sizeof=%d\n", sizeof(t));
    Put(t);
    bt = boost::posix_time::from_time_t(t);
    std::cout << "boost:" << bt << std::endl;

    std::cout << "boost:" << (bt + boost::gregorian::years(29)) << std::endl;

    t += 84600 * 365 * 29;

    Put(t);
    bt = boost::posix_time::from_time_t(t);
    std::cout << "boost:" << bt << std::endl;

}

sizeof=8
time=Mon Nov 02 13:39:33 2009
boost:2009-Nov-02 13:39:33
boost:2038-Nov-02 13:39:33
time=Sat Mar 20 01:09:33 2038
boost:1902-Feb-11 18:41:17 ; おかしい
sizeof(time_t) = 4な環境で駄目なのはしかたないけど, そうじゃない環境では対応しておいて欲しいなあ. 言ってみるか. 
  ptime from_time_t(std::time_t t)
  {
    ptime start(gregorian::date(1970,1,1));
    return start + seconds(static_cast<long>(t));
  }
longにキャストしてるからか? でも64-bit Linuxでも駄目なんだよなあ.
昔
[Back]
御意見は光成滋生<herumi@nifty.com>までお願いします