2020-02-09

近況

輻輳が発生しています……

2020-02-02

64bit符号なし整数を素因数分解するのにかかる時間

64bit符号なし整数を素因数分解するため、事前に32bit符号なし整数で表せる素数をリストアップしておくとします。

方法1: 素数のリストをビットベクトルで保持し、試し割り

32bit符号なし整数で表せる奇数が素数であるかを、0, 1で表した列（ビットベクトル）を事前に作成し、それを素因数分解プログラムのバイナリに埋め込みます。ビットベクトルのエントリ数は2Gi個で、1Byte=8bitであるとして256MiByteのバイナリが完成します。

32bit符号なし整数で表せる素数は203280221個（203M個）存在するので、この方法を使うと最悪で203M回の除算が発生します。

最近のCPUであっても除算器はパイプライン化されていないので、203M回×除算器の占有サイクル数くらいの時間がかかります。

この方法を使ったとき、手元のマシンでは、4.8秒くらいかかりました。ちょっと時間がかかりすぎな気がします。

方法2: 素数の逆元をファイルから読み込み、判定

m % p == 0という剰余演算はpが奇数の定数であれば、pの2⁶⁴における逆元を用いて最適化可能です。これを事前に計算しておく方式です。一つの素数当たり16Byteの情報が必要なため、バイナリに埋め込むことはできませんでした。

事前に計算した定数をファイル（バイナリファイルで3.03GiB、おそらくメモリ上にキャッシュ済み）から読み込んだ場合、素因数分解には5.2秒かかりました。一つの素数で割れるかを確認するのは1cycleで可能で非常に高速なのですが、読み込みに5.0秒くらいかかっているので、この方式をこれ以上改善することはできなさそうです。

mmapで読み込んだら、読み込み1秒、計算1秒くらいまで高速化しました。素因数が小さいとき、読み込みが高速になるのもよい点です。

さすがに事前計算に時間がかかっているだけあって高速な手法です。

方法3: ポラードのρ法

以前に実装した高速なgcd実装を用い、フロイドの循環検出法を用いたオリジナルバージョンを実装しました。

素数であるかわかってない入力が来る場合、ポラードのρ法は素数であるかを確認するのに非常に不向きなので、別途素数判定アルゴリズムを併用することになります。

素数でないとわかっている入力に対し、何らかの因数を見つけるのにかかる時間は、（乱択アルゴリズムなのでぶれが大きいですが）平均5秒、最悪20秒程度でした。

素因数は2³²ほどの大きさであり、これの平方根（65536）くらいの回数の乱数生成で発見できると見積もるようですが、実際には数千万個の乱数を生成しないと素因数が発見できない例が多々ありました。

事前の計算は不要でプログラムサイズも小さいですが、素数のリストを事前に作成した試し割り法よりは低速だということがわかりました。

2020-01-26

ABC152のE問題 Flatten を遅い方法で解いた

以下にはABC152のE問題 Flatten のネタバレが含まれます。

この問題は多倍長整数を使用できる言語であれば愚直解が通ります。Haskellで書かれた愚直解のコードは、例えばSubmission #9632460 - AtCoder Beginner Contest 152が簡潔で読みやすいです。

foldr lcm 1 aの部分は、 $N$ 回「 $O\left(N\log A\right)$ ビットの整数と $O\left(\log A\right)$ ビットの整数の最大公約数（高々 $O\left(\log A\right)$ ビット）を求め、 $O\left(\log A\right)$ ビットの整数をそれで割り、その結果と $O\left(N\log A\right)$ ビットの整数を掛け合わせる」という計算を行うはずです。 $O\left(N\log A\right)$ ビットの整数と $\log A$ ビットの整数の最大公約数をユークリッドの互除法で求めると、多倍長除算は長除法（筆算）を用いるとして、時間計算量が $O\left(\left(N + \log A\right)\left(\log A\right)^2\right)$ となります。

残りの部分は、取るに足らない計算量です。

