Rustで逆三角/三角関数の実行時間を測る

※この記事は，CAMPHOR- Advent Calendar 2023の14日目の記事になります．
CAMPHOR- Advent Calendar 2023についてはこちら：https://advent.camph.net

# この記事について
ロボット制御で三角関数を計算する必要がでてきたので、何回もやる処理なだけにできるだけ実行速度を短くしようと考えました。
また、今の自分の環境では組み込み関数として三角関数が提供されていませんでした。
そこで、今後の開発の足がかりにするため、今勉強しているRustで三角関数と逆三角関数の実装を行い実行時間の計測をやってみようと思いたち、メモ代わりに記事を書きました。

# 三角関数の表現
三角関数$\sin x, \cos x, \tan x$を計算したいわけですが、普通は近似式を使います。第多数の人が真っ先に思いつくのはテイラー展開だと思います。この記事では$\sin x$を計算したいとします。
$$
\sin x = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!} - \cdots
$$
これを単純にGeogebraで描画すると、以下のようになる。

![Screenshot from 2023-12-13 22-07-37.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3628897/38abe1fd-b58a-bad0-9561-e357f075641b.png)

緑色の先は誤差の絶対値を1000倍したものです。
一周気分が計算できれば、$\sin x$の周期性を利用して他の角度での値を求めることができます。
今回は、テイラー展開(9次まで)を愚直に計算することで検証してみたいと思います。
後々テーブルを使った実装も試します。

# 逆三角関数の表現
逆三角関数をベアメタルで使いたいというニッチすぎる需要が一般的に存在するのかはわかりませんが、いざというときに使えたら便利でしょう。
4節リンク機構の一つの角度を入力として別の角度がどうなるかを計算していたとき、なんとか逆三角関数を近似計算できないものかと思っていろいろ調べていたら、Wolfram alphaが面白いことを教えてくれました。

![Screenshot from 2023-12-13 22-41-59.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3628897/b0d575aa-b6e0-f289-2a5b-20dcb21ba497.png)

