melancholic afternoon

SPARC64 SIMD版expf

初めてSPARC64のintrinsicを使ってみる. double x 2はx86/x64と同じ型__m128dのaliasがあるようだ. やりやすい. ただfloat x 4の__m128は無かった. そもそもそういう型が無い?

あと, floatとdoubleの変換関数も見当たらなかった. それが足かせになるのもつまらないのでexp4d(double py[4], double px[4]);というインタフェースで中身はfloatの精度でするものにした. それだけでも30%ほどの速度アップ(twitter).

浮動小数点演算のレイテンシは7～8ぐらいだそうでexp8dを作ったら更に1.8倍ぐらいになった(new_exp8d). まあこのあたりで一段落.

Re: SPARC用fast exp

係数をいじってちょっと精度をよくする. addも1回減らす. 今回はexp(x)を区間[-1/8, 1/8]で5次の多項式で近似してみる. Mapleでやるとこんな感じ.

f:=x->1+A*x+B*x^2+C*x^3+D*x^4+E*x^5;
g:=int((f(x)-exp(x))^2,x=-c..c);
sols:=solve({diff(g,A)=0,diff(g,B)=0,diff(g,C)=0
  ,diff(g,D)=0,diff(g,E)=0},{A,B,C,D,E});
Digits:=10000;
s:=eval(sols,c=1/8);
evalf(s,300);

sollyaだと

> remez(exp(x),5,[-1/8,1/8]);

便利だねえ. ただsollyaの係数だとMapleで求めたのより[-30, 30]での平均二乗誤差が微妙に大きくなった. 私の評価が適当なせいかもしれない. 誤差の話って大昔にほんのちょっとだけ勉強しただけだからちゃんとしたほうがよいかも.

4並列版は自動的にSIMD化されるのかと思っていたけどそうでもない???. intrinsic使ってみるか.

早起き

この1年半ぐらい, 朝4～5時に起きて暗号のpdfを書き続けていたら, 起きたくなくても目が覚めて起きてしまうようになった. もっと寝てたいんだけどなあ.

はっ, これが話題の老人というやつか.

SPARC用fast exp

@syoyoさんが私のfmathを Parallela Epiphanyに移植してくれたら結構速かったらしい. 私はSPARCで速くならないか試してみた.

入力範囲が[-30, 30]で精度が1e-6程度でもよいそうなので, テーブルを使わないバージョンにしてみる. 古典的だがまず
exp(x) = exp(x/2)² = exp(x/4)⁴ = ... = exp(x/128)¹²⁸
という関係式を使って|x| < 1/4にもっていき, それからは単にテイラー展開.

ただやってみるとどうも期待通りの精度がでない. よく考えたらfloatでやると桁落ちするのだった. それで中間変数はdoubleにしたらそこそこになった(expf_test.cpp).

float new_exp(float x)
{
    double t = x / 128;
    // |t| < 1/4
    const double c2 = 1 / 2.0;
    const double c3 = 1 / 3.0;
    const double c4 = 1 / 4.0;
    const double c5 = 1 / 5.0;
    const double c6 = 1 / 6.0;
    double y = 1 + t * (1 + c2 * t * (1 + c3 * t
        * (1 + c4 * t * (1 + t * (c5 + (c5 * c6) * t)))));
    y *= y;
    y *= y;
    y *= y;
    y *= y;
    y *= y;
    y *= y;
    y *= y;
    return (float)y;
}

単純だが手元のSPARCではfast_expの1.7倍で精度もそれよりちょっとよいぐらい. 並列実行もしやすいだろうしこの戦略はよさそうだ. ちなみにx86では逆に遅くなる.

Parallella epiphanyでも遅いそうだ.

『クラウドを支えるこれからの暗号技術』(秀和システム)

本日発売です(amazon). 最近の暗号理論の話題を紹介する入門書です. 関わってくださった方々, 本当にありがとうございました.

I部では駆け足で共通鍵暗号や認証などを紹介します. 不慣れな人は手元に『新版暗号技術入門秘密の国のアリス』(結城浩)があるとよいでしょう. ただそれと同じことをなぞってもつまらないので, POODLE, パディングオラクル攻撃, AEAD, PFS, ビットコイン, ChaCha20など最近の用語を登場させてみました.

