2020-09-19

パソコンの修理がおわった

ここ一か月ほどパソコンが壊れていて修理していたのですが、修理が終わって戻ってきたのでいろいろやっていきます。

2020-09-12

バンド端付近にフェルミエネルギーが来る場合、そこの状態密度が高い構造が選ばれる、について

授業で知ったのですが特に証明がなかったので考えました。

仮定

構造によらずバンドの上端・下端は同じ
構造によらずバンドに入る電子数も同じ（ $N = \int f(E) \, dE$ が一定）
構造によらず電子がバンドに全部入った時のエネルギーも同じ（が一定）
- 状態密度として上下対称なものしか考えない場合、常にこれが満たされます

フェルミエネルギーがバンド下端付近の場合

少数の電子がバンドに入ることになります。

その場合、バンド下端付近に大きな状態密度を持つ構造のほうが、エネルギーの低い状態にたくさんの電子を詰め込めることになります。

すると全エネルギーが低くなるので、そういった構造のほうが安定で選ばれやすいということになります。

フェルミエネルギーがバンド上端付近の場合

一定のエネルギー準位以下の状態に入れられる電子の数を考えてみます。すると、バンド上端付近に大きな状態密度を持つ構造は、そうでない構造と比べて少数の電子しか一定のエネルギー準位以下の状態に入れられません。したがって、バンド上端付近に大きな状態密度を持つ構造は、そうでない構造と比べてフェルミエネルギーが高いです。それにもかかわらず、全エネルギーはバンド上端付近に大きな状態密度を持つ構造のほうが低くなります。

ちょっと不思議ですが、空孔で考えてみると一目瞭然です。

「構造によらず電子がバンドに全部入った時のエネルギーも同じ」と仮定しているので、「電子が 入らなかった 状態のエネルギーが高ければ高いほど、全エネルギーは低い」ことになります。するとフェルミエネルギーがバンド下端の時の議論とほぼ同様になり、「バンド上端付近の高いエネルギー準位にたくさんの空孔を詰め込める構造のほうが有利」ということがわかります。

2020-08-30

浮動小数点数の比較演算に関するメモ

二つの浮動小数点数を比較する演算は、C言語の演算子としては6種類（<、<=、>、>=、==、!=）ありますが、NaNの取り扱いを考えるともっとあります。

一般的なのは、以下の四状況に関してそれぞれtrueとfalseがどうなるかを考えた14種類（2⁴=16のうち、2種類は「常にtrue」「常にfalse」で除外）です*1。

オペランドの少なくとも片方がNaNの場合
両方のオペランドがNaNではなく、左オペランドが右オペランドより小さい場合
両方のオペランドがNaNではなく、左オペランドが右オペランドより等しい場合（ただし、0.0と-0.0は等しいとみなされる）
両方のオペランドがNaNではなく、左オペランドが右オペランドより大きい場合

	どちらかがNaN	左＜右	左＞右	左＝右	備考
OEQ	`false`	`false`	`false`	`true`	C言語の`==`、RISC-VのFEQ
OGT	`false`	`false`	`true`	`false`	C言語の`>`
OGE	`false`	`false`	`true`	`true`	C言語の`>=`
OLT	`false`	`true`	`false`	`false`	C言語の`<`、RISC-VのFLT
OLE	`false`	`true`	`false`	`true`	C言語の`<=`、RISC-VのFLE
ONE	`false`	`true`	`true`	`false`
ORD	`false`	`true`	`true`	`true`
UNO	`true`	`false`	`false`	`false`
UEQ	`true`	`false`	`false`	`true`
UGT	`true`	`false`	`true`	`false`
UGE	`true`	`false`	`true`	`true`
ULT	`true`	`true`	`false`	`false`
ULE	`true`	`true`	`false`	`true`
UNE	`true`	`true`	`true`	`false`	C言語の`!=`

C言語の演算子は、!=以外はorderedな比較を行います。orderedな比較とは、オペランドの少なくとも片方がNaNの場合、falseになるというものです。そのため、整数の場合と異なり、a > bとa <= bの結果がともにfalseであるということが起こります（aもしくはbがNaNの場合に発生します）。

