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

本記事は情報処理学会第63回組込みシステム研究会の同名の発表の補足説明と詳細について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月現在)ので，参照ください．

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

C言語での実装

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

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

RVVでの実装

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

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

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

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

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

プロローグ

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

1. li t0, const_f
   • 浮動小数点数定数のメモリ領域の先頭アドレスをt0にロード
2. flw f0, (t0)
   • 浮動小数点数レジスタf0に単精度32ビット定数$2\pi$をロードする
   • (初期化関数によって(const_f)にあらかじめ$2\pi$が格納しておく．以下同様)
3. flw f1, 4(t0)
   • 浮動小数点数レジスタf1に単精度32ビット定数2.0 * MAX すなわち 2.0 * (1 << 30)をロードする．
   • この値は，浮動小数点数から固定小数点数への変換係数
4. li t0, const_i
   • 固定小数点数定数のメモリ領域の先頭アドレスをt0にロード
5. lw t3, (t0)
   • レジスタt3に$K$(k_value)の32ビット固定小数点数値をロードする
6. lw t4, 4(t0)
   • レジスタt4に$\pi$の32ビット固定小数点数値をロードする
7. 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に，外側のループ中で処理する要素数を格納する
2. vle32.v v0,(a1)
   • 1要素32ビット長として，ベクトルレジスタv0にレジスタa1(valuesへのポインタ，浮動小数点数表記による角度ベクトル$\theta_n$)で示されるメモリ領域を一括ロードする
   • ストライピングにより，実際にはv0, v1, v2, v3の4つのベクトルレジスタにロードされる(以下同様)
3. sub a0, a0, t0
   • レジスタa0(残りのベクトル長)からレジスタt0(外側のループ中で処理する要素数)を引いた値をレジスタa0に格納する
   • これにより，次のループに備える
4. 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に格納する
2. 次の3命令により，ベクトルレジスタv8にベクトルレジスタv0($\theta_n$)の小数部分が格納される
   1. vfcvt.x.f.v v8, v0
      • ベクトルレジスタv0の浮動小数点数値を符号付き整数に変換し，ベクトルレジスタv8に格納する
   2. vfcvt.f.x.v v8, v8
      • ベクトルレジスタv8の符号付き整数値を浮動小数点数に変換し，ベクトルレジスタv8に格納する
   3. vfsubvv v8, v0, v8
      • ベクトルレジスタv0の各要素についてベクトルレジスタv8の各要素の値を引いて，ベクトルレジスタv8に格納する
3. 次の2命令により，ベクトルレジスタv8に，$\frac{\theta_n}{2\pi}$の小数部分を固定小数点数表記したベクトルtheta_iが格納される
   1. vfmul.vf v8, v8, f1
      • ベクトルレジスタv8の各要素に浮動小数点数レジスタf1(浮動小数点数から固定小数点数への変換係数)を掛けた値をベクトルレジスタv8に格納する
   2. 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$が格納される
2. vmslt.vx v0, v8, t5
   • ベクトルレジスタv8(theta_i)の各要素をレジスタt5($-0.5\pi$)と比較し，より小さい場合には1，そうでない場合は0となるようなベクトル値を求めてベクトルレジスタv0に格納する
   • これによりベクトルレジスタv0はベクトルマスクレジスタとして機能するようになる
3. vmmv.m v16, v0
   • ベクトルマスクレジスタv0をベクトルレジスタv16にコピーする
4. vadd.vx v8, v8, t4, v0.t
   • ベクトルマスクレジスタv0の要素が1の要素だけ，ベクトルレジスタv8(theta_i)からレジスタt4($\pi$)を足す
5. sub t5, x0, t5
   • レジスタt5の符号を反転させる
   • これによりレジスタt5には$0.5\pi$が格納される
