20

64 ビット (符号なし) 立方根の高速なコードを探しています。(私は C を使用しており、gcc でコンパイルしていますが、必要な作業のほとんどは言語やコンパイラに依存しないものになると思います。) 64 ビットの符号なし整数を ulong で示します。

入力 n が与えられると、(整数) 戻り値 r が次のようになる必要があります。

r * r * r <= n && n < (r + 1) * (r + 1) * (r + 1)

つまり、切り捨てられた n の立方根が必要です。のような基本的なコード

return (ulong)pow(n, 1.0/3);

範囲の終わりに向かって丸めているため、正しくありません。のような洗練されていないコード

ulong
cuberoot(ulong n)
{
    ulong ret = pow(n + 0.5, 1.0/3);
    if (n < 100000000000001ULL)
        return ret;
    if (n >= 18446724184312856125ULL)
        return 2642245ULL;
    if (ret * ret * ret > n) {
        ret--;
        while (ret * ret * ret > n)
            ret--;
        return ret;
    }
    while ((ret + 1) * (ret + 1) * (ret + 1) <= n)
        ret++;
    return ret;
}

正しい結果が得られますが、必要以上に遅くなります。

このコードは数学ライブラリ用であり、さまざまな関数から何度も呼び出されます。速度は重要ですが、ウォーム キャッシュを期待することはできません (したがって、2,642,245 エントリのバイナリ検索などの提案は適切です)。

比較のために、整数の平方根を正しく計算するコードを次に示します。

ulong squareroot(ulong a) {
    ulong x = (ulong)sqrt((double)a);
    if (x > 0xFFFFFFFF || x*x > a)
        x--;
    return x;
}
4

6 に答える 6

12

「Hacker's Delight」という本には、この問題や他の多くの問題に対するアルゴリズムが掲載されています。コードはオンラインでここにあります。EDIT : そのコードは 64 ビットの int では正しく動作しません。64 ビット用に修正する方法に関する本の指示はやや混乱しています。適切な 64 ビット実装 (テスト ケースを含む) は、こちらでオンラインになっています。

あなたの関数が「正しく」機能するかどうかは疑問です-それは:)ではなく、引数のためsquarerootであるべきです(しかし、すべてのC数学ライブラリに立方根関数があるわけではありませんが、代わりに を使用しても同じアプローチが機能します)。ulong ancbrtsqrt

于 2010-12-02T05:09:24.590 に答える
3

ニュートンの手順を試して、丸め誤差を修正できます。

ulong r = (ulong)pow(n, 1.0/3);
if(r==0) return r; /* avoid divide by 0 later on */
ulong r3 = r*r*r;
ulong slope = 3*r*r;

ulong r1 = r+1;
ulong r13 = r1*r1*r1;

/* making sure to handle unsigned arithmetic correctly */
if(n >= r13) r+= (n - r3)/slope;
if(n < r3)   r-= (r3 - n)/slope;

ニュートンステップは1つで十分ですが、1つずつ（またはそれ以上？）エラーが発生する可能性があります。OQのように、最後のチェックとインクリメントのステップを使用して、それらをチェック/修正できます。

while(r*r*r > n) --r;
while((r+1)*(r+1)*(r+1) <= n) ++r;

またはそのようなもの。

（私は怠惰であることを認めます。それを行う正しい方法は、チェックとインクリメントのどれが実際に必要かを慎重にチェックして決定することです...）

于 2010-12-02T18:31:47.863 に答える
3

powコストが高すぎる場合は、count-leading-zeros 命令を使用して結果の近似値を取得し、次にルックアップ テーブルを使用してから、いくつかのニュートン ステップを使用してそれを終了することができます。

int k = __builtin_clz(n); // counts # of leading zeros (often a single assembly insn)
int b = 64 - k;           // # of bits in n
int top8 = n >> (b - 8);  // top 8 bits of n (top bit is always 1)
int approx = table[b][top8 & 0x7f];

と が与えられるbとtop8、ルックアップ テーブル (私のコードでは 8K エントリ) を使用して の適切な近似値を見つけることができますcuberoot(n)。いくつかのニュートン ステップ (comingstorm の回答を参照) を使用して終了します。

