sqrt(n)의 정수 부분을 얻는 가장 빠른 방법?
있듯이n한 정사각형이 , 「정사각형」이 됩니다.sqrt(n)수가아아 아아아다다에, 발신자가 「」라고 하는 .sqrt(n)부분 계산에도 시간이 걸리기 때문에 그리 빠르지 않을 것입니다.
그래서 제 질문은,
실제 값을 계산하지 않고 sqrt(n)의 정수 부분만 얻을 수 있습니까?sqrt(n)이 더 .sqrt(n):<math.h> ★★★★★★★★★★★★★★★★★」<cmath>
쓸 수 있어요.asm블록도 마찬가지입니다.
나는 빠른 역제곱근 기술을 시도해 볼 것이다.
를 얻을 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수, 수 등.1/sqrt(n)(32비트와 64비트 플랫폼 사이에 존재하는) 일부 비트 트위들링에 기반합니다.
결과를 얻으면 결과를 반전시키면 정수 부분을 얻을 수 있습니다.
물론 좀 더 빠른 방법이 있을 수 있습니다. 이 방법은 좀 둥글기 때문입니다.
편집: 해보자!
먼저 작은 도우미:
// benchmark.h
#include <sys/time.h>
template <typename Func>
double benchmark(Func f, size_t iterations)
{
f();
timeval a, b;
gettimeofday(&a, 0);
for (; iterations --> 0;)
{
f();
}
gettimeofday(&b, 0);
return (b.tv_sec * (unsigned int)1e6 + b.tv_usec) -
(a.tv_sec * (unsigned int)1e6 + a.tv_usec);
}
그리고 본체:
#include <iostream>
#include <cmath>
#include "benchmark.h"
class Sqrt
{
public:
Sqrt(int n): _number(n) {}
int operator()() const
{
double d = _number;
return static_cast<int>(std::sqrt(d) + 0.5);
}
private:
int _number;
};
// http://www.codecodex.com/wiki/Calculate_an_integer_square_root
class IntSqrt
{
public:
IntSqrt(int n): _number(n) {}
int operator()() const
{
int remainder = _number;
if (remainder < 0) { return 0; }
int place = 1 <<(sizeof(int)*8 -2);
while (place > remainder) { place /= 4; }
int root = 0;
while (place)
{
if (remainder >= root + place)
{
remainder -= root + place;
root += place*2;
}
root /= 2;
place /= 4;
}
return root;
}
private:
int _number;
};
// http://en.wikipedia.org/wiki/Fast_inverse_square_root
class FastSqrt
{
public:
FastSqrt(int n): _number(n) {}
int operator()() const
{
float number = _number;
float x2 = number * 0.5F;
float y = number;
long i = *(long*)&y;
//i = (long)0x5fe6ec85e7de30da - (i >> 1);
i = 0x5f3759df - (i >> 1);
y = *(float*)&i;
y = y * (1.5F - (x2*y*y));
y = y * (1.5F - (x2*y*y)); // let's be precise
return static_cast<int>(1/y + 0.5f);
}
private:
int _number;
};
int main(int argc, char* argv[])
{
if (argc != 3) {
std::cerr << "Usage: %prog integer iterations\n";
return 1;
}
int n = atoi(argv[1]);
int it = atoi(argv[2]);
assert(Sqrt(n)() == IntSqrt(n)() &&
Sqrt(n)() == FastSqrt(n)() && "Different Roots!");
std::cout << "sqrt(" << n << ") = " << Sqrt(n)() << "\n";
double time = benchmark(Sqrt(n), it);
double intTime = benchmark(IntSqrt(n), it);
double fastTime = benchmark(FastSqrt(n), it);
std::cout << "Number iterations: " << it << "\n"
"Sqrt computation : " << time << "\n"
"Int computation : " << intTime << "\n"
"Fast computation : " << fastTime << "\n";
return 0;
}
그 결과:
sqrt(82) = 9
Number iterations: 4096
Sqrt computation : 56
Int computation : 217
Fast computation : 119
// Note had to tweak the program here as Int here returns -1 :/
sqrt(2147483647) = 46341 // real answer sqrt(2 147 483 647) = 46 340.95
Number iterations: 4096
Sqrt computation : 57
Int computation : 313
Fast computation : 119
예상대로 고속 연산이 Int 연산보다 훨씬 더 잘 수행됩니다.
그건 그렇고.sqrt 빠릅니다 is is :)
이 . - 이 답변은 바보같다. 사용(int) sqrt(i)
적절한 설정으로 프로파일링 후(-march=native -m64 -O3위의 내용이 훨씬 빨랐습니다.
네, 조금 오래된 질문이지만, "가장 빠른" 대답은 아직 나오지 않았습니다.가장 빠른(내 생각에) 것은 이 Embedded.com 기사에서 자세히 설명한 바이너리 제곱근 알고리즘입니다.
기본적으로는 다음과 같습니다.
unsigned short isqrt(unsigned long a) {
unsigned long rem = 0;
int root = 0;
int i;
for (i = 0; i < 16; i++) {
root <<= 1;
rem <<= 2;
rem += a >> 30;
a <<= 2;
if (root < rem) {
root++;
rem -= root;
root++;
}
}
return (unsigned short) (root >> 1);
}
제 머신(Q6600, Ubuntu 10.10)에서는 1-1000000의 제곱근을 사용하여 프로파일을 작성했습니다.「」를 사용합니다.iqsrt(i)2750ms2,를 사용하였습니다.(unsigned short) sqrt((float) i)3600번으로 , 하다, 하다를 사용해서 한예요.g++ -O3를 합니다.-ffast-math2100밀리초 3100밀리초이것은 조립품 라인 한 개도 사용하지 않기 때문에 훨씬 더 빠를 수 있습니다.
위의 코드는 C와 C++ 모두에서 동작하며 Java에서도 약간의 구문 변경으로 동작합니다.
한정된 범위에서는 바이너리 검색이 더욱 효과적입니다.제 기계에서는 위의 버전이 4배만큼 물 밖으로 날아가 버립니다.안타깝게도 범위가 매우 한정되어 있습니다.
#include <stdint.h>
const uint16_t squares[] = {
0, 1, 4, 9,
16, 25, 36, 49,
64, 81, 100, 121,
144, 169, 196, 225,
256, 289, 324, 361,
400, 441, 484, 529,
576, 625, 676, 729,
784, 841, 900, 961,
1024, 1089, 1156, 1225,
1296, 1369, 1444, 1521,
1600, 1681, 1764, 1849,
1936, 2025, 2116, 2209,
2304, 2401, 2500, 2601,
2704, 2809, 2916, 3025,
3136, 3249, 3364, 3481,
3600, 3721, 3844, 3969,
4096, 4225, 4356, 4489,
4624, 4761, 4900, 5041,
5184, 5329, 5476, 5625,
5776, 5929, 6084, 6241,
6400, 6561, 6724, 6889,
7056, 7225, 7396, 7569,
7744, 7921, 8100, 8281,
8464, 8649, 8836, 9025,
9216, 9409, 9604, 9801,
10000, 10201, 10404, 10609,
10816, 11025, 11236, 11449,
11664, 11881, 12100, 12321,
12544, 12769, 12996, 13225,
13456, 13689, 13924, 14161,
14400, 14641, 14884, 15129,
15376, 15625, 15876, 16129,
16384, 16641, 16900, 17161,
17424, 17689, 17956, 18225,
18496, 18769, 19044, 19321,
19600, 19881, 20164, 20449,
20736, 21025, 21316, 21609,
21904, 22201, 22500, 22801,
23104, 23409, 23716, 24025,
24336, 24649, 24964, 25281,
25600, 25921, 26244, 26569,
26896, 27225, 27556, 27889,
28224, 28561, 28900, 29241,
29584, 29929, 30276, 30625,
30976, 31329, 31684, 32041,
32400, 32761, 33124, 33489,
33856, 34225, 34596, 34969,
35344, 35721, 36100, 36481,
36864, 37249, 37636, 38025,
38416, 38809, 39204, 39601,
40000, 40401, 40804, 41209,
41616, 42025, 42436, 42849,
43264, 43681, 44100, 44521,
44944, 45369, 45796, 46225,
46656, 47089, 47524, 47961,
48400, 48841, 49284, 49729,
50176, 50625, 51076, 51529,
51984, 52441, 52900, 53361,
53824, 54289, 54756, 55225,
55696, 56169, 56644, 57121,
57600, 58081, 58564, 59049,
59536, 60025, 60516, 61009,
61504, 62001, 62500, 63001,
63504, 64009, 64516, 65025
};
inline int isqrt(uint16_t x) {
const uint16_t *p = squares;
if (p[128] <= x) p += 128;
if (p[ 64] <= x) p += 64;
if (p[ 32] <= x) p += 32;
if (p[ 16] <= x) p += 16;
if (p[ 8] <= x) p += 8;
if (p[ 4] <= x) p += 4;
if (p[ 2] <= x) p += 2;
if (p[ 1] <= x) p += 1;
return p - squares;
}
32비트 버전은 https://gist.github.com/3481770 에서 다운로드할 수 있습니다.
정수 합니다. 즉, 정수를 구하다.floor(sqrt(x))반올림 오류 없이 정확하게.
다른 접근법에 관한 문제
- 를 사용합니다.
float★★★★★★★★★★★★★★★★★」double - @ @olp 。
isqrt한 를 낳다isqrt(100) = 15 - 거대한 룩업 테이블을 기반으로 하는 접근법은 32비트 이상으로 실용적이지 않다.
- 빠른 역 sqrt를 사용하는 것은 매우 정확하지 않습니다. 사용하는 것이 좋습니다.
sqrtf - 뉴턴의 접근법은 값비싼 정수 나눗셈과 좋은 초기 추측을 필요로 한다.
마이 어프로치
제 것은 위키피디아에서 제안된 비트 추측 접근법에 기초하고 있습니다.유감스럽게도 Wikipedia에서 제공되는 유사코드에 오류가 있어서 몇 가지 수정을 해야 했습니다.
// C++20 also provides std::bit_width in its <bit> header
unsigned char bit_width(unsigned long long x) {
return x == 0 ? 1 : 64 - __builtin_clzll(x);
}
template <typename Int, std::enable_if_t<std::is_unsigned<Int, int = 0>>
Int sqrt(const Int n) {
unsigned char shift = bit_width(n);
shift += shift & 1; // round up to next multiple of 2
Int result = 0;
do {
shift -= 2;
result <<= 1; // make space for the next guessed bit
result |= 1; // guess that the next bit is 1
result ^= result * result > (n >> shift); // revert if guess too high
} while (shift != 0);
return result;
}
bit_width할 수 됩니다.ceil(bit_width / 2)따라서 64비트 정수의 경우에도 기본 산술 연산과 비트 연산의 최소 32회 반복이 됩니다.
컴파일 출력은 약 20개의 명령어입니다.
성능
에 대해 벤치마킹했다.float -그러면 1번으로 하다..std::numeric_limits<...>::max().
- ★★★★★★에
uint32_t은 약하게 됩니다.25x보다 못하다std::sqrt(float) - ★★★★★★에
uint64_t은 약하게 됩니다.30x보다 못하다std::sqrt(double)
정확성.
이 방법은 부동소수점 수학을 사용하는 접근법과 달리 항상 완벽하게 정확합니다.
- 「」를 사용합니다.
sqrtf는28 [232, 2)] 범위의 반올림을 잘못 제공할 수 있습니다.예를들면,sqrtf(0xffffffff) = 65536이 실제로는 '''일65535.99999. - [260, 264) 범위에서는 두 배의 정밀도가 일관되게 작동하지 않습니다.예를들면,
sqrt(0x3fff...) = 2147483648이 실제로는 '''일2147483647.999999.
를 모두 정밀도64 뿐입니다long double단순히 64비트 정수 전체를 맞출 수 있기 때문입니다.
결론
말씀드렸듯이, 이 솔루션은 모든 입력을 올바르게 처리하고 정수분할을 회피하며 룩업 테이블을 필요로 하지 않는 유일한 솔루션입니다.즉, 정밀도와 무관하고 거대한 룩업 테이블이 필요하지 않은 방법이 필요한 경우 이 방법만이 유일한 옵션입니다.나 '아까', '아까', '아까', '아까'에서 유용하게 쓰일 수 것 요.constexpr성능이 중요하지 않고 100% 정확한 결과를 얻는 것이 훨씬 더 중요한 상황입니다.
뉴턴법을 이용한 대안적 접근법
뉴턴의 방법은 좋은 추측으로 시작할 때 꽤 빠를 수 있다.추측으로는 2의 다음 거듭제곱으로 반올림하여 일정한 시간에 제곱근을 계산합니다.임의의 숫자x 2에 대해 2를 사용하여x/2 제곱근을 구할 수 있습니다.
template <typename Int, std::enable_if_t<std::is_unsigned_v<Int>, int> = 0>
Int sqrt_guess(const Int n)
{
Int log2floor = bit_width(n) - 1;
// sqrt(x) is equivalent to pow(2, x / 2 = x >> 1)
// pow(2, x) is equivalent to 1 << x
return 1 << (log2floor >> 1);
}
오른쪽 시프트 중에 정밀도가 떨어졌기 때문에 이 값은 정확히x/2 2가 아닙니다.대신 2입니다floor(x/2).또, 주의해 주세요.sqrt_guess(0) = 1이는 실제로 첫 번째 반복에서 0으로 분할되는 것을 피하기 위해 필요합니다.
template <typename Int, std::enable_if_t<std::is_unsigned_v<Int>, int> = 0>
Int sqrt_newton(const Int n)
{
Int a = sqrt_guess(n);
Int b = n;
// compute unsigned difference
while (std::max(a, b) - std::min(a, b) > 1) {
b = n / a;
a = (a + b) / 2;
}
// a is now either floor(sqrt(n)) or ceil(sqrt(n))
// we decrement in the latter case
// this is overflow-safe as long as we start with a lower bound guess
return a - (a * a > n);
}
이 대체 접근방식은 첫 번째 제안과 거의 동등하지만 일반적으로 몇 % 포인트 더 빠릅니다.그러나 효율적인 하드웨어 분할에 크게 의존하고 있으며 결과는 크게 다를 수 있습니다.
「 」의 sqrt_guess을 사용하다 .1첫 번째 추측으로
이것은 매우 짧기 때문에 99% 인라인입니다.
static inline int sqrtn(int num) {
int i = 0;
__asm__ (
"pxor %%xmm0, %%xmm0\n\t" // clean xmm0 for cvtsi2ss
"cvtsi2ss %1, %%xmm0\n\t" // convert num to float, put it to xmm0
"sqrtss %%xmm0, %%xmm0\n\t" // square root xmm0
"cvttss2si %%xmm0, %0" // float to int
:"=r"(i):"r"(num):"%xmm0"); // i: result, num: input, xmm0: scratch register
return i;
}
이유xmm0?의 문서
수신처 오퍼랜드는 XMM 레지스터입니다.결과는 행선지 오퍼랜드의 낮은 더블워드에 저장되며, 상위 3개의 더블워드는 변경되지 않습니다.
GCC Intelligent 버전(GCC에서만 실행):
#include <xmmintrin.h>
int sqrtn2(int num) {
register __v4sf xmm0 = {0, 0, 0, 0};
xmm0 = __builtin_ia32_cvtsi2ss(xmm0, num);
xmm0 = __builtin_ia32_sqrtss(xmm0);
return __builtin_ia32_cvttss2si(xmm0);
}
Intel Entrient 버전 (GCC, Clang, ICC에서 테스트 완료):
#include <xmmintrin.h>
int sqrtn2(int num) {
register __m128 xmm0 = _mm_setzero_ps();
xmm0 = _mm_cvt_si2ss(xmm0, num);
xmm0 = _mm_sqrt_ss(xmm0);
return _mm_cvtt_ss2si(xmm0);
}
^^^^ 모두 SSE 1(SSE 2도 아님)이 필요합니다.
주의: GCC의 계산 방법은 다음과 같습니다.(int) sqrt((float) num)-Ofast큰 사이즈로 높은 정밀도를 원하는 경우i 계산하면 (int) sqrt((double) num)(Gumby The Green ( ( ( ( ( ( ( ( ( ( ( ( ( ()
static inline int sqrtn(int num) {
int i = 0;
__asm__ (
"pxor %%xmm0, %%xmm0\n\t"
"cvtsi2sd %1, %%xmm0\n\t"
"sqrtsd %%xmm0, %%xmm0\n\t"
"cvttsd2si %%xmm0, %0"
:"=r"(i):"r"(num):"%xmm0");
return i;
}
또는
#include <xmmintrin.h>
int sqrtn2(int num) {
register __v2df xmm0 = {0, 0};
xmm0 = __builtin_ia32_cvtsi2sd(xmm0, num);
xmm0 = __builtin_ia32_sqrtsd(xmm0);
return __builtin_ia32_cvttsd2si(xmm0);
}
또는 바이너리 검색을 실행하기만 하면 더 간단한 버전 imo를 쓸 수 없습니다.
uint16_t sqrti(uint32_t num)
{
uint16_t ret = 0;
for(int32_t i = 15; i >= 0; i--)
{
uint16_t temp = ret | (1 << i);
if(temp * temp <= num)
{
ret = temp;
}
}
return ret;
}
정수 sqrt를 수행하려면 다음과 같은 newtons 방법을 사용합니다.
Def isqrt(N):
a = 1
b = N
while |a-b| > 1
b = N / a
a = (a + b) / 2
return a
기본적으로 모든 x에 대해 sqrt는 범위(x ... N/x) 내에 있으므로 새로운 추측을 위해 모든 루프에서 해당 간격을 이등분합니다.바이너리 검색과 비슷하지만 컨버전스 속도가 빨라야 합니다.
이것은 매우 빠른 O(loglog(N)로 수렴됩니다.부동소수도 전혀 사용하지 않으며 임의의 정밀도 정수에서도 잘 작동합니다.
근사치라도 괜찮으시다면 제가 만든 이 정수 제곱함수는 어떻습니까?
int sqrti(int x)
{
union { float f; int x; } v;
// convert to float
v.f = (float)x;
// fast aprox sqrt
// assumes float is in IEEE 754 single precision format
// assumes int is 32 bits
// b = exponent bias
// m = number of mantissa bits
v.x -= 1 << 23; // subtract 2^m
v.x >>= 1; // divide by 2
v.x += 1 << 29; // add ((b + 1) / 2) * 2^m
// convert to int
return (int)v.f;
}
이 Wikipedia 문서에 설명된 알고리즘을 사용합니다.제 기계에서는 sqrt보다 거의 2배 빠릅니다. : )
gcc, -fast-math를 사용하는 컴퓨터에서 32비트 정수를 부동으로 변환하고 sqrtf를 사용하는 데 10^9 ops당 1.2초 걸립니다(-fast-math를 사용하지 않는 경우 3.54초 소요).
다음 알고리즘에서는 10^9당 0.87초를 사용하지만 RMS 오류는 0.79에 불과하지만 오류는 -7 또는 +1이 될 수 있습니다.
uint16_t SQRTTAB[65536];
inline uint16_t approxsqrt(uint32_t x) {
const uint32_t m1 = 0xff000000;
const uint32_t m2 = 0x00ff0000;
if (x&m1) {
return SQRTTAB[x>>16];
} else if (x&m2) {
return SQRTTAB[x>>8]>>4;
} else {
return SQRTTAB[x]>>8;
}
}
테이블은 다음을 사용하여 구성됩니다.
void maketable() {
for (int x=0; x<65536; x++) {
double v = x/65535.0;
v = sqrt(v);
int y = int(v*65535.0+0.999);
SQRTTAB[x] = y;
}
}
문장이 정확성을 높인다면 이등분법을 더 잘 할 수 있다는 것을 알게 되었지만, 적어도 -fast-math를 사용하면 sqrtf가 더 빠르다는 점까지 속도가 느려집니다.
제곱근에 대한 성능이 필요하다면 많은 부분을 계산해야 할 것입니다.그럼 왜 답을 캐싱하지 않는 거죠?이 경우 N의 범위는 알 수 없습니다.또한 같은 정수의 제곱근의 몇 배나 계산해도 알 수 없습니다.그렇다면 메서드가 호출될 때마다 결과를 캐시할 수 있습니다(배열이 너무 크지는 않으면 가장 효율적입니다).
「빠른 정수 제곱근」을 검색하면, 많은 옵션을 찾을 수 있을 것 같습니다만, 여기, 잘 될 가능성이 있는 새로운 아이디어가 몇개인가 있습니다(각각의 독립적, 또는 그것들을 조합할 수도 있습니다).
static const지원하는 도메인 내의 모든 완벽한 정사각형 배열과 빠른 분기 없는 이진 검색을 수행합니다.이치노- 숫자를 부동 소수점으로 변환하고 이를 가수와 지수로 나눕니다.지수를 반감하고 가수에 마법의 계수를 곱합니다(찾는 작업).이를 통해 매우 가까운 근사치를 얻을 수 있습니다.정확하지 않은 경우 조정하기 위한 마지막 단계를 포함합니다(또는 위의 이진 검색의 시작점으로 사용).
왜 아무도 가장 빠른 방법을 제안하지 않는 거죠?
다음 경우:
- 숫자의 범위가 한정되어 있다
- 메모리 소비는 중요하지 않다
- 응용 프로그램 시작 시간은 중요하지 않습니다.
을 만듭니다.int[MAX_X]출시 로 채워져 .sqrt(x)sqrt()★★★★★★★★★★★★★★★★★★」
이 모든 조건들이 내 프로그램에 꽤 잘 들어맞습니다. ,,는,int[10000000]됩니다.40MB.
이것에 대해 어떻게 생각하세요?
대부분의 경우 정확한 정수 sqrt 값조차 필요하지 않고 충분한 근사치를 얻을 수 있습니다.(예를 들어 DSP 최적화에서는 32비트 신호가 16비트 또는 16비트~8비트로 압축되어야 하며 0에 가까운 정밀도는 거의 손실되지 않습니다).
유용한 방정식을 찾았습니다.
k = ceil(MSB(n)/2); - MSB(n) is the most significant bit of "n"
sqrt(n) ~= 2^(k-2)+(2^(k-1))*n/(2^(2*k))); - all multiplications and divisions here are very DSP-friendly, as they are only 2^k.
이 방정식은 매끄러운 곡선(n, sqrt(n))을 생성하며, 그 값은 실제 sqrt(n)와 크게 다르지 않으므로 대략적인 정확도가 충분할 때 유용합니다.
언급URL : https://stackoverflow.com/questions/4930307/fastest-way-to-get-the-integer-part-of-sqrtn
'programing' 카테고리의 다른 글
| 액션을 호출하는 방법(NuxtJs (0) | 2022.08.27 |
|---|---|
| Java에서 Map 값을 키별로 정렬하려면 어떻게 해야 합니까? (0) | 2022.08.27 |
| Vue 변경 폭 및 내용 (0) | 2022.08.27 |
| 시리얼라이제이션이란? (0) | 2022.08.27 |
| Java Servlet이란? (0) | 2022.08.27 |