# mpbench — 多倍長計算ライブラリ性能評価のコード一式

記事「多倍長数値計算の性能評価」で使ったベンチマークコードと測定データ。測定仕様は `common/SPEC.md`、共通乱数は `common/xorshift.md`。

```
benchmark_misc/
├── common/     共通仕様 (SPEC.md), 共通乱数 (xorshift.md, xorshift_ref.py)
├── cpp/        C/C++: FLINT (eig, fft), MPLAPACK (eig, gemm), Eigen+mpreal, MPFR+OpenMP FFT, スカラー演算, Makefile
├── julia/      Julia: GenericSchur, GenericLinearAlgebra, MultiFloats, DoubleFloats, Arblib, GenericFFT, FFTW
├── python/     Python: mpmath (+gmpy2), python-flint
├── matlab/     MATLAB + Advanpix: advanpix_bench.m (課題 A/B/C), advanpix_scaling.m (時間・メモリ vs N), advanpix_dense.m (密行列, MPLAPACK と同条件)
├── wolfram/    Mathematica (wolframscript)
├── scripts/    env_info.sh, run_all.sh, analyze.py, run_advanpix_scaling.sh, fit_advanpix.py
└── results/    results.csv (全測定値), env.txt (測定環境)
```

すべてのスクリプトは同一スキーマの CSV（`timestamp,host,library,version,language,task,subtask,digits,n,threads,trial,time_s,check_value,notes`）に追記し、`scripts/analyze.py` が表（`summary.md`）と図（`fig_*.png`）を生成する。

## 準備（Debian / Ubuntu 系の例）

### システムパッケージ

```
apt install build-essential libgmp-dev libmpfr-dev libmpc-dev libflint-dev libeigen3-dev libboost-dev libfftw3-dev python3-pip
```

`libflint-dev` は 3.0 以上（arb 統合後）が必要。QD（double-double / quad-double）は MPLAPACK 同梱のものを使うので `libqd-dev` は不要（下記）。

### mpreal.h（Eigen と scalar_bench で使用）

Pavel Holoborodko 氏の mpreal.h（http://www.holoborodko.com/pavel/mpfr/ 、GPL）を `cpp/mpreal/` に置く。

```
git clone https://github.com/advanpix/mpreal.git cpp/mpreal
```

### MPLAPACK 3.0.0

```
wget https://github.com/nakatamaho/mplapack/releases/download/v3.0.0/mplapack-3.0.0.tar.xz
tar xf mplapack-3.0.0.tar.xz && cd mplapack-3.0.0
./configure --prefix=$HOME/opt/mplapack --enable-gmp=yes --enable-mpfr=yes --enable-dd=yes --enable-qd=yes --enable-binary128=yes --enable-double=yes --enable-optimization=yes --enable-gmpfrxx-mkII --with-system-gmp=yes --with-system-mpfr=yes --with-system-mpc=yes --with-cxx-standard=17
make -j$(nproc) && make install
```

注意点:

- configure は 1 行で実行する（`\` 継続と行末コメントを混ぜると `#` が引数に渡って `invalid host type: #` になる）。
- `--with-system-qd` は付けない。dd/qd の複素型（`dd_complex`, `qd_complex`）は MPLAPACK 同梱の QD にしか無く、システムの libqd では `qd/qd_complex.h` が見つからずビルドが失敗する。
- `--enable-optimization` で OpenMP 版の `libmplapack_<type>_opt.so` ができる。3.0 では BLAS と LAPACK が一つのライブラリにまとまっており（`libmpblas_*` は無い）、`cpp/Makefile` はそれに合わせてある（`MPL_OPT=` を付けると参照版にリンク）。
- gmpfrxx_mkII の GMP バックエンドは `mpf_set_default_prec` を無視する。精度は `gmpxx::set_default_mpf_precision_bits(bits)`（MPFR は `mpfrxx::set_default_precision_bits`）で設定する。`cpp/mplapack_eig.cpp` はそのようにしてある。
- 型名は `mpfr_class` / `mpc_class`（MPFR）、`mpf_class` / `mpfc_class`（GMP）、`mplapack_binary128_t`（binary128、ヘッダは `mpblas_binary128.h`）。

### Julia パッケージ

```
julia -e 'using Pkg; Pkg.add(["GenericSchur","GenericLinearAlgebra","MultiFloats","DoubleFloats","Arblib","GenericFFT","FFTW"])'
```

MultiFloats は `MultiFloat{Float64,2}` と `{Float64,4}` を使う（N=8 は算術が未定義のバージョンがある）。超越関数は `MultiFloats.use_bigfloat_transcendentals()` で MPFR にフォールバックさせている。

### Python パッケージ

```
pip3 install mpmath gmpy2 python-flint pandas matplotlib tabulate
```

（mpmath は gmpy2 があると自動でバックエンドに使う。python-flint は conda より pip の方が確実。tabulate は analyze.py の Markdown 表出力に必要）