制約は $N=10^4, A=10^6$ なので、愚直解は雑な概算でΟの中が $4\times 10^9$ になります。時間制限は2秒であり、 $5\times 10^9$ サイクル以上のCPU時間があるので、これは余裕で間に合います。実際、多くのHaskellでの愚直解の提出は、200msから300msくらいで動作しています。

C++標準には残念ながら多倍長整数がありません。しかし、この問題に限っては位取り記数法を用いた多倍長整数を使う必要はありません。掛け算せず、因数*1の列として保持すればよいです。

それを利用した解法の提出が以下です。

Submission #9637654 - AtCoder Beginner Contest 152

全然間に合っていません。実際、この解法の時間計算量は $O\left(N^2 \left(\log A\right)^3\right)$ であり、愚直解より20倍ほど遅いことがわかります。つまり、3倍程度間に合いません。

手元で作った雑な最大ケースでは間に合ったため大丈夫だろうと思っていたのですが、よくよく考えると真の最大ケースは、「 $N$ 個の $A$ 未満の相異なる素数」です（テストケースの実行時間から推察すると、おそらくmax_02がこれです）。

そのケースを手元で動かした場合からの概算だと、ジャッジサーバーでは6秒くらいかかりそうだという感じだったので、いろいろな計算が大体合います。

そこで、自作の高速化したgcd関数を使えば間に合うのではないかと考えました。

以下は、最大公約数をもっと高速に求める（その4） - よーるで示した高速なgcd関数を、入力が $\frac{2^{64}}5$ を超えない場合に特化したものです。

std::uint64_t gcd_impl( std::uint64_t n, std::uint64_t m ) {
    constexpr std::uint64_t K = 5;
    if( m == 0 ) { return n; }
    for( ; n / 64 < m; ) {
        std::uint64_t t = n - m;
        std::uint64_t s = n - m * K;
        bool q = t < m;
        bool p = t < m * K;
        n = q ? m : t;
        m = q ? t : m;
        if( m == 0 ) { return n; }
        n = p ? n : s;
    }
    return ::gcd_impl( m, n % m );
}
 
std::uint64_t gcd( std::uint64_t n, std::uint64_t m ) {
    return n > m ? ::gcd_impl( n, m ) : ::gcd_impl( m, n );
}

また、64 bit整数には $A^3$ が収まるので、因数を3つまとめて保持することにより、高速化を図ります。これをおしすすめると愚直解になるわけです。

ところで、このgcd関数は、g++でコンパイルすると、mov命令を減らしたいのかわかりませんが低速なcmova命令にコンパイルされてしまいます（最大公約数をもっと高速に求める（その2）【cmova命令は遅い】 - よーる）。clang++だとそのようなことはないため、clang++で提出します。

Submission #9639909 - AtCoder Beginner Contest 152

1825msというぎりぎりですが、通すことができました。

かかった時間から推察すると、自作のgcd関数は一回当たり300サイクル程度で動作しているということがわかりました。

*1:素因数ではありません

2020-01-19

Karatsuba乗算の最適化をした

Karatsuba乗算、あるいはToom–Cook2乗算とは、入力a0、a1、b0、b1が与えられたとき、a0*b0、a0*b1+a1*b0、a1*b1の三つを出力する際、必要な乗算を四回から三回に減らせるアルゴリズムです。ただし、減算ができることを要求します。

長乗法（筆算と同様の手順による多倍長整数の乗算のこと）や数列の畳み込み演算は、愚直に行うとΘ(N²)の計算コストがかかります。

Karatsuba乗算のアルゴリズムを使って多倍長整数の乗算を行う場合、桁数が半分の整数の乗算を3回やることになります。ここで出現した乗算にも再帰的にKaratsuba乗算のアルゴリズムを適用していけば、Ο(N^1.585)の計算コストとなります。ここで、指数の1.585は、正確には $\log3/\log2$ です。

具体的なアルゴリズムは、非常に単純です。

