Квадратный корень в целых числах, Benchmark for ARM Опции

[gladov]

Понадобилось мне в одной задаче приблизительно (но быстро) вычислять квадратный корень. Вся арифметика построена на целых числах, поэтому и корень тоже решено было считать без плавучки. Точность результата +- 1 меня вполне устраивает. На просторах интернета, да и на этом форуме, удалось найти много готовых решений, но их основная масса ИМХО немного устарела, т.к. заточена под АЛУ без умножителя. Меня же интересовал алгоритм для применения на АРМах. Решено было попробовать некоторые из найденных и оценить их скорость работы на АРМах.
Для тестов использовался LPC1758, разогнанный до 100МГц. Отвлекающих от вычислений факторов нет, кроме таймера, отсчитывающего миллисекунды.

CODE

uint32_t sqrt1 (uint32_t L)
{
int32_t temp, div;
uint32_t rslt = L;
if (L <= 0) return 0;
else if (L & 0xFFFF0000L)
if (L & 0xFF000000L)
div = 0x3FFF;
else
div = 0x3FF;
else
if (L & 0x0FF00L) div = 0x3F;
else div = (L > 4) ? 0x7 : L;

while (1)
{
temp = L/div + div;
div = temp >> 1;
div += temp & 1;
if (rslt > div) rslt = (uint32_t)div;
else
{
if (L/rslt == rslt-1 && L%rslt==0) rslt--;
return rslt;
}
}
}

uint32_t sqrt2 (uint32_t src)
{
uint32_t wrk;
uint32_t dst;
int i;
dst = 0x8000;
wrk = 0x8000;
for(i=0; i<16; i++)
{
if(dst*dst>src) dst &= ~wrk;
wrk >>= 1;
dst |= wrk;
}
return dst;
}

uint32_t sqrt3 (uint32_t src)
{
uint32_t mask, sqr = 0, temp;
int j=16;

temp = 0xC0000000;
do {
if( src & temp ) break;
temp>>=2;
} while( --j);
if( j==0 ) return 0;
mask = temp & (temp>>1);
do {
temp = sqr | mask;
sqr >>= 1;
if( temp <= src ) {
sqr |= mask;
src -= temp;
}
mask >>= 2;
} while( --j );

return sqr;
}

uint32_t sqrt4 (uint32_t Val)
{
unsigned int bitSqr = 0x40000000;
unsigned int root = 0;

while (bitSqr != 0)
{
if (Val >= (bitSqr + root))
{
Val = Val - (bitSqr + root);
root = (root >> 1) | bitSqr;
}
else
root = (root >> 1);
bitSqr = (bitSqr >> 2);
}
return(root);
}

int TestSqrt(uint32_t(*func)(uint32_t))
{
tick_count_t starttime = get_tick_count();
for (uint32_t i = 0; i < 10000000; i++) {
uint32_t s = func(i);
//Check the value
uint32_t sq = s * s;
if (!((sq == i) ||
(sq > i) && (s-1)*(s-1) < i ||
(sq < i) && (s+1)*(s+1) > i))
{
while (1);
}
}
return get_tick_count() - starttime;
}

Проверка с бесконечным циклом ни на одном алгоритме не сработала - все вычислялось четко.

Результаты замеров в миллисекундах приведены ниже. В скобках приведено время вычисления без проверки правильности, только вычисление корня в цикле.
sqrt1() - 0x349d (0x296d)
sqrt2() - 0x4286 (0x3986)
sqrt3() - 0x4807 (0x3cca)
sqrt4() - 0x4933 (0x44df)

Явный лидер - sqrt1(). Получается, что среднее время вычисления кв. корня около 1 микросекунды. Алгоритм взят отсюда

Меня результат вполне устраивает, но если у кого-то есть другие интересные алгоритмы, давайте и их проверим. Интересно же найти лучший.

Сообщение отредактировал IgorKossak - Jan 18 2012, 07:55

Причина редактирования: [codebox]

[alex_shevchenko]

Процедура целочисленного извлечения квадратного корня из 16-битного числа методом вычисления "в столбик" и округлением до ближайшего целого. (MPLAB® ASM30 Assembler для семейства PIC24, dsPIC30, dsPIC33.) Время вычисления - max 85 тактов, за исключением значений 0 и 1 - 10 тактов.

Код

;-----------------------------------------------------------------------------------------------------------------------;
;                    An integral calculation of the SQRT16                                    ;
;                                                                                             ;
;    Execution time:    85Tcy (2.88uS @ Tcy=33.9nS)                                                                        ;
;    SFR used:        W0,W1,W2,W3                                                                ;
;    Dependencies:    No dependencies                                                            ;
;    Inputs:            W0 = input 16-bit unsigned int value                                    ;
;    Outputs:        W0 = output SQRT value of input unsigned int                            ;
;                    W1 = remainder                                                            ;
;    Description:    Целочисленное извлечение квадратного корня из 16-битного числа методом вычисления "в столбик";
;                    с округлением до ближайшего целого.                                        ;
;-----------------------------------------------------------------------------------------------------------------------;
_Sqrt16:
        CP            W0,                #0x0001; Сравнение аргумента с '1'.
        BRA            GTU,            .+6; Если аргумент больше - переход к алгоритму извлечения корня, иначе:
        CLR            W1            ; Техническое обнуление регистра остатка.
        RETURN                    ; выход с аргументом в качестве результата.
        MOV            #0x0000,        W1; Предварительная очистка регистра результата.
        MOV            #0x4000,        W3; k = 0x4000
        IOR            W1,                W3,                W2; tmp = k | res
        LSR            W1,                W1; res >>= 1;
        CP             W2,                W0;
        BRA            GTU,            . +6; if( x >= tmp ) {
        SUB            W0,                W2,                W0; x -= tmp;
        IOR            W1,                W3,                W1; res |= k; }
        LSR            W3,                #2,                W3; k >>= 2;
        BRA            NZ,                . -14; пока (k != 0) - переход на очередную итерацию.
        CP            W0,                W1; Сравнение результата и остатка.
        BRA            LEU,            . +4; Если остаток меньше либо равен - пропуск инструкции, иначе:                        
        INC            W1,                W1; инкремент результата.
        EXCH        W0,                W1; Технический своп регистров результата и остатка.
        RETURN                    ; Выход из процедуры.

P.S. После дизассемблирования стало понятно, что Q15sqrt(x) - использует разложение степенной функции в ряд Тейлора, с хорошей оптимизацией алгоритма.

Сообщение отредактировал alex_shevchenko - Sep 18 2012, 12:36