Квадратный корень в целых числах, 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]

[klen]

... В скобках приведено время вычисления ...
sqrt1() - 0x349d (0x296d)

в время в хексах это тру!

[AndyDev]

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

Если написать на чистом асме, то результат явно можно еще улучшить.

[Serj78]

Как-то задавался этим вопросом, скорость библиотечного корня на avr не устраивала, написал свой.
Потом он без изменений переехал на арм.

Позже посмотрю, к какому из ваших алгоритмов он ближе.
На память- к 4-му.

А вы пока проверьте вот что- зависимость длительности от входного числа. Иногда длительность плавает

К слову сказать, целочисленный корень из библиотеки Кейла работает медленнее моего всего на 3-5% примерно, но жрет кода на 500 байт больше.

[ChipKiller]

Вся арифметика построена на целых числах, поэтому и корень тоже решено было считать без плавучки

... подойдет метод Герона (кажется). Гляньте здесь http://algolist.manual.ru/maths/count_fast/intsqrt.php

PS. Результаты замеров обычно пишут в циклах.

[ReAl]

... подойдет метод Герона (кажется). Гляньте здесь http://algolist.manual.ru/maths/count_fast/intsqrt.php

Так «рекордсмен» оттуда же и взят, если я правильно понял.

p.s. offtop:
Блииинн...
http://home.utah.edu/~nahaj/factoring/isqrt.c.html (С) 2003
Чего я в 1998 не прилепил (С) :-) http://groups.google.com/group/fido7.ru.al...c6a27847dfc4a30?
А алгоритм 3 из сравниваемых — это я уже годом позже в embedded
http://groups.google.com/group/fido7.ru.em...d4744f6ff04ac53?

[blackfin]

Чего я в 1998 не прилепил (С) :-)

И тут тоже забыли (С) прилепить: Алгоритм извлечения кубического корня.

[ReAl]

О, за кубический спасибо :-)
А квадратный в столбик у меня в школе в 8 классе по программе был, перед логарифмической линейкой :-) (линейка, кажись, в факультативом курсе)

[alex_shevchenko]

Процедура целочисленного вычисления квадратного корня из 32/16-битного числа методом Ньютона с учетом остатка от деления на 2 и предварительным подбором делителя.
Написана на MPLAB® ASM30 Assembler для семейства dsPIC30 (+PIC24, +dsPIC33) по мотивам статьи Николая Гарбуз "Вычисление квадратного корня из целого числа".


Сообщение отредактировал IgorKossak - Sep 7 2012, 10:15

Причина редактирования: удалил простынь

Прикрепленные файлы

 NewtonSQRT.zip ( 4.36 килобайт ) Кол-во скачиваний: 28

 

[alex_shevchenko]

NewtonSQRT16 считает корень в серднем за 150 машинных тактов (3-6мкс при Tcy=33.9нс). NewtonSQRT32 - в среднем за 3000 тактов при 7-30 итерациях. Основное время сжирает библиотечная функция деления ___udivsi3 - поэтому 88мкс.
Но рекордсменом по скорострельности вычисления квадраных корней является библиотечная функция Q15sqrt из libq.h (см. 16-Bit_Language_Tools_Libraries_51456.pdf), вычисляющая корень из числа в формате с фиксированной точкой Q15 в диапзоне аргумента -2^15...2^15 -1 строго за 79 машинных тактов (2.7мкс! при Tcy=33.9нс).

[Rst7]

Тут есть тонкость - деление не очень производительная операция, так что часто лучше использовать более тупой подход

CODE

unsigned int sqrt(unsigned int v)
{
#define SQRT_ITER(MASK) if (v<r*r) r&=~MASK; r|=(MASK>>1);
unsigned int r=0xC000;
if (v<0x40000000) r=0x4000;
SQRT_ITER(0x4000);
SQRT_ITER(0x2000);
SQRT_ITER(0x1000);
SQRT_ITER(0x0800);
SQRT_ITER(0x0400);
SQRT_ITER(0x0200);
SQRT_ITER(0x0100);
SQRT_ITER(0x0080);
SQRT_ITER(0x0040);
SQRT_ITER(0x0020);
SQRT_ITER(0x0010);
SQRT_ITER(0x0008);
SQRT_ITER(0x0004);
SQRT_ITER(0x0002);
SQRT_ITER(0x0001);
return r;
}


int main()
{
sqrt(1234567890);
return 0;
}

