Rustの3種のべき乗演算の速度差について

# Rustの3種のべき乗演算の速度差について

べき乗（冪乗）とは $b^e$ のように $b$ の $e$ 乗で表される値を差す。ここで $b$ を**底**（base）、$e$ を**冪指数**（exponent）と呼ぶ。Rust の標準ライブラリでは3種類の冪演算が提供されている。

1. **浮動小数点数どうしのべき乗**
   * 底、冪指数共に浮動小数点数型を取り、結果も浮動小数点数型で返す
   * 例：`2f64.powf(8f64)`
2. **浮動小数点数と整数のべき乗**
   * 底は浮動小数点数型、冪指数は整数型を取り、結果は浮動小数点数型で返す
   * 例：`2f64.powi(8i32)`
3. **整数どうしのべき乗**
   * 底、冪指数共に整数型を取り、結果も整数型で返す
   * 例：`2i64.pow(8u32)`

プラットフォーム（OS や CPU アーキテクチャ）によって、それぞれの演算速度に大きな違いがあることに気づいたので、いくつかの環境で試してみた。


## 試そうとしたきっかけ

きっかけとなったのは以下の記事。そこでは各プログラミング言語の計算速度を大まかに比べるために、ライプニッツ級数を用いて円周率の近似値を求めている。

