Cor*_*son 60 c c++ math optimization exponent
我的代码中有热点,我pow()占用了大约10-20%的执行时间.
我的输入pow(x,y)是非常具体的,所以我想知道是否有办法滚动两个pow()近似值(每个指数一个)具有更高的性能:
float向量.如果可以利用平台细节,请立即使用!尽管我对全精度(for float)算法感兴趣,但最大错误率约为0.01%是理想的.
我已经在使用快速pow() 近似,但它没有考虑这些约束.有可能做得更好吗？
Dav*_*men 33
另一个答案,因为这与我以前的答案非常不同,而且这是非常快的.相对误差为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))via(u v)^ a =(u ^ a)(v ^ a)表示正u,v和real 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)积极的你,真正的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],我们差不多了.这qr.rem是r上面步骤3中计算的值.所需要的只是将其乘以2^(5*q)得到所需的结果.ldexp作用.函数逼近
这里的目标是提出一个容易计算的f(x)= x ^(5/12)近似值,它对于手头的问题"足够好".在某种意义上,我们的近似应接近f(x).修辞问题:"接近"是什么意思？两种相互竞争的解释是最小化均方误差与最小绝对误差最小化.
我将用股票市场类比来描述这些之间的差异.假设您想为最终退休储蓄.如果您二十几岁,最好的办法是投资股票或股票市场基金.这是因为在足够长的时间内,股票市场平均击败任何其他投资计划.然而,我们所有人都看到过将资金投入股市是非常糟糕的事情.如果你是五十或六十年代(或四十年代,如果你想退休年轻)你需要更保守地投资.这些下跌可能对您的退休投资组合产生影响.
回到函数逼近:作为某种近似的消费者,您通常担心最坏情况的错误而不是"平均"的性能.使用一些近似构造来"平均"提供最佳性能(例如,最小二乘法),墨菲定律规定您的程序将花费大量时间使用近似值,其性能远低于平均值.你想要的是一个minimax近似值,它可以最大限度地减少某个域上的最大绝对误差.一个好的数学库将采用minimax方法而不是最小二乘方法,因为这可以让数学库的作者给出一些保证其库的性能.
数学库通常使用多项式或有理多项式来近似某个域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近似.第一类的切比雪夫多项式是正交(但不是正交)多项式集合,当涉及函数逼近时具有一些非常好的特征.切比雪夫近似通常非常接近于神奇多项式p(x).(事实上,确实发现圣杯多项式的Remez交换算法通常以Chebyshev近似开始.)
pow512norm(x)
这个函数使用切比雪夫近似来找到近似x ^(5/12)的一些多项式p*(x).这里我用p*(x)来区分这个切比雪夫逼近与上述神奇多项式p(x).Chebyshev近似p*(x)很容易找到; 找到p(x)是一只熊.切比雪夫近似值p*(x)是sum_i Cn [i]*Tn(i,x),其中Cn [i]是切比雪夫系数,并且Tn(i,x)是在x处评估的切比雪夫多项式.
我使用Wolfram alpha Cn为我找到切比雪夫系数.例如,这计算Cn[1].输入框后面的第一个框具有所需答案,在这种情况下为0.166658.这并不像我想的那么多.点击"更多数字",瞧,你会得到更多的数字.Wolfram alpha是免费的; 它将进行多少计算是有限制的.它达到了更高阶条款的限制.(如果您购买或访问mathematica,您将能够高精度地计算这些高阶系数.)
在阵列中计算切比雪夫多项式Tn(x)Tn.除了给出非常接近神奇多项式p(x)的东西之外,使用切比雪夫近似的另一个原因是那些Chebyshev多项式的值很容易计算:从Tn[0]=1和开始Tn[1]=x,然后迭代计算Tn[i]=2*x*Tn[i-1] - Tn[i-2].(我在代码中使用'ii'作为索引变量而不是'i'.我从不使用'i'作为变量名.英语中有多少单词在单词中有'i'？有多少有两个单词连续'我是？)
pow512(x)
pow512是用户代码应该调用的函数.我已经描述了上面这个函数的基础知识.更多细节:数学库函数frexp(x)返回输入的有效数s和指数.(小问题:我希望在1和2之间使用但返回介于0.5和1之间的值.)数学库函数在一个swell foop中返回整数除法的商和余数.最后,我使用数学库函数将这三个部分放在一起形成最终答案.iexpxspow512normfrexpdivldexp
Pot*_*ter 24
在IEEE 754黑客攻击中,这是另一种更快,更少"神奇"的解决方案.它在大约十几个时钟周期内实现了0.08%的误差范围(对于p = 2.4,在Intel Merom CPU上).
浮点数最初是作为对数的近似值发明的,因此您可以使用整数值作为近似值log2.通过将整数转换指令应用于浮点值以获得另一个浮点值,可以在某种程度上实现这一点.
要完成pow计算,可以乘以常数因子并将对数转换回convert-to-integer指令.关于SSE,相关说明是cvtdq2ps和cvtps2dq.
但这并不是那么简单.IEEE 754中的指数字段是有符号的,偏差值为127,表示指数为零.在乘以对数之前必须删除此偏差,并在取幂之前重新添加.此外,通过减法进行偏差调整不会对零起作用.幸运的是,两种调整都可以通过预先乘以常数因子来实现.
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,但是更多的指令会给我们提供大约4个十进制数字的精度,这对于图形应该足够了.
关键是要将高估与低估相匹配,并取平均值.
rsqrtps.(这非常准确,但确实牺牲了零工作的能力.)mulps.rsqrtps.mulps.mulps.mulps.这是高估.mulps.这是低估的.addps,一mulps.指令计数:十四,包括两次延迟= 5的转换和两次倒数平方根估计,吞吐量= 4.
为了正确地取平均值,我们希望通过预期误差对估计值进行加权.低估将误差提高到0.6和0.4的幂,所以我们预计它是错误的1.5倍.加权不添加任何说明; 它可以在前因子中完成.调用系数a:a ^ 0.5 = 1.5 a ^ -0.75,a = 1.38316186.
最终误差约为0.015%,或比初始fastpow结果好2个数量级.对于具有volatile源和目标变量的繁忙循环,运行时大约是十几个循环......虽然它与迭代重叠,但实际使用情况也会看到指令级并行性.考虑到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).
我实现一个小的测试工具和两个X (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操作的限制.
Pen*_*One 20
Ian Stephenson编写了这个他声称优于其中的代码pow().他将这个想法描述如下:
Pow基本上是使用log实现的:
pow(a,b)=x(logx(a)*b).所以我们需要一个快速的日志和快速指数 - 无论什么是x,所以我们使用2.诀窍是浮点数已经是日志样式格式:Run Code Online (Sandbox Code Playgroud)a=M*2E取两边的日志给出:
Run Code Online (Sandbox Code Playgroud)log2(a)=log2(M)+E或更简单地说:
Run Code Online (Sandbox Code Playgroud)log2(a)~=E换句话说,如果我们采用数字的浮点表示,并提取指数,我们得到的东西就像它的日志一样是一个很好的起点.事实证明,当我们通过按摩位模式来实现这一点时,尾数最终给出了错误的良好近似,并且它的效果非常好.
这对于简单的光照计算应该足够好,但是如果你需要更好的东西,你可以提取尾数,并用它来计算非常准确的二次校正因子.
Dav*_*men 17
首先,使用花车不会在现今的大多数机器上买得多.事实上,双打可以更快.你的力量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));
}
Die*_*Epp 15
这可能无法解答您的问题.
在2.4f和1/2.4f让我很怀疑,因为那些正是用来sRGB和线性RGB色彩空间之间进行转换的权力.因此,您可能实际上正在尝试优化它.我不知道,这就是为什么这可能无法回答你的问题.
如果是这种情况,请尝试使用查找表.就像是:
__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位数据,请根据需要进行更改.无论如何我都会将表格设为16位,因此在使用8位数据时,如果需要,可以将结果抖动.如果你的数据是浮点数,这显然不会很好 - 但是将sRGB数据存储在浮点状态并没有用,所以你最好先转换为16位/ 8位数据然后进行从线性到sRGB的转换.
(sRGB没有意义,因为浮点是HDR应该是线性的,sRGB只是方便存储在磁盘上或在屏幕上显示,但不便于操作.)