※この記事は，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で描画すると、以下のようになる。

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

逆三角関数の表現

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

ピュイズー級数なるものがあるとのこと。初めて知りましたが、どうやら分数の指数や負の指数を許した級数の一般化だそうです。

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

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

そして、ピュイズー級数の中には平方根が出てきます。これを、Newton-Raphson法で解いていきましょう。詳しくはWikiを参照してください。

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 main() {
    let ans = newton_raphson(0.0);
    println!("{ans}");
}

結果

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を使います。

[dependencies]
rand = "*"
criterion = { version = "*", features = ["html_reports"] }

記事執筆時点（2023/12/14）ではcriterionクレートのバージョンは0.5.1でした。

そして、Cargo.tomlには上の記述に加えて以下を書き足します。

[[bench]]
name = "main"
path = "<ファイルのディレクトリ>"
harness = false

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

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を参照しました

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

cargo bench

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

テスト結果

結果としては以下のとおりでしょう。

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

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

不足点

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