Jump to content

    

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

Понадобилось мне в одной задаче приблизительно (но быстро) вычислять квадратный корень. Вся арифметика построена на целых числах, поэтому и корень тоже решено было считать без плавучки. Точность результата +- 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 микросекунды. Алгоритм взят отсюда

 

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

Edited by IgorKossak
[codebox]

Share this post


Link to post
Share on other sites
... В скобках приведено время вычисления ...

sqrt1() - 0x349d (0x296d)

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

 

Share this post


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

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

 

Share this post


Link to post
Share on other sites

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

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

 

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

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

 

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

 

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

Share this post


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

 

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

Share this post


Link to post
Share on other sites
... подойдет метод Герона (кажется). Гляньте здесь 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?

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

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

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

NewtonSQRT.zip

Edited by IgorKossak
удалил простынь

Share this post


Link to post
Share on other sites

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нс).

Share this post


Link to post
Share on other sites

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

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. Если известна верхняя граница аргумента, то количество итераций можно уменьшить.

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites
Получилось 82 такта на корень из 32-битного беззнакового числа.

 

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

Share this post


Link to post
Share on other sites
/*
**      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);
       }
}

Edited by IgorKossak
[codebox] для длинного кода!!!

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now
Sign in to follow this