Извлечение квадратного корня. Скорость вычисления на STM8S

@Sovo · Регистрация: 05.10.2007

Студворк — интернет-сервис помощи студентам

В одной программе мне понадобилось вычисление квадратного корня из целого 32-х разрядного числа.
Процессор STM8S103 не имеет такой ассемблерной программы. Использовать библиотечную функцию с плавающей точкой не захотел.
В сети нашёл два алгоритма. Первый с подбором результата и делениями (по-моему, алгоритм Ньютона):

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
/******************************************************************************
* Извлечение квадратного корня из целого числа uint32_t
******************************************************************************/
uint16_t sqrt_newton (uint32_t x)
{
uint32_t temp, div;
uint16_t rslt = (uint16_t)x;
if (x <= 0)
{
return 0;
}
else if (x & 0xFFFF0000L)
{
if (x & 0xFF000000L)
{
div = 0x3FFF;
}
else
{
div = 0x3FF;
}
}
else
{
if (x & 0x0FF00L)
{
div = 0x3F;
}
else
{
div = (x > 4) ? 0x7 : x;
}
}
while (1)
{
temp = x / div + div;
div = temp >> 1;
div += temp & 1;
if (rslt > div)
{
rslt = (uint16_t)div;
}
else
{
if ((x / rslt == rslt - x) && (x % rslt == 0))
{
rslt--;
}
return rslt;
}
}
}

Второй без операций деления, только сдвиги (GLS):

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
uint16_t sqrt_GLS (uint32_t x)  // Hordware algorithm [GLS]
{
uint32_t m, y, b;
m = 0x40000000;
y = 0;
while (m != 0)
{              // Do 16 times.
b = y | m;
y = y >> 1;
if (x >= b)
{
x = x - b;
y = y | m;
}
m = m >> 2;
}
return y;
}

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

Итак, четыре деления плюс другие мат. операции при использовании алгоритма GLS и тактовой частоте 16 МГц заняли 2,5мс, а при использовании алгоритма Ньютона 0,98мс. Разница в два с половиной раза. А на первый взгляд и не скажешь!

@_pv · 22.04.2015, 16:57

что-то там сильно не так.
2.5мс это 40000 тактов, /16 = 2500 тактов на цикл, чтобы пару 32х разрядных циферок сложить/сдвинуть. так не быть не должно

@Sovo · 22.04.2015, 17:05

Поправляюсь - четыре деления, четыре умножения, четыре извлечения корня и четыре сдвига на 10. Плюс прерывание отвлекает раз в мс на 0,2мс.
Если вычесть прерывания, то останется - 0,8 мс против 2,1 мс

@dork_usir · 30.04.2015, 17:51

что-то там сильно не так.

Там с Ньютоном не так... глючная реализация...

@Sovo · 01.05.2015, 11:44

Что, по вашему, не так?
Есть реализация алгоритма ньютона без подбора результата, в два раза дольше. Результаты алгоритма ньютона и GLS совпадают. Что не так?

@wirty · 01.05.2015, 12:30

Результаты алгоритма ньютона и GLS совпадают.

А вы уверены? Возьмите число 45678901 и посмотрите результат.

@Sovo · 01.05.2015, 21:25

Я проверил. На выходе 6770 вместо 6758.
Я думаю, это из-за трёх моментов.
1. Неявное нахождение первого результата rslt.
2. Неправильное условие выхода из программы.
3. Непонятная для меня конструкция коррекции результата. По моему, она никогда не сработает.
Я поправил программу, гляньте:

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
uint16_t sqrt_newton (uint32_t x)
{
int32_t temp, div;
uint16_t rslt;
if (x > 0x0000FFFF)
{
rslt = 0xFFFF;
}
else
{
rslt = (uint16_t)x;
}
//   if (x <= 0)
//   {
//      return 0;
//   }
if (x & 0xFFFF0000L)
{
if (x & 0xFF000000L)
{
div = 0x3FFF;
}
else
{
div = 0x3FF;
}
}
else
{
if (x & 0x0FF00L)
{
div = 0x3F;
}
else
{
div = (x > 4) ? 0x7 : x;
}
}
while (1)
{
temp = x / div + div;
div = temp >> 1;
div += temp & 1;
int16_t error = rslt - div;
if ((error > 1) || (error < -1))
{
rslt = (uint16_t)div;
}
else
{
if ((x / rslt == rslt - 1) && (x % rslt == 0))
{
rslt--;
}
return rslt;
}
}
}

