手元環境と AtCoder 環境で sin の出力が異なるのはな〜んでだ？

人「浮動小数点数の演算なんだから環境によって結果が変わりうるのは当然だろ」

そうか...？

• 前提
  • 補足
• 事象
• 仮説と検証
  • 仮説 1
  • 仮説 2
  • 仮説 3
  • 仮説 4
• 調査
  • 関数定義
  • 実行時の処理
  • おまけ
  • 補足
• あとがき
• おわり

前提

$\gdef\libmsin#1{\operatorname{\texttt{s\hspace{-0.03em}i\hspace{-0.03em}n}}(#1)}$ $\gdef\libmcos#1{\operatorname{\texttt{c\hspace{-0.03em}o\hspace{-0.03em}s}}(#1)}$

次のような前提があります。

• IEEE 754 的には、correctly rounded (c.r.) な値を返すことを要求している
  • 無限精度で計算した値を正確に丸めた場合と同じ値ということ
  • これが各環境で正しく実装されているなら、環境ごとの差異はなくなる
• 実際には「C 言語で提供されている関数」と「IEEE 754 で定められている関数」は異なる
  • 単に名前が一致しているだけで、後者の機能を提供することにはなっていない
  • Annex F の章で一部の関数では対応付けが明記されていて、それらに関しては c.r. が保証されることになる
• C 言語の処理系が IEEE 754 に準拠しているとは限らない
  • 処理系は精度の保証をしなくてもいいという記述が C の規格にある 🤪
  • 広く使われている処理系においては準拠していると思ってよいはず

特に、処理系が IEEE 754 準拠である前提で、四則演算・平方根・融合積和については c.r. が保証されていることを前提にできます。ここで、融合積和 (fused multiply-add; FMA) は浮動小数点数 $x$, $y$, $z$ に対して $\roundcirc{xy+z}$ を返す演算です。$\roundcirc{\roundcirc{xy}+z}$ とは一般に一致しません。また、$\roundcirc{x}$ は実数 $x$ を浮動小数点型に丸めた値を意味します。

ところで、table-maker’s dilemma という名前で知られている問題があり、三角関数を始めとする数学関数たちの c.r. な値を得るのは自明ではありません。

rsk0315.hatenablog.com

そのため、sin のような関数は c.r. であることを期待できません。実際多くの処理系で c.r. でない値が返ってくるケースが存在しますし、たとえば double における sin の TMD は（えびちゃんの知る限り）未解決です*1。

一方で、上記で挙げた一部の演算たちは c.r. であることに注意します。たとえば cos を融合積和を用いて $\libmcos x = \roundcirc{(-x)(x/2) + 1}$ として定義しているライブラリがあったとしましょう*2。使っている演算に可搬性があるので、当然 cos 自体にも（同じライブラリを使っている環境同士なら）可搬性があります。実際のライブラリにおいても（ここまで単純ではないものの）四則演算などを用いて近似をしているため、同様の主張はできます。

なので、共通のライブラリを使っている環境同士で差異があった場合、あまり自然に「当然だろ」とは言えないでしょうという気持ちがあります。

補足

コンパイル時の計算と実行時の計算では精度が異なる場合があります。そのため、実行時の精度について調べようとしているときは、コンパイル時に計算した値が使われていないことを（適宜 disassemble しつつ）確かめておく必要があります。

たとえば、立方根を計算する以下のコードを考えます。

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

int main(void) {
    printf("%.60g\n", cbrt(3375.0));
    printf("%.60g\n", cbrt(strtod("3375.0", NULL)));
}

前者はコンパイル時に計算され、後者は（文字列から変換する部分が実行時に行われて）実行時に計算されて、下記の出力になる場合があります (wandbox)。

15
14.9999999999999982236431605997495353221893310546875

事象

さて、次の入力を考えます。

$$ \begin{aligned} x &= \texttt{1}.\texttt{\small 00000000004AD}_{(16)}\times 2^{0} \\ &= {\small 1.000000000000265}{\footnotesize 787392095262475}{\scriptsize 777417421340942}{\tiny 3828125}. \end{aligned} $$