a0*b0、a1*b1を計算する。乗算がここで二回必要です。
p = a0-a1、q = b0-b1を計算する。減算がここで必要です。
p*qを計算する。乗算がここでもう一回必要です。
a0*b0 + a1*b1 - p*qを計算する。これがa0*b1 + a0*b1になっています。

このように計算結果を再利用することで、乗算の回数を減らすことに成功しています。

以下は、この再帰的にKaratsuba乗算を適用したプログラムを最適化する話です。

問題

Convolution (mod 1,000,000,007)

剰余類環*1 $\mathbb{Z}/1000000007\mathbb{Z}$ における畳み込みを行う問題です。

入力長Nは2¹⁹であり、最も高速な解法は畳み込み定理と高速フーリエ変換（FFT）*2を用いた方法でしょう。

しかし、高速フーリエ変換を用いた手法がKaratsuba乗算より高速になるのはN～5000程度とされており、その百倍程度の入力長であればKaratsuba乗算で強引に押し切ることも可能です*3。

実装

私が書いたもの

#include <cstdio>
#include <cstdint>
#include <cinttypes>

static constexpr int64_t MOD = 1'000'000'007;
static constexpr size_t N_MAX = 1<<19;
static constexpr size_t K = 8;

template<size_t N>
void karatsuba( const int64_t* x, const int64_t* y, int64_t* z ) {
    static int64_t t[N], p[N/2], q[N/2];

    for( size_t i = 0; i < N/2; ++i ) {
        p[i] = (x[i] - x[i+N/2] + MOD) % MOD;
        q[i] = (y[i] - y[i+N/2] + MOD) % MOD;
    }

    karatsuba<N/2>( x, y, z );
    karatsuba<N/2>( x + N/2, y + N/2, z + N );
    karatsuba<N/2>( p, q, t );

    for( size_t i = 0; i < N/2; ++i ) {
        int64_t a = z[i];
        int64_t b = z[i+N/2];
        int64_t c = z[i+N];
        int64_t d = z[i+N+N/2];
        int64_t e = t[i];
        int64_t f = t[i+N/2];

        z[i+N/2] = a + b + c - e;
        z[i+N]   = d + b + c - f;
    }
}

template<>
void karatsuba<K>( const int64_t* x, const int64_t* y, int64_t* z ) {
    uint64_t tmp[K*2] = {};

    for( size_t i = 0; i < K; ++i ) {
        for( size_t j = 0; j < K; ++j ) {
            tmp[i+j] += x[i] * y[j];
        }
    }

    for( size_t i = 0; i < K*2; ++i ) {
        z[i] = tmp[i] % MOD;
    }
}

int main() {
    static int64_t a[N_MAX] = {}, b[N_MAX] = {}, answer[N_MAX*2];

    size_t n, m;
    scanf( "%zd%zd", &n, &m );

    for( size_t i = 0; i < n; ++i ) {
        scanf( "%" SCNd64, &a[i] );
    }

    for( size_t i = 0; i < m; ++i ) {
        scanf( "%" SCNd64, &b[i] );
    }

    karatsuba<N_MAX>( a, b, answer );
    for( size_t i = 0; i < n+m-1; ++i ) {
        if( i ) { putchar(' '); }
        printf( "%" PRId64, (answer[i]%MOD+MOD)%MOD );
    }
    puts("");

    return 0;
}

工夫その1：静的に領域を確保する

Karatsuba乗算に必要な作業領域は、Nが決まれば静的に決定可能で、動的確保は不要です。

一見再帰呼び出し的なコードになっていますが、常にNが半分ずつになっていく再帰呼び出しです。よって、テンプレートで記述すれば毎回異なる関数を呼び出す通常の関数呼び出しとなり、関数内static変数として作業領域を確保することができます。

工夫その2：定数回ループ以外の条件分岐を取り除く

条件分岐は遅くなる原因であるので、取り除きます。

0埋めされた計算とはいえ、常に最大ケースと同じだけの計算量が必要になっているので、小さいケースでも時間がかかります。

