powf を使わずにべき乗を計算

質問

組み込み系で開発しているものです。
どうしても 実数のべき乗 を使わねばならないのですが、
powf があまりに遅くて話になりません。

条件
・sinf, cosf は使ってもよい（高速なライブラリがあるので）
・指数は符号なし整数のみ(0～10000)
・底は 0.0～1.0 のみ
・戻り値も 0.0～1.0 のみ
・精度は 16bit 位で十分
・計算しておいてテーブルから呼んでもよい

自作で powf を作りたいのですが、
３時間いろいろやってみてもうまくいきませんでした。
他のやり方は無いでしょうか？

回答No.7

jacta · ベストアンサー

> 対象ハードは ADSP-BF539 です。

そのプロセッサはFPUがなかったのでは？
そうであれば、主に時間が掛かっているのは浮動小数点演算の部分かと思います（ソフトエミュレーションなので）。

というわけで、#2さんのアルゴリズムで、固定小数点演算を行えばかなり改善するはずです。
固定小数点演算は、加算はそのままで構いませんし、乗算は16ビット×16ビットで計算してから、適当に丸めて16ビット右シフトすればOKです。精度を上げるには、途中の計算を16ビットではなく20ビットにしてもよいでしょう（40ビット整数型が使えるはずなので）。

お礼

回答ありがとうございます。

>乗算は16ビット×16ビットで計算してから、適当に丸めて16ビット右シフト
asuncion さんの解法と組み合わせると、
powf の 12～13倍 は早くなりました。
精度は小数第４位くらいまでなら大体あっています。

まだ高速化できるでしょうか？

// 0.0～1.0 を 0～65536 として計算する。
unsigned int powf2(unsigned int f, unsigned short n)
{
　　unsigned int ret=65536;

　　if(f==0) return 0;
　　if(n==0) return 65536;
　　if(n==1) return f;

　　for(;n;n>>=1)
　　{
　　　　if(n & 1)
　　　　{
　　　　　　ret *= f;
　　　　　　ret >>=16;
　　　　}
　　　　f *= f;
　　　　f >>= 16;
　　}

　　return ret;
}

回答No.11

jacta

> ちなみに inline をつけたら 562 に若干増えました？？？

命令キャッシュが小さいからでしょう。
また、データキャッシュのサイズを考えると、表引きが有効とは思えません。

最適化オプションをどうしているのか分かりませんが、速度を最適化するより、サイズを最適化する方が結果的に高速になる可能性があります。また、案外register指定子が効くことがありますので、これも試してみてください。

お礼

回答ありがとうございます。

register ですが、効果がありませんでした。
（シミュレータで擬似的に動かしているからかもしれません）

最後に。結果として powf より 12倍以上早く計算できました。
開発は現在の方針で進めようと思っています。
長々とお付き合いありがとうございました。

回答No.10

jacta

> まだ高速化できるでしょうか？

微々たるものでしょうが、nの型をunsigned shortではなくunsigned intにしてはどうでしょうか？
汎整数拡張によるオーバーヘッドがそれなりにありそうです。
あとは、for文ではなくdo文で書き直すとわずかに速くなります。

それ以上はアセンブリ言語に踏み込まないと難しそうです。

お礼

回答ありがとうございます。

>unsigned shortではなくunsigned int
消費サイクル数が 584 から 568 になりました。

>for文ではなくdo文で書き直す
消費サイクル数が 584 から 571 になりました。

２つを組み合わせて 584 から 555 になりました。
約 5% の高速化です。

ちなみに inline をつけたら 562 に若干増えました？？？

回答No.9

Yanch

＞例えば 32768 位のテーブルを作ります。
＞0.999 を底に、指数 0～32767 まで計算しておきます。
＞このテーブルを元に、底と指数からテーブルのインデックスを決めて取り出す、
＞ということは可能でしょうか？
＞（散々計算しましたが、インデックスの位置の決め方が分からず、断念しました）

テーブルを使用する場合の例です。
実際に実装する場合は、init関数を使用せずに、
constでテーブルを定義しておくとよいです。
----------------------------------------------------------------------
#include <stdio.h>
#include <math.h>

#define MYPOWF_TABLE_SIZE 32768

double g_mypowf_table999[MYPOWF_TABLE_SIZE];

void init_mypowf_table999()
{
　　int y;

　　for (y = 0; y < MYPOWF_TABLE_SIZE; y++)
　　{
　　　　g_mypowf_table999[y] = powf(0.999, y);
　　}
}

double select_my_powf_table999(int y)
{
　　return g_mypowf_table999[y];
}

int main(int argc, char *argv[])
{
　　double z;
　　int y;

　　/* テーブルを初期化 */
　　init_mypowf_table999();

　　for (y = 0; y < 10; y++)
　　{
　　　　/* y:指数 */
　　　　z = select_my_powf_table999(y);
　　　　printf("select_my_powf_table999(%d){%.15lf}\n", y, z);
　　}

　　return 0;
}

お礼

回答ありがとうございます。

select_my_powf_table999 ですが、
この関数に引数として 底、指数 を持たせて、
0.999 以外の底のべき乗を返すことはできないでしょうか？
（つまり、0.999 のテーブルから 0.999 以外の底の値を取得したいのです）