AtCoder 環境では次のようになりました。 $$ \begin{aligned} \libmsin x &= \texttt{1}.\texttt{\small AED548F0911FB}_{(16)}\times 2^{-1} \\ &= {\small 0.841470984808040}{\footnotesize 056712741261435}{\scriptsize 439810156822204}{\tiny 58984375}. \end{aligned} $$ 手元環境は M2 Mac 上の colima で動かしている Ubuntu です*3。

FROM ubuntu:25.04@sha256:dc4565c7636f006c26d54c988faae576465e825ea349fef6fd3af6bf5100e8b6

% colima start --cpu 2 --memory 2 --disk 30

この環境では次のようになりました。 $$ \begin{aligned} \libmsin x &= \texttt{1}.\texttt{\small AED548F0911FC}_{(16)}\times 2^{-1} \\ &= {\small 0.841470984808040}{\footnotesize 167735043723951}{\scriptsize 093852519989013}{\tiny 671875}. \end{aligned} $$

真の値は $$ \sin(x) = \texttt{1}.\texttt{{\small AED548F0911FB}{\footnotesize 7FEDF452F4881}{\scriptsize \ldots}}_{(16)}\times 2^{-1} $$ なので、AtCoder 側では c.r. な値が得られていることになります。

note: $1.0$ から $2^{-52}$ ずつ増やしていって最初に見つかったケースがこれだったというだけで、一致しないケースはたくさんあります。

note: M2 Mac 側が c.r. になるケースとしては $x = \texttt{1}.\texttt{00000000088BC}_{(16)}\times 2^{0}$ がありました。

う〜〜〜ん。この差異はどのようにして生じるものなのでしょう？

仮説と検証

実際に調べた手順に沿って書いていきます。目次でネタバレになるのを避けるため番号で章立てしています。

仮説 1

安直な予想として、libc のバージョンが違うのではないか？というのが最初に思い当たる部分でした。 実際には sin の実体は libm の中にあるはずですが、コードのバージョンとしては libc と共通だと思われるので、libc のバージョンで調べることにします。

AtCoder のコードテストで下記を実行しました。

cp /lib/x86_64-linux-gnu/libc.so.6 .
chmod u+x libc.so.6
./libc.so.6

元の libc.so.6 の実行権限がなかったことと、それに対しての chmod ができなかったので、コピーしてきたものを実行しました。あるいは下記を実行するのでもよかったかもしれません。

/usr/lib64/ld-linux-x86-64.so.2 /lib/x86_64-linux-gnu/libc.so.6

得られた結果は次の通りです。

GNU C Library (Ubuntu GLIBC 2.36-0ubuntu4) stable release version 2.36.
Copyright (C) 2022 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.
There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.
Compiled by GNU CC version 12.2.0.
libc ABIs: UNIQUE IFUNC ABSOLUTE
Minimum supported kernel: 3.2.0
For bug reporting instructions, please see:
<https://bugs.launchpad.net/ubuntu/+source/glibc/+bugs>.

手元では次の通りでした。

% /lib/x86_64-linux-gnu/libc.so.6

GNU C Library (Ubuntu GLIBC 2.41-6ubuntu1) stable release version 2.41.
Copyright (C) 2025 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.
There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.
Compiled by GNU CC version 14.2.0.
libc ABIs: UNIQUE IFUNC ABSOLUTE
Minimum supported kernel: 3.2.0
For bug reporting instructions, please see:
<https://bugs.launchpad.net/ubuntu/+source/glibc/+bugs>.

たしかにバージョンは違っていますが、このバージョンの違いで挙動に差が出るのかはわかりません。そこで、別のジャッジでも調べてみることにします。

Library Checker ではコードテストがない（と思う）ので、A + B を用いて調べます*4。標準エラー出力が見れるので助かります。

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

int main() {
    int a, b;
    scanf("%d %d", &a, &b);
    printf("%d\n", a + b);

    system("/lib/x86_64-linux-gnu/libc.so.6 >&2");
}

得られた結果は次の通りです。