ジャッジサーバーに負荷をかけてすみません。

工夫その3：小さい部分では再帰しない

Karatsuba乗算は、乗算の回数がオーダーレベルで減る代わりに加減算の回数が増えるという性質があります。

そのため、入力長Nが非常に小さい部分ではΘ(N²)の愚直な方法が高速になります。

定数での剰余演算は、乗算命令二回に最適化されます。これは、a%b == a-a/b*bであること、定数での除算は逆数の乗算に最適化できる*4ことによります。このことを考慮に入れた、切り替え入力長K以下で愚直計算に切り替えた場合の乗算命令数は以下のようになります。

`K`	剰余演算に必要な乗算命令数	愚直計算に必要な乗算回数	合計
`1`	～3¹⁹×8	3¹⁹	104.6G
`2`	～3¹⁸×16	3¹⁸×4	77.5G
`4`	～3¹⁷×32	3¹⁷×16	62.0G
`8`	～3¹⁶×64	3¹⁶×64	55.1G
`16`	～3¹⁵×128	3¹⁵×256	55.1G
`32`	～3¹⁴×256	3¹⁴×1024	61.2G

この結果を見ると、乗算回数としてはK=8またはK=16がよさそうということになります。加減算の数のことまで考えると、N=16の時にKaratsuba法を適用してN=8とするのは乗算回数が減らないのに加減算だけが増える悪手であり、K=16とするのが適切であると結論付けることができます。

しかし、実験の結果、入力長が8のところで方法を切り替えるのが最も高速という結論になりました。

入力長が16のところで切り替えたほうがやや計算量が少なくなるはずですが、コンパイラのインライン展開の限界の影響か、実測では遅くなりました。

切り替え入力長が8であっても、愚直計算のロード命令の数を最小化しようとするとレジスタが20個ほど必要で、x86の持つ汎用レジスタの数16個を超えてしまいます。切り替え入力長を16にしてしまうとさらにロード命令が増えてしまいます。これが切り替え入力長16で遅くなる原因と思われます。

工夫その4：なるべく剰余演算を行わない

定数での剰余演算は、乗算命令二回に最適化されますが、それでも重い演算であることには変わりがなく、なるべく行いたくありません。

プログラムの入力は0以上MOD未満ですが、途中で加減算をするとこの制約が満たされなくなります。そのため、桁あふれを防ぐため、少なくとも乗算の前までには0以上MOD未満の制約を満たすようにしておく必要があります。 karatsuba<K>の最初で剰余をとっても良いですが、もともとの入力は剰余をとらなくてもいいことを考えると、pやqを計算した後に行うのが最も回数が少なくなります（回数の差はほとんどはありません）。

乗算結果はMOD*MOD未満であり、これを8個足すくらいでは*5uint64_tは桁あふれしないため、乗算ごとに剰余をとる必要はありません。 karatsuba<K>が終わるときに剰余をとればよさそうです。

z[i+N/2] = a + b + c - e;の部分で四つ足すため、再帰から一段戻るたびに絶対値が4倍程度になりますが、再帰深度はlog(2¹⁹/8)=16なので、絶対値は2³²倍にしかなりません。 MODの2³²倍はint64_tに収まるため、ここでも剰余をとる必要はありません。最後の出力する部分で剰余をとればよいでしょう。

工夫その5：データのコピー回数を減らす

出力すべき答えのうち前1/4と後1/4は、再帰呼び出しで得られる答えの前1/2や後1/2と正確に一致します。ここのコピーにかかる時間も実行時間にかなり影響してくるので、再帰関数にポインタを渡して直に書き込んでもらうことでコピーコストを減らします。

ポインタで渡すのはエイリアス解析が困難になりそうですが、最近のコンパイラだと意外とやってくれるようです(?)

更なる最適化

うさぎさんが最適化したもの

p[i] = (x[i] - x[i+N/2] + MOD) % MOD;

としていた部分を

