Linux で Arm64 アセンブリプログラミング (09) 浮動小数点型命令

浮動小数点型命令(2017-01-28)

これまで整数型の命令を見てきましたが、ここから ARMv8 の特徴である浮動小数点型命令の解説を進めていきます。ARMv8 では 32bit の ARMと異なり、浮動小数点数演算用の命令がオプションではなく、ARMv8 の標準として組み込まれています。また扱う事のできる最大の整数も浮動小数点数も同じ 64 ビットのため、プログラミングが容易です。

浮動小数点数

まず、浮動小数点数について説明します。コンピュータで浮動小数点数を扱う場合は、浮動小数点数の正負を表す符号(sign)、浮動小数点数の数値自体（有効数字）を表す仮数部(fraction)とべき乗を示す指数部(exponent) に分けて考えます。次の図のように倍精度浮動小数点数は64ビット、単精度浮動小数点数は32ビット、半精度浮動小数点数は16ビットを使用します。数値の精度が決まる仮数部 (fraction) は52ビットですが、正規化した後に仮数部の先頭に来る1を省略可能なため、実質的には有効桁数は53ビットとなります。

倍精度浮動小数点数、単精度浮動小数点数、半精度浮動小数点数として次のような範囲の数値が表現可能です。

形式	最大	最小
倍精度浮動小数点数値	1.79769313486231571e+308	2.22507385850720138e–308
単精度浮動小数点数値	3.40282347e+38	1.17549435e–38
半精度浮動小数点数値(IEEE)	65504	0.00012201070785522461
半精度浮動小数点数値(代替)	131008	0.00012201070785522461

浮動小数点数 (Vector) レジスタ

ARMv8は、浮動小数点数レジスタとして128ビットを使って複数のベクトル の要素を同時に演算する事ができる Vn レジスタを32本持っています。 レジスタの名称を示す文字（V、Q、D、S、H、B）は、大文字と小文字のどちらも使用できます。 ここで n はレジスタの番号で 0 から 31 の値です。

倍精度浮動小数点数(64bit)なら2つ、単精度浮動小数点数(32bit)なら4つを格納して、 同時に演算することができます。

浮動小数点数レジスタとメモリ間の読み書きは、汎用レジスタとメモリ間と 同じように簡単に行なえます。 また浮動小数点数の比較結果がそのまま フラグレジスタ(NZCV)に反映されるため、汎用レジスタと同じ感覚で 扱うことができます。

「Arm64のレジスタ」のページと同じ表をこちらにもにも載せます。 整数用レジスタは64ビットですが、浮動小数点数レジスタは128ビットで64ビットの倍精度浮動小数点数を2つ格納することができます。

また128ビットのVn レジスタ 32本と同じ領域を共有して、 スカラー値(1つの浮動小数点数)を格納することもできます。レジスタの名称が 浮動小数点数の精度毎に決まっています。

ARMv8 では、1つの128ビットレジスタに小さいサイズのレジスタを詰め込むことはせず、 128ビットのうち下位の領域だけを使うため、小さいサイズのレジスタもそれぞれ 32 本までとなっています。

スカラー値用のレジスタとして、128ビットの Qn レジスタ、 倍精度浮動小数点数を格納する64ビットの Dn レジスタ、 単精度浮動小数点数を格納する 32 ビットの Sn レジスタ、 半精度浮動小数点数を格納する 16ビットの Hn レジスタ、 バイトデータを格納する Bn レジスタという名前でアクセスします。 読み出す場合は上位の未使用のビットは無視され、書き込む場合は上位の 未使用のビットには 0 がセットされます。

ここでも n はレジスタの番号で 0 から 31 の値です。

浮動小数点数を複数格納できるベクトルレジスタ (Vn) は、 次のようにデータのサイズと数をレジスタ名の後ろに指定することができます。 データサイズを示す文字は、大文字と小文字のどちらでも構いません。