RSA暗号は安全じゃないよ, 強秘匿性や頑強性がないよ, IND-CCA2最強,
Webのセッション認証にHash(パスワード + メッセージ)を使ったら駄目だよ. SHA-2でもアウト. SHA-3ならセーフ(cf. 伸長攻撃),
Mersenne Twisterを暗号の乱数に使ったら駄目だよ(cf. Mersenne Twisterの次に出す値を推測する),

などの話題もあります. 入門書ではほとんど触れられていない楕円曲線暗号もしっかり説明します.

II部では高機能暗号とか, 組織暗号と呼ばれる再暗号化, 検索可能暗号, 属性ベース暗号, 関数型暗号, 準同型暗号などを紹介します. 本書のメインです. 「こんな暗号があるんだ」と興味を持ってもらえればありがたいです. このあたりは既存の暗号の専門書でもほとんど紹介されていないのでサーベイとして価値があるんじゃないかと思います(網羅的ではありません).

III部では楕円曲線というと天下りに登場するWeierstrassの方程式y² = x³ + ax + bがなぜそれなの?という疑問に答えます. もちろん厳密ではないですが, 穴一つの浮輪の表面上の関数を考えると自然にその式が登場することが伝わればありがたいです. 本当はRiemann面に言及したかったですが.

もともとは高校生向けの講義が発端です. 高校生や専門を決める前の大学生が手にとってみて, 細部は分からなくてもこのジャンルが面白そうだなあと思ってもらえるとありがたいです.

ftrimaddd, ftrismuld, ftrisseld

社内受託開発(?)が一段落してエアポケット的に空きが出来たので久しぶりに調査. 18日のintrinsicではなく, asm出力のコードで探したら見つかった.

どの象限にいるかを選択するためだから2bitだけしか見ないのか. 面白い. こういうトリッキーな命令はIntelには無いなあ(pcmpistriはまた置いといて). 命令数削減の手法は特許も取っているようだ. expの計算に利用できるかあとで試す.

majorとminor

次のコードがgcc/clangでコンパイルエラーになる. -pedanticや-ansiは関係ない.

#include <stdlib.h>

struct A {
    int major;
    A() : major(0) {}
};

int main()
{
}

ちょっとびっくりだがsys/sysmacros.hにmajor, minor, makedevが定義されているから. 大昔, 一度はまったことがあるようなないような.

VCのmin, maxマクロを笑えない(あちらは-DNOMINMAXで回避できる).

MELT upフォーラム参加.

資料はこちらの案内. 久しぶりに笠原先生にお会いする.

宮地先生の国際標準化動向の話. 鍵共有のプロトコルを提案したら, 何度もあるチームからrejectが来るという. そこは形式的証明の手法を使って不備を見つけてるとか. コンピュータが間違いを見つけて人間同士がやりとりしてるというのがちょっと面白かった.

何人かで辻井先生の部屋にお邪魔してエニグマ改良機(NEMA)を見せていただいた(機械式暗号機の傑作～エニグマ登場～).

実際に電源を入れてぽちぽち(カッチン、カッチンという感じ)と触らせていただいた. 歯車が回って, ランプが光って暗号文が一文字ずつ表示される. こんなのを使ってたんだなあ. どうもありがとうございます.

見本

が届いた. 形になるとやっぱりうれしいですね. 関係者の方々, どうもありがとうございました. ただちょっと気持ち, 字が小さいような. もうちょっと大きかった方が読みやすかったかも. ページ数が増えるとコストがかかるので難しいようですが.

_fjsp_trimadd_v2r8

_fjsp_trismul_v2r8, _fjsp_trissel_v2r8のマニュアルを見てもどういう命令なのかいまいち分からない. trimadd, trismul, trisselという命令が発行されるのかと思ったけどasmソースを見てもそれらしいのが見当たらない. あとでもうちょっと調べる.

Effective Modern C++勉強会#6

発表者の方, ありがとうございました. 次回は7/15です. 会場は16階ですのでご注意ください. 飯田橋でするのは最後だ.

Xbyakにmovbe追加

命令があるのは発表されたときから知っていたけど, 使わないかなと思ってほっていたらpull reqが来たので対応. ちなみにcmovbeやfcmovbeは昔からあって一瞬混乱する(意味は全然違う). swapping bytesなんだからmovsbにすればいいのにと思ったがそちらも既にあるのだった. beってbig endianの意味か.