p[i] = x[i] - x[i+N/2];
p[i] += MOD * (p[i] < 0);

と変形することで演算強度を下げる最適化です。x[i]とx[i+N/2]は0以上MOD未満であることが保証されているので、この変形は無問題です。

p[i] += MOD * (p[i] < 0);の部分は、手元のg++-8やg++-9を使った場合、以下のどれを用いてもほとんど同じ速度となりました。

p[i] = p[i] < 0 ? p[i] + MOD : p[i];
p[i] += p[i] < 0 ? MOD : 0;
p[i] <0 && (p[i] += MOD);
if( p[i] < 0 ) p[i] += MOD;

いずれもcmovs命令を用いたコードにコンパイルされます。

ジャッジサーバー*6では、if文を用いたものは低速だったようです。テストケースごとの傾向から、これは分岐予測ミスによるものだと考えられます。ここの分岐は予測不可能な分岐であるため、条件分岐命令にコンパイルされると速度が低下するのでしょう。

ただし、wandboxを使った調査では、g++-6以降やclang++の場合はif文を用いても速度低下はありませんでした。

これを用いた場合、手元のCPU（Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz）では（user時間で）2.0秒程度で実行が完了します。このCPUのターボブースト時クロック周波数は3.50GHzであり、このプログラムに含まれる乗算41.3億回を行うためにはそれだけで最低でも1秒以上かかります*7。 Karatsuba乗算にはこれ以外にも加減算が多く含まれることを考えると、かなり限界に近い性能が出せていることになるでしょう。

*1:環とは、おおざっぱに言って加算・減算・乗算が行えるような数の集合です。

*2:ここでは、複素数環に限らず任意の代数環で行える広義のフーリエ変換を言っています。剰余類環で行うフーリエ変換は、数論変換（NTT）とも呼ばれます。

*3:NlogNとN^1.585に5000と2¹⁹を代入してみると、差はおおよそ10倍程度であるということがわかります

*4:浮動小数点数の除算を逆数乗算に変形した場合はわずかな誤差が発生しますが、整数の場合は正確に求める手法があります

*5:16個でも大丈夫です

*6:g++ 7.4.0を使っているらしいです

*7:乗算命令のスループットは1であり、このCPUでは1秒間には35億回しかできません

2020-01-12

64bit版RISC-Vで即値を生成するときの罠

RISC系の命令セットでは、一命令でレジスタに書き込める即値と、一命令では書き込めない即値が存在する設計になっています。一命令では書き込めない即値の場合、二命令以上かけて即値を生成することになるわけですが、そこに存在する罠に引っかかった記録です。

なんでそんな仕組みになっているのか

RISC系の命令セットでは、一命令の長さを固定（ほとんどの場合、4Byte）にすることが基本です。レジスタ幅が32bit（＝4Byte）だとすると、一命令32bitの中に32bit分の即値を埋め込むことは明らかに不可能です。

一命令の長さをもっと長く、例えば8Byteにすれば入るじゃないか、という感じですが、そうするとプログラムのサイズがほぼ二倍に増えてしまい、受け入れがたいという設計方針のようです。

一命令の長さを可変にすればいいじゃないか、というのはプロセッサの設計が複雑化するため、やはり受け入れがたいというのがRISCの設計方針です。

MIPSの例

MIPSの命令セットでは、LUI命令（命令中の16bit値を16bit左シフトした値をデスティネーションレジスタに書き込む命令）とORI命令（命令中の16bit値をゼロ拡張した値とソースレジスタの値のビット毎論理和をデスティネーションレジスタに書き込む命令）を組み合わせ、任意の32bit即値を作ることができます。

例えば、0x12345678をt0レジスタに書き込みたいときは、

LUI   t0, 0x1234
ORI   t0, t0, 0x5678

のような命令列になります。上位16bitと下位16bitに分ければいいだけなので、非常に簡単です。

32bit版RISC-Vの例

