Integer кубический корень

Ответ 1

В книге "Восторг Хакера" есть алгоритмы для этой и многих других задач. Код здесь онлайн. РЕДАКТИРОВАТЬ: этот код не работает должным образом с 64-разрядными целочисленными, и инструкции в книге о том, как исправить это для 64-разрядных, несколько сбивают с толку. Правильная 64-битная реализация (включая тестовый пример) доступна здесь.

Я сомневаюсь, что ваша функция squareroot работает "правильно" - для аргумента она должна быть ulong a не n :) (но тот же подход будет работать с использованием cbrt вместо sqrt, хотя не все математические библиотеки C имеют корневые функции куба).

Ответ 2

Вы можете попробовать выполнить шаг Ньютона, чтобы исправить ошибки округления:

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;

Один шаг Ньютона должен быть достаточным, но у вас могут быть ошибки "один за другим" (или, возможно, больше?). Вы можете проверить/исправить их с помощью этапа окончательной проверки и увеличения, как в вашем OQ:

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

или некоторые такие.

(Я признаю, что я ленив, правильный способ сделать это - тщательно проверить, чтобы определить, какие (если есть) элементы проверки и приращения действительно необходимы...)

Ответ 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). Используйте несколько шагов Ньютона (см. Подсказку ответа), чтобы закончить его.

Ответ 4

// 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;
} 

Ответ 5

Я адаптировал алгоритм, представленный в 1.5.2 (корень k-й) в современной компьютерной арифметике (Брент и Циммерман). Для случая (k == 3) и с учетом "относительно" точной завышенной оценки первоначального предположения - этот алгоритм, по-видимому, превосходит приведенный выше код "Хакерского восторга".

Не только это, но MCA как текст обеспечивает теоретическое обоснование, а также доказательство правильности и окончательных критериев.

При условии, что мы можем предоставить "относительно" хорошую начальную переоценку, я не смог найти случай, который превышает (7) итераций. (Это эффективно связано с 64-битными значениями, имеющими 2 ^ 6 бит?) В любом случае, это улучшение (21) итераций в коде HacDel!

Первоначальная оценка, которую я использовал, основана на "округлении" количества значащих битов в значении (x). Учитывая (б) значащих бит в (х), мы можем сказать: 2^(b - 1) <= x < 2^b. Я утверждаю без доказательств (хотя это должно быть относительно легко продемонстрировать), что: 2^ceil(b/3) > x^(1/3)

Вот мой код, как это в настоящее время...

static inline uint32_t u64_cbrt (uint64_t x)
{
#if (0) /* an exact IEEE-754 evaluation: */

    if (x <= (UINT64_C(1) << (53)))
        return (uint32_t) cbrt((double) x);
#endif

    int bits_x = (64) - __builtin_clzll(x);

    if (bits_x == 0)
        return (0); /* cbrt(0) */

    int exp_r = (bits_x + 2) / 3;

    /* initial estimate: 2 ^ ceil(b / 3) */
    uint64_t est_r = UINT64_C(1) << exp_r, r;

    do /* quadratic convergence (?) */
    {
        r = est_r;
        est_r = (2 * r + x / (r * r)) / 3;
    }
    while (est_r < r);

    return ((uint32_t) r); /* floor(cbrt(x)) */
}

crbt вероятно, не так уж и полезен - в отличие от вызова sqrt который может быть реализован с максимальной эффективностью на современном оборудовании. Тем не менее, я видел выигрыш для наборов значений до 2^53 (точно представленных в двойных единицах IEEE-754), что меня удивило.

Единственным недостатком является деление на: (r * r) - это может быть медленно, так как задержка целочисленного деления продолжает отставать от других достижений в ALU. Деление на константу: (3) обрабатывается взаимными методами на любом современном оптимизирующем компиляторе.

Ответ 6

Я бы исследовал, как это сделать вручную, а затем перевести это в компьютерный алгоритм, работающий в базе 2, а не в базе 10.

В итоге у нас есть алгоритм, похожий на (псевдокод):

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, наблюдая, что побитовое или эквивалентно добавлению, рефакторинг до result**3 + 3 * i * result ** 2 + 3 * i ** 2 * result + i ** 3, кэширование значений result**3 и result**2 между итерациями и использование сдвигов вместо умножения.

Ответ 7

Я бы настроил справочную таблицу этой формы:

int cbrt_int = {0, 1, 8, 27, 64, 125,...};

И тогда вам просто нужно выполнить бинарный поиск в массиве в поисках n, и индекс будет cbrt (n).