32bit乗算

Blackfinは16bit DSPであるため、ほとんどの演算は16bitで行われます。一方、たまに32bit演算を行いたくなることも当然あります。多くの演算は32bit版も用意されていますが、32bit版固定小数点乗算は用意されていません（整数版はある）。そこで、32bit乗算と32bit積和演算のプログラムを紹介します

乗算

ここでは1.31形式の固定小数点同士を掛け合わせて、1.31形式の積を得るプログラムを考えて見ます。

基本的には16bit乗算を複数回繰り返して筆算のように計算していきます。ただし、中間演算の中には符号付＊符号無しの乗算もありますので注意します。

a0 = r0.H * r1.H;     // H*H
a1 = r0.H * r1.L(M);  // H*L
a1 += r1.H * r0.L(M); // L*H

a1 = a1 >> 15;
a0 += a1;

最初の命令が乗算部です。32bitオペランドを16bitずつに分け、部分積を最初に計算しています。下位＊下位がありませんが、これは省略しています。その結果1LSBだけ誤差が増える可能性がありますが、対応すると非常に速度が落ちるため、目を瞑っています。

部分積のうち、上位＊下位は特殊な演算になります。この部分は符号付＊符号無しなので通常の演算では結果が誤りになります。そこで、このような場合に対応できる(M)オプションが乗算命令に用意されています。

上位＊下位はスケーリングが必要ですのでシフト演算を施し、最後に部分積を足し合わせて32bitの積を得ます。なお、乗算は並列実行可能ですので、束ねることで4サイクルにする事ができます。

a1 = r0.H * r1.L(M), a0 = r0.H * r1.H;  // H*L, H*H
a1 += r1.H * r0.L(M);                   // L*H

a1 = a1 >> 15;
a0 += a1;

乗算には4サイクルかかります。

積和演算

乗算には4サイクルと言う比較的長い時間がかかりますが、積和のサイクルはもっと短くできます。まず、以下のようにループを組んでベクトルの内積プログラムを作ります。

a1 = a0 = 0;

loop mac32 lc0=p0;
loop_begin mac32;
    r0 = [i0++] || r1 = [i1++];
　  a1 += r0.H * r1.L(M), a0 += r0.H * r1.H;
    a1 += r1.H * r0.L(M);
loop_end mac32;
      
a1 = a1 >> 15;
a0 += a1;

積の計算に必要な中間積の和は、最後に回す事ができます。そこでループから追い出すことによって、ループ内部は部分積の計算だけに集中できました。これをさらにパイプライン化すると、

a1 = a0 = 0;
p0 += -1;

    r0 = [i0++] || r1 = [i1++];

loop mac32 lc0=p0;
loop_begin mac32;
　  a1 += r0.H * r1.L(M), a0 += r0.H * r1.H;
    a1 += r1.H * r0.L(M)                     || r0 = [i0++] || r1 = [i1++];
loop_end mac32;

                                             a1 += r0.H * r1.L(M), a0 += r0.H * r1.H;
                                             a1 += r1.H * r0.L(M);
      
a1 = a1 >> 15;
a0 += a1;

となり、32bitの積和あたり2サイクルで済みます。