名称	サイズとデータ数
Vn.8B	8 bit × 8
Vn.16B	8 bit × 16
Vn.4H	16 bit × 4
Vn.8H	16 bit × 8
Vn.2S	32 bit × 2
Vn.4S	32 bit × 4
Vn.1D	64 bit × 1
Vn.2D	64 bit × 2

ベクトルレジスタ (Vn) に格納された複数の浮動小数点数の1つを指定するには、 ゼロから始まるインデックスを付けて配列形式で行います。 「V13 レジスタを倍精度浮動小数点数2つ分に使って、その2つ目」を v13.2d[1] のように表現します。v13.d[1]のようにレジスタ中のデータ数を省略する事もできます。

浮動小数点演算を使ったサンプルコード

浮動小数点演算を使ったサンプルコードを示します。 命令の詳細は今後紹介する予定ですが、まずはアセンブラで行う倍精度の浮動小数点演算が、コンパクトに書ける雰囲気を感じてください。

処理の内容は、当サイトで何度も登場する4行 4列の正方行列の乗算です(32bitARM、各種言語、Java8)。4行 4列の正方行列 M は次のようなものです。

      ┏ m0  m4   m8  m12 ┓
      ┃                  ┃
      ┃ m1  m5   m9  m13 ┃
  M = ┃                  ┃
      ┃ m2  m6  m10  m14 ┃
      ┃                  ┃
      ┗ m3  m7  m11  m15 ┛

この行列 M の16個の要素をメモリ中に { m0, m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15 } の順に格納するとします。同様に、4行 4列の正方行列 A の要素を {a0 .. a15}、 行列 B の要素を {b0 .. b15}、行列 C の要素を {c0 .. c15} のように表すこととします。

行列 A が格納されているメモリの先頭アドレスを x1 レジスタに保存し、 行列 B を格納するメモリの先頭アドレスを x2 レジスタに入っている場合に、 行列 A と 行列 B を乗算して、結果として求められた行列を行列 C として x0 レジスタに格納されている先頭アドレスのメモリ領域に保存するコードを示します。

//----------------------------------------------------------------------
//  file : mul4x4b.s
//  2017/01/27 Jun Mizutani (https://www.mztn.org/)
//----------------------------------------------------------------------

.ifndef __M4x4MUL
__M4x4MUL = 1

