Поиск эффективного целочисленного алгоритма квадратного корня для ARM Thumb2
Самый простой способ чтения из файла и записи в файл:
//Read from a file
string something = File.ReadAllText("C:\\Rfile.txt");
//Write to a file
using (StreamWriter writer = new StreamWriter("Wfile.txt"))
{
writer.WriteLine(something);
}
6 ответов
Я разработал 16-разрядный sqrt для гамма сжатия RGB. Это диспетчеризирует в 3 различных таблицы, на основе более высоких 8 битов. Недостатки: это использует приблизительно килобайт для таблиц, непредсказуемые раунды, если точный sqrt невозможен, и, в худшем случае, использует единственное умножение (но только для очень немногих входных значений).
static uint8_t sqrt_50_256[] = {
114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,
133,134,135,136,137,138,139,140,141,142,143,143,144,145,146,147,148,149,150,
150,151,152,153,154,155,155,156,157,158,159,159,160,161,162,163,163,164,165,
166,167,167,168,169,170,170,171,172,173,173,174,175,175,176,177,178,178,179,
180,181,181,182,183,183,184,185,185,186,187,187,188,189,189,190,191,191,192,
193,193,194,195,195,196,197,197,198,199,199,200,201,201,202,203,203,204,204,
205,206,206,207,207,208,209,209,210,211,211,212,212,213,214,214,215,215,216,
217,217,218,218,219,219,220,221,221,222,222,223,223,224,225,225,226,226,227,
227,228,229,229,230,230,231,231,232,232,233,234,234,235,235,236,236,237,237,
238,238,239,239,240,241,241,242,242,243,243,244,244,245,245,246,246,247,247,
248,248,249,249,250,250,251,251,252,252,253,253,254,254,255,255
};
static uint8_t sqrt_0_10[] = {
1,2,3,3,4,4,5,5,5,6,6,6,7,7,7,7,8,8,8,8,9,9,9,9,9,10,10,10,10,10,11,11,11,
11,11,11,12,12,12,12,12,12,13,13,13,13,13,13,13,14,14,14,14,14,14,14,15,15,
15,15,15,15,15,15,16,16,16,16,16,16,16,16,17,17,17,17,17,17,17,17,17,18,18,
18,18,18,18,18,18,18,19,19,19,19,19,19,19,19,19,19,20,20,20,20,20,20,20,20,
20,20,21,21,21,21,21,21,21,21,21,21,21,22,22,22,22,22,22,22,22,22,22,22,23,
23,23,23,23,23,23,23,23,23,23,23,24,24,24,24,24,24,24,24,24,24,24,24,25,25,
25,25,25,25,25,25,25,25,25,25,25,26,26,26,26,26,26,26,26,26,26,26,26,26,27,
27,27,27,27,27,27,27,27,27,27,27,27,27,28,28,28,28,28,28,28,28,28,28,28,28,
28,28,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,30,30,30,30,30,30,30,30,
30,30,30,30,30,30,30,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,32,32,
32,32,32,32,32,32,32,32,32,32,32,32,32,32,33,33,33,33,33,33,33,33,33,33,33,
33,33,33,33,33,33,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,35,35,
35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,36,36,36,36,36,36,36,36,36,
36,36,36,36,36,36,36,36,36,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,
37,37,37,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,39,39,39,
39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,40,40,40,40,40,40,40,40,
40,40,40,40,40,40,40,40,40,40,40,40,41,41,41,41,41,41,41,41,41,41,41,41,41,
41,41,41,41,41,41,41,41,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,
42,42,42,42,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,
43,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,45,45,
45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,46,46,46,46,
46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,47,47,47,47,47,47,
47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,48,48,48,48,48,48,48,
48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,49,49,49,49,49,49,49,49,
49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,50,50,50,50,50,50,50,50,
50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,51,51,51,51,51,51,51,51,
51,51,51,51,51,51,51,51,51,51,51,51,51,51,51,51,51,51,52,52,52,52,52,52,52,
52,52,52,52,52,52,52,52,52,52,52,52,52,52,52,52,52,52,52,53,53
};
static uint8_t sqrt_11_49[] = {
54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,0,76,77,78,
0,79,80,81,82,83,0,84,85,86,0,87,88,89,0,90,0,91,92,93,0,94,0,95,96,97,0,98,0,
99,100,101,0,102,0,103,0,104,105,106,0,107,0,108,0,109,0,110,0,111,112,113
};
uint16_t isqrt16(uint16_t v) {
uint16_t a, b;
uint16_t h = v>>8;
if (h <= 10) return v ? sqrt_0_10[v>>2] : 0;
if (h >= 50) return sqrt_50_256[h-50];
h = (h-11)<<1;
a = sqrt_11_49[h];
b = sqrt_11_49[h+1];
if (!a) return b;
return b*b > v ? a : b;
}
я сравнил его с log2 базирующийся sqrt, с помощью лязга __ builtin_clz (который должен расшириться до единственного кода операции блока), и sqrtf библиотеки, названный использованием (интервала) sqrtf ((плавание) i). И получил довольно странные результаты:
$ gcc-o3 test.c-o тестируют & &./test isqrt16: 6.498955 sqrtf: 6,981861 log2_sqrt: 61.755873
Лязг скомпилировал вызов в sqrtf к sqrtss инструкции, которая является почти с такой скоростью, как та таблица sqrt. Урок извлечен: на x86 компиляторе может обеспечить достаточно быстро sqrt, который является меньше чем %10 медленнее, чем, что Вы сами можете придумать, тратя впустую много времени, или можете быть в 10 раз быстрее при использовании некоторых ужасных поразрядных взломов. И все еще sqrtss немного медленнее, чем пользовательская функция, поэтому при реальной необходимости в этих 5% можно получить их, и ARM, например, не имеет никакого sqrtss, таким образом, log2_sqrt не должен изолировать это плохо.
На x86, где FPU доступен, старый взлом Quake , http://h14s.p5r.org/2012/09/0x5f3759df.html?mwh=1 , кажется, самый быстрый способ вычислить целое число sqrt. Это в 2 раза быстрее, чем эта таблица или sqrtss FPU.
Если точная точность не требуется, у меня есть быстрое приближение для вас, которое использует 260 байт ОЗУ (вы можете уменьшить это вдвое, но не делайте этого).
int ftbl[33]={0,1,1,2,2,4,5,8,11,16,22,32,45,64,90,128,181,256,362,512,724,1024,1448,2048,2896,4096,5792,8192,11585,16384,23170,32768,46340};
int ftbl2[32]={ 32768,33276,33776,34269,34755,35235,35708,36174,36635,37090,37540,37984,38423,38858,39287,39712,40132,40548,40960,41367,41771,42170,42566,42959,43347,43733,44115,44493,44869,45241,45611,45977};
int fisqrt(int val)
{
int cnt=0;
int t=val;
while (t) {cnt++;t>>=1;}
if (6>=cnt) t=(val<<(6-cnt));
else t=(val>>(cnt-6));
return (ftbl[cnt]*ftbl2[t&31])>>15;
}
Вот код для создания таблиц:
ftbl[0]=0;
for (int i=0;i<32;i++) ftbl[i+1]=sqrt(pow(2.0,i));
printf("int ftbl[33]={0");
for (int i=0;i<32;i++) printf(",%d",ftbl[i+1]);
printf("};\n");
for (int i=0;i<32;i++) ftbl2[i]=sqrt(1.0+i/32.0)*32768;
printf("int ftbl2[32]={");
for (int i=0;i<32;i++) printf("%c%d",(i)?',':' ',ftbl2[i]);
printf("};\n");
В диапазоне 1 → 2 20 максимальная ошибка составляет 11, а в диапазоне 1 → 2 30 - около 256. Вы можете использовать таблицы большего размера и минимизируйте это. Стоит отметить, что ошибка всегда будет отрицательной - то есть, когда это неверно, значение будет МЕНЬШЕ, чем правильное значение.
Возможно, вы сделаете это, выполнив этап уточнения.
Идея достаточно проста: (ab) 0,5 = a 0.b × b 0,5 .
Итак, мы берем входные данные X = A × B, где A = 2 N и 1 ≤ B <2
Тогда у нас есть таблица поиска для sqrt (2 N ), и таблица поиска для sqrt (1 ≤ B <2). Мы сохраняем таблицу поиска для sqrt (2 N ) как целое число, что может быть ошибкой (тестирование не выявило вредных последствий), и мы сохраняем таблицу поиска для sqrt (1 ≤ B <2) как 15 -битовая фиксированная точка.
Мы знаем, что 1 ≤ sqrt (2 N ) <65536, так что это 16-битный, и мы знаем, что действительно можем умножить только 16-битные × 15-битные на ARM, не опасаясь репрессий, вот что мы делаем.
С точки зрения реализации, while (t) {cnt ++; t >> = 1;} фактически является ведущим -bit инструкция (CLB), так что если в вашей версии чипсета она есть, вы выигрываете! Кроме того, инструкцию сдвига можно легко реализовать с помощью двунаправленного переключателя, если он у вас есть?
Здесь есть алгоритм Lg [N] для подсчета самого высокого установленного бита .
Один из распространенных подходов - деление пополам.
hi = number
lo = 0
mid = ( hi + lo ) / 2
mid2 = mid*mid
while( lo < hi-1 and mid2 != number ) {
if( mid2 < number ) {
lo = mid
else
hi = mid
mid = ( hi + lo ) / 2
mid2 = mid*mid
Что-то вроде этого должно работать достаточно хорошо. Он выполняет тесты log2 (число), выполняя
log2 (число) умножается и делится. Поскольку деление - это деление на 2, вы можете заменить его на >> .
Условие завершения может быть неправильным, поэтому обязательно проверьте множество целых чисел, чтобы убедиться, что деление на 2 не колеблется между двумя четными значениями неправильно; они будут отличаться более чем на 1.
Целочисленные квадратные корни Джека В. Креншоу могут быть полезны в качестве еще одной ссылки.
Архив фрагментов C также имеет реализацию целочисленного квадратного корня . Это выходит за рамки просто целочисленного результата и вычисляет дополнительные дробные (с фиксированной точкой) биты ответа. (Обновление: к сожалению, архив фрагментов C теперь не функционирует. Ссылка указывает на веб-архив страницы.) Вот код из архива фрагментов C:
#define BITSPERLONG 32
#define TOP2BITS(x) ((x & (3L << (BITSPERLONG-2))) >> (BITSPERLONG-2))
struct int_sqrt {
unsigned sqrt, frac;
};
/* usqrt:
ENTRY x: unsigned long
EXIT returns floor(sqrt(x) * pow(2, BITSPERLONG/2))
Since the square root never uses more than half the bits
of the input, we use the other half of the bits to contain
extra bits of precision after the binary point.
EXAMPLE
suppose BITSPERLONG = 32
then usqrt(144) = 786432 = 12 * 65536
usqrt(32) = 370727 = 5.66 * 65536
NOTES
(1) change BITSPERLONG to BITSPERLONG/2 if you do not want
the answer scaled. Indeed, if you want n bits of
precision after the binary point, use BITSPERLONG/2+n.
The code assumes that BITSPERLONG is even.
(2) This is really better off being written in assembly.
The line marked below is really a "arithmetic shift left"
on the double-long value with r in the upper half
and x in the lower half. This operation is typically
expressible in only one or two assembly instructions.
(3) Unrolling this loop is probably not a bad idea.
ALGORITHM
The calculations are the base-two analogue of the square
root algorithm we all learned in grammar school. Since we're
in base 2, there is only one nontrivial trial multiplier.
Notice that absolutely no multiplications or divisions are performed.
This means it'll be fast on a wide range of processors.
*/
void usqrt(unsigned long x, struct int_sqrt *q)
{
unsigned long a = 0L; /* accumulator */
unsigned long r = 0L; /* remainder */
unsigned long e = 0L; /* trial product */
int i;
for (i = 0; i < BITSPERLONG; i++) /* NOTE 1 */
{
r = (r << 2) + TOP2BITS(x); x <<= 2; /* NOTE 2 */
a <<= 1;
e = (a << 1) + 1;
if (r >= e)
{
r -= e;
a++;
}
}
memcpy(q, &a, sizeof(long));
}
Я остановился на следующем коде. По сути, это из статьи Википедии о методах вычисления квадратного корня . Но он был изменен для использования stdint.h типов uint32_t и т. Д. Строго говоря, тип возврата можно было изменить на uint16_t .
/**
* \brief Fast Square root algorithm
*
* Fractional parts of the answer are discarded. That is:
* - SquareRoot(3) --> 1
* - SquareRoot(4) --> 2
* - SquareRoot(5) --> 2
* - SquareRoot(8) --> 2
* - SquareRoot(9) --> 3
*
* \param[in] a_nInput - unsigned integer for which to find the square root
*
* \return Integer square root of the input value.
*/
uint32_t SquareRoot(uint32_t a_nInput)
{
uint32_t op = a_nInput;
uint32_t res = 0;
uint32_t one = 1uL << 30; // The second-to-top bit is set: use 1u << 14 for uint16_t type; use 1uL<<30 for uint32_t type
// "one" starts at the highest power of four <= than the argument.
while (one > op)
{
one >>= 2;
}
while (one != 0)
{
if (op >= res + one)
{
op = op - (res + one);
res = res + 2 * one;
}
res >>= 1;
one >>= 2;
}
return res;
}
Хорошая вещь , Я обнаружил, заключается в том, что довольно простая модификация может вернуть «округлый» ответ. Я нашел это полезным в определенном приложении для большей точности. Обратите внимание, что в этом случае возвращаемый тип должен быть uint32_t , потому что округленный квадратный корень из 2 32 - 1 равен 2 16 .
/**
* \brief Fast Square root algorithm, with rounding
*
* This does arithmetic rounding of the result. That is, if the real answer
* would have a fractional part of 0.5 or greater, the result is rounded up to
* the next integer.
* - SquareRootRounded(2) --> 1
* - SquareRootRounded(3) --> 2
* - SquareRootRounded(4) --> 2
* - SquareRootRounded(6) --> 2
* - SquareRootRounded(7) --> 3
* - SquareRootRounded(8) --> 3
* - SquareRootRounded(9) --> 3
*
* \param[in] a_nInput - unsigned integer for which to find the square root
*
* \return Integer square root of the input value.
*/
uint32_t SquareRootRounded(uint32_t a_nInput)
{
uint32_t op = a_nInput;
uint32_t res = 0;
uint32_t one = 1uL << 30; // The second-to-top bit is set: use 1u << 14 for uint16_t type; use 1uL<<30 for uint32_t type
// "one" starts at the highest power of four <= than the argument.
while (one > op)
{
one >>= 2;
}
while (one != 0)
{
if (op >= res + one)
{
op = op - (res + one);
res = res + 2 * one;
}
res >>= 1;
one >>= 2;
}
/* Do arithmetic rounding to nearest integer */
if (op > res)
{
res++;
}
return res;
}
Я остановился на чем-то похожем на алгоритм двоичного числа за цифрой, описанный в этой статье в Википедии .
Это не быстро, но мало и просто:
int isqrt(int n)
{
int b = 0;
while(n >= 0)
{
n = n - b;
b = b + 1;
n = n - b;
}
return b - 1;
}