const 비 정수 지수로 pow () 최적화?
내 코드에 pow()실행 시간의 약 10-20 %를 차지하는 핫스팟이 있습니다.
내 입력 pow(x,y)은 매우 구체적이므로 pow()더 높은 성능으로 두 개의 근사값 (각 지수에 대해 하나씩) 을 롤링하는 방법이 있는지 궁금합니다 .
- 두 개의 상수 지수 : 2.4와 1 / 2.4가 있습니다.
- 지수가 2.4이면 x 는 (0.090473935, 1.0] 범위에 있습니다.
- 지수가 1 / 2.4이면 x 는 (0.0031308, 1.0] 범위에 있습니다.
- SSE / AVX
float벡터를 사용하고 있습니다. 플랫폼 세부 사항을 활용할 수 있다면 당장!
0.01 % 정도의 최대 오류율이 이상적이지만 전체 정밀도 (for float) 알고리즘에도 관심이 있습니다.
나는 이미 빠른 pow() 근사값을 사용하고 있지만 이러한 제약을 고려하지 않습니다. 더 잘할 수 있습니까?
IEEE 754 해킹 분야에서 더 빠르고 덜 "마 법적"인 또 다른 솔루션이 있습니다. 약 12 개의 클럭 사이클에서 0.08 %의 오류 마진을 달성합니다 (Intel Merom CPU에서 p = 2.4의 경우).
부동 소수점 숫자는 원래 로그에 대한 근사값으로 발명되었으므로 정수 값을의 근사값으로 사용할 수 있습니다 log2. 이것은 정수에서 변환 명령어를 부동 소수점 값에 적용하여 다른 부동 소수점 값을 얻음으로써 어느 정도 이식 가능하게 달성 할 수 있습니다.
pow계산 을 완료하기 위해 상수 인수를 곱하고 정수로 변환 명령어를 사용하여 로그를 다시 변환 할 수 있습니다. SSE에서 관련 지침은 cvtdq2ps및 cvtps2dq입니다.
하지만 그렇게 간단하지는 않습니다. IEEE 754의 지수 필드는 부호가 있으며 바이어스 값 127은 지수 0을 나타냅니다. 이 편향은 로그를 곱하기 전에 제거하고 지수화하기 전에 다시 추가해야합니다. 또한 빼기에 의한 편향 조정은 0에서 작동하지 않습니다. 다행히도 두 가지 조정은 미리 상수 계수를 곱하여 수행 할 수 있습니다.
x^p
= exp2( p * log2( x ) )
= exp2( p * ( log2( x ) + 127 - 127 ) - 127 + 127 )
= cvtps2dq( p * ( log2( x ) + 127 - 127 - 127 / p ) )
= cvtps2dq( p * ( log2( x ) + 127 - log2( exp2( 127 - 127 / p ) ) )
= cvtps2dq( p * ( log2( x * exp2( 127 / p - 127 ) ) + 127 ) )
= cvtps2dq( p * ( cvtdq2ps( x * exp2( 127 / p - 127 ) ) ) )
exp2( 127 / p - 127 )상수 요소입니다. 이 함수는 다소 전문화되어 있습니다. 상수 인자가 지수의 역수로 기하 급수적으로 증가하고 오버플로되기 때문에 작은 분수 지수에서는 작동하지 않습니다. 음의 지수에서는 작동하지 않습니다. 큰 지수는 가수 비트가 곱셈에 의해 지수 비트와 혼합되기 때문에 높은 오류로 이어집니다.
그러나 4 개의 빠른 명령에 불과합니다. 사전 곱하기, "정수"에서 (로그로) 변환, 거듭 제곱, "정수"(로그에서)로 변환. 이 SSE 구현에서는 변환이 매우 빠릅니다. 추가 상수 계수를 첫 번째 곱셈으로 짜낼 수도 있습니다.
template< unsigned expnum, unsigned expden, unsigned coeffnum, unsigned coeffden >
__m128 fastpow( __m128 arg ) {
__m128 ret = arg;
// std::printf( "arg = %,vg\n", ret );
// Apply a constant pre-correction factor.
ret = _mm_mul_ps( ret, _mm_set1_ps( exp2( 127. * expden / expnum - 127. )
* pow( 1. * coeffnum / coeffden, 1. * expden / expnum ) ) );
// std::printf( "scaled = %,vg\n", ret );
// Reinterpret arg as integer to obtain logarithm.
asm ( "cvtdq2ps %1, %0" : "=x" (ret) : "x" (ret) );
// std::printf( "log = %,vg\n", ret );
// Multiply logarithm by power.
ret = _mm_mul_ps( ret, _mm_set1_ps( 1. * expnum / expden ) );
// std::printf( "powered = %,vg\n", ret );
// Convert back to "integer" to exponentiate.
asm ( "cvtps2dq %1, %0" : "=x" (ret) : "x" (ret) );
// std::printf( "result = %,vg\n", ret );
return ret;
}
지수가 2.4 인 몇 번의 시행에서는 이것이 약 5 %까지 지속적으로 과대 평가됨을 보여줍니다. (루틴은 항상 과대 평가하는 것이 보장됩니다.) 간단히 0.95를 곱할 수 있지만 몇 가지 명령을 더 사용하면 그래픽에 충분할 정도로 10 진수 4 자리 정도의 정확도를 얻을 수 있습니다.
핵심은 과대 평가와 과소 평가를 일치시키고 평균을 취하는 것입니다.
- x ^ 0.8 계산 : 4 개의 명령어, 오류 ~ + 3 %.
- x ^ -0.4 계산 : 1
rsqrtps. (이것은 충분히 정확하지만 0으로 작업하는 능력을 희생합니다.) - x ^ 0.4 계산 : 1
mulps. - x ^ -0.2 계산 : 1
rsqrtps. - x ^ 2 계산 : 1
mulps. - x ^ 3 계산 : 1
mulps. - x ^ 2.4 = x ^ 2 * x ^ 0.4 : 1
mulps. 이것은 과대 평가입니다. - x ^ 2.4 = x ^ 3 * x ^ -0.4 * x ^ -0.2 : 두
mulps. 이것은 과소 평가입니다. - 위의 평균값 : 하나
addps, 하나mulps.
명령어 집계 : 지연 시간 = 5 인 전환 2 개 및 처리량 = 4 인 2 개의 역 제곱근 추정치를 포함하여 14 개.
평균을 제대로 취하기 위해 예상 오차로 추정값에 가중치를 부여하려고합니다. 과소 추정은 오류를 0.6 대 0.4의 거듭 제곱으로 올립니다. 따라서 오류가 1.5 배가 될 것으로 예상합니다. 가중치는 지침을 추가하지 않습니다. 사전 요소에서 수행 할 수 있습니다. 계수 a를 호출하면 a ^ 0.5 = 1.5 a ^ -0.75, a = 1.38316186입니다.
최종 오차는 약 .015 %로 초기 fastpow결과 보다 2 배 더 좋습니다 . 런타임은 volatile소스 및 대상 변수가 있는 바쁜 루프에 대해 약 12주기입니다. 반복이 겹치는 경우에도 실제 사용에는 명령어 수준 병렬 처리가 표시됩니다. SIMD를 고려하면 3주기 당 스칼라 결과 하나의 처리량입니다!
int main() {
__m128 const x0 = _mm_set_ps( 0.01, 1, 5, 1234.567 );
std::printf( "Input: %,vg\n", x0 );
// Approx 5% accuracy from one call. Always an overestimate.
__m128 x1 = fastpow< 24, 10, 1, 1 >( x0 );
std::printf( "Direct x^2.4: %,vg\n", x1 );
// Lower exponents provide lower initial error, but too low causes overflow.
__m128 xf = fastpow< 8, 10, int( 1.38316186 * 1e9 ), int( 1e9 ) >( x0 );
std::printf( "1.38 x^0.8: %,vg\n", xf );
// Imprecise 4-cycle sqrt is still far better than fastpow, good enough.
__m128 xfm4 = _mm_rsqrt_ps( xf );
__m128 xf4 = _mm_mul_ps( xf, xfm4 );
// Precisely calculate x^2 and x^3
__m128 x2 = _mm_mul_ps( x0, x0 );
__m128 x3 = _mm_mul_ps( x2, x0 );
// Overestimate of x^2 * x^0.4
x2 = _mm_mul_ps( x2, xf4 );
// Get x^-0.2 from x^0.4. Combine with x^-0.4 into x^-0.6 and x^2.4.
__m128 xfm2 = _mm_rsqrt_ps( xf4 );
x3 = _mm_mul_ps( x3, xfm4 );
x3 = _mm_mul_ps( x3, xfm2 );
std::printf( "x^2 * x^0.4: %,vg\n", x2 );
std::printf( "x^3 / x^0.6: %,vg\n", x3 );
x2 = _mm_mul_ps( _mm_add_ps( x2, x3 ), _mm_set1_ps( 1/ 1.960131704207789 ) );
// Final accuracy about 0.015%, 200x better than x^0.8 calculation.
std::printf( "average = %,vg\n", x2 );
}
글쎄요 ...이 글을 더 빨리 게시하지 못해서 죄송합니다. 그리고 x ^ 1 / 2.4로 확장하는 것은 연습으로 남겨집니다; v).
통계로 업데이트
좀 테스트 하니스와 투 엑스 구현 ( 5 / 12 ) 위의 대응 케이스.
#include <cstdio>
#include <xmmintrin.h>
#include <cmath>
#include <cfloat>
#include <algorithm>
using namespace std;
template< unsigned expnum, unsigned expden, unsigned coeffnum, unsigned coeffden >
__m128 fastpow( __m128 arg ) {
__m128 ret = arg;
// std::printf( "arg = %,vg\n", ret );
// Apply a constant pre-correction factor.
ret = _mm_mul_ps( ret, _mm_set1_ps( exp2( 127. * expden / expnum - 127. )
* pow( 1. * coeffnum / coeffden, 1. * expden / expnum ) ) );
// std::printf( "scaled = %,vg\n", ret );
// Reinterpret arg as integer to obtain logarithm.
asm ( "cvtdq2ps %1, %0" : "=x" (ret) : "x" (ret) );
// std::printf( "log = %,vg\n", ret );
// Multiply logarithm by power.
ret = _mm_mul_ps( ret, _mm_set1_ps( 1. * expnum / expden ) );
// std::printf( "powered = %,vg\n", ret );
// Convert back to "integer" to exponentiate.
asm ( "cvtps2dq %1, %0" : "=x" (ret) : "x" (ret) );
// std::printf( "result = %,vg\n", ret );
return ret;
}
__m128 pow125_4( __m128 arg ) {
// Lower exponents provide lower initial error, but too low causes overflow.
__m128 xf = fastpow< 4, 5, int( 1.38316186 * 1e9 ), int( 1e9 ) >( arg );
// Imprecise 4-cycle sqrt is still far better than fastpow, good enough.
__m128 xfm4 = _mm_rsqrt_ps( xf );
__m128 xf4 = _mm_mul_ps( xf, xfm4 );
// Precisely calculate x^2 and x^3
__m128 x2 = _mm_mul_ps( arg, arg );
__m128 x3 = _mm_mul_ps( x2, arg );
// Overestimate of x^2 * x^0.4
x2 = _mm_mul_ps( x2, xf4 );
// Get x^-0.2 from x^0.4, and square it for x^-0.4. Combine into x^-0.6.
__m128 xfm2 = _mm_rsqrt_ps( xf4 );
x3 = _mm_mul_ps( x3, xfm4 );
x3 = _mm_mul_ps( x3, xfm2 );
return _mm_mul_ps( _mm_add_ps( x2, x3 ), _mm_set1_ps( 1/ 1.960131704207789 * 0.9999 ) );
}
__m128 pow512_2( __m128 arg ) {
// 5/12 is too small, so compute the sqrt of 10/12 instead.
__m128 x = fastpow< 5, 6, int( 0.992245 * 1e9 ), int( 1e9 ) >( arg );
return _mm_mul_ps( _mm_rsqrt_ps( x ), x );
}
__m128 pow512_4( __m128 arg ) {
// 5/12 is too small, so compute the 4th root of 20/12 instead.
// 20/12 = 5/3 = 1 + 2/3 = 2 - 1/3. 2/3 is a suitable argument for fastpow.
// weighting coefficient: a^-1/2 = 2 a; a = 2^-2/3
__m128 xf = fastpow< 2, 3, int( 0.629960524947437 * 1e9 ), int( 1e9 ) >( arg );
__m128 xover = _mm_mul_ps( arg, xf );
__m128 xfm1 = _mm_rsqrt_ps( xf );
__m128 x2 = _mm_mul_ps( arg, arg );
__m128 xunder = _mm_mul_ps( x2, xfm1 );
// sqrt2 * over + 2 * sqrt2 * under
__m128 xavg = _mm_mul_ps( _mm_set1_ps( 1/( 3 * 0.629960524947437 ) * 0.999852 ),
_mm_add_ps( xover, xunder ) );
xavg = _mm_mul_ps( xavg, _mm_rsqrt_ps( xavg ) );
xavg = _mm_mul_ps( xavg, _mm_rsqrt_ps( xavg ) );
return xavg;
}
__m128 mm_succ_ps( __m128 arg ) {
return (__m128) _mm_add_epi32( (__m128i) arg, _mm_set1_epi32( 4 ) );
}
void test_pow( double p, __m128 (*f)( __m128 ) ) {
__m128 arg;
for ( arg = _mm_set1_ps( FLT_MIN / FLT_EPSILON );
! isfinite( _mm_cvtss_f32( f( arg ) ) );
arg = mm_succ_ps( arg ) ) ;
for ( ; _mm_cvtss_f32( f( arg ) ) == 0;
arg = mm_succ_ps( arg ) ) ;
std::printf( "Domain from %g\n", _mm_cvtss_f32( arg ) );
int n;
int const bucket_size = 1 << 25;
do {
float max_error = 0;
double total_error = 0, cum_error = 0;
for ( n = 0; n != bucket_size; ++ n ) {
float result = _mm_cvtss_f32( f( arg ) );
if ( ! isfinite( result ) ) break;
float actual = ::powf( _mm_cvtss_f32( arg ), p );
float error = ( result - actual ) / actual;
cum_error += error;
error = std::abs( error );
max_error = std::max( max_error, error );
total_error += error;
arg = mm_succ_ps( arg );
}
std::printf( "error max = %8g\t" "avg = %8g\t" "|avg| = %8g\t" "to %8g\n",
max_error, cum_error / n, total_error / n, _mm_cvtss_f32( arg ) );
} while ( n == bucket_size );
}
int main() {
std::printf( "4 insn x^12/5:\n" );
test_pow( 12./5, & fastpow< 12, 5, 1059, 1000 > );
std::printf( "14 insn x^12/5:\n" );
test_pow( 12./5, & pow125_4 );
std::printf( "6 insn x^5/12:\n" );
test_pow( 5./12, & pow512_2 );
std::printf( "14 insn x^5/12:\n" );
test_pow( 5./12, & pow512_4 );
}
산출:
4 insn x^12/5:
Domain from 1.36909e-23
error max = inf avg = inf |avg| = inf to 8.97249e-19
error max = 2267.14 avg = 139.175 |avg| = 139.193 to 5.88021e-14
error max = 0.123606 avg = -0.000102963 |avg| = 0.0371122 to 3.85365e-09
error max = 0.123607 avg = -0.000108978 |avg| = 0.0368548 to 0.000252553
error max = 0.12361 avg = 7.28909e-05 |avg| = 0.037507 to 16.5513
error max = 0.123612 avg = -0.000258619 |avg| = 0.0365618 to 1.08471e+06
error max = 0.123611 avg = 8.70966e-05 |avg| = 0.0374369 to 7.10874e+10
error max = 0.12361 avg = -0.000103047 |avg| = 0.0371122 to 4.65878e+15
error max = 0.123609 avg = nan |avg| = nan to 1.16469e+16
14 insn x^12/5:
Domain from 1.42795e-19
error max = inf avg = nan |avg| = nan to 9.35823e-15
error max = 0.000936462 avg = 2.0202e-05 |avg| = 0.000133764 to 6.13301e-10
error max = 0.000792752 avg = 1.45717e-05 |avg| = 0.000129936 to 4.01933e-05
error max = 0.000791785 avg = 7.0132e-06 |avg| = 0.000129923 to 2.63411
error max = 0.000787589 avg = 1.20745e-05 |avg| = 0.000129347 to 172629
error max = 0.000786553 avg = 1.62351e-05 |avg| = 0.000132397 to 1.13134e+10
error max = 0.000785586 avg = 8.25205e-06 |avg| = 0.00013037 to 6.98147e+12
6 insn x^5/12:
Domain from 9.86076e-32
error max = 0.0284339 avg = 0.000441158 |avg| = 0.00967327 to 6.46235e-27
error max = 0.0284342 avg = -5.79938e-06 |avg| = 0.00897913 to 4.23516e-22
error max = 0.0284341 avg = -0.000140706 |avg| = 0.00897084 to 2.77556e-17
error max = 0.028434 avg = 0.000440504 |avg| = 0.00967325 to 1.81899e-12
error max = 0.0284339 avg = -6.11153e-06 |avg| = 0.00897915 to 1.19209e-07
error max = 0.0284298 avg = -0.000140597 |avg| = 0.00897084 to 0.0078125
error max = 0.0284371 avg = 0.000439748 |avg| = 0.00967319 to 512
error max = 0.028437 avg = -7.74294e-06 |avg| = 0.00897924 to 3.35544e+07
error max = 0.0284369 avg = -0.000142036 |avg| = 0.00897089 to 2.19902e+12
error max = 0.0284368 avg = 0.000439183 |avg| = 0.0096732 to 1.44115e+17
error max = 0.0284367 avg = -7.41244e-06 |avg| = 0.00897923 to 9.44473e+21
error max = 0.0284366 avg = -0.000141706 |avg| = 0.00897088 to 6.1897e+26
error max = 0.485129 avg = -0.0401671 |avg| = 0.048422 to 4.05648e+31
error max = 0.994932 avg = -0.891494 |avg| = 0.891494 to 2.65846e+36
error max = 0.999329 avg = nan |avg| = nan to -0
14 insn x^5/12:
Domain from 2.64698e-23
error max = 0.13556 avg = 0.00125936 |avg| = 0.00354677 to 1.73472e-18
error max = 0.000564988 avg = 2.51458e-06 |avg| = 0.000113709 to 1.13687e-13
error max = 0.000565065 avg = -1.49258e-06 |avg| = 0.000112553 to 7.45058e-09
error max = 0.000565143 avg = 1.5293e-06 |avg| = 0.000112864 to 0.000488281
error max = 0.000565298 avg = 2.76457e-06 |avg| = 0.000113713 to 32
error max = 0.000565453 avg = -1.61276e-06 |avg| = 0.000112561 to 2.09715e+06
error max = 0.000565531 avg = 1.42628e-06 |avg| = 0.000112866 to 1.37439e+11
error max = 0.000565686 avg = 2.71505e-06 |avg| = 0.000113715 to 9.0072e+15
error max = 0.000565763 avg = -1.56586e-06 |avg| = 0.000112415 to 1.84467e+19
나는 더 정확한 5/12의 정확성이 rsqrt작업에 의해 제한되고 있다고 생각 합니다.
이것은 이전 답변과 매우 다르기 때문에 또 다른 답변입니다. 상대 오차는 3e-8입니다. 더 많은 정확성을 원하십니까? Chebychev 용어를 몇 개 더 추가하십시오. 2 ^ n-epsilon과 2 ^ n + epsilon 사이의 작은 불연속성을 만들기 때문에 순서를 홀수로 유지하는 것이 가장 좋습니다.
#include <stdlib.h>
#include <math.h>
// Returns x^(5/12) for x in [1,2), to within 3e-8 (relative error).
// Want more precision? Add more Chebychev polynomial coefs.
double pow512norm (
double x)
{
static const int N = 8;
// Chebychev polynomial terms.
// Non-zero terms calculated via
// integrate (2/pi)*ChebyshevT[n,u]/sqrt(1-u^2)*((u+3)/2)^(5/12)
// from -1 to 1
// Zeroth term is similar except it uses 1/pi rather than 2/pi.
static const double Cn[N] = {
1.1758200232996901923,
0.16665763094889061230,
-0.0083154894939042125035,
0.00075187976780420279038,
// Wolfram alpha doesn't want to compute the remaining terms
// to more precision (it times out).
-0.0000832402,
0.0000102292,
-1.3401e-6,
1.83334e-7};
double Tn[N];
double u = 2.0*x - 3.0;
Tn[0] = 1.0;
Tn[1] = u;
for (int ii = 2; ii < N; ++ii) {
Tn[ii] = 2*u*Tn[ii-1] - Tn[ii-2];
}
double y = 0.0;
for (int ii = N-1; ii >= 0; --ii) {
y += Cn[ii]*Tn[ii];
}
return y;
}
// Returns x^(5/12) to within 3e-8 (relative error).
double pow512 (
double x)
{
static const double pow2_512[12] = {
1.0,
pow(2.0, 5.0/12.0),
pow(4.0, 5.0/12.0),
pow(8.0, 5.0/12.0),
pow(16.0, 5.0/12.0),
pow(32.0, 5.0/12.0),
pow(64.0, 5.0/12.0),
pow(128.0, 5.0/12.0),
pow(256.0, 5.0/12.0),
pow(512.0, 5.0/12.0),
pow(1024.0, 5.0/12.0),
pow(2048.0, 5.0/12.0)
};
double s;
int iexp;
s = frexp (x, &iexp);
s *= 2.0;
iexp -= 1;
div_t qr = div (iexp, 12);
if (qr.rem < 0) {
qr.quot -= 1;
qr.rem += 12;
}
return ldexp (pow512norm(s)*pow2_512[qr.rem], 5*qr.quot);
}
부록 : 여기서 무슨 일이 일어나고 있습니까?
요청에 따라 다음은 위 코드의 작동 방식을 설명합니다.
개요
두 개의 함수를 정의하는 위의 코드를, double pow512norm (double x)하고 double pow512 (double x). 후자는 스위트의 진입 점입니다. 이것은 사용자 코드가 x ^ (5/12)를 계산하기 위해 호출해야하는 함수입니다. 이 함수 pow512norm(x)는 체비 쇼프 다항식을 사용하여 x ^ (5/12)를 근사하지만 [1,2] 범위의 x에만 해당됩니다. ( pow512norm(x)해당 범위를 벗어난 x 값에 사용 하면 결과는 가비지가됩니다.)
이 함수는 pow512(x)입력 분할 x쌍으로 (double s, int n)되도록 x = s * 2^n하고되도록 1≤ s<2. 의 다른 구획 n으로 (int q, unsigned int r)되도록 n = 12*q + r하고 r(12) 내가 찾는 문제를 분할 수보다 적은 X ^ (5/12) 부분으로 :
x^(5/12)=(s^(5/12))*((2^n)^(5/12))(u v) ^ a = (u ^ a) (v ^ a)를 통해 양의 u, v 및 실수 a.s^(5/12)를 통해 계산됩니다pow512norm(s).(2^n)^(5/12)=(2^(12*q+r))^(5/12)대체를 통해.2^(12*q+r)=(2^(12*q))*(2^r)비아u^(a+b)=(u^a)*(u^b)긍정적 U 실제 A, B에 대한.(2^(12*q+r))^(5/12)=(2^(5*q))*((2^r)^(5/12))더 많은 조작을 통해.(2^r)^(5/12)조회 테이블에 의해 계산됩니다pow2_512.- 계산하면
pow512norm(s)*pow2_512[qr.rem]거의 다되었습니다. 위의 3 단계에서 계산 된 값qr.rem은 다음과 같습니다r. 필요한 것은2^(5*q)원하는 결과를 얻기 위해 이것을 곱하는 것입니다. - 이것이 바로 수학 라이브러리 함수
ldexp가하는 일입니다.
함수 근사
여기서의 목표는 당면한 문제에 대해 '충분히 좋은'f (x) = x ^ (5/12)의 쉽게 계산할 수있는 근사치를 만드는 것입니다. 우리의 근사는 어떤 의미에서 f (x)에 가까워 야합니다. 수사적 질문 : '가까운'이란 무엇을 의미합니까? 두 가지 경쟁 해석은 평균 제곱 오차를 최소화하는 것과 최대 절대 오차를 최소화하는 것입니다.
나는 주식 시장 비유를 사용하여 이들 사이의 차이점을 설명 할 것입니다. 최종 은퇴를 위해 저축을 원한다고 가정하십시오. 20 대라면 가장 좋은 방법은 주식이나 주식 시장 펀드에 투자하는 것입니다. 이는 충분한 시간 동안 주식 시장이 평균적으로 다른 투자 계획을 능가하기 때문입니다. 그러나 우리는 주식에 돈을 투자하는 것이 매우 나쁜 일을 보았습니다. 50 대 또는 60 대 (젊은 은퇴를 원하는 경우 40 대)라면 좀 더 보수적으로 투자해야합니다. 이러한 다운 스윙은 은퇴 포트폴리오에 타격을 줄 수 있습니다.
함수 근사화로 돌아 가기 : 근사치의 소비자로서 일반적으로 "평균"성능보다는 최악의 오류에 대해 걱정합니다. "평균적으로"(예 : 최소 제곱) 최상의 성능을 제공하기 위해 구성된 근사값을 사용하면 Murphy의 법칙에 따라 성능이 평균보다 훨씬 더 나쁜 위치에 정확히 근사값을 사용하여 프로그램이 많은 시간을 소비하게됩니다. 당신이 원하는 것은 일부 영역에서 최대 절대 오차를 최소화하는 최소 근사치입니다. 좋은 수학 라이브러리는 최소 제곱 접근 방식이 아니라 최소 제곱 접근 방식을 사용합니다. 이렇게하면 수학 라이브러리 작성자가 라이브러리의 성능을 보장 할 수 있기 때문입니다.
수학 라이브러리는 일반적으로 다항식 또는 유리 다항식을 사용하여 일부 도메인 a≤x≤b에 대해 일부 함수 f (x)를 근사합니다. 함수 f (x)가이 영역에 대해 분석적이며 함수를 N 차수의 다항식 p (x)로 근사하려고한다고 가정합니다. 주어진 차수 N에 대해 마법적이고 고유 한 다항식 p (x)가 존재하므로 p ( x) -f (x)는 [a, b]에 대해 N + 2 극값을 가지며 이러한 N + 2 극값의 절대 값은 모두 서로 동일합니다. 이 마법의 다항식 p (x)를 찾는 것은 함수 근사 자의 성배입니다.
나는 당신을 위해 성배를 찾지 못했습니다. 대신 Chebyshev 근사를 사용했습니다. 첫 번째 종류의 체비 쇼프 다항식은 함수 근사화와 관련하여 매우 멋진 기능이있는 직교 (직교 정규가 아님) 다항식 집합입니다. Chebyshev 근사는 종종 마법 다항식 p (x)에 매우 가깝습니다. (사실 성배 다항식이 일반적으로 체비 쇼프 근사값으로 시작한다는 것을 발견하는 Remez 교환 알고리즘입니다.)
pow512norm (x)
이 함수는 Chebyshev 근사를 사용하여 x ^ (5/12)에 근사하는 다항식 p * (x)를 찾습니다. 여기서 나는 위에서 설명한 마법 다항식 p (x)와이 Chebyshev 근사를 구별하기 위해 p * (x)를 사용하고 있습니다. 체비 쇼프 근사 p * (x)는 쉽게 찾을 수 있습니다. p (x)를 찾는 것은 곰입니다. 체비 쇼프 근사 p * (x)는 sum_i Cn [i] * Tn (i, x)이며, 여기서 Cn [i]는 체비 쇼프 계수이고 Tn (i, x)는 x에서 평가 된 체비 쇼프 다항식입니다.
저는 Wolfram 알파를 사용하여 체비 쇼프 계수를 찾았 Cn습니다. 예를 들어, 이것은Cn[1] . 입력 상자 다음의 첫 번째 상자에는 원하는 답이 있습니다 (이 경우 0.166658). 내가 원하는만큼의 숫자는 아닙니다. '더 많은 숫자'를 클릭하면 훨씬 더 많은 숫자를 얻을 수 있습니다. Wolfram 알파는 무료입니다. 얼마나 많은 계산을 할 것인지에 대한 제한이 있습니다. 그것은 고차 조건에 대한 제한에 도달합니다. (만약 당신이 mathematica를 구입하거나 접근 할 수 있다면, 높은 정밀도로 이러한 고차 계수를 계산할 수있을 것입니다.)
체비 쇼프 다항식 Tn (x)은 배열에서 계산됩니다 Tn. 시작을 : 매우 근접 마법 다항식 P (x)를 부여하는 것을 넘어서, 체비 셰프 근사를 사용하는 또 다른 이유는 그들 체비 셰프 다항식의 값이 용이하게 계산된다는 것이다 Tn[0]=1하고 Tn[1]=x다음 반복 계산 Tn[i]=2*x*Tn[i-1] - Tn[i-2]. (코드에서 'i'가 아닌 'ii'를 색인 변수로 사용했습니다. 변수 이름으로 'i'를 사용하지 않습니다. 영어에서 단어에 'i'가있는 단어는 몇 개입니까? 연속 'i 's?)
pow512 (x)
pow512 는 사용자 코드가 호출해야하는 함수입니다. 위에서 이미이 함수의 기본 사항을 설명했습니다. 몇 가지 추가 정보 : 수학 라이브러리 함수 frexp(x)는 입력에 대한 유효 s및 지수 iexp를 반환합니다 x. (사소한 문제 : s1에서 2 사이의 값 을 사용하고 pow512norm싶지만 frexp0.5에서 1 사이의 값을 반환합니다.) 수학 라이브러리 함수 div는 한 번의 스웰 퍽에서 정수 나누기에 대한 몫과 나머지를 반환합니다. 마지막으로 수학 라이브러리 함수 ldexp를 사용하여 세 부분을 합쳐 최종 답을 만듭니다.
Ian Stephenson 이이 코드 를 작성 했는데, 이 코드 는 성능이 뛰어납니다 pow(). 그는 아이디어 를 다음과 같이 설명합니다 .
Pow는 기본적으로 log를 사용하여 구현됩니다
pow(a,b)=x(logx(a)*b). 따라서 빠른 로그와 빠른 지수가 필요합니다. x가 무엇인지는 중요하지 않으므로 2를 사용합니다. 트릭은 부동 소수점 숫자가 이미 로그 스타일 형식에 있다는 것입니다.a=M*2E양쪽의 로그를 취하면 다음을 얻을 수 있습니다.
log2(a)=log2(M)+E또는 더 간단하게 :
log2(a)~=E즉, 숫자의 부동 소수점 표현을 취하고 지수를 추출하면 로그로 좋은 시작점이되는 것을 얻게됩니다. 비트 패턴을 마사지하여이 작업을 수행하면 가수가 결국 오류에 대한 좋은 근사치를 제공하고 잘 작동합니다.
이것은 간단한 조명 계산에 충분할 것입니다.하지만 더 나은 것이 필요하다면 가수를 추출하고이를 사용하여 꽤 정확한 2 차 보정 계수를 계산할 수 있습니다.
우선, 플로트를 사용하는 것은 요즘 대부분의 기계에서 많이 구매하지 않을 것입니다. 사실 복식이 더 빠를 수 있습니다. 당신의 힘 1.0 / 2.4는 5/12 또는 1 / 3 * (1 + 1 / 4)입니다. 이것이 cbrt (한 번)와 sqrt (두 번!)를 호출하더라도 여전히 pow ()를 사용하는 것보다 두 배 빠릅니다. (최적화 : -O3, 컴파일러 : i686-apple-darwin10-g ++-4.2.1).
#include <math.h> // cmath does not provide cbrt; C99 does.
double xpow512 (double x) {
double cbrtx = cbrt(x);
return cbrtx*sqrt(sqrt(cbrtx));
}
이것은 귀하의 질문에 대답하지 않을 수 있습니다.
2.4f및 1/2.4f그 정확히 RGB 색 공간 선형의 sRGB와 사이의 변환에 사용되는 힘 때문에, 나를 매우 의심합니다. 그래서 당신은 실제로 최적화하기 위해 시도 할 수있는 , 그것을 구체적으로. 모르겠습니다. 이것이 귀하의 질문에 대한 답변이 아닌 이유입니다.
이 경우 조회 테이블을 사용해보십시오. 다음과 같은 것 :
__attribute__((aligned(64))
static const unsigned short SRGB_TO_LINEAR[256] = { ... };
__attribute__((aligned(64))
static const unsigned short LINEAR_TO_SRGB[256] = { ... };
void apply_lut(const unsigned short lut[256], unsigned char *src, ...
16 비트 데이터를 사용하는 경우 적절하게 변경하십시오. 어쨌든 8 비트 데이터로 작업 할 때 필요한 경우 결과를 디더링 할 수 있도록 표를 16 비트로 만들 것입니다. 데이터가 부동 소수점 인 경우에는 분명히 잘 작동하지 않지만 sRGB 데이터를 부동 소수점에 저장하는 것은 실제로 의미가 없으므로 먼저 16 비트 / 8 비트로 변환하는 것이 좋습니다. 그런 다음 선형에서 sRGB로 변환합니다.
(sRGB가 부동 소수점으로 이해되지 않는 이유는 HDR이 선형이어야하고 sRGB는 디스크에 저장하거나 화면에 표시하는 데만 편리하지만 조작에는 편리하지 않기 때문입니다.)
이항 시리즈 는 상수 지수를 설명하지만 모든 입력을 [1,2) 범위로 정규화 할 수있는 경우에만 사용할 수 있습니다. ((1 + x) ^ a를 계산합니다). 원하는 정확도를 위해 필요한 항 수를 결정하려면 몇 가지 분석을 수행해야합니다.
빠른 sRGB <-> 선형 RGB 변환을 수행하는 방법 인 당신이 정말로 묻고 싶은 질문에 답할 것 입니다. 이를 정확하고 효율적으로 수행하기 위해 다항식 근사를 사용할 수 있습니다. 다음 다항식 근사는 sollya로 생성되었으며 최악의 경우 상대 오차는 0.0144 %입니다.
inline double poly7(double x, double a, double b, double c, double d,
double e, double f, double g, double h) {
double ab, cd, ef, gh, abcd, efgh, x2, x4;
x2 = x*x; x4 = x2*x2;
ab = a*x + b; cd = c*x + d;
ef = e*x + f; gh = g*x + h;
abcd = ab*x2 + cd; efgh = ef*x2 + gh;
return abcd*x4 + efgh;
}
inline double srgb_to_linear(double x) {
if (x <= 0.04045) return x / 12.92;
// Polynomial approximation of ((x+0.055)/1.055)^2.4.
return poly7(x, 0.15237971711927983387,
-0.57235993072870072762,
0.92097986411523535821,
-0.90208229831912012386,
0.88348956209696805075,
0.48110797889132134175,
0.03563925285274562038,
0.00084585397227064120);
}
inline double linear_to_srgb(double x) {
if (x <= 0.0031308) return x * 12.92;
// Piecewise polynomial approximation (divided by x^3)
// of 1.055 * x^(1/2.4) - 0.055.
if (x <= 0.0523) return poly7(x, -6681.49576364495442248881,
1224.97114922729451791383,
-100.23413743425112443219,
6.60361150127077944916,
0.06114808961060447245,
-0.00022244138470139442,
0.00000041231840827815,
-0.00000000035133685895) / (x*x*x);
return poly7(x, -0.18730034115395793881,
0.64677431008037400417,
-0.99032868647877825286,
1.20939072663263713636,
0.33433459165487383613,
-0.01345095746411287783,
0.00044351684288719036,
-0.00000664263587520855) / (x*x*x);
}
그리고 다항식을 생성하는 데 사용 된 sollya 입력 :
suppressmessage(174);
f = ((x+0.055)/1.055)^2.4;
p0 = fpminimax(f, 7, [|D...|], [0.04045;1], relative);
p = fpminimax(f/(p0(1)+1e-18), 7, [|D...|], [0.04045;1], relative);
print("relative:", dirtyinfnorm((f-p)/f, [s;1]));
print("absolute:", dirtyinfnorm((f-p), [s;1]));
print(canonical(p));
s = 0.0523;
z = 3;
f = 1.055 * x^(1/2.4) - 0.055;
p = fpminimax(1.055 * (x^(z+1/2.4) - 0.055*x^z/1.055), 7, [|D...|], [0.0031308;s], relative)/x^z;
print("relative:", dirtyinfnorm((f-p)/f, [0.0031308;s]));
print("absolute:", dirtyinfnorm((f-p), [0.0031308;s]));
print(canonical(p));
p = fpminimax(1.055 * (x^(z+1/2.4) - 0.055*x^z/1.055), 7, [|D...|], [s;1], relative)/x^z;
print("relative:", dirtyinfnorm((f-p)/f, [s;1]));
print("absolute:", dirtyinfnorm((f-p), [s;1]));
print(canonical(p));
For exponents of 2.4, you could either make a lookup table for all your 2.4 values and lirp or perhaps higher-order function to fill in the in-betweem values if the table wasn't accurate enough (basically a huge log table.)
Or, value squared * value to the 2/5s which could take the initial square value from the first half of the function and then 5th root it. For the 5th root, you could Newton it or do some other fast approximator, though honestly once you get to this point, your probably better off just doing the exp and log functions with the appropriate abbreviated series functions yourself.
The following is an idea you can use with any of the fast calculation methods. Whether it helps things go faster depends on how your data arrives. You can use the fact that if you know x and pow(x, n), you can use the rate of change of the power to compute a reasonable approximation of pow(x + delta, n) for small delta, with a single multiply and add (more or less). If successive values you feed your power functions are close enough together, this would amortize the full cost of the accurate calculation over multiple function calls. Note that you don't need an extra pow calculation to get the derivative. You could extend this to use the second derivative so you can use a quadratic, which would increase the delta you could use and still get the same accuracy.
So traditionally the powf(x, p) = x^p is solved by rewriting x as x=2^(log2(x)) making powf(x,p) = 2^(p*log2(x)), which transforms the problem into two approximations exp2() & log2(). This has the advantage of working with larger powers p, however the downside is that this is not the optimal solution for a constant power p and over a specified input bound 0 ≤ x ≤ 1.
When the power p > 1, the answer is a trivial minimax polynomial over the bound 0 ≤ x ≤ 1, which is the case for p = 12/5 = 2.4 as can be seen below:
float pow12_5(float x){
float mp;
// Minimax horner polynomials for x^(5/12), Note: choose the accurarcy required then implement with fma() [Fused Multiply Accumulates]
// mp = 0x4.a84a38p-12 + x * (-0xd.e5648p-8 + x * (0xa.d82fep-4 + x * 0x6.062668p-4)); // 1.13705697e-3
mp = 0x1.117542p-12 + x * (-0x5.91e6ap-8 + x * (0x8.0f50ep-4 + x * (0xa.aa231p-4 + x * (-0x2.62787p-4)))); // 2.6079002e-4
// mp = 0x5.a522ap-16 + x * (-0x2.d997fcp-8 + x * (0x6.8f6d1p-4 + x * (0xf.21285p-4 + x * (-0x7.b5b248p-4 + x * 0x2.32b668p-4)))); // 8.61377e-5
// mp = 0x2.4f5538p-16 + x * (-0x1.abcdecp-8 + x * (0x5.97464p-4 + x * (0x1.399edap0 + x * (-0x1.0d363ap0 + x * (0xa.a54a3p-4 + x * (-0x2.e8a77cp-4)))))); // 3.524655e-5
return(mp);
}
However when p < 1 the minimax approximation over the bound 0 ≤ x ≤ 1 does not appropriately converge to the desired accuracy. One option [not really] is to rewrite the problem y=x^p=x^(p+m)/x^m where m=1,2,3 is a positive integer, making the new power approximation p > 1 but this introduces division which is inherently slower.
There's however another option which is to decompose the input x as its floating point exponent and mantissa form:
x = mx* 2^(ex) where 1 ≤ mx < 2
y = x^(5/12) = mx^(5/12) * 2^((5/12)*ex), let ey = floor(5*ex/12), k = (5*ex) % 12
= mx^(5/12) * 2^(k/12) * 2^(ey)
The minimax approximation of mx^(5/12) over 1 ≤ mx < 2 now converges much faster than before, without division, but requires 12 point LUT for the 2^(k/12). The code is below:
float powk_12LUT[] = {0x1.0p0, 0x1.0f38fap0, 0x1.1f59acp0, 0x1.306fep0, 0x1.428a3p0, 0x1.55b81p0, 0x1.6a09e6p0, 0x1.7f910ep0, 0x1.965feap0, 0x1.ae89fap0, 0x1.c823ep0, 0x1.e3437ep0};
float pow5_12(float x){
union{float f; uint32_t u;} v, e2;
float poff, m, e, ei;
int xe;
v.f = x;
xe = ((v.u >> 23) - 127);
if(xe < -127) return(0.0f);
// Calculate remainder k in 2^(k/12) to find LUT
e = xe * (5.0f/12.0f);
ei = floorf(e);
poff = powk_12LUT[(int)(12.0f * (e - ei))];
e2.u = ((int)ei + 127) << 23; // Calculate the exponent
v.u = (v.u & ~(0xFFuL << 23)) | (0x7FuL << 23); // Normalize exponent to zero
// Approximate mx^(5/12) on [1,2), with appropriate degree minimax
// m = 0x8.87592p-4 + v.f * (0x8.8f056p-4 + v.f * (-0x1.134044p-4)); // 7.6125e-4
// m = 0x7.582138p-4 + v.f * (0xb.1666bp-4 + v.f * (-0x2.d21954p-4 + v.f * 0x6.3ea0cp-8)); // 8.4522726e-5
m = 0x6.9465cp-4 + v.f * (0xd.43015p-4 + v.f * (-0x5.17b2a8p-4 + v.f * (0x1.6cb1f8p-4 + v.f * (-0x2.c5b76p-8)))); // 1.04091259e-5
// m = 0x6.08242p-4 + v.f * (0xf.352bdp-4 + v.f * (-0x7.d0c1bp-4 + v.f * (0x3.4d153p-4 + v.f * (-0xc.f7a42p-8 + v.f * 0x1.5d840cp-8)))); // 1.367401e-6
return(m * poff * e2.f);
}
참고URL : https://stackoverflow.com/questions/6475373/optimizations-for-pow-with-const-non-integer-exponent
'IT TIP' 카테고리의 다른 글
| 페이지의 특정 부분 만 표시하는 iframe (0) | 2020.12.07 |
|---|---|
| 문자열 매개 변수가있는 C # TrimStart (0) | 2020.12.07 |
| 페이지로드시 CSS 전환 실행 중지 (0) | 2020.12.07 |
| Android에서 머티리얼 디자인 원형 진행률 표시 줄을 구현하는 방법 (0) | 2020.12.07 |
| Postgres 9.3에서 현재 시간대 이름을 얻으려면 어떻게해야합니까? (0) | 2020.12.07 |