Re:pairingライブラリの不具合(2)

しかしスタックオーバーフローなんて初歩的なミス, なんでしたんだろうとコードを再確認する. 実際, ちゃんと境界チェックはしているのに. しかもclangとgccでaddress santizerの出す場所が違うのが気になる.

で, よくみたらエラーになっても文句は言えなかった. Intという固定長Lの配列を持つ整数クラスの乗算時にテンポラリをmallocするのがいやで必要な長さをstackに確保して配置newでIntとして使うことをしていたのだった.

Int *t   (クラスのサイズ)
  ｜＿＿＿＿｜
  ｜＿スタック上に領域確保＿｜
  ｜＿→ぶち破る→＿＿＿＿＿｜

そりゃあIntが持つ本来のサイズを超えてアクセスするのでaddress sanitizerのチェックに引っかかる(領域はあるから自分としたはOKのつもり).

さぼらずにちゃんとしたら問題はなくなった. 数値が小さいとちゃんとエラーも表示される.

Re:pairingライブラリの不具合

修正した(fix stack overflow). 実際の所, つまらないバグなのだがこの数値が意図したものじゃなくなってるので本質的に危ない(多分どこかにtypoレベルのミスがある).

当初はWindowsでGMPが無い環境でも動くように汎用的な多倍長クラスを用意しようと思っていたのだが, 実際に使うのは固定長クラスなのでこちらのクラスはとりあえずという感じの性能, 質になってしまってる. 今更そこに力を入れるのもなんなので素直にmpz_classに置き換えていくか.

当時はMPIRが無かったので仕方なかった.

pairingライブラリの不具合

なんかUbuntu 15.04 + g++-4.9.2で落ちるらしい(configurations.txt). clangだと大丈夫. static変数周りで非標準な使い方してたところで罠に落ちたかなあ. 環境作って確かめるか.

C++

以前, tag dispatchを使って書いていた部分を普通に書き下すようにした. 無理に汎用的に書こうとして読みにくくなるデメリットの方が大きいと思ったから. 『Effective Modern C++』のItem 27でも~~tag dispatch~~SFINAE(thanks to @egtra, @kariya_mitsuru)の手法が紹介されているが, 正直そこまでして何がうれしいのかと思ってしまう.

#include <stdio.h>
#include <string>

struct A {
    std::string a;
    template<class T, class = typename std::enable_if<
        !std::is_same<A, typename std::decay<T>::type
        >::value>::type>
    explicit A(T&& a) : a(std::forward<T>(a)) {}
};

struct B {
    std::string a;
    explicit B(const std::string& a) : a(a) {}
    explicit B(std::string&& a) : a(std::move(a)) {}
    explicit B(const char *a) : a(a) {} /* 文字列リテラルのとき一時オブジェクトを作らせないため */
};

int main()
{
    {
        A x("abc");
        A y(x);
    }
    {
        B x("abc");
        B y(x);
    }
}

いや, もちろんstruct Aの方がコードを一つしか書かなくてよいというのはあるし, よくそんな回答を思いついたというのはある. 汎用ライブラリの中では頑張るしかないときもあるのだろう. ただBの方が引数の型が何か明示されてるので分かりやすい. 間違えて異なる型を渡したときのエラーも分かりやすい.

昼休みにC++はなんであんなに複雑怪奇な文法になったんだろうと雑談していて, 「とても頭のいい人がよってたかって考えたから(駄目なんだ)」という結論になった(本気にしないように).

一息

とりあえず一つ作業が終わった. 次のことをしよう. expの計算も再開したいな.

ソースコードで検索するときはackを使っていたのだけどagの方が速くていいよと教えてもらった.

sudo apt-get install silversearcher-ag
alias ag "ag --pager 'less -R'"

として試してみる.

PHP Bug #69619

先々月, 少しPHPの勉強を始めたところでPHPのコードにバグらしきところを見つけたので報告した. memcpyの引数のポインタと実体を取り違えるという初歩的なミス. コメント欄によると, そのパスに入って実際に困るコードを示せ, そうじゃないとclose(正確にはNo Feedback)するとか.

