RISC-Vベクタ拡張を用いたCORDICアルゴリズムによる三角関数の実装

本記事は[情報処理学会第63回組込みシステム研究会の同名の発表](https://researchmap.jp/zacky1972/misc/42825597)の補足説明と詳細についてQiita記事にしたものです．

## 本研究の動機

三角関数は，自動運転車やVR機器に必要な3D処理や，画像処理・信号処理で広く用いられるFFTなどで広く使われています．その典型的な使い方の一つは，複数の角度に対する$\sin$関数と$\cos$関数の値を同時に求める使い方です．

$\sin$関数と$\cos$関数を同時に求めるアルゴリズムとしてCORDICアルゴリズムがあります．これは下記論文に示すように，非常に古くからあるアルゴリズムで，初稿が1956年に遡ります．

* Volder, J. E.: Binary Computation Algorithms for Coordinate Rotation and Function Generation (1956). (internal report), Convair, Aeroelectronics group, IAR-1.148.
* Volder, J. E.: The Birth of CORDIC, Journal of VLSI Signal Processing, Vol. 25, No. 2 (Special issue on CORDIC), pp. 101–105 (online), DOI: 10.1023/A:1008110704586 (2000).  

本研究の動機としては，RISC-Vベクタ拡張(RVV)を用いて，複数の角度に対するCORDICアルゴリズムでの計算をベクトル化してみようというものです．RISC-Vはオープン標準である命令セットアーキテクチャ(ISA)であり，RVVはRISC-Vの目玉機能として規格の策定が進められて，2021年にバージョン1.0として批准されました．

組込みシステム研究会で発表しようと思い立ったのは，近年RISC-Vの組込みシステムでの採用実績が増えてきていることを踏まえてのことです．RISC-Vで標準的にサポートすべき拡張を定めているプロファイルの2023年版では，RVVを必須とする案が審議されている模様です．組み込みシステムにおいて，三角関数を扱う機会は比較的多く，前述のような自動運転車やVR機器は3D処理を行う組込みシステムですし，組込みシステムで信号処理や画像処理を行う機会は多くあります．したがって，表題のような研究は組込みシステム研究会の関心事なのではないかと考えております．

## CORDICアルゴリズム

CORDICアルゴリズムは，三角関数，双曲線関数，平方根，乗算，除算，指数，対数を効率よく計算するアルゴリズムです．三角関数におけるアルゴリズムの導出については，Wikipediaが詳しいです(2023年7月現在)ので，参照ください．

https://en.wikipedia.org/wiki/CORDIC

三角関数を求めるCORDICアルゴリズムは1956年にJack E. Volderによって提案されました．CORDICには，固定小数点数表記で表したときに，加減算とシフト演算のみで実装できるという特徴があります．

## C言語での実装

C言語で実装した固定小数点数表記のCORDICアルゴリズムに基づく$\sin$関数算出のアルゴリズムを次に示します．

```c
#include <stdint.h>
#include <stdbool.h>
#include <math.h>

#define M 32
#define MAX_BIT 30
#define MAX (1 << MAX_BIT)

static bool initialized = false;
static int32_t k_value;
static int32_t angles[M];

void vector_sin_cos_init()
{
    if(!initialized) {
        k_value = MAX;
        for(int i = 0; i < M; i++) {
            angles[i] = (int32_t)(atanf(ldexpf(1.0, -i)) / M_PI * MAX);
            k_value = (int32_t) (k_value / sqrt(1 + ldexp(1.0, -2.0 * i)));
        }
        initialized = true;
    }
}

struct vector_float {
    float x;
    float y;
};

static struct vector_float cordic_cos_sin(float theta)
{
    int32_t theta_i = (int32_t) (theta / M_PI * MAX);
    int32_t x = k_value;
    int32_t y = 0;

    for (int j = 0; j < M; j++) {
        int32_t xx, yy;
        if (theta_i < 0) {
            xx = x + (y >> j);
            yy = y - (x >> j);
            theta_i += angles[j];
        } else {
            xx = x - (y >> j);
            yy = y + (x >> j);
            theta_i -= angles[j];
        }
        x = xx;
        y = yy;
    }

    struct vector_float v;
    v.x = ((float)x / MAX);
    v.y = ((float)y / MAX);
    return v;
}

struct angle_sign_float {
    float angle;
    int sign;
};

static struct angle_sign_float adjust_angle(float theta)
{
    int t = (int)(theta / (2.0 * M_PI));
    theta -= (2.0 * M_PI) * t;
    int sign;
    if (theta < -M_PI * 0.5) {
        theta += M_PI;
        sign = -1;
    } else if (theta <= M_PI * 0.5) {
        sign = 1;
    } else {
        theta -= M_PI;
        sign = -1;
    }
    struct angle_sign_float result;
    result.angle = theta; // 予稿が間違っていました
    result.sign = sign;
    return result;
}

void vector_sin_float(int n, float *values)
{
    vector_sin_cos_init();
    for (int i = 0; i < n; i++) {
        struct angle_sign_float a = adjust_angle(values[i]);
        struct vector_float v = cordic_cos_sin(a.angle);
        values[i] = v.y * a.sign;
    }
}
```

関数は次の4つ定義しています．

* 初期化関数 `vector_sin_cos_init`
* CORDICアルゴリズム本体の関数 `cordic_cos_sin`
* 角度$\theta$を$-\pi \leq \theta \leq \pi$に収まるように変換する関数 `adjust_angle`
* ベクタのsin関数値を求めるドライバ関数`vector_sin_float`

関数`cordic_cos_sin`では$\cos$関数値と$\sin$関数値の両方を返すので，ドライバ関数を差し替えることで，$\cos$関数値を求めることもできますし，$\cos$関数と$\sin$関数の両方を求めることもできます．ここでは単純化のため，$\sin$関数値を求めるものとします．

関数`cordic_cos_sin`では次のような手順で求めます．

1. 浮動小数点数表記の`theta`を，πを最大値とする固定小数点数表記の`theta_i`に変換する．
1. 変数`x`,`y`の初期値をそれぞれ$K$, $0$の固定小数点数表記の値とする．
1. CORDICアルゴリズムにしたがって，加減算とシフト演算によって三角関数値を求める．
1. `x`には$\cos\theta$の固定小数点数値が，`y`には$\sin\theta$の固定小数点数値が，それぞれ格納される．
1. 最後に浮動小数点数表記に改めて，構造体`struct vector_float`に格納し，戻り値として返す．

## RVVでの実装

ではRVVでの実装を示します．この関数の引数は次のとおりです．

* レジスタ`a0`: `n`(角度ベクトル長)
* レジスタ`a1`: `values`へのポインタ(浮動小数点数表記による角度ベクトル$\theta_n$)

結果として，`values`のメモリの値が，$\sin\theta_n$に変換されます．

この実装でポイントになるのは，`if`文を**ベクトルマスクレジスタ**を用いて実装する点が重要です．このように実装することで，ループ中に条件分岐が現れず，制御ハザードが起きにくくなり，実行速度で大きな利点となります．

詳しくは下記Qiita記事を参照ください．

https://qiita.com/zacky1972/items/4ca65dfafbd2a1bc6913

### プロローグ

プロローグでは定数をレジスタにロードしていきます．

1. `li t0, const_f`
    * 浮動小数点数定数のメモリ領域の先頭アドレスを`t0`にロード
1. `flw f0, (t0)`
    * 浮動小数点数レジスタ`f0`に単精度32ビット定数$2\pi$をロードする
    * (初期化関数によって`(const_f)`にあらかじめ$2\pi$が格納しておく．以下同様)
1. `flw f1, 4(t0)` 
    * 浮動小数点数レジスタ`f1`に単精度32ビット定数`2.0 * MAX` すなわち `2.0 * (1 << 30)`をロードする．
    * この値は，浮動小数点数から固定小数点数への変換係数
1. `li t0, const_i`
    * 固定小数点数定数のメモリ領域の先頭アドレスを`t0`にロード
1. `lw t3, (t0)`
    * レジスタ`t3`に$K$(`k_value`)の32ビット固定小数点数値をロードする
1. `lw t4, 4(t0)`
    * レジスタ`t4`に$\pi$の32ビット固定小数点数値をロードする
1. `lw t5, 8(t0)`
    * レジスタ`t5`に$0.5\pi$の32ビット固定小数点数値をロードする

### 外側のループ

1. `loop:`

#### 角度ベクトル読み込み

外側のループの冒頭で角度ベクトルを読み込みます．

1. `vsetvli t0, a0, e32, m4`
    * ベクトルのコンフィギュレーション設定命令
    * `e32`は1要素が32ビット長であることを示す
    * `m4`はストライピングにより，4つのベクトルレジスタをグループ化してまとめることを意味する
        * たとえば`v0, v1, v2, v3`をまとめて1グループで`v0`として表す
        * ベクトルレジスタは32本あるので，`m4`により実質的に4倍の長さの8本のベクトルレジスタが備わっているかのように扱う
    * レジスタ`a0`は残りのベクトル長が格納されている
    * これらの情報をもとにレジスタ`t0`に，外側のループ中で処理する要素数を格納する
1. `vle32.v v0,(a1)`
    * 1要素32ビット長として，ベクトルレジスタ`v0`にレジスタ`a1`(`values`へのポインタ，浮動小数点数表記による角度ベクトル$\theta_n$)で示されるメモリ領域を一括ロードする
    * ストライピングにより，実際には`v0, v1, v2, v3`の4つのベクトルレジスタにロードされる(以下同様)
 1. `sub a0, a0, t0`
     * レジスタ`a0`(残りのベクトル長)からレジスタ`t0`(外側のループ中で処理する要素数)を引いた値をレジスタ`a0`に格納する
     * これにより，次のループに備える
 1. `slli t0, t0, 2`
     * レジスタ`t0`(外側のループ中で処理する要素数)を左に2つシフトした値，すなわち4倍した値をレジスタ`t0`に格納する
     * 1要素が32ビット長であることから4倍となっている
     * これによりレジスタ`t0`は，次のループのために`values`へのポインタに加算する値を格納していることになる

#### 角度の正規化と固定小数点数化

ここまでで次のように各レジスタに格納されています．

* ベクトルレジスタ`v0`: 外側のループで一度に処理する角度ベクトル$\theta_n$(浮動小数点数表記)
* 浮動小数点数レジスタ`f0`: $2\pi$(単精度32ビット値)
* 浮動小数点数レジスタ`f1`: 浮動小数点数から固定小数点数への変換係数`2.0 * MAX`(`2.0 * (1 << 30)`)(単精度32ビット値)

1. `vfdiv.vf v0, v0, f0`
    * ベクトルレジスタ`v0`($\theta_n$)の各要素を浮動小数点数レジスタ`f0`($2\pi$)で割った値をベクトルレジスタ`v0`に格納する
1. 次の3命令により，ベクトルレジスタ`v8`にベクトルレジスタ`v0`($\theta_n$)の小数部分が格納される
    1. `vfcvt.x.f.v v8, v0`
        * ベクトルレジスタ`v0`の浮動小数点数値を符号付き整数に変換し，ベクトルレジスタ`v8`に格納する
    1. `vfcvt.f.x.v v8, v8`
        * ベクトルレジスタ`v8`の符号付き整数値を浮動小数点数に変換し，ベクトルレジスタ`v8`に格納する
    1. `vfsubvv v8, v0, v8`
        * ベクトルレジスタ`v0`の各要素についてベクトルレジスタ`v8`の各要素の値を引いて，ベクトルレジスタ`v8`に格納する
1. 次の2命令により，ベクトルレジスタ`v8`に，$\frac{\theta_n}{2\pi}$の小数部分を固定小数点数表記したベクトル`theta_i`が格納される
    1. `vfmul.vf v8, v8, f1`
        * ベクトルレジスタ`v8`の各要素に浮動小数点数レジスタ`f1`(浮動小数点数から固定小数点数への変換係数)を掛けた値をベクトルレジスタ`v8`に格納する
    1. `vfcvt.x.f.v v8, v8`
        * ベクトルレジスタ`v8`の浮動小数点数値を符号付き整数に変換し，ベクトルレジスタ`v8`に格納する

#### 角度の正規化と結果符号`sign`の生成

ここまでで次のように各レジスタに格納されています．

* レジスタ`t5`: $0.5\pi$(32ビット固定小数点数表記)
* ベクトルレジスタ`v8`: 正規化した角度ベクトル`theta_i`(32ビット固定小数点数表記)
* レジスタ`t4`: $\pi$(32ビット固定小数点数表記)

1. `sub t5, x0, t5`
    * レジスタ`t5`の符号を反転させる
    * レジスタ`x0`は常に0
    * 直訳するとレジスタ`x0`からレジスタ`t5`を引いた値をレジスタ`t5`に格納する
    * これによりレジスタ`t5`には$-0.5\pi$が格納される
1. `vmslt.vx v0, v8, t5`
    * ベクトルレジスタ`v8`(`theta_i`)の各要素をレジスタ`t5`($-0.5\pi$)と比較し，より小さい場合には1，そうでない場合は0となるようなベクトル値を求めてベクトルレジスタ`v0`に格納する
    * これによりベクトルレジスタ`v0`はベクトルマスクレジスタとして機能するようになる
1. `vmmv.m v16, v0`
    * ベクトルマスクレジスタ`v0`をベクトルレジスタ`v16`にコピーする
1. `vadd.vx v8, v8, t4, v0.t`
    * ベクトルマスクレジスタ`v0`の要素が1の要素だけ，ベクトルレジスタ`v8`(`theta_i`)からレジスタ`t4`($\pi$)を足す
1. `sub t5, x0, t5`
    * レジスタ`t5`の符号を反転させる
    * これによりレジスタ`t5`には$0.5\pi$が格納される
1. `vmsgt.vx v0, v8, t5`
    * ベクトルレジスタ`v8`(`theta_i`)の各要素をレジスタ`t5`($0.5\pi$)と比較し，より大きい場合には1，そうでない場合は0となるようなベクトル値を求めてベクトルレジスタ`v0`に格納する
    * これによりベクトルレジスタ`v0`はベクトルマスクレジスタとして機能するようになる
1. `vsub.vx v8, v8, t4, v0.t`
    * ベクトルマスクレジスタ`v0`の要素が1の要素だけ，ベクトルレジスタ`v8`(`theta_i`)からレジスタ`t4`($\pi$)を引く
1. `vmor.mm v16, v0, v16`
    * ベクトルレジスタ`v16`と`v0`をベクトルマスクレジスタとみなしてOR演算を行い，ベクトルレジスタ`v16`に格納する
    * つまり，ベクトルレジスタ`v16`には，最後に$\cos$値と$\sin$値の符号を反転すべき場合には1，そうでない時には0になるようなベクトル値`sign`が格納されている

#### CORDICアルゴリズム本体

ここまでで次のように各レジスタに格納されています．

* ベクトルレジスタ`v8`: 正規化した角度ベクトル`theta_i`(32ビット固定小数点数表記)
* レジスタ`t3`: 定数$K$(`k_value`)

1. `vmv.v.x v24, t3`
    * レジスタ`t3`(定数$K$, `k_value`)の値を各要素に格納したベクトルをベクトルレジスタ`v24`に格納する
1. `vsub.vv v4, v4, v4`
    * 直訳すると，ベクトルレジスタ`v4`から`v4`を引いた値を`v4`に格納する
    * すなわち，ベクトルレジスタ`v4`に0を格納する

以降，次のような値として参照します．

* ベクトルレジスタ`v24`: 変数`x`(内側のループ完了時に$\cos\theta_n$が格納される)
* ベクトルレジスタ`v4`: 変数`y`(内側のループ完了時に$\sin\theta_n$が格納される)

次に内側のループの準備をします．

1. `li t1, 32`
    * レジスタ`t1`に内側のループのループ回数32を格納する
1. `li t2, angles`
    * レジスタ`t2`に定数の静的メモリ領域$\phi_i$(`angles`)の先頭アドレスを格納する
1. `mv t6, x0`
    * レジスタ`t6`に0を格納する

以降，次のような値として参照します．

* レジスタ`t1`: 内側のループのカウンタ値
* レジスタ`t2`: $\phi_i$(`angles`)を参照するときのポインタ値
* レジスタ`t6`: 内側のループの変数`j`

#### 内側のループ

ここまでで次のように各レジスタに格納されています．

* レジスタ`t2`: $\phi_i$(`angles`)を参照するときのポインタ値
* ベクトルレジスタ`v8`: 正規化した角度ベクトル`theta_i`(32ビット固定小数点数表記)
* ベクトルレジスタ`v24`: 変数`x`(内側のループ完了時に$\cos\theta_n$が格納される)
* ベクトルレジスタ`v4`: 変数`y`(内側のループ完了時に$\sin\theta_n$が格納される)

1. `in_loop:`
1. `lw a2, (t2)`
    * レジスタ`a2`にポインタレジスタ`t2`($\phi_i$, `angles`)の値を格納する
1. `vmslt.vx v0, v8, x0`
    * ベクトルレジスタ`v8`の各要素の値が負であるならば1，そうでなければ0であるようなベクトル値をベクトルレジスタ`v0`に格納する
1. `vmv.v.v v20, v24`
    * ベクトルレジスタ`v20`にベクトルレジスタ`v24`(`x`)の値を格納する
1. `vrsub.vx v20, v20, x0, v0.t`
    * ベクトルマスクレジスタ`v0`の各要素の値が1の要素のみ，レジスタ`x0`からベクトルレジスタ`v20`の各要素を引いた値をベクトルレジスタ`v20`に格納する
    * つまり，`theta_i`に各要素の値が負であるような要素のみ，ベクトルレジスタ`v20`の符号を反転する
1. `vmv.v.v v28, v4`
    * ベクトルレジスタ`v28`にベクトルレジスタ`v4`(`y`)の値を格納する

ここまでで次のように各レジスタに格納されています．

* レジスタ`a2`: $\phi_i$(`angles`)の値
* ベクトルレジスタ`v8`: 正規化した角度ベクトル`theta_i`(32ビット固定小数点数表記)
* ベクトルレジスタ`v0`: `theta_i`の各要素の値が負であるならば1，そうでなければ0であるようなベクトル値

1. `vadd.vx v8, v8, a2, v0.t`
    * ベクトルマスクレジスタ`v0`の各要素の値が1の要素のみ，ベクトルレジスタ`v8`にレジスタ`a2`の値を足す
    * すなわち，`theta_i`の各要素が負の要素のみ，ベクトルレジスタ`v8`(`theta_i`)に$\phi_i$(`angles`)の値を足す
1. `vmnot.m v0, v0`
    * ベクトルマスクレジスタ`v0`のNOT値をベクトルマスクレジスタ`v0`に格納する
1. `vsub.vx v8, v8, a2, v0.t`
    * ベクトルマスクレジスタ`v0`の各要素の値が1の要素のみ，ベクトルレジスタ`v8`にレジスタ`a2`の値を引く
    * すなわち，`theta_i`の各要素が正または0の要素のみ，ベクトルレジスタ`v8`(`theta_i`)に$\phi_i$(`angles`)の値を引く

以上で`theta_i`が更新されます．

1. `vrsub.vx v28, v28, x0, v0.t`
    * ベクトルマスクレジスタ`v0`の各要素の値が1の要素のみ，レジスタ`x0`からベクトルレジスタ`v28`を引いた値をベクトルレジスタ`v28`に格納する
    * すなわち，`theta_i`の各要素が正または0の要素のみ，ベクトルレジスタ`v28`の符号を反転する

ここまでで次のように各レジスタに格納されています．

* ベクトルレジスタ`v24`: 変数`x`(内側のループ完了時に$\cos\theta_n$が格納される)
* ベクトルレジスタ`v4`: 変数`y`(内側のループ完了時に$\sin\theta_n$が格納される)
* ベクトルレジスタ`v20`: CORDICアルゴリズムで定める条件に従って変数`x`の符号を反転させた値
* ベクトルレジスタ`v28`: CORDICアルゴリズムで定める条件に従って変数`y`の符号を反転させた値
* レジスタ`t6`: 内側のループの変数`j`

次の2つの命令で，変数`j`の値に従って，ベクトルレジスタ`v20`,`v28`の値を右シフト(`v20 >>= j; v28 >>= j;`)します．
1. `vsra.vx v20, v20, t6`
1. `vsra.vx v28, v28, t6`

次の2つの命令で，変数`x`, `y`を更新します．
1. `vadd.vv v24, v24, v28`
1. `vadd.vv v4, v4, v20`

ここまでで次のように各レジスタに格納されています．

* レジスタ`t1`: 内側のループのカウンタ値
* レジスタ`t2`: $\phi_i$(`angles`)を参照するときのポインタ値
* レジスタ`t6`: 内側のループの変数`j`

次の命令で各値を更新します．

1. `addi t6, t6, 1`
1. `addi t2, t2, 4`
1. `addi t1, t1, -1`

レジスタ`t1`(のループのカウンタ値)が0になるまで内側のループを継続します．

1. `bnez t1, inloop`

#### 結果の符号調整と浮動小数点数表記への変換

ここまでで次のように各レジスタに格納されています．

* ベクトルレジスタ`v16`: 結果の各要素を符号反転させるべき場合に1，そうでない場合に0が入っているベクトル値
* ベクトルレジスタ`v24`: $\cos\theta_n$(固定小数点数表記)
* ベクトルレジスタ`v4`: $\sin\theta_n$(固定小数点数表記)
* 浮動小数点数レジスタ`f1`: 浮動小数点数から固定小数点数への変換係数`2.0 * MAX`(`2.0 * (1 << 30)`)(単精度32ビット値)

1. `vmmv.m v0, v16`
    * ベクトルレジスタ`v16`に格納していたベクトルマスクレジスタ値をベクトルレジスタ`v0`に復旧
1. `vrsub.vx v24, v24, x0, v0.t`
    * ベクトルマスクレジスタ`v0`の各要素の値が1の要素のみ，レジスタ`x0`からベクトルレジスタ`v24`を引いた値をベクトルレジスタ`v24`に格納する
    * すなわち，ベクトルマスクレジスタ`v0`の各要素の値が1の要素のみ，ベクトルレジスタ`v24`の符号を反転する
1. `vrsub.vx v4, v4, x0, v0.t`
    * ベクトルマスクレジスタ`v0`の各要素の値が1の要素のみ，レジスタ`x0`からベクトルレジスタ`v4`を引いた値をベクトルレジスタ`v4`に格納する
    * すなわち，ベクトルマスクレジスタ`v0`の各要素の値が1の要素のみ，ベクトルレジスタ`v4`の符号を反転する
1. `vfcvt.f.x.v v4, v4`
    * ベクトルレジスタ`v4`(符号調整済$\sin\theta_n$, 固定小数点数表記)を浮動小数点数に変換する
1. `addi t1, x0, 2`
    * レジスタ`t1`にレジスタ`x0`に2を加えた値を格納する
    * すなわち，レジスタ`t1`に2を格納する
1. `fcvt.s.w f2, t1`
    * 予稿誤り
    * 浮動小数点数レジスタ`f2`にレジスタ`t1`を浮動小数点数に変換した値を格納する
    * すなわち浮動小数点数レジスタ`f2`に2.0を格納する
1. `fdiv.s f2, f1, f2`
    * 予稿誤り
    * 浮動小数点数レジスタ`f2`に浮動小数点数レジスタ`f1`から`f2`を割った値を格納する
    * すなわち，浮動小数点数レジスタ`f2`に変換係数`MAX`を格納する
1. `vfdiv.vf v4, v4, f2`
    * ベクトルレジスタ`v4`(符号調整済$\sin\theta_n$)を浮動小数点数レジスタ`f2`(浮動小数点数表記への変換係数)で割る

#### 結果の格納と外側のループの継続

ここまでで次のように各レジスタに格納されています．

* ベクトルレジスタ`v4`: $\sin\theta_n$(浮動小数点数表記)
* レジスタ`a1`: `values`へのポインタ
* レジスタ`t0`: 次のループのために`values`へのポインタに加算する値
* レジスタ`a0`: 残りのベクトル長

1. `vse32.v v4, (a1)`
    * ポインタレジスタ`a1`にベクトルレジスタ`v4`の値を書き込む
    * すなわち，`values`に$\sin\theta_n$(浮動小数点数表記)を書き込む
1. `add a1, a1, t0`
    * レジスタ`a1`にレジスタ`t0`に値を足す
    * すなわち`values`のポインタ値を更新
1. `bnez a0, loop`
    * レジスタ`a0`が0でなければ`loop`にジャンプする
    * すなわち，残りのベクトルがある間，外側のループを継続する
1. `ret`
    * リターン

## 評価

詳しくは今後発行される予稿を参照ください(下記でリンクを紹介する予定です)．

https://researchmap.jp/zacky1972/misc/42825597


ベクトルサイズと実行命令数の相関関係と相関係数は次のようになりました．

* 提案手法(固定小数点数表記のRVVハンドコーディング実装): $36.932n + 406.53, 0.9611$
* 固定小数点数表記のC実装，GCC 12.2 (`-O2`): $349.35n + 31.996, 1$
* 固定小数点数表記のC実装，Clang 16.3.2 (Auto-vectorization): $175.16n + 52.681, 0.9887$

提案手法の増分が他に比べてかなり小さいことから，現状ではRVVのアセンブリコードをハンドコーディングすることに意味があると考えています．

提案手法の定数項が想定よりかなり大きいのは，最適化オプションを指定していないためです．最適化オプションを指定しなかった理由は，現状の実装では最適化オプションを指定すると不具合を生じるためです．この原因は，アセンブリコードの生成にGCCのインラインアセンブラ機能を使用したのですが，おそらく申告すべき使用レジスタの指定を誤っているためではないかと考えられます．

## 謝辞

東京大学大学院の塩谷 亮太先生と東京理科大学の福原 淳司先生に助言をいただいた．

本研究はJSPS科研費 22K04657の助成を受けたものである．