憂鬱な午後のひととき

7月31日_

『幻影の手術室: 天久鷹央の事件カルテ』(知念実希人)

面白かった. その真相は予想外だった.

7月30日_

SVE

ちょっとずつ試してたことをまとめてる.

7月29日_

OpenMPでaggregation

サンプルみながら試してみた. custom.cpp. Goから使うのはmacOSのオプションが違う(-fopenmpじゃない)のでだめだな. 素直にC++内部で書く方がよいか.

7月22日_

bls-eth-go-binary

armeabi-v7aでのクロスコンパイル環境を整備した. MIPSでの動作依頼が来てるので次はそれを調べる. 単純なCのコードは動くけどGoのコードは動かない.

% cat a.go
package main
import "fmt"

func main() {
        fmt.Printf("mips\n")
}
///
env GOOS=linux GOARCH=mipsle go build a.go
% ./a
fatal error: sigaction failed

runtime stack:
runtime.throw(0xe4f51, 0x10)
        /home/shigeo/soft/go/src/runtime/panic.go:1114 +0x60
runtime.sysSigaction.func1()
        /home/shigeo/soft/go/src/runtime/os_linux.go:470 +0x48
runtime.sysSigaction(0x41, 0x7fffef64, 0x0)
...

https://golang.org/src/runtime/os_linux.goのsysSigactionでQEMU対策コードが入ってる. でもここで来てるのはsig = 65なので対策が効いてない. 試しにコメントアウトしたら動いた. cgoを使う簡単なプログラムも動いた.

でもblsをbuildしようとするとCALL16 reloc at 0x253f0 not against global symbolというエラー. fix a CALL16で対策のpull reqが来てるけどそれをやってもエラー. この対策はたまたまうまくいった?

もうちょっと調べないとわからない.

7月21日_

xbyak 5.93

対応完了. 次はインタフェースと実装を分離可能にしたいな.

7月20日_

xbyakのnothrowバージョン

以前から対応したいと思っていたのだけど放置していたらpull reqが来てしまった. まじめに対応する.

7月17日_

久しぶりにC#

.netstandard対応のためにはプロジェクトを作り直す必要があるっぽい. ちょっと間が開くと用語が全然わからない.

7月16日_

exp with fexpa(7)

>うーん, 誤差がw/版の方がよければそちらに移行なんだけどなあ. 係数いじって少しよくなった.

どのレジスタに割り当てるかで簡単に速度が倍変わるとか嫌だ. そのうち機械にさせる方法を考えてみる.

7月14日_

exp with fexpa(6)

ループアンロールしたらちゃんと速くなるようになった(fmath-sve.hpp). fexpa無しバージョンとほぼ互角の速度. レジスタ使用数はunrollなしならw/o版が一つ多い. 2回unrollなら同等. 3回unrollならw/o版が一つ少ない. 誤差はw/o版が少しよい. うーん, 誤差がw/版の方がよければそちらに移行なんだけどなあ.

7月13日_

exp with fexpa(5)

なんとなく分かった気がする.

A→B→C→D→Eという依存関係があるときにBとDで同じレジスタを使うとパイプラインがつまる(?) 違うレジスタを割り当てると速くなる. 全部違うレジスタにするともちろん一番よいけど, そこからどこまで減らせるかだろうか. これは機械でやりたいところ.

7月10日_

exp with fexpa(4)

w/ fexpa版が2倍近く遅い. なんとかならんものか.

7月9日_

exp with fexpa(3)

SVEで実装してみた. 命令数は15命令(+1 movprfx)なのでwithout fexpaの13命令より2命令多い. でも必要なレジスタ数は2個減った. 速度は同じっぽい. ループアンロールしたらwithout版は3倍速くなったけどこちらはどうかな. ところで最初のfmin, fmaxっているのだろうか. なんとなくつけていたけどなくてもちゃんと動くなら削除しよう.