GNU C Library (Debian GLIBC 2.36-9+deb12u10) stable release version 2.36.
Copyright (C) 2022 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.
There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.
Compiled by GNU CC version 12.2.0.
libc ABIs: UNIQUE IFUNC ABSOLUTE
Minimum supported kernel: 3.2.0
For bug reporting instructions, please see:
<http://www.debian.org/Bugs/>.

また、同ジャッジサーバの Dockerfile.GCC を手元で動かしてみます。

% /lib/x86_64-linux-gnu/libc.so.6

同じ出力が得られました。

GNU C Library (Debian GLIBC 2.36-9+deb12u10) stable release version 2.36.
Copyright (C) 2022 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.
There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.
Compiled by GNU CC version 12.2.0.
libc ABIs: UNIQUE IFUNC ABSOLUTE
Minimum supported kernel: 3.2.0
For bug reporting instructions, please see:
<http://www.debian.org/Bugs/>.

加えて、各環境で sin(0x1.00000000004ADp+0) の出力も確認しました。 手元のものは最初に見ていたものを「手元 ①」、Dockerfile.GCC を使ったものを「手元 ②」として記載します。

環境	libc のバージョン	sin(0x1.00000000004ADp+0)
AtCoder	Ubuntu GLIBC 2.36-0ubuntu4	0x1.AED548F0911FBp-1
Library Checker	Debian GLIBC 2.36-9+deb12u10	0x1.AED548F0911FBp-1
手元 ①	Ubuntu GLIBC 2.41-6ubuntu1	0x1.AED548F0911FCp-1
手元 ②	Debian GLIBC 2.36-9+deb12u10	0x1.AED548F0911FCp-1

Library Checker と手元 ② は libc のバージョンは同じですが、sin の計算結果は異なっているため、libc のバージョン起因ではないと考えるのが妥当そうです。ということで仮説 1 は棄却です。

仮説 2

コードは sin_0x100000000004adp-52.c として下記を想定しています。

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

int main(void) {
  double x = strtod("1.0000000000002658", NULL);
  printf("%.13a\n", sin(x));
}

これを手元でコンパイルしたときと AtCoder 環境で実行したときとで、別のバイナリができているのではないか？というのを疑いました（各種都合で別のバイナリになることはあり得ると思うが、今回の計算に影響する部分に差異があるのでは？ということ）。そこで、手元でコンパイルしたものを AtCoder 環境で実行してみます。

% gcc -o sin_0x100000000004adp-52{,.c} -O2 -lm
% printf 'cat <<\\EOC | base64 -d | gzip -d > a.out\n%s\nEOC\nchmod u+x a.out\n./a.out\n' "$(cat sin_0x100000000004adp-52 | gzip -9 | base64)"

下記が得られました。

0x1.aed548f0911fbp-1

これは元々 AtCoder 環境でコンパイルしたときと同じ結果ですから、仮説 2 も棄却できます。

仮説 3

libc のバージョンは同じでも、ビルド時の設定（？）（./configure でなんかするやつ）によって差異があるのでは？と思いました。Library Checker と手元 ② で下記を実行します。

% sha256sum /lib/x86_64-linux-gnu/libc.so.6

どちらも同じ値になったので、同じファイルであると思われます*5。

1d25fd63234b59e4c581564c7a6d8f5c6cf36eee757e3d26f4b0808dd36a4896  /lib/x86_64-linux-gnu/libc.so.6

よって、仮説 3 も棄却です。libc に差異があるわけではないということですね。ここを最初に見ていれば仮説 1 を調べる必要はなかったかもですが、まぁそのへんは仕方ないでしょう。

仮説 4

途方に暮れつつ、gdb でステップ実行してみます。

(gdb) file ./sin_0x100000000004adp-52
(gdb) b sin@plt
(gdb) r
(gdb) si 2

__sin_sse2 (x=1.0000000000002658) at ../sysdeps/ieee754/dbl-64/s_sin.c:201

どうやら、__sin_sse2 というのが呼ばれているみたいです。AtCoder 環境では違うものが呼ばれているのかな？ 仮説 3 から libc の実体は共通であろうので、libc の中に __sin_sse2 や __sin_??? などがあり、実行時に何らかの基準によって適切と思われるものを呼んでいるのかな？という気持ちになってきます。

