指数を使用するアルゴリズムの 1 つの漸近的な実行時間を決定しようとしていますが、指数がプログラムでどのように計算されるかわかりません。
特に、倍精度の浮動小数点数に使用される pow() アルゴリズムを探しています。
fdlibm の実装を見る機会がありました。コメントは、使用されたアルゴリズムを説明しています。
* n
* Method: Let x = 2 * (1+f)
* 1. Compute and return log2(x) in two pieces:
* log2(x) = w1 + w2,
* where w1 has 53-24 = 29 bit trailing zeros.
* 2. Perform y*log2(x) = n+y' by simulating muti-precision
* arithmetic, where |y'|<=0.5.
* 3. Return x**y = 2**n*exp(y'*log2)
処理されたすべての特殊なケース (0、1、inf、nan) のリストが続きます。
コードの最も重要なセクションは、すべての特殊なケースの処理の後で、log2
との2**
計算を伴います。そして、それらのいずれにもループはありません。したがって、浮動小数点プリミティブの複雑さにもかかわらず、漸近的に一定時間のアルゴリズムのように見えます。
浮動小数点の専門家 (私はその 1 人ではありません) からのコメントを歓迎します。:-)
彼らがそれを行うためのより良い方法を発見しない限り、三角関数、対数関数、および指数関数 (たとえば、指数関数的な成長と減衰) の近似値は、一般に、正確な近似結果を生成するために算術規則とテイラー級数展開を使用して計算されると思います。要求された精度内に。(べき級数、テイラー級数、マクローリン級数の関数展開の詳細については、微積分の本を参照してください。) これを行ってからしばらく経っているので、たとえば、正確な計算方法などを説明できなかったことに注意してください。含める必要がある級数の項の数は、倍精度計算で無視できるほど小さいエラーを保証します。
たとえば、e^x の Taylor/Maclaurin 級数展開は次のとおりです。
+inf [ x^k ] x^2 x^3 x^4 x^5
e^x = SUM [ --- ] = 1 + x + --- + ----- + ------- + --------- + ....
k=0 [ k! ] 2*1 3*2*1 4*3*2*1 5*4*3*2*1
すべての項 (0 から無限大までの k) を取得すると、この展開は正確かつ完全になります (エラーなし)。
ただし、すべての項を無限大にするのではなく、たとえば 5 項または 50 項などで停止すると、実際の e^x 関数値と余りが異なる近似結果が得られますが、これはかなり簡単です。計算する。
指数関数の良いニュースは、それがうまく収束し、その多項式展開の項が反復的にコーディングするのがかなり簡単であることです。各反復で寄与のサイズをテストし、ゼロに十分近くなったときに停止できるため、エラーが精度よりも小さいことを保証する必要があります。実際には、この戦略が実行可能かどうかはわかりません。試してみる必要があります。私が長い間忘れていた重要な詳細があります。次のようなもの: 機械精度、機械誤差、丸め誤差など。
また、e^x を使用していない場合でも、2^x や 10^x などの別の基数で成長/減衰を行っている場合は、近似多項式関数が変化することに注意してください。
x nの計算には exp(n*ln(x)) を使用できます。x と n はどちらも倍精度の浮動小数点数にすることができます。自然対数と指数関数は、テイラー級数を使用して計算できます。ここで数式を見つけることができます: http://en.wikipedia.org/wiki/Taylor_series
整数指数の場合、a を b に累乗する通常のアプローチは、次のようになります。
result = 1
while b > 0
if b is odd
result *= a
b -= 1
b /= 2
a = a * a
通常、指数のサイズは対数です。このアルゴリズムは、不変式の "a^b*result = a0^b0" に基づいています。ここで、a0 と b0 は a と b の初期値です。
負または非整数の指数については、対数と近似、および数値解析が必要です。実行時間は、使用されるアルゴリズムと、ライブラリが調整されている精度によって異なります。
編集: いくつかの関心があるように見えるので、余分な乗算のないバージョンを次に示します。
result = 1
while b > 0
while b is even
a = a * a
b = b / 2
result = result * a
b = b - 1
Intel を対象とする pow 関数を作成した場合、exp2(log2(x) * y) を返します。Intel の log2 用マイクロコードは、1 年生の微積分と大学院での数値解析を思い出すことができたとしても、私がコーディングできるどのコードよりも確実に高速です。
e^x = (1 + 分数) * (2^指数)、1 <= 1 + 分数 < 2
x * log2(e) = log2(1 + 分数) + 指数、0 <= log2(1 + 分数) < 1
指数 = フロア (x * log2(e))
1 + 分数 = 2^(x * log2(e) - 指数) = e^((x * log2(e) - 指数) * ln2) = e^(x - 指数 * ln2), 0 <= x - 指数* ln2 < ln2