Как извлечь корень в c

Как извлечь корень в c

Разделы

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

18/01/2013 — 01:11 Pavel Bobkov

При программировании микроконтроллеров разработчики иногда сталкиваются с проблемой вычисления квадратного корня. Например, данная операция требуется при выполнении быстрого преобразования Фурье или вычислении среднеквадратического значения сигнала.
В стандартной библиотеке Си – 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);
>

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

где А – фиксированное положительное число, а 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;
>

Related items

• Планировщик для микроконтроллера
• Медианный фильтр
• AVR4027: Трюки и советы по оптимизации Си кода для 8-и разрядных AVR микроконтроллеров. Ч.2
• AVR4027: Трюки и советы по оптимизации Си кода для 8-и разрядных AVR микроконтроллеров. Ч.1
• Что размещать в заголовочном файле .h?

Comments

# Nobody 2013-01-18 02:58
К сожалению, не могу раскрыть подробности его работы, потому что они мне неизвестны.
Это итерационная формула Герона, если добавить ещё одну итерацию, то точность увеличится.
# Nobody 2013-01-18 05:30

Хотя, я был не прав. По итерационной формуле Герона нужно 8 делений для получения приближенного ответа. В описанном вами методе можно увеличить точность добавив ещё одну итерацию:
x = b/x;
a = x = (x+a)>>1;
Но для AVR, данный алгоритм не оптимальный, т.к. деление выполняется долго. Посмотрите на метод бинарного поиска целочисленного квадратного корня (только умножение и сложение), описан в книге «Алгебраические трюки для программиста». Там ещё описан алгоритм без умножения, только сдвиги, сложения и битовые операции. Можно попробовать адаптировать для AVR его, тогда выигрыш во времени будет значительный.

# Pashgan 2013-01-18 06:40

Да нет, все правильно. В википедии приведена формула Герона.
Xn+1 = (Xn + A/Xn)*1/2
A — число из которого требуется вычислить корень, X1 — любое положительное число.

# Pashgan 2013-01-18 07:01

Если написать код так
Code:
unsigned int root(unsigned int a)
unsigned int x;
x = (a/0x3f + 0x3f)>>1;
x = (a/x + x)>>1;
x = (a/x + x)>>1;
return(x);
>

получаются точно такие же результаты. Как по ответам, так и по скорости исполнения кода. А по объему даже небольшой выигрыш. Зачем надо было так мудрить?

# spkik 2016-04-07 13:25
Quoting Nobody:

Хотя, я был не прав. По итерационной формуле Герона нужно 8 делений для получения приближенного ответа. В описанном вами методе можно увеличить точность добавив ещё одну итерацию:
x = b/x;
a = x = (x+a)>>1;
Но для AVR, данный алгоритм не оптимальный, т.к. деление выполняется долго. Посмотрите на метод бинарного поиска целочисленного квадратного корня (только умножение и сложение), описан в книге «Алгебраические трюки для программиста». Там ещё описан алгоритм без умножения, только сдвиги, сложения и битовые операции. Можно попробовать адаптировать для AVR его, тогда выигрыш во времени будет значительный.

подскажите пожалуйста как проделать это с переменной типа long?

# Nikolay 2013-01-18 07:43
Добрый день. Уважаемый, Pashgan, вот ви пишете Quote:
— компактность (~80 байтов), — скорость выполнения (~1000 циклов для AVR).
Можете рассказать как вы определяете сколько занимает и сколько циклов выполняется код?
# Neptun 2013-01-18 09:37

Напишу как делаеться в AVR Studio. Пишеться код — компилим,смотри м сколько занял,добавляем функцию — и смотрим новий размер кода. Разница между новым и старым значение есть размер функции.

Для скорости выполнения. ставим брейкпойнт перед вызовом функции и после,запускаем симуляцию — обнуляем cycle counter,запуска ем симуляцию — и новое значение будеть скоростью выполнения (также можно увидеть сколько время исполняеться функция в мкс или мс).

# Pashgan 2013-01-18 11:19

Для IAR`a . Нужно включить опцию создания листинга программы. Project > Options > C/C++ Compiler > List галочка Output list file. Если включить еще и Assembler mnemonics в lst файле будет ассемблерный код, сгенерированный компилятором из твоей программы. Эта информация полезна для оптимизации сишного кода, конечно, если ты знаешь ассемблер.
Затем запускаешь компиляцию проекта и с левой стороны (в окне отображения структуры проекта) ищешь файлы с расширением *.lst Они будут созданы для каждого программного модуля. В конце этого файла есть табличка со списком функций и значениями занимаемой памяти.

Чтобы прикинуть скорость выполнения какого-нибудь куска кода (обычно функции), я прогоняю этот код в программном симуляторе IAR`a. Включаю опцию Project > Options > Linker > Debug Information . Запускаю компиляцию и отладку с помощью кнопки Debug (Ctrl+D). Устанавливаю брейкпоинты, открываю окно с регистрами микроконтроллер а (меню View > Register) и запускаю код на выполнение по шагам (F11) или непрерывно (f5). В окне регистров в разделе CPU Register есть строка CYCLES. Она отображает число прошедших тактов. По показаниям этого числа можно прикинуть сколько тактов занимает выполнение функции.