RISC系の命令セットでは、必要最小限の比較命令のみ用意して、後は条件を反転、などの方法で他の比較条件を実現することがあります。しかし、こういうことがあるので、コンパイラは安易に条件を反転したりできません。実際にRISC-Vコンパイラ（llc）がどのようにして各条件を実現しているかを調べてみると、以下のようになっていました。

比較	実現方法
OEQ	FEQを使う
OGT	FLTを使う（オペランド順序を逆にする）
OGE	FLEを使う（オペランド順序を逆にする）
OLT	FLTを使う
OLE	FLEを使う
ONE	ORDを計算して、FEQの結果を反転したものとANDをとる
ORD	オペランドごとに自分とのFEQを計算してANDをとる
UNO	ORDを計算して結果を反転
UEQ	UNOを計算して、FEQの結果とORをとる
UGT	FLEの結果を反転
UGE	FLTの結果を反転
ULT	オペランド順序を逆にしたFLEの結果を反転
ULE	オペランド順序を逆にしたFLTの結果を反転
UNE	FEQの結果を反転

いくつかはllvmの自動フォールバック（SelectionDAGLegalize::LegalizeSetCCCondCode関数）で生成されているので一部無駄のあるコード生成になっています。例えば、ONEを計算するにはOGTの結果とOLTの結果をORするのが、命令数の観点から最適と思われます。

実際、(a > b) | (a < b)をコンパイルしてみると以下の残念なコードが出力されます。 C言語フロントエンドは正しく最適化してllvm-IRにONEを出力するものの、RISC-Vバックエンドはその最適化をうまく活用できなかったようです。

        feq.d   a0, ft0, ft0
        feq.d   a1, ft1, ft1
        and     a0, a0, a1
        feq.d   a1, ft1, ft0
        not     a1, a1
        and     a0, a0, a1

ちなみに、UNOは、llvmの自動フォールバックでは「オペランドごとに自分とのUNEを計算してORをとる」ですが、RISC-Vコンパイラは最適な「オペランドごとに自分とのOEQを計算してANDしたものを反転」を出力します（RISCVInstrInfoD.tdに書かれています）。

*1:ほかにもオペランドとしてquiet NaNが与えられた場合に例外を発生させるかを考慮に入れるとさらに状況が二倍に増えます

2020-08-23

gnuplotの等高線でヒートマップを描きたい

二次元データを可視化する方法に、標高線を描くという方法があります。

f:id:lpha_z:20200823223749p:plain — 図1: 標高線を使った可視化

また、ヒートマップを使う方法もあります。

f:id:lpha_z:20200823223248p:plain — 図2: ヒートマップを使った可視化

標高線通りにヒートマップを作ってほしいのですが、gnuplotが標高線を描画するアルゴリズムとヒートマップを作る時の補間アルゴリズムは異なっているため、ずれが発生します。

f:id:lpha_z:20200823224526p:plain — 図3: 標高線とヒートマップを同時に使った可視化

これを一貫するようにしてみたいと思います。

戦略

標高線は、図1や図3だと直線で補間していますが、スプライン曲線で補間することもできます。しかし、スプライン曲線で補間しても、ヒートマップとはずれが生じます。

そこで、gnuplotに標高線データをsvgで出力してもらって、それをスクリプトで加工することを考えます。

実装

まずはこんな感じでgnuplotに標高線データをsvgで出力してもらって、

gnuplot -e "set terminal svg; set output 'contour.svg'; set view map; unset surface; set contour base; set xrange [0:4]; set yrange [0:4]; set xtics 1; set ytics 1; set ticscale 0; set cntrparam bspline; set cntrparam levels incremental 0,2,16; splot 'data' u 1:2:3 w l notitle"

次に以下のRuby スクリプトで加工します。

require 'rexml/document'

palette = ["rgb(96,1,199)", "rgb(136,6,249)", "rgb(167,20,111)", "rgb(193,48,0)", "rgb(216,93,0)", "rgb(236,161,0)", "rgb(255,255,0)"]

xml = REXML::Document.new(File.new("contour.svg"))

contour = xml.root.get_elements('g/g/g/path').map { |e| e.attribute("d") }.reverse

xml.root.elements.each_with_index('g/g/g/path') { |e, i|
    e.add_attribute("d", contour[i])
    e.add_attribute("fill", palette[i])
    e.add_attribute("stroke", palette[i])
}