.text
//----------------------------------------------------------------------
// [c:x0] = [a:x1] * [b:x2]
// [a]  v1= m0, m1   v2= m2, m3   v3= m4, m5   v4= m6, m7
//      v5= m8, m9   v6=m10,m11   v7=m12,m13   v8=m14,m15
// [b] v11= m0, m1  v12= m2, m3  v13= m4, m5  v14= m6, m7
//     v15= m8, m9  v16=m10,m11  v17=m12,m13  v18=m14,m15
// [c] v21= m0, m1  v22= m2, m3  v23= m4, m5  v24= m6, m7
//     v25= m8, m9  v26=m10,m11  v27=m12,m13  v28=m14,m15
// M4x4= m0, m1, m2, m3, m4, m5, m6, m7, m8, m9,m10,m11,m12,m13,m14,m15
//        x,  0,  0,  0,  0,  y,  0,  0,  0,  0,  z,  0, tx, ty, tz,  1
//----------------------------------------------------------------------
mul4x4b:

        stp     x0, x30, [sp, #-16]!
        stp     x1, x2,  [sp, #-16]!
        stp     q1, q2,  [sp, #-32]!    // save v1.2d, v2.2d
        stp     q3, q4,  [sp, #-32]!    // save v3.2d, v4.2d
        stp     q5, q6,  [sp, #-32]!    // save v5.2d, v6.2d
        stp     q7, q8,  [sp, #-32]!    // save v7.2d, v8.2d
        stp     q11,q12, [sp, #-32]!    // save v11.2d, v12.2d
        stp     q13,q14, [sp, #-32]!    // save v13.2d, v14.2d
        stp     q15,q16, [sp, #-32]!    // save v15.2d, v16.2d
        stp     q17,q18, [sp, #-32]!    // save v17.2d, v18.2d
        stp     q21,q22, [sp, #-32]!    // save v21.2d, v22.2d
        stp     q23,q24, [sp, #-32]!    // save v23.2d, v24.2d
        stp     q25,q26, [sp, #-32]!    // save v25.2d, v26.2d
        stp     q27,q28, [sp, #-32]!    // save v27.2d, v28.2d

        ld1 {v1.2d, v2.2d, v3.2d, v4.2d},[x1],#64  // a0 .. a7
        ld1 {v5.2d, v6.2d, v7.2d, v8.2d},[x1]      // a8 .. a15
        ld1 {v11.2d,v12.2d,v13.2d,v14.2d},[x2],#64 // b0 .. b7
        ld1 {v15.2d,v16.2d,v17.2d,v18.2d},[x2]     // b8 .. b15

        fmul v21.2d, v1.2d, v11.2d[0] //  c.m0 =a.m0 * b.m0,  c.m1 =a.m1 *b.m0
        fmul v22.2d, v2.2d, v11.2d[0] //  c.m2 =a.m2 * b.m0,  c.m3 =a.m3 *b.m0
        fmul v23.2d, v1.2d, v13.2d[0] //  c.m4 = a.m0 *b.m4,  c.m5 =a.m1 * b.m4
        fmul v24.2d, v2.2d, v13.2d[0] //  c.m6 =a.m2 * b.m4,  c.m7 =a.m3 * b.m4
        fmul v25.2d, v1.2d, v15.2d[0] //  c.m8 =a.m0 * b.m8,  c.m9 =a.m1 * b.m8
        fmul v26.2d, v2.2d, v15.2d[0] // c.m10 =a.m2 *b.m12, c.m11 =a.m3 *b.m12
        fmul v27.2d, v1.2d, v17.2d[0] // c.m12 =a.m0 * b.m8, c.m13 =a.m1 * b.m8
        fmul v28.2d, v2.2d, v17.2d[0] // c.m14 =a.m2 *b.m12, c.m15 =a.m3 *b.m12

        fmla v21.2d, v3.2d, v11.2d[1] //  c.m0+=a.m4 * b.m1,  c.m1+=a.m5 *b.m1
        fmla v22.2d, v4.2d, v11.2d[1] //  c.m2+=a.m6 * b.m1,  c.m3+=a.m7 *b.m1
        fmla v23.2d, v3.2d, v13.2d[1] //  c.m4+= a.m4 *b.m5,  c.m5+=a.m5 * b.m5
        fmla v24.2d, v4.2d, v13.2d[1] //  c.m6+=a.m6 * b.m5,  c.m7+=a.m7 * b.m5
        fmla v25.2d, v3.2d, v15.2d[1] //  c.m8+=a.m4 * b.m9,  c.m9+=a.m5 * b.m9
        fmla v26.2d, v4.2d, v15.2d[1] // c.m10+=a.m6 *b.m13, c.m11+=a.m7 *b.m13
        fmla v27.2d, v3.2d, v17.2d[1] // c.m12+=a.m4 * b.m9, c.m13+=a.m5 * b.m9
        fmla v28.2d, v4.2d, v17.2d[1] // c.m14+=a.m6 *b.m13, c.m15+=a.m7 *b.m13

        fmla v21.2d, v5.2d, v12.2d[0] //  c.m0+=a.m8 * b.m2,  c.m1+=a.m9 *b.m2
        fmla v22.2d, v6.2d, v12.2d[0] //  c.m2+=a.m10* b.m2,  c.m3+=a.m11*b.m2
        fmla v23.2d, v5.2d, v14.2d[0] //  c.m4+= a.m8 *b.m6,  c.m5+=a.m9 * b.m6
        fmla v24.2d, v6.2d, v14.2d[0] //  c.m6+=a.m10* b.m6,  c.m7+=a.m11* b.m6
        fmla v25.2d, v5.2d, v16.2d[0] //  c.m8+=a.m8 *b.m10,  c.m9+=a.m9 *b.m10
        fmla v26.2d, v6.2d, v16.2d[0] // c.m10+=a.m10*b.m14, c.m11+=a.m11*b.m14
        fmla v27.2d, v5.2d, v18.2d[0] // c.m12+=a.m8 *b.m10, c.m13+=a.m9 *b.m10
        fmla v28.2d, v6.2d, v18.2d[0] // c.m14+=a.m10*b.m14, c.m15+=a.m11*b.m14

        fmla v21.2d, v7.2d, v12.2d[1] //  c.m0+=a.m12* b.m3,  c.m1+=a.m13*b.m3
        fmla v22.2d, v8.2d, v12.2d[1] //  c.m2+=a.m14* b.m3,  c.m3+=a.m15*b.m3
        fmla v23.2d, v7.2d, v14.2d[1] //  c.m4+= a.m12*b.m7,  c.m5+=a.m13* b.m7
        fmla v24.2d, v8.2d, v14.2d[1] //  c.m6+=a.m14* b.m7,  c.m7+=a.m15* b.m7
        fmla v25.2d, v7.2d, v16.2d[1] //  c.m8+=a.m12*b.m11,  c.m9+=a.m13*b.m11
        fmla v26.2d, v8.2d, v16.2d[1] // c.m10+=a.m14*b.m15, c.m11+=a.m15*b.m15
        fmla v27.2d, v7.2d, v18.2d[1] // c.m12+=a.m12*b.m11, c.m13+=a.m13*b.m11
        fmla v28.2d, v8.2d, v18.2d[1] // c.m14+=a.m14*b.m15, c.m15+=a.m15*b.m15

        st1     {v21.2d,v22.2d,v23.2d,v24.2d},[x0],#64
        st1     {v25.2d,v26.2d,v27.2d,v28.2d},[x0]
        ldp     q27, q28,  [sp], #32  // restore v27.2d, v28.2d
        ldp     q25, q26,  [sp], #32  // restore v25.2d, v26.2d
        ldp     q23, q24,  [sp], #32  // restore v23.2d, v24.2d
        ldp     q21, q22,  [sp], #32  // restore v21.2d, v22.2d
        ldp     q17, q18,  [sp], #32  // restore v17.2d, v18.2d
        ldp     q15, q16,  [sp], #32  // restore v15.2d, v16.2d
        ldp     q13, q14,  [sp], #32  // restore v13.2d, v14.2d
        ldp     q11, q12,  [sp], #32  // restore v11.2d, v12.2d
        ldp     q7, q8,  [sp], #32    // restore v7.2d, v8.2d
        ldp     q5, q6,  [sp], #32    // restore v5.2d, v6.2d
        ldp     q3, q4,  [sp], #32    // restore v3.2d, v4.2d
        ldp     q1, q2,  [sp], #32    // restore v1.2d, v2.2d
        ldp     x1, x2,  [sp], #16
        ldp     x0, x30, [sp], #16
        ret

.endif

コードの先頭2行で整数レジスタのX0とX1をスタック領域に退避して、 最後の3行でスタック領域に退避していた整数レジスタのX0とX1を復帰して 呼び出し元に戻ります。コードの先頭2行と最後の3行以外は、 すべて浮動小数点型の命令となります。

使用する浮動小数点レジスタを stp 命令でスタック領域に退避して、 ld1 命令でメモリから行列の要素を読み込み、fmul 命令で乗算、 fmla 命令で前の結果と乗算結果を加算、st1 命令で計算結果を メモリに書き出し、 ldp 命令でスタック領域に退避していたレジスタの 内容をもとに戻した後に ret 命令で呼び出し元に戻ります。

浮動小数点レジスタを stp 命令で退避して、ldp 命令で復帰していますが、 浮動小数点レジスタの内容を保存する必要がなければ、これらの24行は不要で、 さらに短くなります。

上記の行列の乗算ルーチンをテストしてみましょう。 D001 と D002 というラベルの位置に行列の値を置いて、乗算して結果を ラベル D003 の位置に格納します。 4x4行列の乗算を1000万回行った後に、結果となる行列の16個の要素の 倍精度浮動小数点数を print_d0 で文字列に変換して表示しています。 倍精度浮動小数点数を文字列に変換するサブルーチンのコードは double2string.s ですが、長くなるので、次回 紹介します。

青字はコメントです。

//-------------------------------------------------------------------------
//  file : multestb.s
//  2017/01/27 Jun Mizutani (https://www.mztn.org/)
//-------------------------------------------------------------------------

.global _start

.include        "mul4x4b.s"
.include        "double2string.s"

.text

_start:
        adr     x1, D001
        adr     x2, D002
        adr     x0, D003

        ldr     x4, NUM       // repetition 
1:      bl      mul4x4b
        subs    x4, x4, #1
        bne     1b

        mov     x3, #16       // print 4x4 matrix
2:      ldr     d0, [x0], #8
        bl      print_d0
        subs    x3, x3, #1
        bne     2b
        bl      Exit

print_d0:
        stp     x0, x30, [sp, #-16]!
        bl      Double2String
        bl      OutAsciiZ
        bl      NewLine
        ldp     x0, x30, [sp], #16
        ret

  //------------------------------------
NUM:    .quad    10000000      // number of repetition
D001:
        .double  1.111111111111E1
        .double  2.222222222222E2
        .double  3.333333333333E3
        .double  4.444444444444E4
        .double  5.555555555555E5
        .double  6.666666666666E6
        .double  7.777777777777E7
        .double  8.888888888888E8
        .double  9.999999999999E9
        .double  1.000000000000E10
        .double  1.111111111111E11
        .double  2.222222222222E12
        .double  3.333333333333E13
        .double  4.444444444444E14
        .double  5.555555555555E15
        .double  6.666666666666E16
D002:
        .double  1.111111111111E1
        .double  2.222222222222E2
        .double  3.333333333333E3
        .double  4.444444444444E4
        .double  5.555555555555E5
        .double  6.666666666666E6
        .double  7.777777777777E7
        .double  8.888888888888E8
        .double  9.999999999999E9
        .double  1.000000000000E10
        .double  1.111111111111E11
        .double  2.222222222222E12
        .double  3.333333333333E13
        .double  4.444444444444E14
        .double  5.555555555555E15
        .double  6.666666666666E16

.bss
        .align  2
D003:   .skip   128

以上のコード (multestb.s) を Orange Pi PC2 でアセンブル、リンクして実行しました。

jun@Orangepi:~/arm64asm0$ as -o multestb.o multestb.s
jun@Orangepi:~/arm64asm0$ ld -o multestb multestb.o
jun@Orangepi:~/arm64asm0$ time ./multestb
 1.481514814937975E18
 1.975311975456395E19
 2.469139506345185E20
 2.962970370567309E21
 2.963040741110519E22
 3.950625062172050E23
 4.938280247431114E24
 5.925943210467954E25
 7.407518519072603E25
 9.876554321652368E26
 1.234569135880003E28
 1.481483950705881E29
 2.222277778024247E30
 2.962968518814223E31
 3.703709876888149E32
 4.444456790517631E33

real  0m1.835s
user  0m1.800s
sys 0m0.000s

2500円のシングルボードコンピュータで 11億2千万回（16 x (4 + 3) x 10,000,000）の倍精度浮動小数点演算が 2秒弱でできています。 浮動小数点レジスタの退避と復帰を省略すると、さらに速くなり、１秒で実行できます。 以上のコードは次のページからダウンロードできます。

続く...

ご意見・ご感想は、こちらの掲示板まで。