То же самое можно делать и в AVR Studio. Там это даже лучше получается, потому что студия моделирует прерывания, а IAR нет.

# Nobody 2013-01-18 09:44

У меня нет компилятора для AVR под рукой.
Можете проверить функцию:
Code:
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;
>

Сколько занимает и как долго выполняется?

# Neptun 2013-01-18 11:04
Quoting Nobody:

У меня нет компилятора для AVR под рукой.
Можете проверить функцию:
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;
>
Сколько занимает и как долго выполняется?

F = 8 MHz, ATmega8, optimization O0 (none):

размер 14 байт.
скорость 540 циклов — 67.5uS

F = 8 MHz, ATmega8, optimization Os (none):

размер 14 байт.
скорость 2 циклf — 0.25uS

# Neptun 2013-01-18 11:07

Код которий тестировал:

unsigned int value = 0;

unsigned int isqrt(unsigned int x)
unsigned int m, y, b;
m = 0x4000;
y = 0;
while (m != 0)
b = y | m;
y = y >> 1;

int main(void)
asm(«nop»);
value = isqrt(4096);
asm(«nop»);

# Pashgan 2013-01-18 11:33
У меня в IAR` получилось 52 байта, 180 тактов при LOW оптимизации по размеру кода
# spkik 2016-04-07 13:27
Quoting Nobody:

У меня нет компилятора для AVR под рукой.
Можете проверить функцию:
Code:
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;
>

Сколько занимает и как долго выполняется?

подскажите пожалуйста как проделать это с переменной типа long?
# Neptun 2016-04-07 15:49

тоже самое но с типом long
unsigned long isqrt(unsigned long x)
unsigned long 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;
>

# Васьок 2014-10-09 12:34

Пользовался детским алгоритмом, на мой взгляд достаточно быстр и достаточно компактный. Идея в том что от числа последовательно отнимаются все нечётные числа, и сколько вычитаний удалось сделать, таков и корень числа. Пример, число 49;
1) 49 — 1 = 48
2) 48 — 3 = 45
3) 45 — 5 = 40
4) 40 — 7 = 33
5) 33 — 9 = 24
6) 24 — 11 = 13
7) 13 — 13 = 0

7 циклов, корень числа 49 — 7.

И кстати при работе с МК типа AVR-ки лучше избегать делений, т.к. у AVR ядра нет аппаратного деления, а программное занимает дофига тактов. Другое дело ARM Cortex-M3 и выше, у которых деление выполняется за 2. 12 тактов.

# Петгосян 2016-11-20 13:37

У функции корня есть некоторые свойства симметрии, которые позволяют вычислять ее только на некотором отрезке, а потом решение распространить на всю ось. Например,
sqrt(a*2^16)=2^ 8*sqrt(a).

Удобно в качестве такого отрезка взять значения [2^30-2^31), потому что остальные значения можно свести к нему побитовым сдвигом и при этом не будет происходить потеря точности. Сначала вычисляем первый значащий бит (программно половинным делением или процессорной инструкцией, например на ARM это __clz). Потом сдвигаем входное число на это кличество бит и вычисляем корень, полученное значение сдвигаем обратно на в два раза меньшее количество).
Для вычисления корня на отрезке интерполируем его многочленом Лагранжа (параболой). Например, возьмем в качестве точек многочлена 2^30, 1,5 * 2^30, 2^31. Можно воспользоваться сторонним сервисом, и не возиться с вычислением коэффициентов. У меня получилась такая формула:
-x^2/499100218444523 + x/52370 + 14575
Очевидно, напрямую её использовать нельзя, потому что значения не влазят даже в диапазон целых. Но надо учесть, что нам важны только 16 бит результата, поэтому можно немного схитрить и вынести что-то за скобки.
(-x/9530269590 + 1) * x/52370 + 14575
(-x/145420 + 65536) * (x/65536) / 52370 + 14575
Ну и последнее — заменить деление на умножение. Допустим, у нас в резерве 30 бит числа. Мы хотим поделить некое число x, например, на 543. Вычисляем, в числе 543 есть 10 бит, в х 16 бит.
x / 543 * 2^26 / 2^26
x * (2^26 / 543) / 2^26
x * 123589 / 2^26
Теперь эти знания применяем к своему многочлену.
(-x/2^14 * 7384 / 2^16 + 2^16) * (x/2^16) / 2^16 * 20503 / 2^14 + 14575
Не ручаюсь за правильность коэффициентов, надо внимательно проверить.
Когда писал, не учел одну штуку, число бит может быть нечетным, отрезок надо брать больше.

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

# Петгосян 2016-11-20 13:38

Если умножение делается за один такт, можно сделать вычисление корня побитовым подбором. На первой итерации выставляем 16 бит в 1, возводим в квадрат, сравниваем с входным числом. Если больше, сбрасываем бит. Потом с 15 битом повторяем то же самое и т.д. Как в АЦП.

# nebelwerfer 2017-01-22 14:19

А что за магическое число:
Code: m = 0x4000; ?
Это половина от максимума int ?
А вот 2й вариант работает отлично, спасибо!
Мне нужно считать корни из больших чисел до 250 000 000, поэтому увеличил количество:
Code: x = (a/x + x)>>1;
до 7 и точность приемлима.

У вас недостаточно прав для комментирования.

Как извлечь корень в c

Сообщения: 467
Благодарности: 25

Конфигурация компьютера
Процессор: Intel(R) Core(TM)2 Quad CPU Q9550 @ 2.83GHz
Материнская плата: Gigabyte GA-EP45-UD3L
Память: Corsair XMS2 CM2X1024-6400 X 2
HDD: WDC WD1001FALS-00J7B0
Видеокарта: Asus EN9800GTX
Звук: Creative SB X-Fi XtremeGamer
Монитор: Samsung SyncMaster 960B(F)/MagicSyncMaster CX915T (Digital)

Кстати, может быть ктонить кинет ссыль на описание всех функций cmath, люблю математику я (кстати пой плагин для проверки орфографии в браузере предложил вариант этого слова как маразматику )

——-
Мне надо не так много «Полезных сообщений», чтоб сровнять их количество с моими постами :).

Сообщения: 3806
Благодарности: 824

crashtuak, на кой вам плагин проверки орфографии, если всё равно игнорируете его подсказки?

Сообщения: 467
Благодарности: 25

Конфигурация компьютера
Процессор: Intel(R) Core(TM)2 Quad CPU Q9550 @ 2.83GHz
Материнская плата: Gigabyte GA-EP45-UD3L
Память: Corsair XMS2 CM2X1024-6400 X 2
HDD: WDC WD1001FALS-00J7B0
Видеокарта: Asus EN9800GTX
Звук: Creative SB X-Fi XtremeGamer
Монитор: Samsung SyncMaster 960B(F)/MagicSyncMaster CX915T (Digital)

Цитата Busla:

crashtuak, на кой вам плагин проверки орфографии, если всё равно игнорируете его подсказки? »

У моего плагина словарь не безграничный и он не такой умный как россиянин, тоесть мой плагин немного улучшает правописание украинца (тоесть меня), но полюбому, даже с плагином я не смогу написать так, как это сделает носитель языка!

——-
Мне надо не так много «Полезных сообщений», чтоб сровнять их количество с моими постами :).

Сообщения: 7318
Благодарности: 1204

Конфигурация компьютера
Процессор: Intel Core i5-2500K (-> 4.2 ГГц)
Материнская плата: AsRock P67 Pro3 (Socket 1155)
Память: 2×4096 МБ DDR3-1333 Kingston (-> 1600 МГц)
HDD: Hitachi 2 ТБ + 0,75 ТБ (SATA2), SSD: OCZ Vertex3 (120 ГБ, SATA3)
Видеокарта: VTX3D Radeon 7950
Звук: Asus Xonar DS (PCI)
Блок питания: Chieftec CTG-550-80P (550 Вт, октябрь 2011)
CD/DVD: DVD-RW LiteOn iHAS124 (SATA2)
Монитор: Dell U2412M (24″, 1920×1200, e-IPS)
ОС: Windows 7 Ultimate x64
Индекс производительности Windows: 7.6
Прочее: Кулер Zalman 10X Performa

Цитата AlexTNT:

А для возведения в квадрат можно использовать функцию pow(x,y) [из той же cmath] , где x- нужное число, а у — степень (в нашем случае 2). »

sqrt , sqrtf , sqrtl

Поскольку C++ допускает перегрузку, можно вызывать перегрузки sqrt , которые принимают типы float или long double . В программе C, если вы не используете макрос для вызова этой функции, sqrt всегда принимает и возвращает. double

При использовании sqrt() макроса тип аргумента определяет, какая версия функции выбрана. Дополнительные сведения см . в разделе «Математика с универсальным типом».

По умолчанию глобальное состояние этой функции ограничивается приложением. Чтобы изменить это поведение, см . статью «Глобальное состояние» в CRT.

Возвращаемое значение

Функции sqrt возвращают квадратный корень x . По умолчанию, если x отрицательный, sqrt возвращает неопределенное время NaN .

Входные данные	Исключение SEH	Исключение _matherr
± QNaN, IND	нет	_DOMAIN
-INF	нет	_DOMAIN
x < 0	нет	_DOMAIN

Требования

Функция	Заголовок C	Заголовок C++
sqrt , sqrtf , sqrtl
sqrt Макрос

Дополнительные сведения о совместимости см. в разделе Совместимость.

Пример

// crt_sqrt.c // This program calculates a square root. #include #include #include int main( void )

The square root of 45.35 is 6.73 

Math. Sqrt(Double) Метод

Некоторые сведения относятся к предварительной версии продукта, в которую до выпуска могут быть внесены существенные изменения. Майкрософт не предоставляет никаких гарантий, явных или подразумеваемых, относительно приведенных здесь сведений.

Возвращает квадратный корень из указанного числа.

public: static double Sqrt(double d);

public static double Sqrt (double d);

static member Sqrt : double -> double

Public Shared Function Sqrt (d As Double) As Double

Параметры

Число, квадратный корень которого требуется найти.

Возвращаемое значение

Одно из значений, перечисленных в следующей таблице.

Параметр d	Возвращаемое значение
Нуль или положительное число	Положительный квадратный корень из d .
Отрицательное число	NaN
Равняется NaN	NaN
Равняется PositiveInfinity	PositiveInfinity

Примеры

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

// Create an array containing the area of some squares. Tuple[] areas = < Tuple.Create("Sitka, Alaska", 2870.3), Tuple.Create("New York City", 302.6), Tuple.Create("Los Angeles", 468.7), Tuple.Create("Detroit", 138.8), Tuple.Create("Chicago", 227.1), Tuple.Create("San Diego", 325.2) >; Console.WriteLine("  \n", "City", "Area (mi.)", "Equivalent to a square with:"); foreach (var area in areas) Console.WriteLine("   miles per side", area.Item1, area.Item2, Math.Round(Math.Sqrt(area.Item2), 2)); // The example displays the following output: // City Area (mi.) Equivalent to a square with: // // Sitka, Alaska 2,870.3 53.58 miles per side // New York City 302.6 17.40 miles per side // Los Angeles 468.7 21.65 miles per side // Detroit 138.8 11.78 miles per side // Chicago 227.1 15.07 miles per side // San Diego 325.2 18.03 miles per side 

open System // Create a list containing the area of some cities. let areas = [ "Sitka, Alaska", 2870.3 "New York City", 302.6 "Los Angeles", 468.7 "Detroit", 138.8 "Chicago", 227.1 "San Diego", 325.2 ] printfn "%-18s %14s> %2s\n" "City" "Area (mi.)" "Equivalent to a square with:" for city, area in areas do printfn $"   miles per side" // The example displays the following output: // City Area (mi.) Equivalent to a square with: // // Sitka, Alaska 2,870.3 53.58 miles per side // New York City 302.6 17.40 miles per side // Los Angeles 468.7 21.65 miles per side // Detroit 138.8 11.78 miles per side // Chicago 227.1 15.07 miles per side // San Diego 325.2 18.03 miles per side 

Module Example Public Sub Main() ' Create an array containing the area of some squares. Dim areas() As Tuple(Of String, Double) = < Tuple.Create("Sitka, Alaska", 2870.3), Tuple.Create("New York City", 302.6), Tuple.Create("Los Angeles", 468.7), Tuple.Create("Detroit", 138.8), Tuple.Create("Chicago", 227.1), Tuple.Create("San Diego", 325.2) >Console.WriteLine("  ", "City", "Area (mi.)", "Equivalent to a square with:") Console.WriteLine() For Each area In areas Console.WriteLine("   miles per side", area.Item1, area.Item2, Math.Round(Math.Sqrt(area.Item2), 2)) Next End Sub End Module ' The example displays the following output: ' City Area (mi.) Equivalent to a square with: ' ' Sitka, Alaska 2,870.3 53.58 miles per side ' New York City 302.6 17.40 miles per side ' Los Angeles 468.7 21.65 miles per side ' Detroit 138.8 11.78 miles per side ' Chicago 227.1 15.07 miles per side ' San Diego 325.2 18.03 miles per side 

Применяется к

См. также раздел