Вид одной итерации:

CODE

   \   0000000E   01FB01F2           MUL      R2,R1,R1
   \   00000012   9042               CMP      R0,R2
   \   00000014   38BF               IT       CC
   \   00000016   21F48041           BICCC    R1,R1,#0x4000
   \   0000001A   41F40051           ORR      R1,R1,#0x2000

Итого 5 тактов на итерацию, а итераций нужно разрядность_аргумента/2. Если известна верхняя граница аргумента, то количество итераций можно уменьшить.

[Alex11]

Я делал методом последовательных приближений, но под DSP от TI. Получилось 82 такта на корень из 32-битного беззнакового числа. Для меньшей разрядности почти пропорционально меньше. Но он на ассемблере, напрямую перенести на АРМ не получится.

[Rst7]

Получилось 82 такта на корень из 32-битного беззнакового числа.

Да вот мой код не более 80 (5*16) тактов занимает для 32х бит. Там еще IT не очень понятно когда 1 такт, а когда 0, так что бывает и меньше.

[Genadi Zawidowsk...]

CODE

/*
** ISQRT.C
**
** Calculate integer sqare root.
**
** Copyright © MocroGenSf 1992
**
** Created 23-Sep-1992.
**
*/

#include <stdio.h>
#include <stdlib.h>

typedef unsigned int uint;
typedef signed int sint;

typedef unsigned long ulong;
typedef signed long slong;

/* Calculate integer value of sqare root */
uint
isqrt(ulong a)
{ auto ulong x0, x1;
auto slong delta0, delta1;

if (a < 2)
return (a);
delta0 = 0;
x1 = a / 2; /* Initial approximation. */
for (;;)
{ x0 = x1;
x1 = (a / x0 + x0) >> 1;
if ((delta1 = x1 - x0) == 0) return ((uint) x1);
if ((delta0 + delta1) == 0) return ((uint) x0);
delta0 = delta1;
}
}

uint
ihypot(int dx, int dy)
{ return isqrt((slong)dx * (slong)dx + (slong)dy * (slong)dy);
}

void main(void)
{
sint dx, dy;
uint hyp;
char buff[128];
for (;;)
{ printf("dx = ");
//gets(buff, 128);
dx = atoi(buff);
printf("dy = ");
//gets(buff, 128);
dy = atoi(buff);
hyp = ihypot(dx, dy);
printf("ihypot(%d,%d) = %u\n", dx, dy, hyp);
}
}

Сообщение отредактировал IgorKossak - Sep 11 2012, 07:19

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

[zhz]

Andrew N. Sloss, Dominic Symes, Chris Wright
ARM System Developer’s Guide. Designing and Optimizing System Software.

7.4 Square Roots
7.4.1 Square Root by Trial Subtraction
... following optimized assembly to implement the preceding algorithm in
only 50 cycles including the return. ...

7.4.2 Square Root by Newton-Raphson Iteration
... It uses a table lookup followed by two Newton-Raphson iterations and is accurate
to a maximum error of 2**-29. On an ARM9E the code takes 34 cycles including the return....