Вот что значит копипастить, даже из сорсфоржа :(

P.S. Я нашёл описание этого алгоритма:
http://www.codenet.ru/progr/alg/sqrtint.php
и понял, что корень там приблизительный, а условие коррекции должно быть другое. Немного поправил код под спойлером.

@wirty · 02.05.2015, 01:52

При 4026318845 опять лажа.

@_pv · 02.05.2015, 02:27

а тупо бинарный перебор с умножением не быстрее будет чем через деление на stm8?
по ссылке выше алгоритмы проверяются на х86 он делить числа умеет куда лучше чем stm8

Code
1
2
3
4
5
6
7
8
9
10
u16 sqrt(u32 x){
u16 result = 0xFFFF;
u16 d = 0x8000;
u8 n = 16;
while (n--){
if (result * result < x) result += d; else result -= d;
d >>= 1;
}
return result;
}

@Sovo · 02.05.2015, 11:43

Сообщение от wirty

При 4026318845 опять лажа.

Здесь вы испытываете переменную div на переполнение. Изначально код не рассчитан на такие большие значения.
Но я поправил ещё раз. Теперь, конечно, код гораздо более громоздкий, и нужно заново тестировать его производительность. Но это я смогу сделать только после всех майских праздников.

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
uint16_t sqrt_newton (uint32_t x)
{
uint32_t temp, div;
uint32_t rslt;
if (x > 0x0000FFFFL)
{
rslt = 0x0000FFFFL;
}
else
{
rslt = x;
}
if (x & 0xFFFF0000L)
{
if (x & 0xFF000000L)
{
div = 0x00003FFFL;
}
else
{
div = 0x000003FFL;
}
}
else
{
if (x & 0x0FF00L)
{
div = 0x0000003FL;
}
else
{
div = (x > 4) ? 0x00000007L : x;
}
}
while (1)
{
temp = x / div + div;
div = temp >> 1;
div += temp & 1;
int32_t error = rslt - div;
if ((error > 1) || (error < -1))
{
rslt = div;
}
else
{
if ((x / rslt == rslt - 1) && (x % rslt == 0))
{
rslt--;
}
return (uint16_t)rslt;
}
}
}

_pv, я проверю и этот вариант.

P.S. Если при определении предполагаемого результата учесть старший разряд и присвоить div = 0x7FFF, то число итераций уменьшается с семи до пяти.

@dork_usir · 02.05.2015, 12:16

Изначально код не рассчитан на такие большие значения.

Значит он тоже фигня... раз не выдерживается во всём диапазоне значений...

@_pv · 02.05.2015, 12:35

Сообщение от dork_usir

Сообщение от pv

а тупо бинарный перебор с умножением не быстрее будет чем через деление на stm8?

не быстрее
16х16 что-то около 70 тактов.
32/16 должно быть раз в 5-7 медленнее.
соответственно 16 умножений пожалуй всё-таки побыстрее будет чем 5 делений.

@dork_usir · 02.05.2015, 13:04

Под Кейл СТМ32 медленнее, но верный результат...
Под ИАР СТМ8 быстрее, но неверный результат...
.
Если вместо u16 result = 0xFFFF; написать u32 result = 0xFFFF; ... то результат под СТМ8 верный... но медленнее...
Так что... не быстрее...

@Sovo · 02.05.2015, 14:04

Благодарю за подсказки и указание на ошибки. Как только сравню на реальном устройстве, результаты отпишу сюда.

@Sovo · 02.05.2015, 14:08

Сообщение от _pv

16х16 что-то около 70 тактов.

ИМХО тут обязательно нужно 32*16, иначе результат останется в u16, и будет переполнение.
И мой IAR не даёт заглянуть в код подпрограммы умножения, только в режиме пошаговой отладки. Поэтому я не знаю, какой используется алгоритм умножения и деления 32 бит.

@_pv · 02.05.2015, 15:26

Сообщение от dork_usir

Под Кейл СТМ32 медленнее, но верный результат...
Под ИАР СТМ8 быстрее, но неверный результат...
Если вместо u16 result = 0xFFFF; написать u32 result = 0xFFFF; ... то результат под СТМ8 верный... но медленнее...
Так что... не быстрее...

ну так stm32 числа делить умеет одной инструкцией за 2-12 тактов, а stm8 - нет, так что не удивительно,
а вот то что компилятор не может два 16ти битных числа перемножить и результат сложить в 32 это косяк компилятора.

@wirty · 02.05.2015, 16:17

Можно компилятору подсказать
if ((u32)result * result < x) result += d; else result -= d;
Но всё равно, на малых числах Ньютон в 2 раза быстрее умножения. На больших числах - наоборот.

@Sovo · 02.05.2015, 16:48

Сообщение от _pv

а вот то что компилятор не может два 16ти битных числа перемножить и результат сложить в 32 это косяк компилятора.

Это не косяк. Я давно об этом читал. По правилам си разрядность произведения равна большей из разрядности множителей. Компилятор не должен своевольничать с типами переменных.

@_pv · 02.05.2015, 19:41

Сообщение от SOVO

Это не косяк. Я давно об этом читал. По правилам си разрядность произведения равна большей из разрядности множителей. Компилятор не должен своевольничать с типами переменных.

своевольничать не должен, но если ему известно, что множители 16бит, то делать 32х32->32 вместо 16х16->32, особенно на 8ми битном МК, тупо кучу времени умножая/складывая нули из старших байт чтобы потом их выкинуть - это косяк компилятора по части оптимизации.

http://www.ti.com/lit/an/spra683/spra683.pdf

@VyvotzorD · 02.05.2015, 19:57

Можно реализовать ряд Тейлора для функции sqrt(1+x)
Педивикия знает формулу :-)

