立方根計算の高速化‐VC#2010EE

Question

立方根を求めるメソッドを定義したのですが、ニュートン法によるものですがSystem.Mat.Sqrt(double)より、約１５倍以上かかります。

よりより最適化は可能でしょうか？

よろしくお願いします。

以下ソース

namespace User
{
    public static class Math
    {
        public static double Cbrt(double d)
        {
            double rv = 1, rt = 1, b;
            //符号ルーチン
            bool sgn = d >= 0 ? true : false;
            if (sgn) b = d;
            else b = -d;
            //桁調整ルーチン
            while (b >= 64)
            {
                b /= 64;
                rt *= 4;
            }
            while (b < 1)
            {
                b *= 64;
                rt /= 4;
            }
            //ニュートン法による近似ルーチン
            for (int i = 0; i < 10; ++i)
            {
                rv = rv * 2 + b / (rv * rv);
                rv /= 3;
            }
            //出力ルーチン
            if (sgn)
                return rv * rt;
            else
                return -rv * rt;
        }
    }
}

Answer 1

値の1/3乗で計算してはダメ？

public static double Cbrt(double d)
{
  return System.Math.Pow(d , 1 / 3d);
}

Answer 2

外池と申します。

取り組んでおられのは、、立方根を求めるアルゴリズムの考察とお見受けします。お示し頂いているコードをどうチューンアップするかですよね？

Math.Sqrtの１５倍程度の時間で計算できるというのは、なかなか良い結果と思います。と言うのも、Math.Sqrtを使うと、結局のところ、FPUで平方根を求める１つの命令（IA-32命令セットなら「fsqrt」）に展開されるからで・・・。お示し頂いているメソッドをさらに速くするのは、かなり大変だと思います。たぶん、いじるとかえって遅くなるのではないかと。

ソースをざっと眺めて思うことは、桁調整のところが、改善の余地があるかも・・・、です。「rt」は浮動小数点ではなく、int型でも良いのではないかと思います。そうすると「rt*=4」や「rt/=4」はシフト命令にできて、浮動小数点命令で実行するより若干速くなるかなぁ。あと、「b」の指数部だけを取り出せれば、ループせずに桁調整できる可能性があります。大幅な桁調整が必要な数値の場合は速くなるでしょうが、桁調整不要な数値の場合はかえって遅くなりますね。

gekkaさんの仰るMath.Powは・・・、別途メソッドを呼び出すようで、さて、どれぐらい速いかは、試してみないとわからないですね。

---

（ホームページを再開しました）

Answer 3

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

計算速度は向上しました。Math.Sqrtの５倍程度になりました。

しかしながら、若干誤差が気になりました。※

と言いますと、元々立方根に特化したメソッドを作りたかったのが、動機でありましたので。

せっかく答えていただいたのに恐縮です。

※

引数の1/3がまず丸め誤差が入り込み、及びMath.Powのアルゴリズムはそれよりかなり複雑なため、ひょっとしたら自作したほうが速くなるかもと思いました。ちなみに、Powだと凡そ相対誤差は10^-15～10^-14程度あり、計算精度は、上で挙げたアルゴリズムのほうが良かったです。

Answer 4

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

>結局のところ、FPUで平方根を求める１つの命令（IA-32命令セットなら「fsqrt」）に展開されるからで・・・。

なるほど、そういう事情でしたか。

>「rt」は浮動小数点ではなく、int型でも良いのではないかと思います。そうすると「rt*=4」や「rt/=4」はシフト命令にできて、浮動小数点命令で実行するより若干速くなるかなぁ。

引数が10^100オーダーも考慮に入れると、結局のところ

>指数部だけを取り出せれば、ループせずに桁調整できる可能性があります。

が出来ればよいのですが、System.Doubleのメソッドでそれを実現できそうもないので、System.Int64へ強引に（Objectを経由して）キャストして、指数部のビットを取り出すような方法でしょうか。

何れにせよ、ハードの仕様により、期待している速度は得られそうにないですね。

Answer 5

外池です。

速さを、Double型で表現可能な数値の範囲全体にわたって、まんべんなく求めるのか、極端な値は多少遅くても良いのかによりますが・・・、

Double型のビット列をInt64型で取り出す一番速い方法は、私は、構造体でStructLayout属性やFieldOffset属性を使ってDouble型とInt64型の２つのメンバー変数を同じメモリアドレスに割り当てる方法ではないかと思います。C++で言うところのＵＮＩＯＮと同じようなことができますので、.Net Framework組み込みのメソッドを呼び出すことなく、かなり高速にビット列を取り出せます。

VBでやった例は拙ＨＰにありますので参考にしてください。
http://homepage3.nifty.com/numericworld/computer/vb/dotnetmemo.htm

FPUで処理できる平方根と、独自のアルゴリズムで計算する立法根と、まぁ、速さの違いは致し方ないかと。

Math.Powのメソッドは、おそらく底の変換をFPUの対数関数（FYL2XやFYL2XP1）でやって、その上で、FPUの指数関数（F2XM1）を呼び出していると想像しますが、ホント、まったく想像です。丸め誤差を気にされるようであれば、なおさら、ご自身でアルゴリズムを書き下されることが得策と思います。

---

（ホームページを再開しました）

Answer 6

外池さん

