实现基于表查找的三角函数

Question

对于我在业余时间实现的视频游戏，我尝试使用查找表实现我自己的 sinf()、cosf() 和 atan2f() 版本。目的是实现更快，但准确性较低。

我的初始实现如下。这些函数工作，并返回良好的近似值。唯一的问题是它们比调用标准 sinf()、cosf() 和 atan2f() 函数慢。

那么，我做错了什么？

// Geometry.h includes definitions of PI, TWO_PI, etc., as
// well as the prototypes for the public functions
#include "Geometry.h"

namespace {
    // Number of entries in the sin/cos lookup table
    const int SinTableCount = 512;

    // Angle covered by each table entry
    const float SinTableDelta = TWO_PI / (float)SinTableCount;

    // Lookup table for Sin() results
    float SinTable[SinTableCount];

    // This object initializes the contents of the SinTable array exactly once
    class SinTableInitializer {
    public:
        SinTableInitializer() {
            for (int i = 0; i < SinTableCount; ++i) {
                SinTable[i] = sinf((float)i * SinTableDelta);
            }
        }
    };
    static SinTableInitializer sinTableInitializer;

    // Number of entries in the atan lookup table
    const int AtanTableCount = 512;

    // Interval covered by each Atan table entry
    const float AtanTableDelta = 1.0f / (float)AtanTableCount;

    // Lookup table for Atan() results
    float AtanTable[AtanTableCount];

    // This object initializes the contents of the AtanTable array exactly once
    class AtanTableInitializer {
    public:
        AtanTableInitializer() {
            for (int i = 0; i < AtanTableCount; ++i) {
                AtanTable[i] = atanf((float)i * AtanTableDelta);
            }
        }
    };
    static AtanTableInitializer atanTableInitializer;

    // Lookup result in table.
    // Preconditions: y > 0, x > 0, y < x
    static float AtanLookup2(float y, float x) {
        assert(y > 0.0f);
        assert(x > 0.0f);
        assert(y < x);

        const float ratio = y / x;
        const int index = (int)(ratio / AtanTableDelta);
        return AtanTable[index];    
    }

}

float Sin(float angle) {
    // If angle is negative, reflect around X-axis and negate result
    bool mustNegateResult = false;
    if (angle < 0.0f) {
        mustNegateResult = true;
        angle = -angle;
    }

    // Normalize angle so that it is in the interval (0.0, PI)
    while (angle >= TWO_PI) {
        angle -= TWO_PI;
    }

    const int index = (int)(angle / SinTableDelta);
    const float result = SinTable[index];

    return mustNegateResult? (-result) : result;
}

float Cos(float angle) {
    return Sin(angle + PI_2);
}

float Atan2(float y, float x) {
    // Handle x == 0 or x == -0
    // (See atan2(3) for specification of sign-bit handling.)
    if (x == 0.0f) {
        if (y > 0.0f) {
            return PI_2;
        }
        else if (y < 0.0f) {
            return -PI_2;
        }
        else if (signbit(x)) {
            return signbit(y)? -PI : PI;
        }
        else {
            return signbit(y)? -0.0f : 0.0f;
        }
    }

    // Handle y == 0, x != 0
    if (y == 0.0f) {
        return (x > 0.0f)? 0.0f : PI;
    }

    // Handle y == x
    if (y == x) {
        return (x > 0.0f)? PI_4 : -(3.0f * PI_4);
    }

    // Handle y == -x
    if (y == -x) {
        return (x > 0.0f)? -PI_4 : (3.0f * PI_4);
    }

    // For other cases, determine quadrant and do appropriate lookup and calculation
    bool right = (x > 0.0f);
    bool top = (y > 0.0f);
    if (right && top) {
        // First quadrant
        if (y < x) {
            return AtanLookup2(y, x);
        }
        else {
            return PI_2 - AtanLookup2(x, y);
        }
    }
    else if (!right && top) {
        // Second quadrant
        const float posx = fabsf(x);
        if (y < posx) {
            return PI - AtanLookup2(y, posx);
        }
        else {
            return PI_2 + AtanLookup2(posx, y);
        }
    }
    else if (!right && !top) {
        // Third quadrant
        const float posx = fabsf(x);
        const float posy = fabsf(y);
        if (posy < posx) {
            return -PI + AtanLookup2(posy, posx);
        }
        else {
            return -PI_2 - AtanLookup2(posx, posy);
        }
    }
    else { // right && !top
        // Fourth quadrant
        const float posy = fabsf(y);
        if (posy < x) {
            return -AtanLookup2(posy, x);
        }
        else {
            return -PI_2 + AtanLookup2(x, posy);
        }
    }

    return 0.0f;
}

Answer 1

内置函数已经优化得很好，所以要打败它们真的很难。就个人而言，我会寻找其他地方来获得性能。

也就是说，我可以在您的代码中看到一项优化：

// Normalize angle so that it is in the interval (0.0, PI)
while (angle >= TWO_PI) {
    angle -= TWO_PI;
}

可以替换为：

angle = fmod(angle, TWO_PI);

Answer 2

请记住，您有一个协处理器...如果是 1993 年，您会看到速度有所提高...但是今天您将难以击败本机内在函数。

尝试查看 disassebly to sinf。

Answer 3

有人已经对此进行了基准测试，看起来Trig.Math 函数似乎已经过优化，并且比您能想出的任何查找表都要快：

http://www.tommti-systems.de/go.html?http://www.tommti-systems.de/main-Dateien/reviews/languages/benchmarks.html

（他们没有在页面上使用锚点，因此您必须向下滚动大约 1/3）

Answer 4

“过早的优化是万恶之源”——Donald Knuth

现在的编译器为三角函数提供了非常有效的内在函数，可以充分利用现代处理器（SSE 等），这解释了为什么您几乎无法击败内置函数。不要在这些部分上浪费太多时间，而是专注于您可以使用分析器发现的真正瓶颈。

Answer 5

我很担心这个地方：

// Normalize angle so that it is in the interval (0.0, PI)
while (angle >= TWO_PI) {
    angle -= TWO_PI;
}

但您可以： 为所有功能添加计时器，编写特殊性能测试，运行性能测试，打印时间测试报告。我想你会在这个测试之后知道答案。

您还可以使用一些分析工具，例如 AQTime。