29

固定小数点クラス (10.22) があり、pow、sqrt、exp、log 関数が必要です。

残念ながら、これをどこから始めればよいかわかりません。有用な記事へのリンクを提供してくれる人、またはコードを提供してくれる人はいますか?

exp 関数があれば、 pow と sqrt を実装するのが比較的簡単になると思います。

pow( x, y ) => exp( y * log( x ) )
sqrt( x )   => pow( x, 0.5 )

私が難しいと感じているのは、それらの exp 関数と log 関数だけです (ログ規則のいくつかを覚えているかのように、それ以外のことはあまり覚えていません)。

おそらく、sqrt と pow のためのより高速な方法もあるでしょう。

注意: これはクロス プラットフォームであり、純粋な C/C++ コードである必要があるため、アセンブラーの最適化は使用できません。

4

4 に答える 4

32

非常に簡単な解決策は、まともなテーブル駆動型の近似を使用することです。入力を正しく削減すれば、実際には多くのデータは必要ありません。exp(a)==exp(a/2)*exp(a/2)、つまり、実際に計算する必要があるのは だけexp(x)です1 < x < 2。その範囲で、runga-kutta 近似は ~16 エントリ IIRC で妥当な結果をもたらします。

同様に、sqrt(a) == 2 * sqrt(a/4) == sqrt(4*a) / 2これは のテーブル エントリのみが必要であることを意味します1 < a < 4。Log(a) は少し難しいです: log(a) == 1 + log(a/e). これはかなり遅い反復ですが、log(1024) はわずか 6.9 であるため、多くの反復はありません。

pow: にも同様の「整数優先」アルゴリズムを使用しますpow(x,y)==pow(x, floor(y)) * pow(x, frac(y))。これpow(double, int)は些細なこと (分割統治) であるため機能します。

[編集] の整数コンポーネントのlog(a)場合、テーブルを格納すると便利な場合があるため、そのテーブル内の a の単純なハードコードされたバイナリ検索によって1, e, e^2, e^3, e^4, e^5, e^6, e^7削減できます。log(a) == n + log(a/e^n)7 ステップから 3 ステップへの改善はそれほど大きくありませんが、 で割るe^n代わりに で1n回割ればよいことを意味しますe。

[編集 2] そして、その最後のlog(a/e^n)項については、使用できます- 各反復は、テーブル ルックアップによってlog(a/e^n) = log((a/e^n)^8)/8さらに 3 ビットを生成します。これにより、コードとテーブルのサイズが小さく保たれます。これは通常、組み込みシステム用のコードであり、大きなキャッシュはありません。

[編集 3] それは私の側ではまだ賢明ではありません。log(a) = log(2) + log(a/2). 固定小数点値を保存し、log2=0.6931471805599先頭のゼロの数を数え、aルックアップ テーブルに使用される範囲にシフトし、そのシフト (整数) を固定小数点定数で乗算するだけlog2です。命令数は 3 つまでです。

削減ステップに使用eすると、「良い」log(e)=1.0定数が得られますが、それは誤った最適化です。0.6931471805599 は 1.0 と同じくらい良い定数です。どちらも 10.22 固定小数点の 32 ビット定数です。範囲削減の定数として 2 を使用すると、除算にビット シフトを使用できます。

[編集 5] Q10.22 に保存しているので、log(65536)=11.09035488 の方が適切に保存できます。(16 x log(2))。「x16」は、さらに 4 ビットの精度が利用可能であることを意味します。

編集 2 からのトリックを引き続き取得しますlog(a/2^n) = log((a/2^n)^8)/8。基本的に、これにより結果が得られます(a + b/8 + c/64 + d/512) * 0.6931471805599-範囲[0,7]のb、c、d。a.bcd実際には 8 進数です。パワーとして 8 を使用したので、驚くことではありません。(このトリックは、べき乗が 2、4、または 16 のいずれでも同じように機能します。)