そういえば colima で AVX2 とかを有効にすることができたような...？というのを思い出したので、有効にしてみます。

% colima stop
% colima start --cpu 2 --memory 2 --disk 30 --cpu-type max,+avx,+avx2

からの、実行です。

% ./sin_0x100000000004adp-52
#> 0x1.aed548f0911fbp-1

おっと......？？ これは AtCoder 環境での出力と同じものですね。さっきと同様に gdb でも確認してみると、下記のものが呼ばれているようです。

__sin_fma (x=1.0000000000002658) at ../sysdeps/ieee754/dbl-64/s_sin.c:201

どうやら、__sin_fma と __sin_sse2 というのがあり、ハードウェアの機能に応じて別のものが呼ばれていそうです。FMA を使う場合とそうでない場合で計算結果は一般に一致しないので、計算結果が異なるのは納得できます。

ということで、どういう関数が用意されていて、どういう条件で誰が呼び分けているのか？というところが気になってきます。

調査

関数定義

glibc のコードを見ていきます。タグは最新の glibc-2.41.9000 を見ましたが、他もたぶんあまり変わらないでしょう*6。手元 ① の環境と合わせたバージョンですね。

まず、__sin_??? は下記のファイルで定義されています。

• __sin_sse2 in sysdeps/x86_64/fpu/multiarch/s_sin.c
• __sin_fma in sysdeps/x86_64/fpu/multiarch/s_sin-fma.c
• __sin_fma4 in sysdeps/x86_64/fpu/multiarch/s_sin-fma4.c
• __sin_avx in sysdeps/x86_64/fpu/multiarch/s_sin-avx.c

それぞれ、下記のような記述があります。

#define __sin __sin_fma
#define SECTION __attribute__ ((section (".text.fma")))
#include <sysdeps/ieee754/dbl-64/s_sin.c>

適宜端折りますが、sysdeps/ieee754/dbl-64/s_sin.c では次のようになっています。

double SECTION __sin (double x) { /* ... */ }
libm_alias_double (__sin, sin)

マクロを展開すると次のような感じですね。

double __attribute__ ((section (".text.fma"))) __sin_fma (double x) { /* ... */ }
libm_alias_double (__sin_fma, sin)

sysdeps/x86_64/fpu/multiarch/Makefile 内で下記のような記述があり、それぞれに合わせた最適化をしているのだろうと思っています。

CFLAGS-s_sin-fma.c = -mfma -mavx2
CFLAGS-s_sin-fma4.c = -mfma4
CFLAGS-s_sin-avx.c = -msse2avx -DSSE2AVX

また、sysdeps/x86_64/fpu/multiarch/s_sin.c に下記の記述があります。

extern double __redirect_sin (double);
# define SYMBOL_NAME sin
# include "ifunc-avx-fma4.h"
libc_ifunc_redirected (__redirect_sin, __sin, IFUNC_SELECTOR ());
libm_alias_double (__sin, sin)

libm_alias_double (__sin, sin) に関しては、sysdeps/generic/libm-alias-double.h などを見つつ、結局 weak_alias(__sin, sin) などが呼ばれていそうです。 libc_ifunc_redirected (...) に関しては、include/libc-symbols.h に説明のコメントが書かれています。“The following macros are used for indirect function symbols in libc.so.” から始まる、サンプルつきで 80 行程度のものです。

int foo (int __bar) という関数に対して、libc.so が使われる全部の環境で使えるとは限らない（実行時には判定できる）ハードウェア機能を用いた実装 int __foo_special (int __bar) と、フォールバックの実装 int __foo_default (int __bar) を提供する場合のことが書かれています。まさに〜？という感じがします。

sysdeps/generic/ifunc-init.h に、REDIRECT_NAME, OPTIMIZE, IFUNC_SELECTOR の #define があります。また、IFUNC_SELECTOR (sin_ifunc_selector？) の定義については sysdeps/x86_64/fpu/multiarch/ifunc-avx-fma4.h に記述があります。

