Перейти к содержанию
Поиск
    

Квадратный корень в целых числах

Ответить

gladov    0

gladov
Опубликовано (изменено) · Жалоба

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

Для тестов использовался LPC1758, разогнанный до 100МГц. Отвлекающих от вычислений факторов нет, кроме таймера, отсчитывающего миллисекунды.

 

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
[codebox]

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

klen    1

klen
Опубликовано · Жалоба
... В скобках приведено время вычисления ...

sqrt1() - 0x349d (0x296d)

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

 

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

AndyDev    0

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

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

 

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

Serj78    0

Serj78
Опубликовано · Жалоба

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

Потом он без изменений переехал на арм.

 

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

На память- к 4-му.

 

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

 

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

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

ChipKiller    0

ChipKiller
Опубликовано · Жалоба
Вся арифметика построена на целых числах, поэтому и корень тоже решено было считать без плавучки
... подойдет метод Герона (кажется). Гляньте здесь http://algolist.manual.ru/maths/count_fast/intsqrt.php

 

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

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

ReAl    0

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    28

blackfin
Опубликовано · Жалоба
Чего я в 1998 не прилепил (С) :-)

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

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

ReAl    0

ReAl
Опубликовано · Жалоба

О, за кубический спасибо :-)

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

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

alex_shevchenko    0

alex_shevchenko
Опубликовано (изменено) · Жалоба

Процедура целочисленного вычисления квадратного корня из 32/16-битного числа методом Ньютона с учетом остатка от деления на 2 и предварительным подбором делителя.

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

NewtonSQRT.zip

Изменено пользователем IgorKossak
удалил простынь

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

alex_shevchenko    0

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    5

Rst7
Опубликовано · Жалоба

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

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;
}

 

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

   \   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    6

Alex11
Опубликовано · Жалоба

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

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

Rst7    5

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

 

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

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

GenaSPB    11

GenaSPB
Опубликовано (изменено) · Жалоба 
/*
**      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
[codebox] для длинного кода!!!

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

zhz    0

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....

Поделиться сообщением

Ссылка на сообщение
Поделиться на другие сайты

Присоединяйтесь к обсуждению

Вы можете написать сейчас и зарегистрироваться позже. Если у вас есть аккаунт, авторизуйтесь, чтобы опубликовать от имени своего аккаунта.

Гость
К сожалению, ваш контент содержит запрещённые слова. Пожалуйста, отредактируйте контент, чтобы удалить выделенные ниже слова.
Ответить в этой теме...

×   Вставлено с форматированием.   Вставить как обычный текст

  Разрешено использовать не более 75 эмодзи.

×   Ваша ссылка была автоматически встроена.   Отображать как обычную ссылку

×   Ваш предыдущий контент был восстановлен.   Очистить редактор

×   Вы не можете вставлять изображения напрямую. Загружайте или вставляйте изображения по ссылке.

×
×
×
  • Создать...