コード上明らかにおかしいのを指摘しただけで, もちろんそれを修正するかどうかは向こうの判断だろうけど, 再現コードをこちら待ちというのはどうなんだろう. PHPにそこまでコミットする気はないので調べないけど. 重要じゃないからという判断なら分からないでも無いが.

あれ, 今見たらだれかOpenし直している.

将棋と数独とpext

数独のパターンマッチを大量にしている人がいたのでコンピュータ将棋における Magic Bitboard の提案と実装が何かヒントになるかなと読んでみた.

チェスでは駒の効き目のビットパターンにマッチしてるかどうかを出来るだけ高速に判定したい. 局面を表すデータ(8x8bit)にマジックナンバーをかけて小さい値にデータを押し込む手法がある(magic bitboard). それを将棋の9x9にも応用しようというもの.

チェスなら64bitの一つのレジスタに入るけど将棋や数独は9x9なので少なくとも二つの64bitレジスタが必要. 二つのレジスタへの持ち方をうまくして似たようなことができたよという話.

以前128bitレジスタに入れて実験したこともあるけどそのときは整数とSIMDレジスタの変換コストを考えると整数のままでもよさげだった. ところで今ならHaswell以降にあるpextを使うのが一番なんだろうな(cf. フラグを愛でるの23ページ目).

本当は試したいところなんだが, こういうのは実際に使われているコードで試さないと効果の程が分からない. 加えてデータ構造に直結するので大規模な書き換えが必要. そして一番問題なのが速くしたいマシンはHaswellじゃない(w.

ポンピング補題が正規言語の必要条件であること

@sinya8282さんに反復補題の不十分性を教えていただいた.

Lを回文の文字列最後に何か文字列vをつけた集合とする. w ∈ Lはw=uu^Rvの形. u^Rはuのビット列を反転させたもの. uの先頭ビットをu₀, 残りのビットをu'とする. u=u₀.u', v=v₀.v'. |w|>=4とする. w=xyzと分解する.

n=1なら|v|>=2でx=uu^R, y=v₀, z=v'でxyⁱz =uu^Rと何かの文字列 ∈ L.

n>=2ならx={}, y=u₀, z=u'u^Rvとする.

xz=u' u'^Ru₀ vでOK.
i>=2ならxyⁱz=u₀.u₀と残りでOK.

Lが正規でないことは, Myhill-Nerodeの定理を括弧の対応が正規言語でないのと同じように右同値類が無限個あることを示すのだろうな. なるほど.

『正規表現技術入門』(新屋良磨, 鈴木勇介, 高田謙)

前回の続き. 今回は残りのAppendixを読んだ. 注釈多すぎ….

ポンピング補題, 最初読んだときは「おまえは何を言ってるんだ」だったけど, 何度か読んでると雰囲気は分かった. Myhill-Nerodeの定理はきれいだね, 証明読んでみたいね(でも今は気力がないけど).

再帰や後方参照が非正規であることの証明をそれぞれ2通り書いてあるのがすばらしい. 比較することでよく分かる. 後方参照があると合成数個ならぶ文字列ができるので駄目という論法がおもしろい. Myhill-Nerode定理による後方参照が非正規であることの証明はまだ分かってない.

欲を言えば, ポンピング補題の仮定を満たすけど正規言語じゃない例が載ってるとよかった.

*を含まない正規表現に書き換えるときにΦを使うのはよいけどnot(.*)なのだから間接的に*を使ってるのでは. 直感的にはnot(.*0)が「0で終わる」をチェックすればよいだけなので効率がよいのは分かるけど.

syntactic monoidでx^N!=u^N!+1を示すのはどれぐらいのコストなんだろう. オートマトンの状態数がnならNは最大nⁿぐらいいくのか? 結構でかい?

以下はeditorial error.

p.291 補題A.1の「次の条件を(2)満たす」→「次の条件(2)を満たす」
p.293 「対象律」→「対称律」
p.294 [ε]の左の括弧の位置がおかしい.
p.298 一つ目のw=((...の下の説明「l個の閉じ括弧」→「l個の開き括弧」
p.298 「つまりxyⁿz(n ≠ 0) は」→「つまりxyⁿz(n ≠ 1) は」
p.313 図A.4の「(1|01)*」の縦棒が無い.

ささいな表記の揺れ.

「一つ」と「1つ」
「変え得る」と「変わり得る」

憂鬱な午後のひととき