CPU_FEATURE_USABLE_P (cpu_features, FMA) については、sysdeps/x86/include/cpu-features.h に定義があり、諸々あって

((ptr->features[index_cpu_FMA].active.reg_FMA & bit_cpu_FMA) != 0)

からの

((ptr->features[0].active.ecx & (1u << 12)) != 0)

に展開されていそうです。

__get_cpu_features () については、ざっくり飛ばしつつ、関連するファイルは下記たちです。

sysdeps/x86/include/cpu-features.h

#define __get_cpu_features() _dl_x86_get_cpu_features()

sysdeps/x86/dl-get-cpu-features.c

__ifunc (__x86_cpu_features, __x86_cpu_features, NULL, void,
     _dl_x86_init_cpu_features);

const struct cpu_features *
_dl_x86_get_cpu_features (void)
{
  return &GLRO(dl_x86_cpu_features);
}

sysdeps/generic/ldsodefs.h

# define GLRO(name) _##name

sysdeps/x86/cpu-features.c

static inline void
init_cpu_features (struct cpu_features *cpu_features)

実行時の処理

実行ファイルを実行するとき、実際には program interpreter と呼ばれるプログラムが実行されていたりします。./foo.sh を実行するときに shebang に書いたプログラムが実行されるのと似ている気がします？*7 なんのプログラムが実行されるかは、readelf -l の出力の Program Headers の INTERP の部分で確認することができます。

% readelf -l sin_0x100000000004adp-52
#:
#> Program Headers:
#>   Type           Offset             VirtAddr           PhysAddr
#>                  FileSiz            MemSiz              Flags  Align
#:
#>   INTERP         0x00000000000003a4 0x00000000000003a4 0x00000000000003a4
#>                  0x000000000000001c 0x000000000000001c  R      0x1
#>       [Requesting program interpreter: /lib64/ld-linux-x86-64.so.2]
#:

これはシンボリックリンクになっており、/usr/lib/x86_64-linux-gnu/ld-linux-x86-64.so.2 を指しています。ld.so も同様のファイルを指しているようなので、名前の短いこちらを使っていくことにします。実際、下記でも同じ出力が得られることが確かめられますね。

% ld.so ./sin_0x100000000004adp-52
#> 0x1.aed548f0911fbp-1

ということで、ld.so にあたる部分のコードを追っていきましょう（gdb を使いつつ調べています）。これも、コードとしては glibc と一緒に提供されています。

sysdeps/x86_64/dl-machine.h の RTLD_START の部分（関数としては _start にあたる）に下記のコメントがあります。

Initial entry point code for the dynamic linker. The C function ‘_dl_start’ is the real entry point; its return value is the user program's entry point.

_dl_start は elf/rtld.c にあります。ファイル冒頭に “Run time dynamic linker.” とあるので、名前はそういう意味だと思います。下記のような呼び出しがあり、CPU の各機能が使えるかが調べられます。

• _dl_start (elf/rtld.c:581)
• _dl_start_final (elf/rtld.c:496)
• _dl_sysdep_start (sysdeps/unix/sysv/linux/dl-sysdep.c:119)
• dl_platform_init (sysdeps/x86_64/dl-machine.h:209)
• _dl_x86_init_cpu_features (sysdeps/x86/dl-get-cpu-features.c:41)
• init_cpu_features (sysdeps/x86/cpu-features.c:754)

note: dl_start 内部で dl_start_final を呼び、その内部で dl_sysdep_start を呼び、... という感じです。

こうして前半パートで見た _dl_x86_cpu_features を初期化します。dl_sysdep_start 内部から、下記の呼び出しが続きます。

• _dl_start (elf/rtld.c:581)
• _dl_start_final (elf/rtld.c:496)
• _dl_sysdep_start (sysdeps/unix/sysv/linux/dl-sysdep.c:141)
• dl_main (elf/rtld.c:2267)
• _dl_relocate_object (elf/dl-reloc.c:346)
• _dl_relocate_object_no_relro (elf/dl-reloc.c:296)
• elf_dynamic_do_Rela (elf/do-rel.h:138)
• elf_machine_rela (sysdeps/x86_64/dl-machine.h:309)
  • value = ((ElfW(Addr) (*) (void)) value) ();