于 2010-12-02T21:38:42.210 に答える
1 
// On my pc: Math.Sqrt 35 ns, cbrt64 <70ns, cbrt32 <25 ns, (cbrt12 < 10ns)

// cbrt64(ulong x) is a C# version of:
// http://www.hackersdelight.org/hdcodetxt/acbrt.c.txt     (acbrt1)

// cbrt32(uint x) is a C# version of:
// http://www.hackersdelight.org/hdcodetxt/icbrt.c.txt     (icbrt1)

// Union in C#:
// http://www.hanselman.com/blog/UnionsOrAnEquivalentInCSairamasTipOfTheDay.aspx

using System.Runtime.InteropServices;  
[StructLayout(LayoutKind.Explicit)]  
public struct fu_32   // float <==> uint
{
[FieldOffset(0)]
public float f;
[FieldOffset(0)]
public uint u;
}

private static uint cbrt64(ulong x)
{
    if (x >= 18446724184312856125) return 2642245;
    float fx = (float)x;
    fu_32 fu32 = new fu_32();
    fu32.f = fx;
    uint uy = fu32.u / 4;
    uy += uy / 4;
    uy += uy / 16;
    uy += uy / 256;
    uy += 0x2a5137a0;
    fu32.u = uy;
    float fy = fu32.f;
    fy = 0.33333333f * (fx / (fy * fy) + 2.0f * fy);
    int y0 = (int)                                      
        (0.33333333f * (fx / (fy * fy) + 2.0f * fy));    
    uint y1 = (uint)y0;                                 

    ulong y2, y3;
    if (y1 >= 2642245)
    {
        y1 = 2642245;
        y2 = 6981458640025;
        y3 = 18446724184312856125;
    }
    else
    {
        y2 = (ulong)y1 * y1;
        y3 = y2 * y1;
    }
    if (y3 > x)
    {
        y1 -= 1;
        y2 -= 2 * y1 + 1;
        y3 -= 3 * y2 + 3 * y1 + 1;
        while (y3 > x)
        {
            y1 -= 1;
            y2 -= 2 * y1 + 1;
            y3 -= 3 * y2 + 3 * y1 + 1;
        }
        return y1;
    }
    do
    {
        y3 += 3 * y2 + 3 * y1 + 1;
        y2 += 2 * y1 + 1;
        y1 += 1;
    }
    while (y3 <= x);
    return y1 - 1;
}

private static uint cbrt32(uint x)
{
    uint y = 0, z = 0, b = 0;
    int s = x < 1u << 24 ? x < 1u << 12 ? x < 1u << 06 ? x < 1u << 03 ? 00 : 03 :
                                                         x < 1u << 09 ? 06 : 09 :
                                          x < 1u << 18 ? x < 1u << 15 ? 12 : 15 :
                                                         x < 1u << 21 ? 18 : 21 :
                           x >= 1u << 30 ? 30 : x < 1u << 27 ? 24 : 27;
    do
    {
        y *= 2;
        z *= 4;
        b = 3 * y + 3 * z + 1 << s;
        if (x >= b)
        {
            x -= b;
            z += 2 * y + 1;
            y += 1;
        }
        s -= 3;
    }
    while (s >= 0);
    return y;
}

private static uint cbrt12(uint x) // x < ~255
{
    uint y = 0, a = 0, b = 1, c = 0;
    while (a < x)
    {
        y++;
        b += c;
        a += b;
        c += 6;
    }
    if (a != x) y--;
    return y;
}
于 2013-09-24T15:11:21.817 に答える
0

手作業でそれを行う方法を研究し、それを基数 10 ではなく基数 2 で動作するコンピューター アルゴリズムに変換します。

(疑似コード) のようなアルゴリズムになります。

Find the largest n such that (1 << 3n) < input.
result = 1 << n.
For i in (n-1)..0:
    if ((result | 1 << i)**3) < input:
        result |= 1 << i.

(result | 1 << i)**3ビットごとの OR が加算と同等であることを観察し、 にリファクタリングし、 と の値を反復間result**3 + 3 * i * result ** 2 + 3 * i ** 2 * result + i ** 3でキャッシュし、乗算の代わりにシフトを使用することで、 の計算を最適化できます。result**3result**2

于 2010-12-02T06:25:59.560 に答える