[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

[alex_shevchenko]

Вариант адаптивной процедуры целочисленного извлечения квадратного корня из 16-битного числа методом вычисления "в столбик" и округлением до ближайшего целого. Время вычисления сокращается в зависимости от разрядности числа и составляет 24-90 машинных циклов.

CODE

;-----------------------------------------------------------------------------------------------------------------------;
; An integral calculation of the SQRT16 ;
; ;
; Execution time: max 90Tcy ;
; 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;
; Description: Целочисленное извлечение квадратного корня из 16-битного числа методом вычисления "в столбик";
; с округлением до ближайшего целого. ;
;-----------------------------------------------------------------------------------------------------------------------;
_Sqrt16:
CP W0, #0x0001; Сравнение аргумента с '1'.
BRA GTU, .+4; Если аргумент больше - переход к алгоритму извлечения корня, иначе:
RETURN ; выход с аргументом в качестве результата.
MUL.UU W2, #0, W2; Предварительная очистка регистра накопления результата/остатка.
FF1L W0, W1; Определение номера старшего разряда.
DEC W1, W1; Технический декремент номера разряда.
BCLR W1, #0; Проверка номера разряда на четность +
SUBR W1, #14, W1; разбиение разрядов на пары.
BSW.C W3, W1; Инициализация начального приближения корня.
IOR W2, W3, W1; tmp = W2 = k | res
LSR W2, W2; res >>= 1;
CP W1, W0;
BRA GTU, .+6; if( x >= tmp ) {
SUB W0, W1, W0; x -= tmp;
IOR W2, W3, W2; res |= k; }
LSR W3, #2, W3; k >>= 2;
BRA NZ, .-14
CP W0, W2; Срвнение результата и остатка.
BRA LEU, .+4; Если остаток меньше либо равен - пропуск инструкции, иначе
INC W2, W2; инкремент результата.
EXCH W0, W2; Технический своп регистров результата и остатка.
RETURN ; Выход из процедуры.

Сообщение отредактировал IgorKossak - Sep 22 2012, 16:44

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

[Rst7]

Вариант адаптивной процедуры ...

Это все хорошо, а название топика читали?

[alex_shevchenko]

Все это легко перепихивается и на другие платформа. В данный момент пишу под PICи, а в АРМ влазить нет времени.

Еще пару десятков тактов в среднем можно сэкономить, если развернуть цикл в ущерб объему (для тех, кому важна скорость):

CODE

_Sqrt16:
SUB W0, #0x0000, W2; Сравнение аргумента с '0' + техническое перемещение.
BRA GTU, .+4 ; Если аргумент больше - переход к алгоритму извлечения корня, иначе:
RETURN ; выход с аргументом в качестве результата.
FF1L W0, W0 ; Определение номера старшего разряда.
DEC W0, W0 ; Технический декремент номера разряда.
LSR W0, W0 ; Разбиение разрядов числа на пары.
MUL.UU W0, #7, W0; Вычисление длины прыжка с обнулением регистра результата W1.
BRA W0 ; Переход к итерации с необходимым номером.
; Итерация 1
MOV #0x4000, W0 ; k = 0x4000.
BRA .+4 ; Сокращение одного машинного такта.
NOP ; Техническое заполнение тела итерации.
CP W0, W2 ;
BRA GTU, .+6 ; if(x >= tmp) {
SUB W2, W0, W2; x -= tmp;
IOR W1, W0, W1; res |= k; }
; Итерация 2
MOV #0x1000, W0 ;
IOR W0, W1, W3; tmp = k | res
LSR W1, W1 ; res >>= 1;
CP W3, W2 ;
BRA GTU, .+6 ; if(x >= tmp) {
SUB W2, W3, W2; x -= tmp;
IOR W1, W0, W1; res |= k; }
; Итерация 3
MOV #0x0400, W0 ;
IOR W0, W1, W3; tmp = k | res
LSR W1, W1 ; res >>= 1;
CP W3, W2 ;
BRA GTU, .+6 ; if(x >= tmp) {
SUB W2, W3, W2; x -= tmp;
IOR W1, W0, W1; res |= k; }
; Итерация 4
MOV #0x0100, W0 ;
IOR W0, W1, W3; tmp = k | res
LSR W1, W1 ; res >>= 1;
CP W3, W2 ;
BRA GTU, .+6 ; if(x >= tmp) {
SUB W2, W3, W2; x -= tmp;
IOR W1, W0, W1; res |= k; }
; Итерация 5
MOV #0x0040, W0 ;
IOR W0, W1, W3; tmp = k | res
LSR W1, W1 ; res >>= 1;
CP W3, W2 ;
BRA GTU, .+6 ; if(x >= tmp) {
SUB W2, W3, W2; x -= tmp;
IOR W1, W0, W1; res |= k; }
; Итерация 6
MOV #0x0010, W0 ;
IOR W0, W1, W3; tmp = k | res
LSR W1, W1 ; res >>= 1;
CP W3, W2 ;
BRA GTU, .+6 ; if(x >= tmp) {
SUB W2, W3, W2; x -= tmp;
IOR W1, W0, W1; res |= k; }
; Итерация 7
MOV #0x0004, W0 ;
IOR W0, W1, W3; tmp = k | res
LSR W1, W1 ; res >>= 1;
CP W3, W2 ;
BRA GTU, .+6 ; if(x >= tmp) {
SUB W2, W3, W2; x -= tmp;
IOR W1, W0, W1; res |= k; }
; Итерация 8
IOR W1, #0x0001, W3; tmp = k | res
LSR W1, W0 ; res >>= 1;
CP W3, W2 ;
BRA GTU, .+6 ; if(x >= tmp) {
SUB W2, W3, W2; x -= tmp;
IOR W0, #0x0001, W0; res |= k; }
; Округление до ближайшего целого
CP W2, W0 ; Сравнение результата с остатком.
BRA LEU, .+4 ; Если остаток меньше, либо равен - пропуск инструкции, иначе:
INC W0, W0 ; инкремент результата.
RETURN ; Выход из процедуры.

Сообщение отредактировал IgorKossak - Sep 22 2012, 16:42

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

[alex_shevchenko]

Маленькая оговорка:
В адаптивной процедуре _Sqrt16 для инициализации начального приближения корня применена пара инструкций, использующих особенности работы процессора.
Если где-то между ними:

SUBR W1, #14, W1 ; разбиение разрядов на пары.
BSW.C W3, W1 ; Инициализация начального приближения корня

произойдет прерывание без сохранения контекста (LB статусного регистра SR), то результат вычисления может быть потерян.
Если нет желания возиться с контекстом, то лучше использовать что-то вроде:

SUBR W1, #14, W1 ; разбиение разрядов на пары.
MOV #0x0001, W3 ;
SL W3, W1, W3 ; Инициализация начального приближения корня.

или

DISI #2 ; выкл. прерываний на две инструкции.
SUBR W1, #14, W1 ; разбиение разрядов на пары.
BSW.C W3, W1 ; Инициализация начального приближения корня

Т.е. на один такт больше.

Сообщение отредактировал alex_shevchenko - Sep 21 2012, 14:42

[Genadi Zawidowsk...]

прерывание без сохранения контекста

Ой... поясните, где в жизни на такое можно напороться?

[ReAl]

Еще пару десятков тактов в среднем можно сэкономить, если развернуть цикл в ущерб объему (для тех, кому важна скорость):

А чтобы меньше людям мешать -- в тег codebox завернуть этот оффтопик нельзя было? Раз не выходит не лезть в ARM-ветку со своим пиком?

(кстати, я уж если пишу на ассемблере, то свой инструмент знаю и в таких местах для разворачивания циклов .rept использую или подряд вызовы по месту склёпанного макроса одной итерации — и читается легче, и на экране короче… а так тупо вручную писать… тьху…)

Ладно бы ещё было на С, но с какой-то пик-фичой типа атрибутов, без которой код работоспособен, но просто с ними он на пике эффективнее.
Если считаете, что это легко спортить на арм, то спортьте и покажите. Нет -- зачем это тут?

[IgorKossak]

alex_shevchenko, длинные простыни Вашего кода позаворачивал в тэги codebox там где это было возможно.
Вклинивание в ветку ARM с пиками пока терплю.
Но на этом мои предупреждения заканчиваются, в следующий раз будут удаления сообщений с баном на неделю и больше.
Модератор.

[alex_shevchenko]

To IgorKossak,

Банами меня не напугаешь. Зависание в форумах для меня не есть жизненно необходимым, а так - от нехрен делать в перерывах между делом. Обойдусь легко. И если мало интересует содержание и наполненность ресурса - можешь банить и удалять. Не имею возражений.
Но все-таки, наверное, будет правильно, если перенести это с одноименным названием в ветку к пикам. А то некоторые армисты при слове PIC аж из себя выходят и бьются в конвульсиях. Сердечко может не выдержать.
Please.


P.S. Крайне неудачно настроен интерфейс движка. Одна минута на редактирование - это издевательство над пользователями. Редактор - совсем не wysiwyg. И ту ерунду, которую он высыпает в и тоге на экран, я успеваю только просмотреть, но изменить что-то уже нет возможности. Поэтому - или миритесь с конечным результатом, или СДЕЛАЙТЕ УДОБНЫЙ ИНТЕРФЕЙС ПОЛЬЗОВАТЕЛЯ.

[IgorKossak]

наверное, будет правильно, если перенести это с одноименным названием в ветку к пикам.

Какой же Вы умный задним числом. Отчего же сразу не разместить тему где положено?

или миритесь с конечным результатом, или СДЕЛАЙТЕ УДОБНЫЙ ИНТЕРФЕЙС ПОЛЬЗОВАТЕЛЯ.

Я поставлю вопрос иначе - или соблюдайте правила форума на том, что есть или на появляйтесь здесь вообще. Ка Вы без форума, так и форум без Вас вполне способен обойтись.

И последнее - публичное обсуждение действий модератора запрещено правилами, так что получИте обещанные 15 суток.
Модератор.