おそらく、libm.so.6 から取得できる sin は、

libc_ifunc_redirected (__redirect_sin, __sin, sin_ifunc_selector ());
libm_alias_double (__sin, sin)

のおかげで

extern __typeof (__redirect_sin) __sin;
__typeof (__redirect_sin) *__sin_ifunc (void) __asm__ ("__sin");
__attribute__ ((__optimize__ ("-fno-stack-protector")))
__typeof (__redirect_sin) *__sin_ifunc (void)
  {
    INIT_ARCH ();
    __typeof (__redirect_sin) *res = sin_ifunc_selector ();
    return res;
  }
__asm__ (".type " "__sin" ", %gnu_indirect_function");

extern __typeof (__sin) sin __attribute__ ((weak, alias ("__sin"))) __attribute_copy__ (__sin);

のようになっており、sin（という名前で公開されているやつ）を呼んで返ってきたアドレスを sin@got.plt に入れるという感じになっているように見えます (cf. 6.12.5 Controlling Names Used in Assembler Code)。

実際に sin@got.plt に __sin_fma が入るのは下記です。RELRO の設定によるかも。今回は full RELRO になっています。

• _dl_start (elf/rtld.c:581)
• _dl_start_final (elf/rtld.c:496)
• _dl_sysdep_start (sysdeps/unix/sysv/linux/dl-sysdep.c:141)
• dl_main (elf/rtld.c:2267)
• _dl_relocate_object (elf/dl-reloc.c:346)
• _dl_relocate_object_no_relro (elf/dl-reloc.c:296)
• elf_dynamic_do_Rela (elf/do-rel.h:138)
• elf_machine_rela (sysdeps/x86_64/dl-machine.h:418)
  • *(Elf64_Addr *) reloc_addr = (Elf64_Addr) value + reloc->r_addend;

こうして、条件に応じて __sin_sse2 や __sin_fma などが sin@got.plt にお届けされるわけですね。

おまけ

ld.so(8) や 38.6 Hardware Capability Tunables を眺めていると、glibc.cpu.hwcaps というのの存在に気づきます。環境変数 GLIBC_TUNABLES によって指定できるようです。

% ./sin_0x100000000004adp-52 
#> 0x1.aed548f0911fbp-1

% GLIBC_TUNABLES=glibc.cpu.hwcaps=-FMA,-AVX ./sin_0x100000000004adp-52
#> 0x1.aed548f0911fcp-1

あらあら、結果が変わることが確かめられましたね。オプション --list-diagnostics によって情報を確認できるので、差分も見てみましょう。

% diff <(ld.so --list-diagnostics) <(GLIBC_TUNABLES=glibc.cpu.hwcaps=-FMA,-AVX ld.so --list-diagnostics)
#:
#> 93c94
#> < x86.cpu_features.features[0x0].active[0x2]=0x7ed83203
#> ---
#> > x86.cpu_features.features[0x0].active[0x2]=0x6ed82203

たしかに、該当するビット (cf. sysdeps/x86/include/cpu-features.h) に影響があることがわかりました。

補足

GCC でも Clang でも、FMA が使える環境向けには（-ffast-math なしでも）x * y + z の形の式に対して FMA を使うような最適化をするようです。 今回触れた __sin_sse2 と __sin_fma で値が変わる理由として、本質的な部分はこの最適化のせいですね。 そんなことしていいのか...？

あとがき

いや〜満足満足。調査から記事まで、土日一回分で済んでうれしいです。他の進捗はありません。 GOT, PLT, RELRO などの説明を省きましたが、興味がある人は適宜調べてくださればと思います。CTF の pwn とかをすると詳しさが増すかも？という気もします。

C で呼び出す場合に限らず Rust や Python などの言語からでも（あるいはもちろん、GCC でコンパイルする場合に限らず Clang を使う場合でも）libm に依存していると影響を受ける部分なので、気をつけなきゃだなと思いました。

おわり

おわりです。