File.write("contour.svg", xml)

reverseでひっくり返しているのは、描画順を入れ替えるためです。 gnuplotは標高線データを値の大きいほうから順番に出力してきます。

この方法でうまくいくためには、少なくとも以下の条件が要求されます。

等高線に谷のような部分がなく山だけがある（複数の頂上があるのは大丈夫）
- 標高の低い部分の塗りつぶしから順にキャンバスに書くことで、間が塗りつぶされているように見せているため
一周しない等高線の始点と終点は同じ辺に属している
- もともと線だったものを面に変換しているので、一周しない線の場合は始点と終点が結ばれる
- 標高線が切れる部分が同じ辺にある場合は始点と終点が結ばれて問題なくなる
- 標高線が切れる部分が違う辺にある場合は斜めに結ばれて正しくなくなる

結果

f:id:lpha_z:20200823231713p:plain — 図4: 標高線と一貫したヒートマップによる可視化

作図に用いたコマンド

元データ

使ったgnuplotコマンド

図1

gnuplot -e "set terminal png; set output 'contour.png'; set view map; unset surface; set xrange [0:4]; set yrange [0:4]; set xtics 1; set ytics 1; set ticscale 0; set pm3d at b; set pm3d interpolate 100,100; set cbrange [0:16]; set palette maxcolors 8; splot 'data' u 1:2:3 w l notitle"

図2

gnuplot -e "set terminal png; set output 'contour.png'; set view map; set contour base; unset surface; set xrange [0:4]; set yrange [0:4]; set xtics 1; set ytics 1; set ticscale 0; set cntrparam linear; set cntrparam levels incremental 0,2,16; splot 'data' u 1:2:3 w l notitle"

図3

gnuplot -e "set terminal png; set output 'contour.png'; set view map; set contour base; set xrange [0:4]; set yrange [0:4]; set xtics 1; set ytics 1; set ticscale 0; set cntrparam linear; set cntrparam levels incremental 0,2,16; set pm3d at b; set pm3d interpolate 100,100; set cbrange [0:16]; set palette maxcolors 8; splot 'data' u 1:2:3 w l notitle"

2020-08-16

完全精度sincosf関数の作り方（その1）

前に完全精度expf関数を作った（高速な完全精度 expf 関数の作り方 - よーる）ので、今度は完全精度sincosf関数を作っていきます。

作戦

まず、引数の絶対値が小さい時のsin関数とcos関数を（taylor展開などを利用して）作ります。引数の絶対値が大きい場合、引数 $x$ を $x = t +\frac\pi 2 k$ と分解します。その後 $k$ の値に応じて先に作ったsin関数やcos関数を使い分ければ答えが求まります。

ただし、この分解は簡単にできることではありません。今回は、ある程度簡単に実行可能な、引数の絶対値がそれほど大きくない（4×10⁸程度以下）場合のみを作ることにします。

ベースとなるTaylor展開ベースの実装

まず、sin関数とcos関数をtaylor展開で求めるコードを書きます。double精度で求めて最後にfloatに丸める作戦なので、そこまで精度は必要ではありません。

double cos_taylor(double x) {
        double x2 = x * x;
        double a = -1./87178291200;
        a = std::fma(a, x2, 1./479001600);
        a = std::fma(a, x2, -1./3628800);
        a = std::fma(a, x2, 1./40320);
        a = std::fma(a, x2, -1./720);
        a = std::fma(a, x2, 1./24);
        a = std::fma(a, x2, -1./2);
        return std::fma(a, x2, 1);
}

double sin_taylor(double x) {
        double x2 = x * x;
        double a = -1./1307674368000;
        a = std::fma(a, x2, 1./6227020800);
        a = std::fma(a, x2, -1./39916800);
        a = std::fma(a, x2, 1./362880);
        a = std::fma(a, x2, -1./5040);
        a = std::fma(a, x2, 1./120);
        a = std::fma(a, x2, -1./6);
        return std::fma(a, x2 * x, x);
}

$\frac\pi 2$ で割ったあまりを求める

引数 $x$ を $x = t + \frac\pi 2 k$ と分解するためには、以下のような計算を行います。これは、完全精度expf関数を作る時に行った、引数 $x$ を $x = t + \frac{\log 2}{2^N}k$ と分解するのと同様の手順です。

