|
Квадратный корень в целых числах, Benchmark for ARM |
|
|
|
Jan 18 2012, 06:41
|
Частый гость
 
Группа: Свой
Сообщений: 169
Регистрация: 10-11-05
Из: Воронеж
Пользователь №: 10 687

|
Понадобилось мне в одной задаче приблизительно (но быстро) вычислять квадратный корень. Вся арифметика построена на целых числах, поэтому и корень тоже решено было считать без плавучки. Точность результата +- 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]
|
|
|
|
2 страниц
< 1 2
|
 |
Ответов
(15 - 24)
|
Sep 18 2012, 12:33
|
Группа: Новичок
Сообщений: 8
Регистрация: 21-01-09
Пользователь №: 43 758

|
Процедура целочисленного извлечения квадратного корня из 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
|
|
|
|
|
Sep 21 2012, 12:30
|
Группа: Новичок
Сообщений: 8
Регистрация: 21-01-09
Пользователь №: 43 758

|
Вариант адаптивной процедуры целочисленного извлечения квадратного корня из 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] для длинного кода!!!!
|
|
|
|
|
Sep 21 2012, 12:58
|
Группа: Новичок
Сообщений: 8
Регистрация: 21-01-09
Пользователь №: 43 758

|
Все это легко перепихивается и на другие платформа. В данный момент пишу под 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] для длинного кода!!!!!!!
|
|
|
|
|
Sep 21 2012, 14:39
|
Группа: Новичок
Сообщений: 8
Регистрация: 21-01-09
Пользователь №: 43 758

|
Маленькая оговорка: В адаптивной процедуре _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
|
|
|
|
|
Sep 22 2012, 12:49
|

Нечётный пользователь.
     
Группа: Свой
Сообщений: 2 033
Регистрация: 26-05-05
Из: Бровари, Україна
Пользователь №: 5 417

|
Цитата(alex_shevchenko @ Sep 21 2012, 15:58)  Еще пару десятков тактов в среднем можно сэкономить, если развернуть цикл в ущерб объему (для тех, кому важна скорость): А чтобы меньше людям мешать -- в тег codebox завернуть этот оффтопик нельзя было? Раз не выходит не лезть в ARM-ветку со своим пиком? (кстати, я уж если пишу на ассемблере, то свой инструмент знаю и в таких местах для разворачивания циклов .rept использую или подряд вызовы по месту склёпанного макроса одной итерации — и читается легче, и на экране короче… а так тупо вручную писать… тьху…) Ладно бы ещё было на С, но с какой-то пик-фичой типа атрибутов, без которой код работоспособен, но просто с ними он на пике эффективнее. Если считаете, что это легко спортить на арм, то спортьте и покажите. Нет -- зачем это тут?
--------------------
Ну, я пошёл… Если что – звоните…
|
|
|
|
|
Sep 22 2012, 19:32
|
Группа: Новичок
Сообщений: 8
Регистрация: 21-01-09
Пользователь №: 43 758

|
To IgorKossak, Банами меня не напугаешь. Зависание в форумах для меня не есть жизненно необходимым, а так - от нехрен делать в перерывах между делом. Обойдусь легко. И если мало интересует содержание и наполненность ресурса - можешь банить и удалять. Не имею возражений. Но все-таки, наверное, будет правильно, если перенести это с одноименным названием в ветку к пикам. А то некоторые армисты при слове PIC аж из себя выходят и бьются в конвульсиях. Сердечко может не выдержать. Please.  P.S. Крайне неудачно настроен интерфейс движка. Одна минута на редактирование - это издевательство над пользователями. Редактор - совсем не wysiwyg. И ту ерунду, которую он высыпает в и тоге на экран, я успеваю только просмотреть, но изменить что-то уже нет возможности. Поэтому - или миритесь с конечным результатом, или СДЕЛАЙТЕ УДОБНЫЙ ИНТЕРФЕЙС ПОЛЬЗОВАТЕЛЯ.
|
|
|
|
|
Sep 22 2012, 19:44
|

Шаман
     
Группа: Модераторы
Сообщений: 3 064
Регистрация: 30-06-04
Из: Киев, Украина
Пользователь №: 221

|
Цитата(alex_shevchenko @ Sep 22 2012, 22:32)  наверное, будет правильно, если перенести это с одноименным названием в ветку к пикам. Какой же Вы умный задним числом. Отчего же сразу не разместить тему где положено? Цитата(alex_shevchenko @ Sep 22 2012, 22:32)  или миритесь с конечным результатом, или СДЕЛАЙТЕ УДОБНЫЙ ИНТЕРФЕЙС ПОЛЬЗОВАТЕЛЯ. Я поставлю вопрос иначе - или соблюдайте правила форума на том, что есть или на появляйтесь здесь вообще. Ка Вы без форума, так и форум без Вас вполне способен обойтись. И последнее - публичное обсуждение действий модератора запрещено правилами, так что получИте обещанные 15 суток. Модератор.
|
|
|
|
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0
|
|
|