Print this page

Быстрое вычисление квадратного корня на Си

18/01/2013 - 01:11

  При программировании микроконтроллеров разработчики иногда сталкиваются с проблемой вычисления квадратного корня. Например, данная операция требуется при выполнении быстрого преобразования Фурье или вычислении среднеквадратического значения сигнала.
   В стандартной библиотеке Си – math.h, есть функция для вычисления квадратного корня sqrt(), которой при желании можно воспользоваться. Она работает с числами типа float, обеспечивает высокую точность результата, но требует для своей работы длительного времени. Для микроконтроллера AVR это порядка 3000 циклов тактовой частоты (проверено в компиляторе IAR на разных уровнях оптимизации).
  Если к точности вычисления корня не предъявляются высокие требования, можно воспользоваться упрощенным алгоритмом, занимающим меньше места в памяти и выполняющим вычисления в несколько раз быстрее.

Алгоритм выглядит так. 


unsigned int root(unsigned int x)
{

   unsigned int a,b;
   b = x;
   a = x = 0x3f;
   x = b/x;
   a = x = (x+a)>>1;
   x = b/x;
   a = x = (x+a)>>1;
   x = b/x;
   x = (x+a)>>1;
   return(x); 
}

 Как мне подсказали умные люди, алгоритм основан на итерационной формуле Герона.

Xn+1 = (A/Xn + Xn)*1/2

где А – фиксированное положительное число, а X1 – любое положительное число.
Итерационная формула задаёт убывающую (начиная со 2-го элемента) последовательность, которая при любом выборе X1 быстро сходится к квадратному корню из числа А.

   Ради интереса я переписал алгоритм в явном виде. Скомпилированный, он ничуть не потерял ни в быстродействии, ни в объеме. Объем даже на пару байтов уменьшился.


unsigned int root1(unsigned int a)
{
   unsigned int x;
   x = (a/0x3f + 0x3f)>>1;
   x = (a/x + x)>>1;
   x = (a/x + x)>>1;
   return(x);
}


   Недостатки приведенного кода в том, что он работает только с целыми 16-ти разрядными числами и при больших значениях аргумента вычисления становятся не точными. Правда, точность вычислений можно повысить, добавив еще несколько итераций, но за это естественно придется платить быстродействием.

   Код занимает прядка 70 байт и выполняется ~ за 700 циклов. Данные получены в компиляторе IAR AVR при medium оптимизация по скорости.

   Точность вычисления данного алгоритма можно оценить по приведенному ниже графику. Синий график построен по значениям, полученным c помощью библиотечной функции sqrt(), красный график по значениям функции root().


Сравнение двух функций извлечения квадратного корня

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


unsigned int isqrt(unsigned int x)
{
   unsigned int m, y, b;
   m = 0x4000;
   y = 0;
   while (m != 0){
      b = y | m;
      y = y >> 1;
      if (x >= b) {
         x = x - b;
         y = y | m;
      }
      m = m >> 2;
   }
   return y;
}


~50 байт, <200 циклов для IAR AVR с medium оптимизацией по скорости. Точность не оценивал.

Related items