constexpr double harf_pi_hi = 0x1.921fb5p+0;
constexpr double harf_pi_lo = 0x1.110b4611a6263p-26;
constexpr double inv_harf_pi = 0x1.45f306dc9c883p-1;

double k = std::fma(x, inv_harf_pi, 0x3.p+51) - 0x3.p+51;
assert( -0x1.p+28 <= k && k <= 0x1.p+28 );
double t = std::fma(-harf_pi_lo, k, std::fma(-harf_pi_hi, k, x));

ここで、harf_pi_hi、harf_pi_lo、inv_harf_piは、以下の条件を満たすdoubleの定数です。

harf_pi_hiの下位Nビットは0
harf_pi_hiの値とharf_pi_loの値を足し合わせると、 $\frac\pi 2$ に近い
inv_harf_piの値は $\frac 2\pi$ に近い

harf_pi_hiの下位Nビットを0にするのは、k倍したときに誤差が生じないようにするためです。ここでN=28を選べば、先のコードになります。 N=28の時、xの絶対値が4.2e+8程度以下であればassertが満たされます。

なぜN=28を選んだか

なぜN=28を選ぶ必要があるのでしょうか。例えば、Nをもっと大きくすれば、対応できるxの範囲が広がりよさそうに思われます。しかし、実際には、Nを大きくすると $\frac\pi 2$ の有効桁数が小さくなり、それに起因する誤差が増加し、完全精度を達成できなくなります。このことを説明してみます。

doubleの精度は53ビットですが、二つのdoubleの和で表された数の精度は107ビットです。これは、下位を担当するdoubleの符号ビットも有効活用されるためです。ここで、Nビットを強制的に0にしなければならないので、 $\frac\pi 2$ の実効的な精度は107-Nビットです。

$t = x - \frac\pi 2 k$ を計算するとき、桁落ちが発生します。 $\frac\pi 2 k$ と $x$ とで打ち消しあう部分の長さはおおよそN+1ビットなので、最終結果には106-2Nビットしか残りません。

N=29の場合、小数部には48ビットしか使われないということになります。これはdoubleの精度を下回っています。実際、0x1.8db252p+25を $\frac\pi 2$ で割ったあまりを求めてみると、正確には0x1.5292b563fb139p-5程度になるべきところ、0x1.5292b563fb120p-5と求まってしまいます。なんと25ULPも誤差が発生しています。運の悪いことにcos(0x1.8db252p+25) = 0x1.527a08ffffff0p-5とfloatに丸めるときの丸め境界から16ULP程度しか離れていないため、ここの計算を間違えてしまいます。

よって、このような問題を発生させないためにも、N=28よりも大きくすることはできません（N=28であっても、運よくこのような問題が発生していないだけであり、やや精度不足であることには変わりありません）。

絶対値がそれほど大きくない場合の実装

void sincos_reduction(double x, double* s, double* c) {
        constexpr double harf_pi_hi = 0x1.921fb5p+0;
        constexpr double harf_pi_lo = 0x1.110b4611a6263p-26;
        constexpr double inv_harf_pi = 0x1.45f306dc9c883p-1;

        double k = std::fma(x, inv_harf_pi, 0x3.p+51) - 0x3.p+51;
        assert( -0x1.p+28 <= k && k <= 0x1.p+28 );
        double t = std::fma(-harf_pi_lo, k, std::fma(-harf_pi_hi, k, x));

        std::int32_t k_ = static_cast<std::int32_t>(k);

        switch(k_ & 3) {
                case 0: *s = sin_taylor(t), *c = cos_taylor(t); break;
                case 1: *s = cos_taylor(t), *c = -sin_taylor(t); break;
                case 2: *s = -sin_taylor(t), *c = -cos_taylor(t); break;
                case 3: *s =-cos_taylor(t), *c = sin_taylor(t); break;
        }
}

8/28 追記

入力が±0x1.33333p+13の時、sin(x) ~ ∓0x1.63f4bafffffffa5cp-2であり、完全精度なら∓0x1.63f4bap-2を返すべきですが、上記コードは∓0x1.63f4bcp-2を返すので完全精度になっていませんでした。

