Comments 29
микроконтроллерам stm32
Вы бы конкретизировали хотя бы процессорное ядро, а то у STM32 есть Cortex M0, Cortex M0+, Cortex M3, Cortex M4F, Cortex M7, Cortex A57. Мне пришлось листать до конца документа, чтобы понять, что речь о Cortex M0.
Микроконтроллеры работают, преимущественно, с целыми числами; регистров для действительных чисел у них, как правило, нет.
У M4F, M7 FPU есть. Как и операция быстрого вычисления квадратного корня (VSQRT.F32).
double sqrt (double num);
Еще и двойной точности. Ее стоит использовать только на Cortex M7, оборудованных DP-FPU.
У M4F, M7 FPU есть. Как и операция быстрого вычисления квадратного корня (VSQRT.F32).
Это из всей-то линейки? Выходит все же, что
регистров для действительных чисел у них, как правило, нет.
У Cortex A-серии FPU, как правило, есть. Просто у ST есть единственный STM32MP15x с Cortex A7 на борту.
Если разработчику не критично микропотребление, но при этом нужна математика — он берет ядро с FPU, а не пишет монографию об оптимизации расчета. Стоимость даже Cortex M4F довольно низка, как правило обвязка вокруг микросхемы (собственно сам прибор) стоит дороже
есть единственный STM32MP15x с Cortex A7))
Еще стоит заметить, что кроме разработчиков (скорее даже в большей степени) здесь присутствуют хоббисты/любители/студенты и подход к выбору компонентов у которых несколько иной — как правило, что есть в наличии и не дорогое то и используется. Гляньте стоимость приведенного вами STM32MP15x если брать не партию.
Стараюсь понять Вашу мысль, Dima_Sharihin.
Это очевидно, что хорошая программа будет работать даже на бюджетном железе.
Однако сомнительно, что премиальное железо сделает плохую программу лучше.
Статья о приёме программирования, сокращающем число итераций через увеличение размера кода до цикла. В данном контексте железо служит фоном, подсвечивающим данный метод.
Ну, вообще есть готовые библиотеки математики с фиксированной точкой. А-ля IQmath.
Согласен. Существует много хороших программ.
Вопрос выбора.
Однако, вводит в некоторое заблуждение отсутствие конкретики в замечаниях.
Думаю, многим было бы интересно узнать, во сколько раз алгоритм sqrt_fpu опережает sqrt_env на премиальных платформах Cortex?
Или, во сколько раз IQmath быстрее на бюджетном микроконтроллере, чем sqrt_env?
Сколько требуются вычислительных ресурсов и какие появляются риски в том или ином случае при производстве устройств для работы на ответственных участках?
Возвращаясь к прикладной задаче.
Извлечение квадратного корня из целого
с округлением до ближайшего целого.Задача была решена в установленных технических условиях. Метод и порядок решения раскрыты в публикации. Результат можно подвергать сомнениям, проверять и улучшать.
Допускаю, что существуют отличные решения для отличных технических условий. Буду рад случаю узнать о них что-нибудь конкретно.
rslt += L < 0 ? 0 : 1;
L объявлен как uint32_t, компилятор выкинет эту проверку, а ещё возможно переполнение, если div > L.
Для целого аргумента без знака функция sqrt_odd может быть, например, такой:
uint16_t sqrt_odd ( uint32_t L )
{
if ( L < 2 )
return ( uint16_t ) L;
uint16_t div = 1, rslt = 1;
while ( 1 )
{
div += 2;
if ( ( uint32_t ) div >= L )
return rslt;
L -= div, rslt++;
}
}Спасибо за внимание к моей статье. Ошибку в тексте исправил.
Как вариант можно чуть упростить и на первой итерации выкинуть деление:
//tested on VC Express, Win7
//max error: 0.999989
int isqrt(long x)
{
if (x == 0) return 0;
long xout = 1;//начальное приближение
int por = 0;//порядок начального приближения
//уточняем его в зависимости от двоичного порядка входного аргумента:
if (x & 0xFFFF0000)
{
int nx = x >> 16;
xout = (1 << 8);
por = 8;
while (nx != 0)
{
xout = xout << 1;
nx = nx >> 2;
por++;
}
}
else
{
int nx = x;
while (nx != 0)
{
xout = xout << 1;
nx = nx >> 2;
por++;
}
}
//далее выполняем итерации методом Ньютона:
xout = ((xout + (x >> por)) >> 1);//первая итерация - деление заменяем сдвигом, т.к. начальное приближение xout = 2^por
xout = ((xout + x / xout) >> 1);//вторая
xout = ((xout + x / xout) >> 1);//третья
xout = ((xout + x / xout) >> 1);//четвертая
//xout = ((xout + x / xout) >> 1);//пятая - лишняя, для 32-х бит хватает 4-х
return xout;
}
p.s. стмок под рукой нету, поэтому тестировал и писал в vc express edition, при копировании на стм надо будет изменить типы на соответствующие int -16bit, long -32bit.
FGV, искренне рад Вашему интересу к прикладной задаче.
Ускорение работы есть, но можно ли считать задачу решённой?
Если следовать строго букве задачи:
Извлечение корня из целого
с округлением до ближайшего целого
На платформе, обозначенной в статье, фиксируются случайные ошибки округления на всём множестве значений аргумента.
test for quality with in range[0..1000]
isqrt(7)=2 <=> 3 (2.65)
isqrt(8)=2 <=> 3 (2.83)
isqrt(13)=3 <=> 4 (3.61)
isqrt(14)=3 <=> 4 (3.74)
too many mistakes ... test failed.
test for quality with in range[1000..10000]
isqrt(1000)=31 <=> 32 (31.62)
isqrt(1001)=31 <=> 32 (31.64)
isqrt(1002)=31 <=> 32 (31.65)
isqrt(1003)=31 <=> 32 (31.67)
too many mistakes ... test failed.
test for quality with in range[10000..100000]
isqrt(10101)=100 <=> 101 (100.50)
isqrt(10102)=100 <=> 101 (100.51)
isqrt(10103)=100 <=> 101 (100.51)
isqrt(10104)=100 <=> 101 (100.52)
too many mistakes ... test failed.
***********
test range=[0..1000], repeat=100
_fpu ... , _evn ... , isqrt ... ,
test range=[0..10000], repeat=10
_fpu ... , _evn ... , isqrt ... ,
test range=[0..100000], repeat=1
_fpu ... , _evn ... , isqrt ... ,
***********
range _fpu _evn isqrt
1000 5120 680 619
10000 5142 782 693
100000 5148 922 716
done.
на стм надо будет изменить типы на соответствующие int -16bit, long -32bit
Типы переменных в коде оставил как есть. Если правильная работа Вашего кода зависит от размерности и знака целых, назначьте нужные типы, пользуясь определениями из <stdint.h>.
Извлечение корня из целого
с округлением до ближайшего целого
Тогда для прохождения тестов чуть допилить надо:
int isqrt(long x)
{
x=x<<4;//mul x4
...
return (xout+2)>>2;//div 2
}
в этом случае максимальная ошибка для чисел от 0 до 100000 не превышает 0.5.
int isqrt(long x)
{
if (x == 0) return 0;
...
int er0=x-xout*xout;
int er1=x-(xout+1)*(xout+1);
if (er0<0) er0=-er0;
if (er1<0) er1=-er1;
if (er0>er1) return xout+1;
return xout;
}
Публикуйте, пожалуйста код функции без пропусков (многоточий). Теряется контекст.
Относительно округления результата до ближайшего целого — без вариантов, на всём множестве значений аргумента [ 0… 0xFFFFFFFFU ].
int isqrt(long x)
{
if (x == 0) return 0;
long xout = 1;//начальное приближение
int por = 0;//порядок начального приближения
//уточняем его в зависимости от двоичного порядка входного аргумента:
if (x & 0xFFFF0000)
{
int nx = x >> 16;
xout = (1 << 8);
por = 8;
while (nx != 0)
{
xout = xout << 1;
nx = nx >> 2;
por++;
}
}
else
{
int nx = x;
while (nx != 0)
{
xout = xout << 1;
nx = nx >> 2;
por++;
}
}
//далее выполняем итерации методом Ньютона:
xout = ((xout + (x >> por)) >> 1);//первая итерация - деление заменяем сдвигом, т.к. начальное приближение xout = 2^por
xout = ((xout + x / xout) >> 1);//вторая
xout = ((xout + x / xout) >> 1);//третья
xout = ((xout + x / xout) >> 1);//четвертая
//xout = ((xout + x / xout) >> 1);//пятая - лишняя, для 32-х бит хватает 4-х
//коррекция результата до ближайшего целого:
int er0=x-xout*xout;
int er1=x-(xout+1)*(xout+1);
if (er0<0) er0=-er0;
if (er1<0) er1=-er1;
if (er0>er1) return xout+1;
return xout;
}//Возвращаемый результат - 32 бита, т.к. для
//больших входных чисел вроде 65536*65536-1,-2,...
//с учетом максимальной ошибки по модулю не более 0.5
//возвращаемый результат должен быть равен 65536,
//что не укладывается в 16 бит.
uint32_t isqrt(uint32_t x)
{
if (x == 0) return 0;
uint32_t xout = 1;//начальное приближение
uint16_t por = 0;//порядок начального приближения
//уточняем начальное приближение и порядок в зависимости от порядка входного аргумента:
if (x & 0xFFFF0000)
{
uint16_t nx = x >> 16;
xout = (1 << 8);
por = 8;
while (nx != 0)
{
xout = xout << 1;
nx = nx >> 2;
por++;
}
}
else
{
uint16_t nx = x;
while (nx != 0)
{
xout = xout << 1;
nx = nx >> 2;
por++;
}
}
//далее выполняем итерации методом Ньютона:
xout = ((xout + (x >> por)) >> 1);//первая итерация - деление заменяем сдвигом, т.к. начальное приближение xout = 2^por
xout = ((xout + x / xout) >> 1);//вторая
xout = ((xout + x / xout) >> 1);//третья
xout = ((xout + x / xout) >> 1);//четвертая
//xout = ((xout + x / xout) >> 1);//пятая - лишняя, для 32-х бит хватает 4-х
//коррекция результата:
uint32_t er0 = x - xout*xout;//ошибка для корня выч. с отбрасыванием дробной части
uint32_t er1 = er0-(xout<<1)-1;//ошибка для корня выч. с отбрасыванием дробной части + 1
if (er0 & 0x80000000) er0=(~er0)+1;//приведение результата к положительному
if (er1 & 0x80000000) er1 = (~er1) + 1;//
if (er0<er1) return xout;
return xout+1;
}
Когда-то давно писал код для вычисления целочисленного квадратного корня
http://cdeblog.ru/user-sqrt
uint16_t sqrxi32(uint32_t y)
{
uint32_t xh = y>0xffffu?0xffff:y;
uint32_t xl = 0;
uint32_t xc;
for(int k=0;k<17;k++)
{
xc = (xh+xl)>>1;
if(xc*xc>y)
{
xh = xc;
}
else
{
xl = xc;
}
}
return xc;
}
repeat 9000
uint16_t sqrxi32(uint32_t y)
{
uint32_t xh = y>0xffffu?0xffffu:y;
uint32_t xl = 0;
uint32_t xc;
for(int k=0;k<16;k++)
{
xc = (xh+xl)>>1;
if(xc*xc>y)
{
xh = xc;
}
else
{
xl = xc;
}
}
return (xh+xl)>>1;
}
sqrt_new :336 ms
sqrt_evn :133 ms
isqrt: 93 ms
sqrxi32: 28 ms
Замечательный, быстрый код.
Вы проверяли алгоритм на отсутствие ошибок при округлении?
В качестве эталона для сверки лучше брать стандартную функцию из "math.h". В статье, как эталон, применялась sqrt_fpu(...).
Эта немного медленнее, но с округлением к nearest.
Тест точность на PC.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>
uint16_t sqrxi32(uint32_t y)
{
if(y==1) return 1;
uint32_t xh = y>0x10000ul?0x10000ul:y;
uint32_t xl = 0;
uint32_t xc ;
for(int k=0;k<16;k++)
{
xc = (xh+xl)>>1ul;
if(xc*xc-xc>=y)
{
xh = xc;
}
else
{
xl = xc;
}
}
return (xh+xl)>>1ul;
}
int main()
{
printf("%x %x %x\n",sqrxi32(0),sqrxi32(1),sqrxi32(0xfffffffful));
for(uint32_t y=0;y<0xfffffffful;y++)
{
uint32_t xUint = sqrxi32(y);
double xDoub = sqrt(y);
if(fabs(xDoub-xUint)>0.5)
{
printf("error in y = %x xUint =%x xDoub =%f \n",y,xUint,xDoub);
return 1;
}
}
printf("end OK\n");
return 0;
}Output:
0 1 ffff
error in y = ffff0001 xUint =ffff xDoub =65535.500006
— как и ожидалось!
Вообще-то на библиотеки в таких случаях лучше не полагаться.
А классический метод нахождения корней функции(любой функции!) методом деления отрезка пополам нужно знать.
sqrxi32 3107 ms
sqrxi32 nearest 3109 ms
isqrt 10431 ms
sqrt_evn 14731 ms
sqrt_new 37283 ms
Вот замеры с stm32f407 ( with fp32 fpu )168MHz на миллион чисел:
sqrxi32 = 1217 ms
sqrxi32 nearest = 1331 ms
isqrt = 776 ms
sqrt_evn = 625 ms
sqrt_new = 1494 ms
math library sqrtf = 382 ms
//sqrtf — 24 bit mantissa.
stm32. Смотрим в корень