コードにホットスポットがpow()あり、実行時間の約10〜20％を占めています。
私の入力は非常に具体的であるため、より高いパフォーマンスで2つの近似（各指数に1つ）pow(x,y)をロールする方法があるかどうか疑問に思っています。pow()
floatベクトルを使用しています。プラットフォームの詳細を利用できる場合は、すぐに！float完全な精度（の）アルゴリズムにも関心がありますが、最大エラー率は約0.01％が理想的です。
私はすでに高速pow() 近似を使用していますが、これらの制約は考慮されていません。もっとうまくやることは可能ですか？
これは私の前の答えとは非常に異なっており、これは非常に速いので、別の答え。相対誤差は3e-8です。より正確にしたいですか？Chebychevの用語をさらにいくつか追加します。2^n-イプシロンと2^n +イプシロンの間に小さな不連続性が生じるため、順序を奇数に保つのが最善です。
#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);
}
補遺：ここで何が起こっているのですか？
リクエストごとに、以下は上記のコードがどのように機能するかを説明しています。
概要
上記のコードは、2つの関数とを定義しdouble pow512norm (double x)てdouble pow512 (double x)います。後者はスイートへのエントリポイントです。これは、x ^（5/12）を計算するためにユーザーコードが呼び出す必要のある関数です。この関数pow512norm(x)は、チェビシェフ多項式を使用してx ^（5/12）を近似しますが、範囲[1,2]のxに対してのみです。（pow512norm(x)その範囲外のxの値に使用すると、結果はガベージになります。）
この関数は、1≤<code>s<2となるようにpow512(x)着信xをペアに分割します。とが12未満になるようにさらに分割すると、x ^（5/12）を見つける問題を部分に分割できます。(double s, int n)x = s * 2^nn(int q, unsigned int r)n = 12*q + rr
x^(5/12)=(s^(5/12))*((2^n)^(5/12))正のu、vおよび実数aの場合は（u v）^ a =（u ^ a） （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）に近いはずです。修辞的な質問：「近い」とはどういう意味ですか？2つの競合する解釈は、平均二乗誤差を最小化することと最大絶対誤差を最小化することです。
これらの違いを説明するために、株式市場の例えを使用します。あなたがあなたの最終的な引退のために貯金したいとします。20代の場合、最善の方法は株式や株式市場のファンドに投資することです。これは、十分に長い期間にわたって、株式市場は平均して他の投資スキームを上回っているためです。しかし、私たちは皆、お金を株に入れることが非常に悪いことであるのを見てきました。50代または60代（または若くして引退したい場合は40代）の場合は、もう少し保守的に投資する必要があります。これらの下降傾向は、退職後のポートフォリオに悪影響を与える可能性があります。
関数近似に戻る：ある近似の消費者として、あなたは通常、「平均的な」パフォーマンスではなく、最悪の場合のエラーについて心配しています。「平均」で最高のパフォーマンス（最小二乗法など）を実現するように構築された近似を使用します。マーフィーの法則では、パフォーマンスが平均よりはるかに悪い場合に、プログラムは近似を使用して多くの時間を費やします。必要なのは、ミニマックス近似です。これは、一部のドメインで最大絶対誤差を最小化するものです。優れた数学ライブラリは、最小二乗アプローチではなくミニマックスアプローチを採用します。これにより、数学ライブラリの作成者は、ライブラリのパフォーマンスを保証できるからです。
数学ライブラリは通常、多項式または有理多項式を使用して、ある定義域a≤x≤bである関数f（x）を近似します。関数f（x）がこの定義域で分析的であり、関数を次数Nの多項式p（x）で近似したいとします。与えられた次数Nに対して、p（ x）-f（x）は[a、b]上にN + 2極値を持ち、これらのN+2極値の絶対値はすべて互いに等しくなります。この魔法の多項式p（x）を見つけることは、関数近似器の聖杯です。
私はあなたのためにその聖杯を見つけられませんでした。代わりに、チェビシェフ近似を使用しました。第1種のチェビシェフ多項式は、関数近似に関していくつかの非常に優れた機能を備えた直交（正規直交ではない）多項式のセットです。チェビシェフ近似は、多くの場合、その魔法の多項式p（x）に非常に近いものです。（実際、聖杯多項式が通常チェビシェフ近似で始まることを検出するRemez交換アルゴリズム。）
pow512norm（x）
この関数は、チェビシェフ近似を使用して、x ^（5/12）を近似する多項式p *（x）を見つけます。ここでは、p *（x）を使用して、このチェビシェフ近似を上記の魔法の多項式p（x）と区別しています。チェビシェフ近似p*（x）は簡単に見つけることができます。p（x）を見つけるのはクマです。チェビシェフ近似p*（x）はsum_i Cn [i] * Tn（i、x）です。ここで、Cn [i]はチェビシェフ係数であり、Tn（i、x）はxで評価されたチェビシェフ多項式です。
Wolframalphaを使用してチェビシェフ係数を見つけCnました。たとえば、これはを計算しCn[1]ます。入力ボックスの後の最初のボックスには、目的の答え（この場合は0.166658）があります。それは私が望むほど多くの桁ではありません。「その他の桁」をクリックすると、出来上がり、さらに多くの桁が表示されます。Wolframalphaは無料です。実行する計算量には制限があります。高階項でその制限に達します。（mathematicaを購入するか、数学にアクセスできる場合は、これらの高次係数を高精度で計算できます。）
チェビシェフ多項式Tn（x）は、配列で計算されTnます。魔法の多項式p（x）に非常に近いものを与える以外に、チェビシェフ近似を使用するもう1つの理由は、これらのチェビシェフ多項式の値が簡単に計算Tn[0]=1されるTn[1]=xことTn[i]=2*x*Tn[i-1] - Tn[i-2]です。（コードでは「i」ではなく「ii」をインデックス変数として使用しました。変数名として「i」を使用することはありません。英語で単語に「i」が含まれる単語はいくつありますか？2つある単語はいくつありますか？連続した'i's？）
pow512（x）
pow512は、ユーザーコードが呼び出す必要のある関数です。この関数の基本については、すでに説明しました。さらにいくつかの詳細：数学ライブラリ関数は、入力の仮数と指数をfrexp(x)返します。（マイナーな問題： 1から2の間で使用したいのですが、0.5から1の間の値を返します。）数学ライブラリ関数は、1つのスウェルフープで整数除算の商と剰余を返します。最後に、数学ライブラリ関数を使用して3つの部分を組み合わせ、最終的な答えを作成します。siexpxspow512normfrexpdivldexp
IEEE 754ハッキングの流れの中で、より高速で「魔法」の少ない別のソリューションがあります。約12クロックサイクルで.08％の許容誤差を達成します（IntelMeromCPUでp=2.4の場合）。
浮動小数点数は元々対数の近似値として考案されたものであるため、整数値をの近似値として使用できますlog2。これは、convert-from-integer命令を浮動小数点値に適用して、別の浮動小数点値を取得することで、ある程度移植可能になります。
計算を完了するにはpow、定数係数を掛けて、整数に変換する命令を使用して対数を変換し直します。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を掛けることができますが、さらにいくつかの命令を実行すると、10進数で約4桁の精度が得られます。これは、グラフィックスには十分です。
重要なのは、過大評価と過小評価を一致させ、平均を取ることです。
rsqrtps。（これは十分に正確ですが、ゼロで動作する能力を犠牲にします。）mulps。rsqrtps。mulps。mulps。mulps。これは過大評価です。mulps。これは過小評価です。addps1、1 mulps。命令の集計：レイテンシー= 5の2つの変換と、スループット=4の2つの逆数平方根推定値を含む14。
平均を適切にとるために、予想される誤差によって推定値に重みを付けます。過小評価すると、誤差が0.6対0.4の累乗になるため、誤差の1.5倍になると予想されます。重み付けは命令を追加しません。プレファクターで行うことができます。係数をaと呼ぶ：a ^ 0.5 = 1.5 a ^ -0.75、およびa=1.38316186。
fastpow最終的な誤差は約.015％であり、最初の結果より2桁優れています。実行時間は、ソース変数と宛先変数を使用するビジーループの場合、約12サイクルですvolatile。反復と重複していますが、実際の使用では、命令レベルの並列性も見られます。SIMDを考慮すると、これは3サイクルあたり1つのスカラー結果のスループットです。
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）。
小さなテストハーネスと、上記に対応する2つの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は、基本的にログを使用して実装されます
pow(a,b)=x(logx(a)*b)。したがって、高速の対数と高速の指数が必要です。xが何であるかは関係ないため、2を使用します。トリックは、浮動小数点数がすでに対数形式になっていることです。a=M*2E両側のログを取ると、次のようになります。
log2(a)=log2(M)+Eまたはもっと簡単に：
log2(a)~=E言い換えると、数値の浮動小数点表現を取得し、指数を抽出すると、対数として適切な開始点となるものが得られます。ビットパターンをマッサージしてこれを行うと、仮数はエラーに適切に近似することになり、かなりうまく機能することがわかります。
これは単純な照明計算には十分なはずですが、より良いものが必要な場合は、仮数を抽出し、それを使用してかなり正確な2次補正係数を計算できます。
まず、フロートを使用しても、最近のほとんどのマシンではあまり購入されません。実際、doubleはより速くなる可能性があります。あなたのパワー、1.0 /2.4は5/12または1/3*（1 + 1/4）です。これはcbrt（1回）とsqrt（2回！）を呼び出していますが、それでもpow（）を使用する場合の2倍の速度です。（最適化：-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));
}
これはあなたの質問に答えないかもしれません。
これらは、sRGB2.4fと1/2.4f線形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を計算することに注意してください）。必要な精度を得るために必要な用語の数を決定するには、いくつかの分析を行う必要があります。
2.4の指数の場合、すべての2.4値とlirpのルックアップテーブルを作成するか、テーブルの精度が十分でない場合は高階関数を使用して中間値を入力できます（基本的には巨大なログテーブル）。
または、値の2乗*値を2/5に設定します。これにより、関数の前半から最初の2乗値を取得し、5番目にルートすることができます。5乗根については、ニュートンまたは他の高速近似を行うことができますが、正直なところ、このポイントに到達したら、適切な省略された系列関数を使用してexp関数とlog関数を自分で実行する方がよいでしょう。
以下は、高速計算方法のいずれかで使用できるアイデアです。それが物事をより速く進めるのに役立つかどうかは、データがどのように到着するかによって異なります。とを知っている場合は、累乗の変化率を使用して、1回の乗算と加算（多かれ少なかれ）で、の妥当な近似値を計算できるという事実をx使用pow(x, n)できますpow(x + delta, n)。delta電力関数に供給する連続する値が十分に近い場合、これにより、複数の関数呼び出しにわたる正確な計算の全コストが償却されます。導関数を取得するために追加のパワー計算は必要ないことに注意してください。これを拡張して2次導関数を使用できるため、2次式を使用できます。これにより、使用できる量が増え、delta同じ精度が得られます。
したがって、伝統的には、を作成するようにpowf(x, p) = x^p書き直すことで解決されます。これにより、問題が2つの近似値＆に変換されます。これには、より大きな電力で動作するという利点がありますが、欠点は、これが一定の電力および指定された入力境界を超える場合の最適なソリューションではないことです。xx=2^(log2(x))powf(x,p) = 2^(p*log2(x))exp2()log2()pp0 ≤ 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適切に目的の精度に収束しない場合。[実際には]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のミニマックス近似は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);
}