6. vmsgt.vx v0, v8, t5
   • ベクトルレジスタv8(theta_i)の各要素をレジスタt5($0.5\pi$)と比較し，より大きい場合には1，そうでない場合は0となるようなベクトル値を求めてベクトルレジスタv0に格納する
   • これによりベクトルレジスタv0はベクトルマスクレジスタとして機能するようになる
7. vsub.vx v8, v8, t4, v0.t
   • ベクトルマスクレジスタv0の要素が1の要素だけ，ベクトルレジスタv8(theta_i)からレジスタt4($\pi$)を引く
8. 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に格納する
2. 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を格納する
2. li t2, angles
   • レジスタt2に定数の静的メモリ領域$\phi_i$(angles)の先頭アドレスを格納する
3. 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:
2. lw a2, (t2)
   • レジスタa2にポインタレジスタt2($\phi_i$, angles)の値を格納する
3. vmslt.vx v0, v8, x0
   • ベクトルレジスタv8の各要素の値が負であるならば1，そうでなければ0であるようなベクトル値をベクトルレジスタv0に格納する
4. vmv.v.v v20, v24
   • ベクトルレジスタv20にベクトルレジスタv24(x)の値を格納する
5. vrsub.vx v20, v20, x0, v0.t
   • ベクトルマスクレジスタv0の各要素の値が1の要素のみ，レジスタx0からベクトルレジスタv20の各要素を引いた値をベクトルレジスタv20に格納する
   • つまり，theta_iに各要素の値が負であるような要素のみ，ベクトルレジスタv20の符号を反転する
6. 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)の値を足す
2. vmnot.m v0, v0
   • ベクトルマスクレジスタv0のNOT値をベクトルマスクレジスタv0に格納する
3. 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
2. vsra.vx v28, v28, t6

次の2つの命令で，変数x, yを更新します．

1. vadd.vv v24, v24, v28
2. vadd.vv v4, v4, v20

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

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

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

1. addi t6, t6, 1
2. addi t2, t2, 4
3. 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に復旧
2. vrsub.vx v24, v24, x0, v0.t
   • ベクトルマスクレジスタv0の各要素の値が1の要素のみ，レジスタx0からベクトルレジスタv24を引いた値をベクトルレジスタv24に格納する
   • すなわち，ベクトルマスクレジスタv0の各要素の値が1の要素のみ，ベクトルレジスタv24の符号を反転する
3. vrsub.vx v4, v4, x0, v0.t
   • ベクトルマスクレジスタv0の各要素の値が1の要素のみ，レジスタx0からベクトルレジスタv4を引いた値をベクトルレジスタv4に格納する
   • すなわち，ベクトルマスクレジスタv0の各要素の値が1の要素のみ，ベクトルレジスタv4の符号を反転する
4. vfcvt.f.x.v v4, v4
   • ベクトルレジスタv4(符号調整済$\sin\theta_n$, 固定小数点数表記)を浮動小数点数に変換する
5. addi t1, x0, 2
   • レジスタt1にレジスタx0に2を加えた値を格納する
   • すなわち，レジスタt1に2を格納する
6. fcvt.s.w f2, t1
   • 予稿誤り
   • 浮動小数点数レジスタf2にレジスタt1を浮動小数点数に変換した値を格納する
   • すなわち浮動小数点数レジスタf2に2.0を格納する
7. fdiv.s f2, f1, f2
   • 予稿誤り
   • 浮動小数点数レジスタf2に浮動小数点数レジスタf1からf2を割った値を格納する
   • すなわち，浮動小数点数レジスタf2に変換係数MAXを格納する
8. 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$(浮動小数点数表記)を書き込む
2. add a1, a1, t0
   • レジスタa1にレジスタt0に値を足す
   • すなわちvaluesのポインタ値を更新
3. bnez a0, loop
   • レジスタa0が0でなければloopにジャンプする
   • すなわち，残りのベクトルがある間，外側のループを継続する
4. ret
   • リターン

評価

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

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

• 提案手法(固定小数点数表記の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の助成を受けたものである．