void genExp1SVE(const PReg& p, const ZReg& tz0, const ZReg& tz1, const ZReg& tz2, const ExpParam& para)
{
    fmin(tz0.s, p, para.expMax.s);
    fmax(tz0.s, p, para.expMin.s);
    fmul(tz0.s, tz0.s, para.log2_e.s);
    frintm(tz2.s, p, tz0.s); // floor : float -> float
    fcvtzs(tz1.s, p, tz2.s); // n = float -> int
    fsub(tz2.s, tz0.s, tz2.s); // a
    fadd(tz2.s, tz2.s, para.one.s); // b = 1 + a
    lsr(tz0.s, tz2.s, 17); // bL
    fexpa(tz0.s, tz0.s); // c = fexpa(bL)
    fscale(tz0.s, p, tz1.s); // tz0 *= 2^n
    and_(tz1.d, tz2.d, para.not_mask17.d);
    fsub(tz1.s, tz2.s, tz1.s); // z
    movprfx(tz2.s, p, para.coeff2.s);
    fmad(tz2.s, p, tz1.s, para.log2.s);
    fmad(tz2.s, p, tz1.s, para.one.s);
    fmul(tz0.s, tz0.s, tz2.s);
}
7月8日_

exp with fexpa(2)

y = n + bにしたあと、bの仮数部の上位6bitを使ったあとは、それをマスクした値を引けばよいのか. そうすると2次の項でも十分な精度がでるようだ. あと, 2^z = exp(z log(2))に変換して近似計算と思ったけど, log(2)をかけるのを近似計算に組み込んでしまえばよいのか.

float expC2(float x)
{
    static const float log2 = (float)log(2.0);
    static const float log2_e = float(1 / log2);
    static const float c2 = 0.5 * log2 * log2;
    float y = x * log2_e;
    int n = (int)floor(y);
    float a = y - n; // 0 <= a < 1
    float b = 1 + a; // 1 <= b < 2, y = (n-1) + b
    uint32_t bu = f2u(b);


    float bL = u2f(bu >> 17);
    float z = b - u2f(bu & ~FexpaTbl::mask(17));
    /*
        split b into bL and z where bL is for fexpa and z is remain
    */
    float c = fexpaEmu(bL);
    float d = 1 + (log2 + z * c2) * z;
    return powf(2, (float)n) * c * d;
}
7月7日_

exp with fexpa

y = log_2(e) xとしてfexpaを使うにはy = n + b, n∈Z, 1<=b<2にしないといけない. テーブル引きしたあと残りをどうやって計算するのかな. 下記のようだと3次の項まで必要だから余計な計算コストを考えるとfexpaを使わないよりメリットがなさそう.

e^x = 2^(log_2(e) x)
y = log_2(e) x
z = y - 1
n = floor(z)
a = z - n ; 0 <= a < 1
b = 1 + a ; 1 <= b < 2, y = z + 1 = (n + a) + 1 = n + b
b1 = (b & (mask(32) << 17)
b2 = (b & ~(mask(6) << 17)) - 1 ; b = b1 + b2
d = fexpa[b1] ; d = 2^b1
f = b2 * loge_2
g = 1 + f + f^2/2 + f^3/6

e^x = 2^y = 2^(n + b1 + b2)
= (2^n) (2^b1) (2^b2)
= (2^n) (2^b1) e^f
7月6日_

fexpa

SVEのfexpaの仕様がややこしい. 入力値の下6bit x=[0:5]に対してexp(x/64)の下位23bitと元の[31:23]のorになる. 下位23bitはともかく残りはどう使うのだろう.

struct FexpaTbl {
    static const size_t N = 64;
    uint32_t tbl[N];
    union fi {
        float f;
        uint32_t i;
    };

    static inline uint32_t mask(int n)
    {
        if (n == 32) return uint32_t(-1);
        return (1u << n) - 1;
    }
    FexpaTbl()
    {
        for (size_t i = 0; i < N; i++) {
            fi fi;
            fi.f = (float)pow(2.0, i / float(N));
            tbl[i] = fi.i & mask(23);
        }
    }
    float calc(float f) const
    {
        fi fi;
        fi.f = f;
        uint32_t v = fi.i;
        fi.i = ((v >> 6) & mask(8)) << 23;
        fi.i |= tbl[v & 63];
        return fi.f;
    }
};

float fexpaEmu(float x)
{
    static FexpaTbl tbl;
    return tbl.calc(x);
}
7月3日_

mcl/blsのコード整理

Ethereumの古いハッシュ関数の実装を削除. すっきりした.

7月2日_

有給

疲れがたまってるので休みにした. いろいろやることはあるので普段としてることはあまり変わらない気もするが気分的にまし(なのか?).

Microsoft MVP受賞. いろいろありがたいです.

7月1日_

xbyakのAMXサポート

Intelからのパッチのリファクタリングと細かいエラーチェック追加. 実用上は使われないであろう命令エンコーディングは未対応(なのでエラーにした). エミュレータが公開されたらちゃんと対応する.


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