По моему ряды Тейлора для разных функций - это то что надо, ибо итеративно и для каждой конкретной задачи можно ограничиться необходимой точностью в угоду быстродействию.

Новые блоги и статьи Все статьи Все блоги /
SDL3 для Web (WebAssembly): Реализация движения на Box2D v3 - трение и коллизии с повёрнутыми стенами 8Observer8 20.02.2026 Содержание блога Box2D позволяет легко создать главного героя, который не проходит сквозь стены и перемещается с заданным трением о препятствия, которые можно располагать под углом, как верхнее. . .	Конвертировать закладки radiotray-ng в m3u-плейлист damix 19.02.2026 Это можно сделать скриптом для PowerShell. Использование . \СonvertRadiotrayToM3U. ps1 <path_to_bookmarks. json> Рядом с файлом bookmarks. json появится файл bookmarks. m3u с результатом. # Check if. . .	Семь CDC на одном интерфейсе: 5 U[S]ARTов, 1 CAN и 1 SSI Eddy_Em 18.02.2026 Постепенно допиливаю свою "многоинтерфейсную плату". Выглядит вот так: https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11617&stc=1&d=1771445347 Основана на STM32F303RBT6. На борту пять. . .	Камера Toupcam IUA500KMA Eddy_Em 12.02.2026 Т. к. у всяких "хикроботов" слишком уж мелкий пиксель, для подсмотра в ESPriF они вообще плохо годятся: уже 14 величину можно рассмотреть еле-еле лишь на экспозициях под 3 секунды (а то и больше),. . .
И ясному Солнцу zbw 12.02.2026 И ясному Солнцу, и светлой Луне. В мире покоя нет и люди не могут жить в тишине. А жить им немного лет.	«Знание-Сила» zbw 12.02.2026 «Знание-Сила» «Время-Деньги» «Деньги -Пуля»	SDL3 для Web (WebAssembly): Подключение Box2D v3, физика и отрисовка коллайдеров 8Observer8 12.02.2026 Содержание блога Box2D - это библиотека для 2D физики для анимаций и игр. С её помощью можно определять были ли коллизии между конкретными объектами и вызывать обработчики событий столкновения. . . .	SDL3 для Web (WebAssembly): Загрузка PNG с прозрачным фоном с помощью SDL_LoadPNG (без SDL3_image) 8Observer8 11.02.2026 Содержание блога Библиотека SDL3 содержит встроенные инструменты для базовой работы с изображениями - без использования библиотеки SDL3_image. Пошагово создадим проект для загрузки изображения. . .

@_pv 0 / 0 / 0 Регистрация: 06.06.2011 Сообщений: 2,514
	22.04.2015, 16:57
	что-то там сильно не так. 2.5мс это 40000 тактов, /16 = 2500 тактов на цикл, чтобы пару 32х разрядных циферок сложить/сдвинуть. так не быть не должно 0

@Sovo 0 / 0 / 0 Регистрация: 05.10.2007 Сообщений: 498
	22.04.2015, 17:05
	Поправляюсь - четыре деления, четыре умножения, четыре извлечения корня и четыре сдвига на 10. Плюс прерывание отвлекает раз в мс на 0,2мс. Если вычесть прерывания, то останется - 0,8 мс против 2,1 мс 0

@Sovo 0 / 0 / 0 Регистрация: 05.10.2007 Сообщений: 498
	01.05.2015, 11:44
	Что, по вашему, не так? Есть реализация алгоритма ньютона без подбора результата, в два раза дольше. Результаты алгоритма ньютона и GLS совпадают. Что не так? 0

@wirty 1 / 1 / 0 Регистрация: 14.02.2013 Сообщений: 408
	02.05.2015, 01:52
	При 4026318845 опять лажа. 0

@dork_usir 0 / 0 / 0 Регистрация: 03.10.2012 Сообщений: 701
	02.05.2015, 13:04
	Под Кейл СТМ32 медленнее, но верный результат... Под ИАР СТМ8 быстрее, но неверный результат... . Если вместо u16 result = 0xFFFF; написать u32 result = 0xFFFF; ... то результат под СТМ8 верный... но медленнее... Так что... не быстрее... 0

@Sovo 0 / 0 / 0 Регистрация: 05.10.2007 Сообщений: 498
	02.05.2015, 14:04
	Благодарю за подсказки и указание на ошибки. Как только сравню на реальном устройстве, результаты отпишу сюда. 0

@wirty 1 / 1 / 0 Регистрация: 14.02.2013 Сообщений: 408
	02.05.2015, 16:17
	Можно компилятору подсказать if ((u32)result * result < x) result += d; else result -= d; Но всё равно, на малых числах Ньютон в 2 раза быстрее умножения. На больших числах - наоборот. 0

@VyvotzorD 0 / 0 / 0 Регистрация: 06.12.2016 Сообщений: 2,309
	02.05.2015, 19:57
	Можно реализовать ряд Тейлора для функции sqrt(1+x) Педивикия знает формулу :-) По моему ряды Тейлора для разных функций - это то что надо, ибо итеративно и для каждой конкретной задачи можно ограничиться необходимой точностью в угоду быстродействию. 0