### MATLAB / Mathematica

Advanpix のパスは `~/Documents/MATLAB/startup.m` で `addpath` しておく（ベンチマークの関数ファイルの先頭に書くと「Local function name must be different from the script name」になる）。`matlab` と `wolframscript` は PATH に無ければ `MATLAB=...`、`WOLFRAMSCRIPT=...` の環境変数で `run_all.sh` に渡す。

## ビルド

```
make -C cpp flint mpfrfft scalar eigen MPREAL_DIR=./mpreal HAVE_BOOST=1
make -C cpp mplapack mplapack_gemm MPLAPACK=$HOME/opt/mplapack
```

## 実行

```
bash scripts/env_info.sh          # 環境記録 → results/env_host>.txt
bash scripts/run_all.sh smoke     # 全スクリプトの疎通確認 (数分)
bash scripts/run_all.sh scalar    # 課題 A
bash scripts/run_all.sh fft       # 課題 C
bash scripts/run_all.sh eig       # 課題 B (長時間。下記の環境変数で分割する)
python3 scripts/analyze.py        # results/summary.md と fig_*.png
```

`run_all.sh` の環境変数（既定値）:

| 変数 | 意味 |
|---|---|
| `THREADS=1,2,4,8,16,32,64,<nproc>` | FLINT / Julia / Advanpix のスレッド数スイープ（Julia の eig と fft はスレッドを使わないので最初の段だけ） |
| `MPL_THREADS=1,8,<nproc>` | MPLAPACK のスレッド数（Rgeev/Cheev は並列 BLAS の恩恵をほぼ受けないので段数を絞ってある） |
| `N_EIG=10,20,50,...,1000` / `N_FFT=...` / `DIGITS=100,1000` / `TRIALS=3` | サイズ・桁数・試行回数 |
| `ONLY=mplapack,flint,julia,advanpix,seq` | 回すライブラリ群を絞る（`seq` = Eigen, GenericLinearAlgebra, mpmath, Mathematica） |
| `VERIFIED=` | FLINT の保証付き固有値（`eig_multiple`）を省く（1000 桁 N=300 で 3 時間かかる） |

例: `ONLY=advanpix THREADS=32,64 N_EIG=500,700,1000 bash scripts/run_all.sh eig`

### Advanpix の時間・メモリ vs N と見積り式

```
MATLAB=... HERM=1 bash scripts/run_advanpix_scaling.sh      # 1 点 1 プロセスで N=100〜2000 (既定 64 スレッド)
python3 scripts/fit_advanpix.py results/advanpix_scaling_<host>.csv   # フィットと見積り表 → results/advanpix_fit.md, fig_advanpix_*.png
```

### MPLAPACK の Rgemm / Rgehrd / Rgeev（密行列）と Advanpix の同条件比較

```
for t in 1 8 64; do OMP_NUM_THREADS=$t cpp/mplapack_gemm --digits 100,1000 --n 300,1000 --threads $t --trials 1 --out results/results_<host>.csv; done
matlab -batch "addpath('matlab'); advanpix_dense('digits',[100 1000],'n',[300 1000],'threads',[1 8 64],'out','results/results_<host>.csv')"
```

### 個別実行の例

```
OMP_NUM_THREADS=8 cpp/mplapack_eig_mpfr --digits 100 --n 100,200 --threads 8 --check --herm
cpp/flint_eig --digits 100,1000 --n 10,50,100 --threads 1 --verified
cpp/mpfr_fft_omp --digits 1000 --n 1048576 --threads 16
julia -t 1 julia/eig_bench.jl --libs genericschur,multifloats --digits 100 --n 100 --herm
python3 python/eig_bench.py --digits 100 --n 100 --verified --herm
matlab -batch "addpath('matlab'); advanpix_bench('task','eig','digits',100,'n',[100 200],'threads',8,'herm',true)"
wolframscript -file wolfram/bench.wls task=fft digits=100 n=65536
```

すべてのスクリプトは `--smoke`（MATLAB は `'smoke',true`、Wolfram は `smoke=1`）で最小設定の動作確認ができる。

## 測定上の注意（記事本文の「落とし穴」も参照）

- スカラー演算と gemm の入力は sqrt(u) で仮数を目標精度いっぱいに埋めている。倍精度からの変換値をそのまま使うと arb や MPFR が短い仮数で速く見える。
- Circus 行列は上 Hessenberg 形なので、対角化の時間に Hessenberg 化は含まれない。
- Advanpix は物理コア数を超えるスレッド（SMT）で大きく遅くなる。`maxNumCompThreads` を物理コア数にする。
- MPLAPACK の `_opt` 版は小さい行列でスレッド起動のオーバーヘッドが支配的になるので、`OMP_NUM_THREADS` を明示する（未設定だと全論理コアを使う）。
- Julia の GenericSchur と GenericLinearAlgebra は同時に読み込むと衝突するので別プロセスで走らせる（`run_all.sh` は分けてある）。