[ピュイズー級数](https://en.wikipedia.org/wiki/Puiseux_series "Puiseux series")なるものがあるとのこと。初めて知りましたが、どうやら分数の指数や負の指数を許した級数の一般化だそうです。

Geogebraで描画してみるとこのような感じになります。

![Screenshot from 2023-12-13 22-44-50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3628897/66f6d30f-e9d9-cb25-1833-8729d3cfef78.png)

例によって赤い線は誤差の100倍です。三角関数と比べると少し精度は劣りますがとても面白い数式ではないでしょうか。

そして、ピュイズー級数の中には平方根が出てきます。これを、[Newton-Raphson法](https://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95 "ニュートン法")で解いていきましょう。詳しくはWikiを参照してください。
```Rust
fn newton_raphson(x: f64) -> f64 {
    let eps = 0.0000001;
    let mut res = 1.0;
    let mut tmp = 0.0;
    let mut error = 1.0;
    while error > eps{
        // println!("{res}");
        tmp = res;
        res = res / 2.0 + (1.0 - x) / res;
        error = res - tmp;
        if error < 0.0 {
            error = -error;
        }
    };
    res
}
```

ちゃんと動いているかテストしてみましょう。
```Rust
fn main() {
    let ans = newton_raphson(0.0);
    println!("{ans}");
}
```

結果
```shell
1
1.5
1.4166666666666665
1.4142156862745097
1.4142135623746899
1.414213562373095
```
$\sqrt{2}$をいい感じに計算できています。


# 環境・前提
今回は、`std`にある組み込み関数との速度比較をやってみたいので、`#![no_std]`環境ではなく`std`ありでやっていきます。
時間があるときに記事の続きとして`#![no_std]`環境でやってみようと思います。

環境は以下のとおりです。
- Ubuntu 20.04
- rustup 1.26.0
- rustc 1.76.0-nightly
- criterion 0.5.1

余談ですが、rustcがnightlyなのは組み込み環境でのコンパイルに必要でインストールした名残です。

# テストのやりかた
`criterion`クレートと`cargo bench`を使います。
```toml
[dependencies]
rand = "*"
criterion = { version = "*", features = ["html_reports"] }
```
記事執筆時点（2023/12/14）では`criterion`クレートのバージョンは0.5.1でした。

そして、`Cargo.toml`には上の記述に加えて以下を書き足します。
```toml
[[bench]]
name = "main"
path = "<ファイルのディレクトリ>"
harness = false
```

コード本体は以下のとおりです。なお、ピュイズー級数については全体を-1倍して$x$を$-x$に置き換えて青い線を点対称に変換しました。
```Rust
extern crate criterion;

const PI: f64 = 3.14159265359;
use criterion::*;

fn newton_raphson(x: f64) -> f64 {
    let eps = 0.0000001;
    let mut res = 1.0;
    let mut tmp = 0.0;
    let mut error = 1.0;
    while error > eps{
        // println!("{res}");
        tmp = res;
        res = res / 2.0 + (1.0 - x) / res;
        error = res - tmp;
        if error < 0.0 {
            error = -error;
        }
    };
    res
}

fn sin_embedded(x: f64) -> f64 {
    f64::sin(x)
}

fn sin_taylor(x: f64) -> f64 {
    x * (1.0 - (x*x / 6.0) * (1.0 - (x*x / 20.0) * (1.0 - (x*x / 42.0) * (1.0 - x*x / 72.0))))
}

fn sin_taylor_2(x: f64) -> f64 {
    let x_sq = x * x;
    x * (1.0 - (x_sq / 6.0) * (1.0 - (x_sq / 20.0) * (1.0 - (x_sq / 42.0) * (1.0 - x_sq / 72.0))))
}

fn trig_arcsin(x: f64) -> f64 {
    f64::asin(x)
}

fn puiseux_powf(x: f64) -> f64 {
    PI / 2.0 - f64::sqrt(2.0 - 2.0 * x) * (1.0 - (x - 1.0) / 12.0 + 3.0 * (x - 1.0).powf(2.0) / 160.0 - 5.0 * (x - 1.0).powf(3.0) / 896.0 + (x - 1.0).powf(4.0) / 18432.0)
}

fn puiseux_powi(x: f64) -> f64 {
    PI / 2.0 - f64::sqrt(2.0 - 2.0 * x) * (1.0 - (x - 1.0) / 12.0 + 3.0 * (x - 1.0).powi(2) / 160.0 - 5.0 * (x - 1.0).powi(3) / 896.0 + (x - 1.0).powi(4) / 18432.0)
}

fn puiseux_horner(x: f64) -> f64 {
    PI / 2.0 - f64::sqrt(2.0 - 2.0 * x) * (1.0 - (x - 1.0) * (1.0 / 12.0 - (x - 1.0) * (3.0 / 160.0 - (x - 1.0) * (5.0 / 896.0 - (x - 1.0) * 35.0 / 18432.0)) ))
}

fn puiseux_horner_2(x: f64) -> f64 {
    let c = 1.0 - x;
    PI / 2.0 - f64::sqrt(2.0 * c) * (1.0 + c * (1.0 / 12.0 + c * (3.0 / 160.0 + c * (5.0 / 896.0 + c * 35.0 / 18432.0)) ))
}

fn bench_sin_embedded(c: &mut Criterion) {
    c.bench_function("sin_embedded", |b| b.iter(|| sin_embedded(black_box(0.5))));
}

fn bench_sin_taylor(c: &mut Criterion) {
    c.bench_function("sin_taylor", |b| b.iter(|| sin_taylor(black_box(0.5))));
}

fn bench_sin_taylor_2(c: &mut Criterion) {
    c.bench_function("sin_taylor_2", |b| b.iter(|| sin_taylor_2(black_box(0.5))));
}

fn bench_trig_arcsin(c: &mut Criterion) {
    c.bench_function("trig_arcsin", |b| b.iter(|| trig_arcsin(black_box(0.5))));
}

fn bench_puiseux_powf(c: &mut Criterion) {
    c.bench_function("puiseux_powf", |b| b.iter(|| puiseux_powf(black_box(0.5))));
}

fn bench_puiseux_powi(c: &mut Criterion) {
    c.bench_function("puiseux_powi", |b| b.iter(|| puiseux_powi(black_box(0.5))));
}

fn bench_puiseux_horner(c: &mut Criterion) {
    c.bench_function("puiseux_horner", |b| b.iter(|| puiseux_horner(black_box(0.5))));
}

fn bench_puiseux_horner_2(c: &mut Criterion) {
    c.bench_function("puiseux_horner_2", |b| b.iter(|| puiseux_horner_2(black_box(0.5))));
}


criterion_group!(benches, 
    bench_sin_embedded, 
    bench_sin_taylor, 
    bench_sin_taylor_2, 
    bench_trig_arcsin, 
    bench_puiseux_powf, 
    bench_puiseux_powi, 
    bench_puiseux_horner, 
    bench_puiseux_horner_2);

criterion_main!(benches);
```

コードのポイントは、
- 組み込み前提なので、`std::f64::PI`を用いず円周率を自前で定義した
- 多項式計算にHorner法を採用
- 何度も使う値を予め計算しない場合（無印）とした場合（`_2`を付加）を作り比較した
- `black_box`は、コンパイラが定数値を代入して予め結果を計算してしまうことを防ぐ
    - 参考：https://docs.rs/criterion/latest/criterion/fn.black_box.html
- 実装はcriterionの[user guide](https://bheisler.github.io/criterion.rs/book/getting_started.html "criterion user guide")を参照しました




コードが書けましたら、以下を実行します。
```shell
cargo bench
```

しばらく待つと、`target/criterion/report/index.html`にレポートが生成されています。


# テスト結果


<table>
<tr>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3628897/bf129e33-f2c4-e276-6c0b-5344bb51d940.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3628897/73da7c89-9116-d52e-d9ae-78b07774fe23.png"></td>
</tr>
<tr>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3628897/c3b67255-161d-c72b-833f-5a058b6ad39d.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3628897/857fa726-e742-a29b-10d2-d442fbce3d25.png"></td>
</tr>
<tr>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3628897/931122eb-82bf-9fa0-e5ed-609e948a817d.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3628897/6e21d815-db37-5eea-de04-5a308a62614a.png"></td>
</tr>
<tr>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3628897/7f0fd4d7-6e74-cbef-4a62-b9e8863b9e4e.png"></td>
<td><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3628897/07377b0c-a42b-8b5c-7212-64d88b08b6ec.png"></td>
</tr>
</table>

結果としては以下のとおりでしょう。
- 三角関数についても逆三角関数についても半分の実行時間で計算できている
- `powf`を使うと極端に遅い
- 自前の逆三角関数について、Honer法による実装と`powi`による実装であまり実行時間に差がない

なお、`cargo bench`を複数回実行すると、前回との差分をご丁寧にも用意してくれます。下の画像は、`trig_arcsin`について前回と今回の結果を比較しているものです。

![trig_arcsin_diff.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3628897/8dfd0a21-bb48-c8e9-557f-25931e389792.png)



# 不足点
- 精度に注目していないので、今後は精度を確保する実装を考えていきたい
    - 精度を担保すると今より遅くなるだろう
    - 場合分けで精度良く計算する方法がある
- テーブルを利用する方法・Newton-Raphson法を利用するやり方もあるのでこちらも要検討