Побитовый целочисленный алгоритм корня куба

Вот простой способ вычисления целочисленного квадратного корня:

int isqrt(int num)
{
    int root=0;
    int b = 0x8000;
    int a=0, c=0;

    while (b) {
        c = a|b;

        if (c*c <= num)
            a |= b;   

        b >>= 1;
    }       
}

Ingeniously (спасибо Wikipedia), это можно оптимизировать следующим образом:

int sqrt(short num)
{
    int op = num;
    int res = 0;
    int one = 1 << 30;

    while (one > op)
        one >>= 2;

    while (one != 0) {
        if (op >= res + one) {
            op -= res + one;
            res = (res >> 1) + one;
        }
        else
          res >>= 1;
        one >>= 2;
    }
    return res;
}

Мой вопрос: Можно ли написать аналогично оптимизированный алгоритм для целочисленного кубического корня? (Это должно выполняться на небольшом микроконтроллере, который предпочитает не делать умножения)

Ответы

Ответ 1

В соответствии с этим вопросом SO и ответом, помеченным, из книги Hacker Delight вы можете найти эту реализацию:

int icbrt2(unsigned x) {
   int s;
   unsigned y, b, y2;

   y2 = 0;
   y = 0;
   for (s = 30; s >= 0; s = s - 3) {
      y2 = 4*y2;
      y = 2*y;
      b = (3*(y2 + y) + 1) << s;
      if (x >= b) {
         x = x - b;
         y2 = y2 + 2*y + 1;
         y = y + 1;
      }
   }
   return y;
}

Ответ 2

Да, алгоритм может быть расширен до корня куба, даже без умножений.

Смотрите этот код: http://www.hackersdelight.org/HDcode/icbrt.c.txt

И подумайте о покупке книги Hackers Delight, откуда он приходит. Если вам приходится решать такие проблемы чаще, чем раз в год, вы обязательно должны это прочитать!

Ответ 3

Это (экстремальная) оптимизированная версия кода Hacker Delight, оптимизированная на С#, как упоминалось другими.

Для справки (на моем компьютере): Math.Sqrt занимает около 35 нс, cbrt < 15 нс.

Используются умножения по малым числам, но их легко заменить на сдвиги и добавляет. Например, наибольшее значение multipication (последняя строка): "12 * z" == > "(z < 3) + (z < 2)"

Трудно судить, допустим ли размер кода для небольшого микроконтроллера.

Первый шаг: бинарный поиск, операторы if, большие значения ( >= 1u < 24) найдены относительно быстрее, при выполнении поиска обрабатываются небольшие значения (< 64).

Второй шаг: прыжок в развернутый цикл, "метки".

private static uint cbrt(uint x)
{
    uint y = 2, z = 4, b = 0;
    if (x < 1u << 24)
        if (x < 1u << 12)
            if (x < 1u << 06)
                if (x < 1u << 03)
                    return x == 0u ? 0u : 1u;
                else
                    return x < 27u ? 2u : 3u;
            else
                if (x < 1u << 09) goto L8; else goto L7;
        else
            if (x < 1u << 18)
                if (x < 1u << 15) goto L6; else goto L5;
            else
                if (x < 1u << 21) goto L4; else goto L3;
    else
        if (x >= 1u << 30) goto L0;
        else
            if (x < 1u << 27) goto L2; else goto L1;

L0: x -= 1u << 30; if (x >= 19u << 27)
    { x -= 19u << 27; z = 9; y = 3; } goto M0;
L1: x -= 1u << 27; if (x >= 19u << 24)
    { x -= 19u << 24; z = 9; y = 3; } goto M1;
L2: x -= 1u << 24; if (x >= 19u << 21)
    { x -= 19u << 21; z = 9; y = 3; } goto M2;
L3: x -= 1u << 21; if (x >= 19u << 18)
    { x -= 19u << 18; z = 9; y = 3; } goto M3;
L4: x -= 1u << 18; if (x >= 19u << 15)
    { x -= 19u << 15; z = 9; y = 3; } goto M4;
L5: x -= 1u << 15; if (x >= 19u << 12)
    { x -= 19u << 12; z = 9; y = 3; } goto M5;
L6: x -= 1u << 12; if (x >= 19u << 09)
    { x -= 19u << 09; z = 9; y = 3; } goto M6;
L7: x -= 1u << 09; if (x >= 19u << 06)
    { x -= 19u << 06; z = 9; y = 3; } goto M7;
L8: x -= 1u << 06; if (x >= 19u << 03)
    { x -= 19u << 03; z = 9; y = 3; } goto M8;

M0: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 24; 
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M1: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 21; 
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M2: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 18;
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M3: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 15;
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M4: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 12; 
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M5: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 09; 
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M6: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 06; 
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M7: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 03; 
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M8: y *= 2; return x <= 3 * y + 12 * z ? y : y + 1;
}