[編集 4] まだオープンエンドがありました。pow(x, frac(y)はちょうどpow(sqrt(x), 2 * frac(y))であり、私たちはまともな1/sqrt(x). これにより、はるかに効率的なアプローチが可能になります。frac(y)=0.101バイナリ、つまり 1/2 プラス 1/8 と言ってください。ということx^0.101は、 です(x^1/2 * x^1/8)。でも、ただいまx^1/2です。もう 1 つの演算を保存すると、ニュートン ラフソンが与えるので、 を計算します。最終結果を反転するだけで、sqrt 関数を直接使用しないでください。sqrt(x)x^1/8(sqrt(sqrt(sqrt(x)))NR(x)1/sqrt(x)1.0/(NR(x)*NR((NR(NR(x)))

于 2011-01-11T12:57:19.557 に答える
19

以下は、Clay S. Turner の固定小数点対数基数 2 アルゴリズムの C 実装の例です [ 1 ]。このアルゴリズムには、ルックアップ テーブルは一切必要ありません。これは、多くのマイクロコントローラの場合のように、メモリの制約が厳しく、プロセッサに FPU がないシステムで役立ちます。対数の基数eと対数の基数 10 は、対数のプロパティを使用することによってもサポートされます。

          logₘ(x)
logₙ(x) = ───────
          logₘ(n)

ここで、このアルゴリズムでは、mは 2 です。

この実装の優れた機能は、可変精度をサポートしていることです。精度は実行時に決定できますが、範囲が犠牲になります。私が実装した方法では、プロセッサ (またはコンパイラ) は、中間結果を保持するために 64 ビット演算を実行できる必要があります。64 ビットのサポートを必要としないように簡単に適合させることができますが、範囲は縮小されます。

これらの関数を使用する場合、xは、指定された に従ってスケーリングされた固定小数点値であることが期待されますprecision。たとえば、precision16 の場合、 x2^16 (65536) でスケーリングする必要があります。結果は、入力と同じスケール係数を持つ固定小数点値です。の戻り値はINT32_MIN負の無限大を表します。の戻り値はINT32_MAXエラーを示し、入力精度が無効であることを示すerrnoに設定されます。EINVAL

#include <errno.h>
#include <stddef.h>

#include "log2fix.h"

#define INV_LOG2_E_Q1DOT31  UINT64_C(0x58b90bfc) // Inverse log base 2 of e
#define INV_LOG2_10_Q1DOT31 UINT64_C(0x268826a1) // Inverse log base 2 of 10

int32_t log2fix (uint32_t x, size_t precision)
{
    int32_t b = 1U << (precision - 1);
    int32_t y = 0;

    if (precision < 1 || precision > 31) {
        errno = EINVAL;
        return INT32_MAX; // indicates an error
    }

    if (x == 0) {
        return INT32_MIN; // represents negative infinity
    }

    while (x < 1U << precision) {
        x <<= 1;
        y -= 1U << precision;
    }

    while (x >= 2U << precision) {
        x >>= 1;
        y += 1U << precision;
    }

    uint64_t z = x;

    for (size_t i = 0; i < precision; i++) {
        z = z * z >> precision;
        if (z >= 2U << (uint64_t)precision) {
            z >>= 1;
            y += b;
        }
        b >>= 1;
    }

    return y;
}

int32_t logfix (uint32_t x, size_t precision)
{
    uint64_t t;

    t = log2fix(x, precision) * INV_LOG2_E_Q1DOT31;

    return t >> 31;
}

int32_t log10fix (uint32_t x, size_t precision)
{
    uint64_t t;

    t = log2fix(x, precision) * INV_LOG2_10_Q1DOT31;

    return t >> 31;
}

この実装のコードは、標準入力から読み取った数値から対数を計算して表示するためにこの関数を使用する方法を示すサンプル/テスト プログラムと共に、Githubにもあります。

[ 1 ] CS Turner、「A Fast Binary Logarithm Algorithm」、IEEE Signal Processing Mag. 、pp。124,140、2010 年 9 月。

于 2013-02-14T21:59:31.653 に答える
5

良い出発点は、Jack Crenshaw の著書「Math Toolkit for Real-Time Programming」です。さまざまな超越関数のアルゴリズムと実装についての良い議論があります。

于 2011-01-11T12:16:04.010 に答える
4

整数演算のみを使用して固定小数点の sqrt の実装を確認してください。発明するのは楽しかったです。今ではかなり古い。

https://groups.google.com/forum/?hl=fr%05aacf5997b615c37&fromgroups#!topic/comp.lang.c/IpwKbw0MAxw/discussion

それ以外の場合は、アルゴリズムのCORDICセットを確認してください。これが、リストしたすべての関数と三角関数を実装する方法です。

編集：レビューされたソースをここのGitHubで公開しました

于 2012-05-20T15:37:54.827 に答える