4

私は現在、余弦の近似に取り組んでいます。最終的なターゲット デバイスは 32 ビット浮動小数点 ALU/LU で動作する自己開発であり、C に特化したコンパイラがあるため、C ライブラリの数学関数 (cosf など) を使用できません。精度と命令/サイクル数の点で異なるさまざまな方法をコーディングすることを目指しています。

fdlibm、テイラー展開、パデ近似、maple を使用した remez アルゴリズムなど、さまざまな近似アルゴリズムを既に試しました。

しかし、浮動小数点精度のみを使用してそれらを実装するとすぐに、精度が大幅に失われます。そして、確かに:倍精度では、はるかに高い精度がまったく問題にならないことを私は知っています...

現在、pi/2 (最大のエラーが発生する範囲) 付近で数千 ulp までの正確な近似値がいくつかありますが、単精度変換によって制限されていると感じています。

トピックの引数削減に対処するには: 入力はラジアンです。引数を減らすと、除算/乗算により精度がさらに低下すると思います....私の全体的な入力範囲は0..piしかないので、引数を0..pi/2に減らすことにしました。

したがって、私の質問は次のとおりです。コサイン関数の単精度近似を高精度 (および最良の場合は高効率) で知っている人はいますか? 単精度の近似を最適化するアルゴリズムはありますか? 組み込みの cosf 関数が内部的に単精度または倍精度で値を計算するかどうか知っていますか? 〜

float ua_cos_v2(float x)
{
    float output;
    float myPi = 3.1415927410125732421875f;
    if (x < 0) x = -x;
    int quad = (int32_t)(x*0.63661977236f);//quad = x/(pi/2) = x*2/pi
    if (x<1.58f && x> 1.57f) //exclude approximation around pi/2
    {
        output = -(x - 1.57079637050628662109375f) - 2.0e-12f*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f) + 0.16666667163372039794921875f*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f) + 2.0e-13f*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)+ 0.000198412701138295233249664306640625f*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f);
        output -= 4.37E-08f;
    }
    else {
        float param_x;
        int param_quad = -1;
        switch (quad)
        {
        case 0:
            param_x = x;
            break;
        case 1:
            param_x = myPi - x;
            param_quad = 1;
            break;
        case 2:
            param_x = x - myPi;
            break;
        case 3:
            param_x = 2 * myPi - x;
            break;
        }
        float c1 = 1.0f,
            c2 = -0.5f,
            c3 = 0.0416666679084300994873046875f,
            c4 = -0.001388888922519981861114501953125f,
            c5 = 0.00002480158218531869351863861083984375f,
            c6 = -2.75569362884198199026286602020263671875E-7f,
            c7 = 2.08583283978214240050874650478363037109375E-9f,
            c8 = -1.10807162057025010426514199934899806976318359375E-11f;
        float _x2 = param_x * param_x;
        output = c1 + _x2*(c2 + _x2*(c3 + _x2*(c4 + _x2*(c5 + _x2*(c6 + _x2*(c7 
        + _x2* c8))))));
        if (param_quad == 1 || param_quad == 0)
            output = -output;
    }
    return output;
}

情報を忘れた場合は、お気軽にお問い合わせください。

前もって感謝します

4

2 に答える 2