我的代码中有热点，我正在pow()占用大约 10-20% 的执行时间。
我的输入pow(x,y)非常具体，所以我想知道是否有办法以pow()更高的性能滚动两个近似值（每个指数一个）：
float向量。如果可以利用平台细节，那就对了！大约 0.01% 的最大错误率是理想的，尽管我也对全精度 (for float) 算法感兴趣。
我已经在使用快速pow() 近似，但它没有考虑到这些约束。有没有可能做得更好？
另一个答案，因为这与我之前的答案非常不同，而且速度非常快。相对误差为 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^n1≤<code>s<2。进一步划分n为小于 12 让我将寻找 x^(5/12) 的问题分成几部分(int q, unsigned int r)：n = 12*q + rr
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)viau^(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]一下，我们就快到了。这qr.rem是r上面步骤 3 中计算的值。所需要的只是将其乘以2^(5*q)产生所需的结果。ldexp所做的。函数逼近
这里的目标是提出一个易于计算的 f(x)=x^(5/12) 逼近，它对于手头的问题“足够好”。在某种意义上，我们的近似值应该接近 f(x)。反问：“接近”是什么意思？两种相互竞争的解释是最小化均方误差与最小化最大绝对误差。
我将使用股票市场的类比来描述它们之间的区别。假设您想为最终的退休储蓄。如果你二十多岁，最好的办法是投资股票或股票市场基金。这是因为在足够长的时间跨度内，股市平均优于任何其他投资计划。然而，我们都见过把钱投入股票是一件非常糟糕的事情的时候。如果你是五六十岁（或者四十岁，如果你想年轻退休），你需要更保守地投资。这些下滑可能会对您的退休投资组合产生影响。
回到函数逼近：作为某种逼近的消费者，您通常担心最坏情况的错误，而不是“平均”性能。使用一些构造来“平均”提供最佳性能（例如最小二乘）的近似值，墨菲定律规定您的程序将花费大量时间在性能远低于平均水平的情况下使用近似值。你想要的是一个极小极大近似值，它可以最小化某个域上的最大绝对误差。一个好的数学库将采用极小值方法而不是最小二乘方法，因为这可以让数学库的作者为他们的库提供一些有保证的性能。
数学库通常使用多项式或有理多项式在某个域 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 交换算法通常以切比雪夫近似开始。）
pow512norm(x)
此函数使用切比雪夫近似来找到一些近似 x^(5/12) 的多项式 p*(x)。在这里，我使用 p*(x) 来区分这个 Chebyshev 近似和上面描述的神奇多项式 p(x)。Chebyshev 近似 p*(x) 很容易找到；发现 p(x) 是熊。Chebyshev 近似 p*(x) 是 sum_i Cn[i]*Tn(i,x)，其中 Cn[i] 是 Chebyshev 系数，Tn(i,x) 是在 x 处计算的 Chebyshev 多项式。
我使用 Wolfram alpha 为我找到切比雪夫系数Cn。例如，这计算Cn[1]. 输入框之后的第一个框有所需的答案，在本例中为 0.166658。这不是我想要的那么多数字。单击“更多数字”，瞧，您将获得更多数字。Wolfram alpha 是免费的；它会做多少计算是有限制的。它达到了高阶项的限制。（如果您购买或使用了mathematica，您将能够以高精度计算这些高阶系数。）
Chebyshev 多项式 Tn(x) 在数组 中计算Tn。除了给出非常接近神奇多项式 p(x) 的值之外，使用切比雪夫近似的另一个原因是这些切比雪夫多项式的值很容易计算：从 和 开始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 之间的值。）数学库函数在一个膨胀循环中返回整数除法的商和余数。最后，我使用数学库函数将这三个部分放在一起，形成最终答案。iexpxspow512normfrexpdivldexp
在 IEEE 754 hacking 中，这是另一种更快且不那么“神奇”的解决方案。它在大约十几个时钟周期内实现了 0.08% 的误差范围（对于 p=2.4 的情况，在 Intel Merom CPU 上）。
浮点数最初是作为对数的近似值而发明的，因此您可以将整数值用作 的近似值log2。通过将整数转换指令应用于浮点值以获得另一个浮点值，这在某种程度上是可以实现的。
要完成pow计算，您可以乘以一个常数因子，然后使用转换为整数指令将对数转换回来。在上交所，相关指令为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 ⁄<sub>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操作的限制。
Ian Stephenson 编写了这段代码，他声称它的性能优于其他代码pow()。他将这个想法描述如下：
Pow 基本上是使用 log's: 来实现的
pow(a,b)=x(logx(a)*b)。所以我们需要一个快速的对数和快速的指数——x 是什么并不重要，所以我们使用 2。诀窍是浮点数已经是对数样式的格式：a=M*2E取双方的日志给出：
log2(a)=log2(M)+E或更简单地说：
log2(a)~=E换句话说，如果我们采用一个数字的浮点表示，并提取指数，我们就会得到一个很好的起点作为它的对数。事实证明，当我们通过按摩位模式来做到这一点时，尾数最终给出了一个很好的误差近似值，而且效果很好。
这对于简单的照明计算应该足够好，但如果您需要更好的东西，您可以提取尾数，并使用它来计算非常准确的二次校正因子。
首先，现在在大多数机器上使用浮点数不会买太多。事实上，双打可以更快。你的力量，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因为这些正是用于在 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只方便存储在磁盘上或显示在屏幕上，不方便操作。）
我将回答您真正想问的问题，即如何进行快速 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));
二项式级数确实考虑了一个常数指数，但只有当您可以将所有输入标准化为范围 [1,2) 时，您才能使用它。（请注意，它计算 (1+x)^a）。您必须进行一些分析才能确定需要多少项才能达到所需的准确性。
以下是您可以与任何快速计算方法一起使用的想法。它是否有助于加快速度取决于您的数据是如何到达的。您可以使用这样一个事实，即如果您知道x和pow(x, n)，您可以使用幂的变化率来计算pow(x + delta, n)小的合理近似值delta，只需一次乘法和加法（或多或少）。如果你提供给你的幂函数的连续值足够接近，这将分摊多个函数调用的准确计算的全部成本。请注意，您不需要额外的 pow 计算来获得导数。您可以将其扩展为使用二阶导数，这样您就可以使用二次方，这将增加delta您可以使用的并且仍然获得相同的精度。
对于 2.4 的指数，您可以为所有 2.4 值和 lirp 创建一个查找表，或者如果表不够准确（基本上是一个巨大的日志表），则可以使用高阶函数来填充中间值。
或者，值平方 * 值到 2/5 秒，它可以从函数的前半部分获取初始平方值，然后对其进行第 5 次根。对于第 5 个根，您可以使用 Newton 或其他一些快速逼近器，但老实说，一旦您到达这一点，您最好自己使用适当的缩写系列函数执行 exp 和 log 函数。
所以传统上powf(x, p) = x^p是通过重写x为make 来解决x=2^(log2(x))的powf(x,p) = 2^(p*log2(x))，这将问题转换为两个近似值exp2()& log2()。这具有处理更大功率的优点p，但缺点是这不是恒定功率p和超过指定输入范围的最佳解决方案0 ≤ x ≤ 1。
当幂为 时p > 1，答案是在边界 上的一个平凡的极小极大多项式0 ≤ x ≤ 1，p = 12/5 = 2.4如下所示：
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);
}
但是，当p < 1边界上的极小极大近似值0 ≤ x ≤ 1不能适当地收敛到所需的精度时。一种选择[不是真的]是重写正整数的问题，进行新的幂近似y=x^p=x^(p+m)/x^m，但这会引入本质上较慢的除法。m=1,2,3p > 1
然而，还有另一种选择是将输入分解x为其浮点指数和尾数形式：
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)
mx^(5/12)over的 minimax 逼近1 ≤ mx < 2现在收敛速度比以前快得多，没有除法，但需要 12 点 LUT 用于2^(k/12). 代码如下：
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);
}