Извлечение квадратного корня. Скорость вычисления на 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)
Педивикия знает формулу :-)

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

Новые блоги и статьи Все статьи Все блоги /
Новый CodeBlocs. Версия 25.03 palva 04.01.2026 Оказывается, недавно вышла новая версия CodeBlocks за номером 25. 03. Когда-то давно я возился с только что вышедшей тогда версией 20. 03. С тех пор я давно снёс всё с компьютера и забыл. Теперь. . .	Модель микоризы: классовый агентный подход anaschu 02.01.2026 Раньше это было два гриба и бактерия. Теперь три гриба, растение. И на уровне агентов добавится между грибами или бактериями взаимодействий. До того я пробовал подход через многомерные массивы,. . .	Учёным и волонтёрам проекта «Einstein@home» удалось обнаружить четыре гамма-лучевых пульсара в джете Млечного Пути Programma_Boinc 01.01.2026 Учёным и волонтёрам проекта «Einstein@home» удалось обнаружить четыре гамма-лучевых пульсара в джете Млечного Пути Сочетание глобально распределённой вычислительной мощности и инновационных. . .	Советы по крайней бережливости. Внимание, это ОЧЕНЬ длинный пост. Programma_Boinc 28.12.2025 Советы по крайней бережливости. Внимание, это ОЧЕНЬ длинный пост. Налог на собак: https:/ / ********/ gallery/ V06K53e Финансовый отчет в Excel: https:/ / ********/ gallery/ bKBkQFf Пост отсюда. . .	Кто-нибудь знает, где можно бесплатно получить настольный компьютер или ноутбук? США. Programma_Boinc 26.12.2025 Нашел на реддите интересную статью под названием Anyone know where to get a free Desktop or Laptop? Ниже её машинный перевод. После долгих разбирательств я наконец-то вернула себе. . .
Thinkpad X220 Tablet — это лучший бюджетный ноутбук для учёбы, точка. Programma_Boinc 23.12.2025 Рецензия / Мнение/ Перевод Нашел на реддите интересную статью под названием The Thinkpad X220 Tablet is the best budget school laptop period . Ниже её машинный перевод. Thinkpad X220 Tablet —. . .	PhpStorm 2025.3: WSL Terminal всегда стартует в ~ and_y87 14.12.2025 PhpStorm 2025. 3: WSL Terminal всегда стартует в ~ (home), игнорируя директорию проекта Симптом: После обновления до PhpStorm 2025. 3 встроенный терминал WSL открывается в домашней директории. . .	Как объединить две одинаковые БД Access с разными данными VikBal 11.12.2025 Помогите пожалуйста !! Как объединить 2 одинаковые БД Access с разными данными.	Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro	Музыка, написанная Искусственным Интеллектом volvo 04.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .

@_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