回答No.8

Yanch

＞・計算しておいてテーブルから呼んでもよい
テーブルを使用する場合、
必要になるテーブル領域を単純に計算すると、
精度は16bit×0x4000×10000＝327,680,000バイト
となりますから、
＞ROM/RAM あわせて100KBくらいまで使えます。
の方法では、少し無理がありそうです。

お礼

回答ありがとうございます。

例えば 32768 位のテーブルを作ります。
0.999 を底に、指数 0～32767 まで計算しておきます。
このテーブルを元に、底と指数からテーブルのインデックスを決めて取り出す、
ということは可能でしょうか？
（散々計算しましたが、インデックスの位置の決め方が分からず、断念しました）

回答No.6

asuncion

＞#3さん
＞（２）方法その２
＞底や、指数の範囲も決まっているみたいですので、
＞固定小数点で計算してから、浮動小数点に変換する方法も、よさそうですね。

＞・底は 0.0～1.0 のみ
という情報だけからは、底（の範囲）が決まっているとは言えなそうです。
きざみ幅などの情報がないからです。
0.9の場合があるかもしれませんし、0.99999999の場合があるかもしれません。
小数点以下の桁数がもっと多い場合のことも考える必要があるのかもしれません。

よって、固定小数点で計算するのはむずかしそうです。
範囲が決められないからです。

お礼

きざみ幅などの情報がないからです。
>ご指摘ありがとうございます。

刻み幅ですが、0.0～1.0 と書きましたが実際には
渡されるのは符号なし整数で 0x0000～0x4000 が範囲になります。
0x0000 を 0.0、0x4000 を 1.0 に直してから計算しています。

なので刻み幅は 1.0 / 0x4000 = 0.00006103515625 です。

回答No.5

jacta

#1です。

先ほどの補足要求の際に一緒にたずねるべきでしたが、まずはCPUとコンパイラを明確にしてください。メモリサイズに関する情報も必須です。
FPUの有無によっても、当然どうすべきかが変わってきます。

お礼

コンパイラは VisualDSP4.5 です。

対象ハードは ADSP-BF539 です。
http://www.analog.com/jp/embedded-processing-dsp/blackfin/processors/index.html
ROM/RAM あわせて100KBくらいまで使えます。

現在は設計段階で、冶具もないために
VisualDSPのシミュレータでサイクル数を比較しています。
最終目標は作ろうとしているモジュールの MIPS の概算を出して、
ソフトウェアとして破綻しないことを提示することです。
（現在サイクル数からMIPSを計算しています。あくまで概算なので）

回答No.4

asuncion

どれくらい性能がよいかはわかりません。

#include <stdio.h>

double mypowf(double a, int n)
{
double p;

if (a == 0.0) return 0.0;
if (a == 1.0) return 1.0;
if (n == 1) return a;
for (p = 1; n; n /= 2) {
if (n & 1) {
p *= a;
}
a *= a;
}
return p;
}

int main(void)
{
int i;

for (i = 0; i < 33; ++i) {
printf("%.15f\n", mypowf(.5, i));
}
return 0;
}

お礼

回答ありがとうございます。
powf より２倍早いです。

回答No.3

Yanch

（１）方法その１
＞・計算しておいてテーブルから呼んでもよい
とあるので、そのように実装してみては、如何でしょう。
速度的には速くなると思いますよ。

（２）方法その２
底や、指数の範囲も決まっているみたいですので、
固定小数点で計算してから、浮動小数点に変換する方法も、よさそうですね。
こちらの方法だと、（１）よりは若干遅くなりますが、消費メモリを大幅に
抑えられそうです。

お礼

一番早いのは多分テーブルから呼び出すことだと思いますが、
テーブルのどこから持ってくればいいのかがわかりません。

実際に表を組んで見ました。
底 0.95 で大体１になるのが 150乗 として テーブルサイズ 150 で、
あらかじめ表を作ってみました。
指数をインデックスとして、底 0.95、指数 35 が渡された場合、
テーブルの 35 番目を返します。

では、底 0.837、指数 23 だった場合、
どのインデックスを参照するればよいのでしょうか？

回答No.2

asuncion

＞・指数は符号なし整数のみ(0～10000)

この条件から、「指数の回数」だけ掛け算を繰り返すのがベースになろうかと思います。
ただし、単純に繰り返すのではなく、例えば指数が16ならば
（（（底の2乗）の2乗）の2乗）の2乗
とすることで、掛け算の回数を減らすことができます。
この考え方を拡張して、指数が2のべき乗（1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192）の
場合を記憶しておきます。
そして、1～10000の指数（0乗は自明なので省略）を2のべき乗数の和で表わすと、
劇的とはいかないかもしれませんが、単純なループで掛け算を繰り返すよりは
いくぶんかでも効率がよくなるかもしれません。

お礼

回答ありがとうございます。
組んでみましたが、powf より２倍は早い結果を得られました。

回答No.1

jacta

> ３時間いろいろやってみてもうまくいきませんでした。
> 他のやり方は無いでしょうか？

他のやり方といわれても、まずはそのうまくいかなかったものを出していただかないとどうしようもありません。