これはテイラー展開の誤差由来やリダクション部分の誤差由来ではなく、doubleを経由していること、つまり二回丸めていることが原因でした。sin(x) ~ ∓0x1.63f4bafffffffa5cp-2はdoubleに丸めると∓0x1.63f4bb0000000p-2に丸められ、floatの丸め境界上の値に丸まってしまいます。これをfloatに丸めると偶数丸めになるため最終桁が0になる∓0x1.63f4bcp-2に丸められ、正しくなくなります。

2020-08-09

RISC-Vのビット操作系拡張（B拡張）のまとめ（その5）

ブール環（ $\mathbb{F}_2$ 、 $\mathbb{Z}/2\mathbb{Z}$ 、“足し算”がxorで“掛け算”がandであるような演算体系、“2で割ったあまり”）に関係する命令群です。

`clmul`, `clmulh`, `clmulr`

繰上りなしの掛け算です。ビット列の畳み込みと考えることもできます。

clmulは掛け算の結果の下32/64ビットを返します。 clmulhは掛け算の結果の上32/64ビット（下から32～63/64～127ビット目）を返します。ここで、繰上りがないため、clmulhの最上位ビットは常に0です。

clmulrはその再上位ビットを除いた上32/64ビット（下から31～62/63～126ビット目）を返します。名前のrは、reverseに由来します。clmulrの結果は、ビット順序逆転した結果同士のclmulの結果のビット順序逆転と同じになるからです。

巡回冗長検査（Cyclic Redundancy Check, CRC）やハッシュ関数の実装などに有用だそうです。自分自身とのclmul/clmulhを使うと、0を00に、1を01に変換する（32/64ビットのビット列の各ビットの間に0を挿入して64/128ビットのビット列にする）操作を実現します。また、-1とのclmulはprefix sum（ここで言うsumはxor）を求めるのに使えます（グレイコードのデコードに使えます）。

`crc32`, `crc32c`

巡回冗長検査（CRC）をまさにそのまま計算する命令です。

`bmatxor`, `bmator`, `bmatflip`

64bitのレジスタを、8×8のビット行列だとみなして、行列演算を行う命令です。

bmatxorは、“足し算”がxorであるような行列乗算、bmatorは“足し算”がorであるような凝結乗算、bmatflipは行列転置です。

bmatflipは、zip三回分に相当します。

2020-08-08

RISC-Vのビット操作系拡張（B拡張）のまとめ（その4）

ビットフィールドやビットベクトルの操作系命令をまとめます。

`bfp`

rs1のoffビット目からlenビットをrs2の下lenビットに置き換えます。ここで、offとlenは（オペランド数の都合上）rs2の上位ビットを使って指定します。

offは、RV32であればrs2の16ビット目から20ビット目、RV64であればrs2の32ビット目から37ビット目または48ビット目から53ビット目を使います。 lenは、RV32であればrs2の24ビット目から27ビット目、RV64であればrs2の40ビット目から44ビット目または56ビット目から60ビット目を使います。

「または」と書いてある部分は、rs2の62ビット目から63ビット目が2であれば後者を、そうでなければ前者を使います。

lenは1から16/32になるはずです。そこで、4/5bitのフィールドは、すべてが0の時に16/32と解釈する工夫がなされます。

offやlenを指定する定数をrs2の上位ビットに送り込むには、pack命令が有用です。 RV32の場合、addiで作った12ビットの定数をpack命令で上位ビットに送り込めば、任意のlenとoffを実現できます。 RV64の場合、luiで作った32ビットの定数をpack命令で上位ビットに送り込めば、任意のlenとoffを実現できます（「または」と書いてある部分を後者で使うモードを利用します）。

`bext`, `bdep`

bext命令は、rs1のうち、rs2の中で1がたっているビットの部分だけを下位ビットに集約します。 bdep命令は、その逆操作を行います。つまり、rs2の中で1がたっているビットの位置に、rs1の各ビットを展開します。

bext命令を使うと、Magic Bitboardみたいなことが簡単にできます。

また、bdep命令を使うと簡潔データ構造で用いられるようなselect1(rank)演算が実現可能です。先頭から10ビット目の1の位置を知りたい場合、bdep a0, 1ull<<10, a0; ctz a0, a0のようにすることで求まります。

仮定