Статья будет очень короткой, но надеюсь весьма полезной.

В одном из проектов необходимо было достаточно быстро вычислять некоторую формулу которая в себя включает вычисление целочисленного квадратного корня 32-х разрядного целого числа. Как показал опыт функция из стандартной библиотеки выполняет вычисления достаточно долго. В итоге пришли к выводу, что необходимо сделать всё самим.

Так и родился этот достаточно короткий и быстрый код:

/**
@function MATH_sqrt(x) − Вычисление квадратного корня.
@param uint32_t value − число квадратный корень которого нужно вычислять
@retval uint32_t значение квадратного корня
 */
uint32_t MATH_sqrt(register uint32_t value)
{
    //lint -save -e40 -e10 -e845 -e522
    register uint32_t d = 0, ax = 0;

    // Находим количество битов дополняющих до 32, по сути a + b = 32, где и номер старшего бита
    // Затем из 32 вычитаем значение, получаем номер старшего бита b = 32 - a и кладем в 'а'

    #if (__ARMCC_VERSION >= 6010050)
        __asm volatile (
        "CLZ    %[A], %[VALUE]"           "\n\t"
        "RSB    %[A], %[A], 0x1F"
                : [A]"=&r" (ax)
                : [VALUE]"r" (value), [A]"r" (ax)
        );
    #else
        __asm {
            CLZ    ax, value
            RSB    ax, ax, 0x1F
        };
    #endif

    /* Нахождение квадратного корня методом Ньютона. Основная формула для вычисления:
    *     sqrt = (sqrt' + value / sqrt') / 2
    *  где:
    *     value - значение корень которого находим
    *     sqrt - текущее значение приближения
    *     sqrt' - значение приближения из предыдущего шага
    *  Для ускорения первым приближением выбирается значение равное value сдвинутому на половину
    */
    d = value >> (ax >> 1);    // первое приближение
    ax = (d + value / d) >> 1;
    d = (ax + value / ax) >> 1;
    ax = (d + value / d) >> 1;
    value = (ax + value / ax) >> 1;
    //lint -restore
    return value;
}

Код достаточно хорошо прокомментирован и лишний раз думаю объяснять его не стоит. Но если кратко, то на ассемблере находится номер старшего единичного бита числа, это означает, что:

// 41646465 = 00000010 01111011 01111001 10000001
// номер старшего единичного бита будет 25, необходимо учесть, счёт начинается с 0

Далее в качестве первого приближения используется исходное значение, но сдвинутое на половину номера максимального бита и как в классике (кто не особо знаком с данным алгоритмом прошу обратиться за помощью в Википедию).

Здесь стоит остановиться, изначально алгоритм предполагал использование цикла с необходимым для нахождения количеством повторений, но от этого пришлось отказаться ввиду больших накладных расходов. Поэтому было решено выстроить линейную цепочку из звеньев вычисления, число которых гарантированно даст адекватный результат.

Описанная функция входит в мой пакет для keil, который можно скачать по ссылке - devprodest.Lib.pack

Код будет работать только на ядре Cortex M3 и старше, на M0 лучше использовать другие методы.

Буду рад любым комментариям по поводу ускорения данного кода.

• Немного истории. K-line адаптер.