- [色々な言語で計算速度を比較してみた](https://qiita.com/hanada/items/c91788bcac2a40f1bb05)

ライプニッツ級数は以下の式で示され、収束が遅いという特徴がある。

$$
\sum_{n=0}^{\infty}\frac{(-1)^n}{2n+1} = \frac{\pi}{4}
$$

元記事ではこの特徴を利用して各言語の計算速度を計測している。

ところが、元記事のコメントにあるように、OS だったり、プログラムのちょっとした書き方の違いによって選ばれる冪演算のタイプによって、プログラムの実行時間に意外に大きな差が現れることがわかった。

たとえば、いくつかの言語は、Linux x86_64 では整数どうしの冪演算が速い傾向にあり、macOS では浮動小数点数どうしの冪演算が速い傾向にあった。ライプニッツ級数では $(-1)^n$ を繰り返し計算するので、冪演算の速度の違いが、言語の速さの違いとして大きく現れてしまう。そこで、Rust ではどうなっているのか調べたくなったわけだ。

ところで、もしライプニッツ級数の計算を高速化したいなら、演算コストの高い冪演算を使用するべきではない。$(-1)^n$ は $n$ が偶数なら $1$ となり、奇数なら $-1$ となるので、以下のような低コストな演算で置き換えられる。

```rust
if n % 2 == 0 { 1.0 } else { -1.0 }
```

しかし本記事ではこのような高速化については考慮せず、元記事と同じ冪演算を用いたプログラムを用いて、冪演算のタイプごとの速度の違いを調べる。


## プログラムを書く

Rustでライプニッツ級数を計算するプログラムを実装し、ベンチマーク機能を使用して計算にかかった時間の中央値を測定する。なお現時点ではベンチマークに必要な test クレートが安定化されていないため、Rust Nightly を用いる。

```
$ rustup run nightly rustc --version --verbose
rustc 1.23.0-nightly (6160040d8 2017-11-18)
binary: rustc
commit-hash: 6160040d8547222e761ad876cbe3a48c9c90a5bf
commit-date: 2017-11-18
host: x86_64-apple-darwin
release: 1.23.0-nightly
LLVM version: 4.0
```

適当なディレクトリに移動して、`cargo new` で新しいバイナリプロジェクトを作る。

```
$ cargo new --bin leibniz_formula
$ cd leibniz_formula
```


### ライプニッツ級数を求める関数群

ライプニッツ級数を求めるための式を `leibniz_common()` 関数として定義する。どの冪演算を用いるかは `pow` 引数を通してクロージャで指定する。このクロージャは u32 型の引数を $n$ として取り、$(-1)^n$ を計算して答えを f64 型で返す。

```rust:src/main.rs
fn leibniz_common<F>(terms: u32, pow: &F) -> f64
    where F: Fn(u32) -> f64
{
    (0..(terms + 1))
        .map(|n| pow(n) / (2 * n + 1) as f64)
        .sum::<f64>() * 4.0
}
```

ここでは`map()` と `sun()` メソッドを使ってループを回したが、以下のように `for` を使って書いてもいい。Rustのゼロコスト抽象化により、どちらを使っても性能に違いはない。

```rust:src/main.rs
fn leibniz_common<F>(terms: u32, pow: &F) -> f64
    where F: Fn(u32) -> f64
{
    let mut s = 0.0;
    for n in 0..(terms + 1) {
        s+= pow(n) / (2 * n + 1) as f64;
    }
    s * 4.0
}
```

次に $(-1)^n$ を計算するクロージャを用意して `leibniz_common()` を呼び出す関数を定義する。Rustの標準ライブラリは3種類の冪演算を提供しているので、それに対応した3つの関数を定義した。

```rust:src/main.rs
// f64::powf - 底: f64, 冪指数: f64, 戻り値: f64
// LLVMが提供するllvm.pow.f64を使用
// https://github.com/rust-lang/rust/blob/f1ea23e2cc72/src/libstd/f64.rs#L353
pub fn leibniz_f64powf(terms: u32) -> f64 {
    leibniz_common(terms, &|n| (-1f64).powf(n as f64))
}

// f64::powi - 底: f64, 冪指数: i32, 戻り値: f64
// LLVMが提供するllvm.powi.f64を使用
// https://github.com/rust-lang/rust/blob/917260ebc28a/src/libcore/num/f64.rs#L231
pub fn leibniz_f64powi(terms: u32) -> f64 {
    leibniz_common(terms, &|n| (-1f64).powi(n as i32))
}

// i64::pow - 底: i64, 冪指数: u32, 戻り値: i64
// Rustプログラムとして実装されている
// https://github.com/rust-lang/rust/blob/dbeb5bf89021/src/libcore/num/mod.rs#L1090
pub fn leibniz_i64pow(terms: u32) -> f64 {
    leibniz_common(terms, &|n| (-1i64).pow(n) as f64)
}
```

- 1つ目の関数は、浮動小数点数どうしの冪演算を用いる。LLVM が提供する `llvm.pow.f64` が使われる。
- 2つ目の関数は、浮動小数点数と整数の冪演算を用いる。LLVM が提供する `llvm.powi.f64` が使われる。
- 3つ目の関数は、整数どうしの冪演算を用いる。Rust で書かれた関数として実装されている。


### ユニットテスト

念のためテストを書いておく。

```rust:src/rust.rs
#[cfg(test)]
mod test {
    use super::*;

    static TERMS: u32 = 1e8 as u32;
    static EXPECTED: &str = "3.1415926635893259";

    #[test]
    fn test_leibniz_f64powf() {
        assert_approx_eq_f64(EXPECTED, leibniz_f64powf(TERMS));
    }

    #[test]
    fn test_leibniz_f64powi() {
        assert_approx_eq_f64(EXPECTED, leibniz_f64powi(TERMS));
    }

    #[test]
    fn test_leibniz_i64pow() {
        assert_approx_eq_f64(EXPECTED, leibniz_i64pow(TERMS));
    }

    fn assert_approx_eq_f64(expected: &str, actual: f64) {
        assert_eq!(expected, format!("{:.16}", actual));
    }

}
```

試しにテストを実行しよう。ここまでならRustのStable版で実行できる。デバッグモードでは実行に時間がかかるので、リリースモードでビルドする。

```
$ cargo test --release
running 3 tests
test test::test_leibniz_f64powf ... ok
test test::test_leibniz_i64pow ... ok
test test::test_leibniz_f64powi ... ok

test result: ok. 3 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out
```

問題ないようだ。


### ベンチマーク

Rustのライブラリが提供するベンチマーク機能を使うために非安定な test クレートを使用する。feature でベンチマーク機能をオン・オフできるようにしておこう。以下の行を Cargo.toml に追加する。

```toml:Cargo.toml
[features]
benches = []
```

main.rs にベンチマークを行うモジュールを追加する。

```rust:src/main.rs
#![cfg_attr(feature = "benches", allow(unstable_features))]
#![cfg_attr(feature = "benches", feature(test))]

#[cfg(feature = "benches")]
extern crate test as bench;

#[cfg(feature = "benches")]
mod benches {
    use bench::Bencher;
    use super::*;

    // static TERMS: u32 = 1e8 as u32;
    static TERMS: u32 = 1e6 as u32;

    #[bench]
    fn bench_leibniz_f64powf(b: &mut Bencher) {
        b.iter(|| leibniz_f64powf(TERMS));
    }

    #[bench]
    fn bench_leibniz_f64powi(b: &mut Bencher) {
        b.iter(|| leibniz_f64powi(TERMS));
    }

    #[bench]
    fn bench_leibniz_i64pow(b: &mut Bencher) {
        b.iter(|| leibniz_i64pow(TERMS));
    }

}
```

このように速度を測定したい式をクロージャで包み、`Bencher` オブジェクトの `iter()` 関数に与えればよい。ベンチマークではこの式が繰り返し評価され、かかった時間の中央値とばらつきがナノ秒で表示される。

なお元記事ではループを $10^8$ 回まわしていたが、それだとベンチマークによる繰り返しの評価で時間がかかりすぎるので $10^6$ 回に減らしている。（現状ではベンチマークが評価する回数は変更できないようだ）

試しに実行してみよう。以下のように `RUSTFLAGS` で追加の最適化オプションを指定し、さらに、Rust の Nightly 版でベンチマークを実行する。

```
$ RUSTFLAGS='-C opt-level=3 -C target-cpu=native' \
  cargo +nightly bench --features benches
...
test benches::bench_leibniz_f64powf ... bench:   9,173,296 ns/iter (+/- 951,528)
test benches::bench_leibniz_f64powi ... bench:  50,389,351 ns/iter (+/- 2,325,417)
test benches::bench_leibniz_i64pow  ... bench:  19,426,510 ns/iter (+/- 1,753,125)
```

時間（ナノ秒）が短い方が高速。これは macOS での実行結果だが、ご覧の通り `f64::powf` が最速だった。この環境では浮動小数点数型の結果がほしい時だけでなく、整数型の結果がほしい時にも `f64::powf` を使った方がいいかもしれない。


### main関数

ベンチマークやユニットテストだけでなく、コマンドラインからもプログラムを実行できるようにする。コマンドライン引数を使って冪演算のタイプを指定できるようにした。

```rust:src/main.rs
fn main() {
    // 繰り返しの回数
    let terms = 1e8f64 as u32;

    // 第1引数を取得して、対応する関数を評価する。
    let mut pi = None;
    if let Some(t) = std::env::args().nth(1) {
        match t.as_str() {
            "f64powf" => pi = Some(leibniz_f64powf(terms)),
            "f64powi" => pi = Some(leibniz_f64powi(terms)),
            "i64pow" => pi = Some(leibniz_i64pow(terms)),
            _ => (),
        }        
    }

    // 答えが得られていたら表示する。そうでなければusageを表示する。
    if let Some(p) = pi {
        println!("{:.16}", p);
    } else {
        println!("Usage: {} f64powf | f64powi | i64pow",
                 std::env::args().nth(0).unwrap_or("leibniz_formula".to_string()));
        std::process::exit(4);
    }
}
```

以下のように追加の最適化オプションを指定し、さらに、Rust の Stable 版でビルドする。

```
$ RUSTFLAGS='-C opt-level=3 -C target-cpu=native' cargo build --release
```

実行してみよう。

```
$ ./target/release/leibniz_formula
Usage: ./target/release/leibniz_formula f64powf | f64powi | i64pow

$ ./target/release/leibniz_formula f64powf
3.1415926635893259
```


## 速度の計測

`cargo bench`コマンドでベンチマークを実行して、プラットフォームごとにどの冪演算が速いか調べてみよう。以下の3台のマシンを使用した。

- 第3世代 Core i7 搭載の Mac mini
- 第7世代 Core i5 搭載の PC
- 32ビット ARM Cortex-A8 搭載の Beagle Board xM


### Intel Core i7 3720QM @ 2.60GHz

Core i7-3720QM は2012年に発売された第3世代（Ivy Bridge）のCore i7 プロセッサ。まずこのマシン（Mac mini）で macOS、Linux、Windows、FreeBSD を稼働させて、速度を計測してみた。

- Rust バージョン：1.23.0-nightly (6160040d8 2017-11-18)
- 計算回数： $10^6$ 回

|         | macOS<br/>High Sierra | Arch Linux | Windows 10<br/>(MSVC) | FreeBSD 11   |
| :------ | ----------------: | ------------: | ------------: | ------------: |
| f64powf |    9,173,296 ns   | 44,584,977 ns |  7,984,275 ns | 12,565,797 ns |
| f64powi |   50,389,351 ns   | 50,300,094 ns | 50,241,521 ns | 52,415,875 ns |
| i64pow  |   19,426,510 ns   | 19,359,836 ns | 20,374,670 ns | 20,855,907 ns |

並べてみると Linux の f64powf のみ異様に遅くなっているようだ。なお Arch Linux 環境のカーネルバージョンは `4.13.12-1-ARCH`。それ以外の環境では、たとえ結果が整数でほしくても f64powf で計算した方がいいかもしれない。


### Intel Core i5 7200U @ 2.50GHz

2016年に発売された第7世代（Kaby Lake）の Core i5 プロセッサ。

- Rust バージョン：1.23.0-nightly (6160040d8 2017-11-18)
- 計算回数： $10^6$ 回

|         | Arch Linux    | FreeBSD 11    |
| :------ | ------------: | ------------: |
| f64powf | 47,014,258 ns |  8,572,597 ns |
| f64powi | 53,550,763 ns | 53,445,031 ns |
| i64pow  | 17,449,901 ns | 17,516,663 ns |

Core i7 3720QM と同じ傾向だ。


### 32ビット ARM Cortex-A8 @ 1.00GHz

他のアーキテクチャの CPU を試してみようと思ったが、手元には 32ビットの ARM と 64ビットの PowerPC G4 しかなかった。PowerPC は使っている人はもういなさそうなので、ARM のみ計測した。本当は 64ビットの ARM があるとよかったのだが。

使用したのは 2011 年に購入した [Beagle Board xM](http://beagleboard.org/beagleboard-xm) の Rev C モデル。テキサス・インスツルメンツ（TI）製のSoCを採用し、シングルコアの ARM Cortex-A8 プロセッサが搭載されている。OS は Arch Linux ARM を使用した。

- Rust バージョン：1.23.0-nightly (6160040d8 2017-11-18)
  * armv7-unknown-linux-gnueabihf
- 計算回数： $10^6$ 回

|         | Arch Linux ARM   |
| :------ | ------------:    |
| f64powf | 1,802,574,320 ns |
| f64powi |   794,286,023 ns |
| i64pow  |   646,267,843 ns |

f64powi が健闘している。Core i プロセッサでは f64powi は使いどころがない雰囲気だったが、ARM Cortex-A8 では活用できそう。

以下、環境についての情報。

```
$ rustup run nightly rustc --version --verbose
rustc 1.23.0-nightly (6160040d8 2017-11-18)
binary: rustc
commit-hash: 6160040d8547222e761ad876cbe3a48c9c90a5bf
commit-date: 2017-11-18
host: armv7-unknown-linux-gnueabihf
release: 1.23.0-nightly
LLVM version: 4.0

$ uname -srm
Linux 4.14.0-1-ARCH armv7l

$ cat /proc/cpuinfo 
processor	: 0
model name	: ARMv7 Processor rev 2 (v7l)
BogoMIPS	: 794.62
Features	: half thumb fastmult vfp edsp thumbee neon vfpv3 tls vfpd32 
CPU implementer	: 0x41
CPU architecture: 7
CPU variant	: 0x3
CPU part	: 0xc08
CPU revision	: 2

Hardware	: Generic OMAP36xx (Flattened Device Tree)
Revision	: 0000
Serial		: 0000000000000000
```


## まとめ

ライプニッツ級数を求めるプログラムを使用して Rust の標準ライブラリが提供している3種類の冪演算の性能を測定した。

- x86_64 系プロセッサでは、テストした OS の大半で浮動小数点数どうしの冪演算が最も高速だったが、Linuxだけは結果が異なった。
- 32ビット ARM Cortex-A8 では、整数どうしの冪演算と、浮動小数点数と整数の演算が速かった。

今回のプログラムでは、採用した冪演算が最適なものでなかった場合、最適なものを選択した時よりも実行時間が約2倍から5倍長くなってしまった。速度を重視するプログラムでは、本記事で使用した Rust のベンチマーク機能を使って対象のプラットフォームで実際に計測することをおすすめしたい。