32bit版RISC-Vの命令セットでは、LUI命令（命令中の20bit値を12bit左シフトした値をデスティネーションレジスタに書き込む命令）とADDI命令（命令中の12bit値を符号拡張した値とソースレジスタの値のビット毎論理和をデスティネーションレジスタに書き込む命令）を組み合わせ、任意の32bit即値を作ることができます。

MIPSの命令セットの場合と異なる点は、

上位16bitと下位16bitに分割するのではなく、上位20bitと下位12bitに分割する
ADDI命令の即値は符号拡張される

の二点です。

符号拡張される点が厄介で、例えば0x89abcdefをt0レジスタに書き込みたいとき、

LUI   t0, 0x89abd
ADDI  t0, t0, 0xdef

のような命令列で作る必要があります。LUIの即値が、作りたい即値0x89abcdefの上位20bitをそのまま切り出した値ではなく、それより1だけ大きい値になっている点に注意してください。 ADDI命令の即値0xdefは、符号拡張すると0xfffffdefになるので、余計な0xfffffを打ち消すために1だけずれた値を作る必要があります。

64bit版RISC-Vの例

64bit版RISC-Vの命令セットでは、LUI命令（命令中の20bit値を符号拡張し12bit左シフトした値をデスティネーションレジスタに書き込む命令）とADDI命令（命令中の12bit値を符号拡張した値とソースレジスタの値の加算結果をデスティネーションレジスタに書き込む命令）を組み合わせ、任意の符号付き32bit即値を作ることができ……ません。

この二命令で作れる値の範囲は、[-2³¹-2048, 2³¹-2048)の範囲の数です。符号付き32bit値の範囲は[-2³¹, 2³¹)なので、微妙にずれています。

なぜこんなことになってしまうのでしょうか？

符号付き32bit値の範囲に収まっており、かつLUI命令とADDI命令を組み合わせても作れない範囲の数である0x000000007ffffabcをt0レジスタに書き込みたいときのことを考えてみます。

まず、LUI命令では下位12bitに触ることができませんので、ADDI命令の即値は必然的に0xabcになります。この値は符号拡張すると0xfffffffffffffabcになるので、余分な0xfffffffffffffを打ち消すために、LUI命令では1だけずれた値を作る必要があります（32bit版の時と同じ論理です）。よって、LUI命令で0000000080000000をt0レジスタに書き込めばよさそうです。

しかし、これは命令セットの都合上不可能です。LUI命令は命令中の20bit値を符号拡張するため、その20bitの中の最上位ビットが1となる正の数（0x0000000080000000）をレジスタに書き込むことはできません。

解決法

32bit符号付き整数に収まる値を生成したい場合、LUI命令とADDIW命令を使うのが正解です。ADDIW命令は、ADDI命令と似ていますが、加算結果の下32bitを符号拡張した値をデスティネーションレジスタに書き込む命令です。これを使うと0x000000007ffffabcを以下のように作成できます。

LUI   t0, 0x80000
ADDIW t0, t0, 0xabc

まず、最初のLUI命令で、t0レジスタには0xffffffff80000000が書き込まれます。次のADDIW命令では、t0レジスタに入っている値である0xffffffff80000000と、即値を符号拡張した値である0xfffffffffffffabcが加算され、0xffffffff7ffffabcが途中結果として得られます。そして、下32bitが符号拡張され、0x000000007ffffabcがt0レジスタに書き込まれ、目的が達成されます。

gccでは

gccでは、極力64bit命令を使いたいらしく、問題となる[2³¹-2048, 2³¹)の区間の即値を作る場合、LUI命令とXORI命令を組み合わせたコードが出力されます（gcc7.1.1で確認）。

LUI   t0, 0x80000
XORI  t0, t0, -1348

これはなるほどという感じですが、この命令列に特にメリットがあるようには思えません。逆に、C.ADDIW命令に圧縮する機会を失っているという明確なデメリットがあります*1。

余談：ほかの命令ではこのような問題は発生しないのか

