浮動小数点数オタクがM1 Macを触ってみた

どうも、浮動小数点数オタクのmod_poppoです。

去年はMacのCPUがx86からArm (Apple Silicon)へ移行することが発表され、実際にApple M1搭載のハードウェアが発売されました。Apple Silicon MacにはRosetta 2というエミュレーターが搭載され、x86系のアプリケーションが（仮想化等の一部を除いて）そのまま動作します。

浮動小数点数オタクとして気になるのは、

- Armv8で規定された浮動小数点数関連の命令セット拡張のうち、Apple M1がどれを実装しているのか
- Rosetta 2がx86_64とAArch64の浮動小数点数に関する違いをどのように吸収しているのか

という点です。

# Armv8のあれこれ

## long doubleについて

long doubleは環境によって精度がまちまちなのは以前の記事に書いた通りです。

* [long doubleの話 - Qiita](https://qiita.com/mod_poppo/items/8860505f38e2997cd021)

Armの [Procedure Call Standard for the Arm 64-bit Architecture](https://developer.arm.com/documentation/ihi0055/d/?lang=en) ではlong doubleはbinary128（四倍精度）と定められていますが、Appleはこの規定を上書きして、long doubleはbinary64（倍精度、doubleと同じ）としました。

* [Writing ARM64 Code for Apple Platforms | Apple Developer Documentation](https://developer.apple.com/documentation/xcode/writing_arm64_code_for_apple_platforms)

M1 Mac上での筆者が作成したテストプログラム [floating-point-test/eval-method.c at master · minoki/floating-point-test](https://github.com/minoki/floating-point-test/blob/master/eval-method.c) の出力結果は以下のようになりました：

```
FLT_RADIX = 2
FLT_ROUNDS = 1; to nearest
FLT_EVAL_METHOD = 0; just to the range and precision of the type
__STDC_IEC_559__ is not defined
---
sizeof(float) = 4
FLT_MANT_DIG = 24
FLT_MIN_EXP = -125
FLT_MAX_EXP = 128
FLT_MAX = 0x1.fffffep+127
FLT_MIN = 0x1p-126
FLT_TRUE_MIN = 0x1p-149
float_t is float
float is probably IEEE binary32 format
---
sizeof(double) = 8
DBL_MANT_DIG = 53
DBL_MIN_EXP = -1021
DBL_MAX_EXP = 1024
DBL_MAX = 0x1.fffffffffffffp+1023
DBL_MIN = 0x1p-1022
DBL_TRUE_MIN = 0x1p-1074
double_t is double
double is probably IEEE binary64 format
---
sizeof(long double) = 8
LDBL_MANT_DIG = 53
LDBL_MIN_EXP = -1021
LDBL_MAX_EXP = 1024
LDBL_MAX = 0x1.fffffffffffffp+1023
LDBL_MIN = 0x1p-1022
LDBL_TRUE_MIN = 0x1p-1074
long double is probably IEEE binary64 format
---
0x1.00002fffp+0 * 0x1.000000008p+0 = 0x1.00002fff80001p+0; correctly rounded
0x1.fffe0effffffep-51 * 0x1.0000000000001p-1000 = 0x1.fffe0ep-1051; correctly rounded
three_mul_1(0x1p+1000, 0x1p+1000, 0x1p-1000) = inf
three_mul_2(0x1p+1000, 0x1p+1000, 0x1p-1000) = inf (variable)
three_mul_3(0x1p+1000, 0x1p+1000, 0x1p-1000) = inf (cast)
```

特に、 `__STDC_IEC_559__` は定義されません。

## 浮動小数点例外のトラップ

Armでは浮動小数点例外のトラップ（浮動小数点例外が発生した時に実行を止める機能）の実装は任意で、例えばRaspberry Pi 4のCPUでは浮動小数点例外のトラップは実装されていません。

一方、Apple M1は浮動小数点例外のトラップを実装しているようです。テスト用のコードを置いておきます：

```c
#include <stdio.h>
#include <fenv.h>
#include <stdint.h>

#pragma STDC FENV_ACCESS ON

// #if defined(__clang__)
// __builtin_aarch64_get_fpcr and __builtin_aarch64_set_fpcr are not available on clang
__attribute__((always_inline))
uint64_t get_fpcr()
{
    uint64_t fpcr;
    asm volatile("mrs %0, fpcr" : "=r"(fpcr));
    return fpcr;
}
__attribute__((always_inline))
void set_fpcr(uint64_t x)
{
    asm volatile("msr fpcr, %0" : : "r"(x));
}
// #endif

int main(void)
{
    uint64_t fpcr = get_fpcr();
    set_fpcr(fpcr | 0x9f00u); // IDE + IXE + UFE + OFE + DZE + IOE
    fpcr = get_fpcr();
    if (fpcr & 0x9f00u) {
        puts("Traps are supported.");
    } else {
        puts("Traps are not supported.");
    }
    volatile double zero = 0.0, one = 1.0;
    printf("1.0 / 0.0 = %g\n", one / zero);
}
```

実行例（Apple M1 / macOS）：

```
% clang trap-aarch64.c
trap-aarch64.c:5:14: warning: pragma STDC FENV_ACCESS ON is not supported, ignoring pragma [-Wunknown-pragmas]
#pragma STDC FENV_ACCESS ON
             ^
1 warning generated.
% ./a.out 
Traps are supported.
zsh: illegal hardware instruction  ./a.out
```

実行例（Raspberry Pi 4 / Linux）：

```
$ gcc trap-aarch64.c 
trap-aarch64.c:17:6: warning: always_inline function might not be inlinable [-Wattributes]
   17 | void set_fpcr(uint64_t x)
      |      ^~~~~~~~
trap-aarch64.c:10:10: warning: always_inline function might not be inlinable [-Wattributes]
   10 | uint64_t get_fpcr()
      |          ^~~~~~~~
$ ./a.out
Traps are not supported.
1.0 / 0.0 = inf
```

Apple M1ではFPCRを使ってトラップを有効にでき、実際にゼロ除算でプログラムが落ちています（SIGFPEではなくSIGILLなのが気になりますが）。一方、Raspberry Pi 4ではトラップを有効にできず、ゼロ除算では無限大が返されます。

ちなみに、GCCでは `__builtin_aarch64_get_fpcr` / `__builtin_aarch64_set_fpcr` という組み込み関数でFPCRにアクセスできますが、Clangではこれらはサポートされていないようなのでinline assemblyを使っています。

## Arm拡張に関して

Armv8/AArch64と言っても一枚岩ではなく、種々の拡張が存在します。先日紹介した「JavaScript専用命令」もその一つです。

* [ArmにあるというJavaScript専用命令とは何か、あるいは浮動小数点数を整数に変換する方法について - Qiita](https://qiita.com/mod_poppo/items/66663d9a05790579b5e4)

事前情報として、LLVMの定義ファイル [AArch64.td](https://github.com/llvm/llvm-project/blob/main/llvm/lib/Target/AArch64/AArch64.td) によるとApple A13はArmv8.4を実装しています。 ~~Apple A14/Apple M1の情報は現在のところAArch64.tdには載っていませんが、Armv8.4は実装していると考えるのが妥当でしょう。~~ 記事の仕上げをサボっている間に[Apple A14 に対する定義](https://github.com/llvm/llvm-project/blob/02bc4c95f0729cc819776f73ec94a25405579183/llvm/lib/Target/AArch64/AArch64.td#L871-L901)が追加されていました。Apple M1はA14がベースらしいのでA14で使える命令がそのままApple M1でも使えると考えて良いでしょう。

### sysctl -a の結果

まずは `sysctl -a` を実行してみます。

```shell
$ sysctl -a | grep -E '(cpu|hw.optional)'
kern.cpu_checkin_interval: 4000
hw.ncpu: 8
hw.activecpu: 8
hw.physicalcpu: 8
hw.physicalcpu_max: 8
hw.logicalcpu: 8
hw.logicalcpu_max: 8
hw.cputype: 16777228
hw.cpusubtype: 2
hw.cpu64bit_capable: 1
hw.cpufamily: 458787763
hw.cpusubfamily: 2
hw.optional.floatingpoint: 1
hw.optional.watchpoint: 4
hw.optional.breakpoint: 6
hw.optional.neon: 1
hw.optional.neon_hpfp: 1
hw.optional.neon_fp16: 1
hw.optional.armv8_1_atomics: 1
hw.optional.armv8_crc32: 1
hw.optional.armv8_2_fhm: 1
hw.optional.armv8_2_sha512: 1
hw.optional.armv8_2_sha3: 1
hw.optional.amx_version: 2
hw.optional.ucnormal_mem: 1
hw.optional.arm64: 1
machdep.cpu.cores_per_package: 8
machdep.cpu.core_count: 8
machdep.cpu.logical_per_package: 8
machdep.cpu.thread_count: 8
machdep.cpu.brand_string: Apple M1
```

SHA512やSHA3などを実装していることは読み取れますが、他はよくわかりません。

浮動小数点数関連の面白そうなArm拡張をいくつか挙げてみます：

- JCVT (Armv8.3)
	- 以前[記事](https://qiita.com/mod_poppo/items/66663d9a05790579b5e4)にした「JavaScript専用命令」です。
- FRINTTS (Armv8.5)
	- AArch64の浮動小数点数命令には浮動小数点数を固定長整数に変換する命令がすでに含まれますが、それでは飽き足らずArmv8.5では別系統の命令が追加されました。
	- 結果は浮動小数点数レジスターに格納されます。
	- 範囲外の場合の挙動が異なります。
	- 「何をする」命令なのかは調べれば色々出てきますが、「何のために」追加された命令なのかは全然出てきません。筆者的には、「x86の浮動小数点数→固定長整数の変換命令を模倣するため」に追加された可能性を考えています。
- SVE/SVE2
    - 対応してたら面白いですね。まあ現状では対応してないんですけど。
- BFloat16
	- 機械学習で使われるやつです。機械学習に関してはApple的にはCPUよりもGPUやNeural Engine等のアクセラレーターを使って欲しそうな気がする？

これらについて、実際に対応の有無を検証してみましょう。

### JCVT: 対応

詳しいことは

* [ArmにあるというJavaScript専用命令とは何か、あるいは浮動小数点数を整数に変換する方法について](https://qiita.com/mod_poppo/items/66663d9a05790579b5e4)

を参照してください。。特別なコンパイルオプションをつけなくても、組み込み関数 `__jcvt` が普通に使えます。

```
$ clang arm-jcvt.c 
$ ./a.out
(int32_t)(-2.8) = -2
(int32_t)1.99 = 1
(int32_t)(-Infinity) = -2147483648
(int32_t)Infinity = 2147483647
(int32_t)NaN = 0
(int32_t)(0x1p50 + 123.0) = 2147483647
__jcvt(-2.8) = -2
__jcvt(1.99) = 1
__jcvt(-Infinity) = 0
__jcvt(Infinity) = 0
__jcvt(NaN) = 0
__jcvt(0x1p50 + 123.0) = 123
```

### FRINT{32,64}[XZ]: 対応

本来ならFRINT系の命令に対応する組み込み関数 (ACLE) として `__rint{32,64}{x,z}` が使えるはずなのですが、現行のClangは対応していません。なので、インラインアセンブリーで呼び出します。

```c
#include <stdio.h>
#include <math.h>
#include <arm_acle.h>

double rint32z(double x)
{
  double result;
  asm("frint32z %d0, %d1" : "=w"(result) : "w"(x));
  return result;
}
double rint64x(double x)
{
  double result;
  asm("frint64x %d0, %d1" : "=w"(result) : "w"(x));
  return result;
}

int main(void)
{
#if defined(__ARM_FEATURE_FRINT)
  puts("__ARM_FEATURE_FRINT is defined.");
  printf("__rint32z(%g) = %a\n", INFINITY, __rint32z(INFINITY));
  printf("__rint64x(%g) = %a\n", NAN, __rint64x(NAN));
#else
  puts("__ARM_FEATURE_FRINT is not defined.");
#endif
  printf("rint32z(%g) = %a\n", INFINITY, rint32z(INFINITY));
  printf("rint64x(%g) = %a\n", NAN, rint64x(NAN));
}
```

ビルドの際には

```
$ clang -march=armv8.5-a arm-frint.c
```

と指定してやらないと「非対応の命令を使っているぞ」と怒られます。実行結果は

```
__ARM_FEATURE_FRINT is not defined.
rint32z(inf) = -0x1p+31
rint64x(nan) = -0x1p+63
```

でした。つまり、コンパイラーは対応していないけどCPU (M1) は対応している、という結果です。

### SVE: 非対応

AArch64の従来のSIMD (NEON) は幅が128ビットですが、SVEでは実装次第でもっと広い幅のレジスターを使えます。

Apple Siliconが実装したら面白いな〜〜と思いますが、Apple M1の時点では実装されていないようです。

テストプログラム：

```c
#include <stdio.h>
#include <arm_sve.h>

#if defined(__GNUC__)
#define NOINLINE __attribute__((noinline))
#else
#define NOINLINE
#endif

NOINLINE
void add(size_t n, double result[restrict n], const double x[restrict n], const double y[restrict n])
{
    for (size_t i = 0; i < n; ++i) {
        result[i] = x[i] + y[i];
    }
}

NOINLINE
void add_sve(size_t n, double result[restrict n], const double x[restrict n], const double y[restrict n])
{
    for (size_t i = 0; i < n; i += svcntd()) {
        svbool_t p = svwhilelt_b64((uint64_t)i, (uint64_t)n);
        svfloat64_t xx = svld1(p, &x[i]);
        svfloat64_t yy = svld1(p, &y[i]);
        svfloat64_t zz = svadd_z(p, xx, yy);
        svst1(p, &result[i], zz);
    }
}

int main()
{
    double x[] = {1.0, 2.0, 3.0};
    double y[] = {-0.5, 2.1, 3.14};
    double result[3];
    add(3, result, x, y);
    for (size_t i = 0; i < 3; ++i) {
        printf("result[%zu]=%g\n", i, result[i]);
        result[i] = 0.0;
    }
    puts("==SVE==");
    add_sve(3, result, x, y);
    for (size_t i = 0; i < 3; ++i) {
        printf("result[%zu]=%g\n", i, result[i]);
    }
}
```

実行結果（MacPortsのclang-11を利用）：

```shell
$ clang-mp-11 -O2 -march=armv8.4-a+sve arm-sve.c
$ ./a.out                                            
result[0]=0.5
result[1]=4.1
result[2]=6.14
==SVE==
zsh: illegal hardware instruction  ./a.out
```

というわけで、Apple M1はSVEに対応していません。

なお、Apple版のClang

```shell
$ clang --version
Apple clang version 12.0.0 (clang-1200.0.32.28)
Target: arm64-apple-darwin20.2.0
Thread model: posix
InstalledDir: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin
```

はSVEに非対応（`<arm_sve.h>` が用意されていない）でした。

### BFloat16: 非対応

Arm拡張としてのBF16はArmv8.2に対するもので、Armv8.6で必須となっています。

ここでは単精度との変換を行うBFCVT命令を試してみましょう。

```c
#include <stdio.h>
#include <arm_acle.h>
#include <arm_neon.h>

int main(void)
{
#if defined(__ARM_FEATURE_BF16)
    puts("__ARM_FEATURE_BF16 is defined");
#else
    puts("__ARM_FEATURE_BF16 is not defined");
#endif
#if defined(__ARM_BF16_FORMAT_ALTERNATIVE)
    puts("__ARM_BF16_FORMAT_ALTERNATIVE is defined");
#else
    puts("__ARM_BF16_FORMAT_ALTERNATIVE is not defined");
#endif
#if defined(__ARM_BF16_FORMAT_ALTERNATIVE)
    float32_t x = 3.14f;
    bfloat16_t y = vcvth_bf16_f32(x);
    float32_t z = vcvtah_f32_bf16(y);
    printf("%a -> %a\n", x, z);
#endif
}
```

次のようにコンパイル・実行します。

```shell
$ clang-mp-11 -march=armv8.6-a arm-bf16.c 
$ ./a.out
__ARM_FEATURE_BF16 is defined
__ARM_BF16_FORMAT_ALTERNATIVE is defined
zsh: illegal hardware instruction  ./a.out
```

というわけで、Apple M1（のCPU）はBF16に対応していないことがわかりました。

ちなみに上記のコードはQEMU (qemu-aarch64 version 4.2.1 (Debian 1:4.2-3ubuntu6.11)) でも動作しませんでした。

# Rosseta 2について

## Rosetta 2が実装している命令セット

すでにあちこちで報告されていますが、改めて書いておきます。

まずは `sysctl -a` で自己申告させてみましょう。

```shell
$ arch -arch x86_64 sysctl -a | grep -E '(cpu|hw.optional)'
kern.cpu_checkin_interval: 4000
hw.ncpu: 8
hw.activecpu: 8
hw.physicalcpu: 8
hw.physicalcpu_max: 8
hw.logicalcpu: 8
hw.logicalcpu_max: 8
hw.cputype: 7
hw.cpusubtype: 4
hw.cpu64bit_capable: 1
hw.cpufamily: 1463508716
hw.cpusubfamily: 0
hw.cpufrequency: 2400000000
hw.cpufrequency_min: 2400000000
hw.cpufrequency_max: 2400000000
hw.optional.floatingpoint: 1
hw.optional.mmx: 1
hw.optional.sse: 1
hw.optional.sse2: 1
hw.optional.sse3: 1
hw.optional.supplementalsse3: 1
hw.optional.sse4_1: 1
hw.optional.sse4_2: 1
hw.optional.x86_64: 1
hw.optional.aes: 1
hw.optional.avx1_0: 0
hw.optional.rdrand: 0
hw.optional.f16c: 0
hw.optional.enfstrg: 0
hw.optional.fma: 0
hw.optional.avx2_0: 0
hw.optional.bmi1: 0
hw.optional.bmi2: 0
hw.optional.rtm: 0
hw.optional.hle: 0
hw.optional.adx: 0
hw.optional.mpx: 0
hw.optional.sgx: 0
hw.optional.avx512f: 0
hw.optional.avx512cd: 0
hw.optional.avx512dq: 0
hw.optional.avx512bw: 0
hw.optional.avx512vl: 0
hw.optional.avx512ifma: 0
hw.optional.avx512vbmi: 0
hw.optional.watchpoint: 4
hw.optional.breakpoint: 6
hw.optional.neon: 1
hw.optional.neon_hpfp: 1
hw.optional.neon_fp16: 1
hw.optional.armv8_1_atomics: 1
hw.optional.armv8_crc32: 1
hw.optional.armv8_2_fhm: 1
hw.optional.armv8_2_sha512: 1
hw.optional.armv8_2_sha3: 1
hw.optional.amx_version: 2
hw.optional.ucnormal_mem: 1
hw.optional.arm64: 1
machdep.cpu.cores_per_package: 8
machdep.cpu.core_count: 8
machdep.cpu.logical_per_package: 8
machdep.cpu.thread_count: 8
machdep.cpu.brand_string: Apple M1
machdep.cpu.features: FPU VME DE PSE TSC MSR PAE MCE CX8 APIC SEP MTRR PGE MCA CMOV PAT PSE36 CLFSH DS ACPI MMX FXSR SSE SSE2 SS HTT TM PBE SSE3 PCLMULQDQ DTSE64 MON DSCPL VMX EST TM2 SSSE3 CX16 TPR PDCM SSE4.1 SSE4.2 AES SEGLIM64
machdep.cpu.feature_bits: 151121000215084031
machdep.cpu.family: 6
```

x86版の `sysctl -a` の出力にはArm拡張の情報も混ざっています。

SIMD幅が128ビットなのでAVX系が実装されていないのは良いとして、VEXエンコーディングを使うFMAも実装されていません。実際にFMAを含むコードを実行してみるとillegal instructionで落ちます。

## x86_64とAArch64の違いとRosetta 2

x86_64（のSSE系命令）とAArch64の浮動小数点数演算はいずれもIEEE 754という浮動小数点数の標準に基づいています。ただ、だからといって動作が完全に同じとは限りません。具体的には、ある種の浮動小数点例外の発生条件についてIEEEは実装による選択の余地を許しています。

その一つがアンダーフロー例外の発生条件、別の一つがFMA(0, inf, qNaN)によるinvalid例外の有無です。詳しくは筆者の同人誌「浮動小数点数小話」を読んで欲しいのですが、アンダーフロー例外の判定についてはx86は丸めの後に行い、AArch64は丸めの前に行います。

確認するためには以下のコードを実行します。

```c
#include <stdio.h>
#include <math.h>
#include <fenv.h>
#include <stdbool.h>

#pragma STDC FENV_ACCESS ON

void raise_invalid_on_fma_with_nan(void)
{
    volatile double x = 0.0;
    volatile double y = INFINITY;
    volatile double z = NAN;
    feclearexcept(FE_INVALID);
    volatile double w = fma(x, y, z);
    bool invalid = fetestexcept(FE_INVALID) != 0;
    printf("fma(%g, %g, %g), which yields %g, %s\n", x, y, z, w, invalid ? "raises INVALID" : "does not raise INVALID");
}

void underflow_after_rounding(void)
{
    volatile double x = 0x1.0000001p-1022;
    volatile double y = 0x0.fffffffp0;
    // x * y = 0x0.ffffffffffffffp-1022 -(rounding)-> 0x1.0000000000000p-1022
    feclearexcept(FE_UNDERFLOW);
    volatile double z = x * y;
    bool underflow = fetestexcept(FE_UNDERFLOW) != 0;
    printf("%a * %a, which yields %a, %s\n", x, y, z, underflow ? "raises UNDERFLOW; underflow is detected before rounding" : "does not raise UNDERFLOW; underflow is detected after rounding");
}

int main(void)
{
#if defined(FP_FAST_FMA)
    puts("FP_FAST_FMA is defined");
#else
    puts("FP_FAST_FMA is not defined");
#endif
    raise_invalid_on_fma_with_nan();
    underflow_after_rounding();
}
```

実行結果：

```shell
$ clang exception.c
exception.c:6:14: warning: pragma STDC FENV_ACCESS ON is not supported, ignoring pragma [-Wunknown-pragmas]
#pragma STDC FENV_ACCESS ON
             ^
1 warning generated.
$ ./a.out
FP_FAST_FMA is defined
fma(0, inf, nan), which yields nan, raises INVALID
0x1.0000001p-1022 * 0x1.ffffffep-1, which yields 0x1p-1022, raises UNDERFLOW; underflow is detected before rounding
$ clang -arch x86_64 exception.c
exception.c:6:14: warning: pragma STDC FENV_ACCESS ON is not supported, ignoring pragma [-Wunknown-pragmas]
#pragma STDC FENV_ACCESS ON
             ^
1 warning generated.
$ ./a.out                       
FP_FAST_FMA is not defined
fma(0, inf, nan), which yields nan, does not raise INVALID
0x1.0000001p-1022 * 0x1.ffffffep-1, which yields 0x1p-1022, does not raise UNDERFLOW; underflow is detected after rounding
```

というわけで、アンダーフローの発生条件についてはRosetta 2はきちんとx86の動作をエミュレートしているようです。

「Rosetta 2はAVXに対応しない」との事前情報から予測されたことですが、Rosetta 2ではx86_64のFMA命令は動きません。なのでFMAの0, inf, NaNの件は関係ありません。Rosetta 2上でmacOSのfma関数を呼び出してみた感じ、x86_64の挙動が真似られているようです。

FMA命令の実行結果：

```shell
$ clang -arch x86_64 -mfma -O2 exception.c
exception.c:6:14: warning: pragma STDC FENV_ACCESS ON is not supported, ignoring pragma [-Wunknown-pragmas]
#pragma STDC FENV_ACCESS ON
             ^
1 warning generated.
$ ./a.out                                 
FP_FAST_FMA is defined
zsh: illegal hardware instruction  ./a.out
```

という風に `-mfma` でFMA命令を使うようにするとSIGILLで落ちます。