再度のご返事ありがとうございます。

>構造体でStructLayout属性やFieldOffset属性を使ってDouble型とInt64型の２つのメンバー変数を同じメモリアドレスに割り当てる方法ではないかと思います。

を参考にいたしました。

あと、ニュートン法を行う前に線形近似もいたしました。その結果、計算時間は最初に示したアルゴリズムに比べ８割程度に改善いたしました。

以上、ご回答ありがとうございました。

p.s.

もし、さらにご意見が在る場合の参考として、改良したコードを挙げます。

namespace User
{
    [StructLayout(LayoutKind.Explicit)]
    struct DoubleField
    {
        [FieldOffset(0)]
        public double r;
        [FieldOffset(0)]
        public uint i0;
        [FieldOffset(4)]
        public uint i1;
        public DoubleField(double d) { i0 = 0; i1 = 0; r = d; }
        public DoubleField(int i) { r = 0.0; i0 = (uint)i; i1 = 0; }
        public DoubleField(uint i) { r = 0.0; i0 = i; i1 = 0; }
        public DoubleField(long i) { r = 0.0; i0 = (uint)(i & 0xFFFFFFFF); i1 = (uint)((ulong)i >> 32); }
        public static implicit operator DoubleField(double d) { return new DoubleField(d); }
        public static implicit operator DoubleField(int i) { return new DoubleField(i); }
        public static implicit operator DoubleField(long i) { return new DoubleField(i); }
    }
    //--------------------------------------------------------------
    public static class Math
    {
        public static double Cbrt(double d)
        {
            if (d == 0 || d == -0) return 0;
            DoubleField rv = d;
            uint sgn  = rv.i1 & 0x80000000U, //符号フラグ
                 sexp = (rv.i1 & 0x7FF00000U) >> 20, //指数フラグ
                 exp  = (sexp / 0x3U + 0x2AAU) << 20; //演算後の指数フラグ
            rv.i1 = (rv.i1 & 0x000FFFFFU) | (((sexp % 0x3U) + 0x3FFU) << 20); //1 ≦ rv ＜ 8 の算出
            double b = rv.r;
            //線形化ルーチン
            rv.r = rv.r * 0.14285714285714286 + 0.9571428571428571;
            //ニュートン法による近似ルーチン
            for (int i = 0; i < 4; ++i)
            {
                rv.r = rv.r * 2.0 + b / (rv.r * rv.r);
                rv.r /= 3.0;
            }
            //出力ルーチン
            rv.i1 &= 0x000FFFFFU;
            rv.i1 = rv.i1 | sgn | exp;
            return rv.r;
        }
    }
}

Answer 7

このコードで実行時間が削減できたのは、ニュートン法の段数を減らしたのが一番効いてるんじゃないかなと感じました。

どうせならですが、Pow(1/3D)を初期値にしてニュートン法を１段だけにするのはいかがでしょうか。

1段だけにする根拠ですが、a^xのグラフの0.33...付近での傾きを考えてのことです。

まず、aは桁数調整の結果8未満になっています。a^xの傾きは ln a * a^x ですからｘ＝0.33...付近ではたかだか2.809...未満。

1/3Dの丸め誤差（1ビット）は2.8倍程度までしか拡大せず、それは仮数部での3ビット以下ということになります。

もともとのニュートン法によるコードでは仮数部52ビットの精度を10段で確保していました。ニュートン法の誤差桁数は段数に比例しますので、1段で最低5ビットは精度を上げられていたわけです。

ですから、3ビットまで拡大した可能性のある丸め誤差も1段で十分、と。

 

# ただ、桁数調整を自力でやる必要ってあるんでしょうか？　仮数部は変わらないので精度も変わらないはずかと・・・　共用体つかって、CPU可搬性をなくしているだけのように見えます。

Answer 8

外池です。

miuras_netさんのご提案、大変良いと思います。漸化式の初期値がなんとかならないか、とも思っていたのですが、Math.Powを使えば、非常にいいですよね。

これは、基本、Math.Powを使って立法根を求めて、解を改善する、ということですね。

「CPU可搬性をなくしているとのご指摘」　ごもっともです。ただ、IEEEに定める倍精度浮動小数点数のフォーマットであれば、使えるとは思います。

---

（ホームページを再開しました）

Answer 9

miuras_net さん

ご意見ご尤もです。

Math.Powを近似値として、ニュートン法を一段行えば、目的のものが得られ、Math.Sqrtの６～７倍の処理時間でした。恐らく現状ではベストと思います。さらに、この場合桁調整をする必要すらありません。ニュートン法は、凡そn段処理を行えば、2^n桁も止まります。ですので、８桁ほどの精度さえあれば、１段で求まります。欲を言えば、その程度の精度でかまわないので、MathPowより速いものが有ればいいとは思います。

それが故に、Powより速く計算できるのかな思ってしまったのが事の顛末です。しかし、今回の質疑応答で、ほとんどの場合FPUに任せるのが一番得策であるということでした。

外池 さん

ただ、前回のコードに比べ、桁が大きい数値でも、桁調整なしのときの処理時間と同等に出来たのは、良い成果と思います。さらに、共用体の宣言や、IEEEの仕様を実技によって勉強になったことは、大きな成果です。

ご意見ありがとうございます。