先ほどの例では、64bit用の命令であるADDI命令ではなく、32bit用の命令であるADDIWを使うことで問題が解決されました。

しかし、32bit用の対応する命令がない命令も多く存在します。

例えば、BLT命令（二つのソースレジスタの値rs1、rs2を二の補数表現で表された符号付き64bit整数として比較し、rs1 < rs2なら分岐する命令）や、BLTU命令（二つのソースレジスタの値rs1、rs2を符号なし64bit整数として比較し、rs1 < rs2なら分岐する命令）には、32bit版がありません。

これらは大丈夫なのでしょうか？

実は、不思議なことに大丈夫なようになっています。

まず、RISC-Vでは、32bit値は64bitレジスタ中に符号拡張された表現で入ることになっています。よって64bit値とみなして符号付き比較を行ってもまったく問題ありません。

これだと符号なし比較の方はまずそうですが、実は大丈夫です。なぜなら、符号拡張された表現は以下のようになるため、これを符号なし比較すれば大小関係は一致するからです。

32bit値	符号拡張された64bit値
`0x00000000`	`0x0000000000000000`
`0x00000001`	`0x0000000000000001`
`0x00000002`	`0x0000000000000002`
:	:
`0x7ffffffe`	`0x000000007ffffffe`
`0x7fffffff`	`0x000000007fffffff`
`0x80000000`	`0xffffffff80000000`
`0x80000001`	`0xffffffff80000001`
:	:
`0xfffffffe`	`0xfffffffffffffffe`
`0xffffffff`	`0xffffffffffffffff`

他にも、ビット毎論理演算命令は32bit版がありませんが、符号拡張していれば上位33bitには全く同じビットが入っているため、結果は32bit版も64bit版も同じになります。

まとめ

64bit版RISC-Vにおいて32bit符号付き整数をレジスタに書き込みたいとき、LUI命令とADDIW命令を使うのが最も楽です。LUI命令とADDI命令の組み合わせでは[2³¹-2048, 2³¹)の範囲の値が生成できないという罠があります。

*1:このコード例では即値が大きすぎて圧縮することは不可能ですが、即値が非常に小さい場合でもXORI命令が使われることを確認しています

2020-01-05

新年になりました

今年もよろしくお願いします。

2020ってちょっと不思議な感じのする数字です。

今年は、きっかけを待つのではなくて、自分で行動を起こせる年にしたいです。

あと、お酒飲むのやめます。心身の健康に悪いので……。

2019-12-29

write関数の変な使い方

ゴルフの話ではないです（ゴルフに役立たないとは言っていません）。

C言語のwrite関数は、unistd.hで宣言された、ファイルディスクリプタを通じて書き込む関数です。その宣言は以下のようになっています。

ssize_t write(int fd, void* buf, size_t len)

コードゴルフであればfdに変なファイルディスクリプタを渡したりするところですが*1、今回はそういう話ではないです。

次に変なものを渡せるのは、bufです。ここに変なポインタを渡すと、lenが0でなければエラーになり、write関数の返り値は-1になります。 errnoは、14、つまりEFAULTにセットされます。

「lenが0でなければ」なので、lenが0なら正当です。write関数の返り値は0になります。

一見意味のない行為に思えますが、どうもfflushの実装に使えるようです。実際、軽量libc実装のmuslでは、このテクニックが使われています。

git.musl-libc.org

（該当箇所引用）

 /* If writing, flush output */
    if (f->wpos != f->wbase) {
        f->write(f, 0, 0);
        if (!f->wpos) {
            FUNLOCK(f);
            return EOF;
        }
    }

エミュレータでこのコードを動かした時、おかしな挙動になっていて苦労したという話でした。

システムコールエミュレータが「ポインタbufが未割当の領域を指していておかしい！」と早合点する実装になっていたのです。

*1:read関数の第一引数に0(stdin)ではなく1(stdout)を渡して読み込まないというのは常套手段ですが、write関数でこういうことをやるのはあまり見たことがない気がします。不勉強なだけかもしれません。