libmのpow()関数を高速化改造する - SPRの日記（*BSDとかLinuxとかの覚え書き）

FreeBSDのpow()関数とpowf()関数（べき乗）がとてつもなく動作が遅い事に気が付いた。
ひたすらpow()するだけのプログラムのFreeBSD版とLinux版を用意してFreeBSDで両方を動かしてみると、何故かLinux版の方が速い。
しかもFreeBSD版は倍の時間が掛かる。

/*
 * コンパイルは「gcc hoge.c -lm」とやる。
 */
#include 
#include 

int main()
{
    double i, j, k = 0;

    for (i = 0.0001; i <= 0.2; i += 0.0001)
        for (j = 0.0001; j <= 0.2; j += 0.0001)
            k += pow(i, j);

    printf("%f\n", k);
}

何で遅いのかソースファイルを探してみると、該当ソースは/usr/src/lib/msun/srcにあるe_pow.cとe_powf.cと分かる。
FreeBSDのpow()関数はC言語で書いてあるが、glibcのソースを見てみると、こっちはpow()関数をアセンブラで（浮動小数演算命令を使って）書いてある。
ここでかなり差が出てると思われるので、glibcのソースをFreeBSDにマージして高速化してみる。

以下作業手順

具体的なやり方

★環境

浮動小数演算機能付きx86系のFreeBSD。
FreeBSDのバージョン6.2、glibcはバージョン2.3.6のソースで試した。

★ファイルのコピー

glibcにあるe_pow.Sとe_powf.Sを/usr/src/lib/msun/i387にコピーする。
e_pow.Sとe_powf.Sはglibc-2.3.6/sysdeps/i386/fpuある。

★ソースを編集

コピーしたe_pow.Sとe_powf.Sをエディタで編集する。

「#include 」の下辺りに次の３行を付け加える。

#define ASM_TYPE_DIRECTIVE(name,typearg) .type name,typearg;
#define ASM_SIZE_DIRECTIVE(name) .size name,.-name
#define ALIGNARG(log2) log2

「ENTRY(__ieee754_pow)」、「ENTRY(__ieee754_powf)」をそれぞれ「ENTRY(pow)」、「ENTRY(powf)」にする。
これは必要なのか検証してないが、何となく回りのソースに合わせた。

一番下の「END(__ieee754_pow)」とEND「(__ieee754_powf)」を削る。

★Makefileの編集

/usr/src/lib/msun/i387/Makefile.incをエディタで編集する。
下の方に、「ARCH_SRCS+= e_pow.S e_powf.S」を加える。

# long double counterparts
ARCH_SRCS+= s_ceill.S s_copysignl.S s_floorl.S s_scalbnl.S s_truncl.S

# ここに付け加える
ARCH_SRCS+= e_pow.S e_powf.S

LDBL_PREC = 64  # XXX 64-bit format, but truncated to 53 bits

★libmをコンパイル＆インストール

下記の通りにコマンドを入れるとインストールされる。

# cd /usr/src/lib/msun
# make obj && make depend && make && make install

元に戻したいなら、Makefile.incを編集し直して、/usr/obj/usr/src/lib/msunを消してからコンパイルし直す。

★比較

Celeron 2.66GHzマシンで最初の方に書いてあるのプログラムを動かしたタイム。
上が改造前、下が改造後のタイム。
倍以上のスピードアップに成功してる。