C++ log2 使用位黑客实现双精度？

Question

我知道我可以在 C++ 中使用 cmath 中的 log2，但我想实现一种用于学习目的的有效方法。

我看到了很多相关的问题，但没有一个给出快速实现，将双精度浮点作为输入并将输入的 log2 近似值输出为双精度。

我想使用位运算来做到这一点，我从这个answer移植了代码：

#include <cmath>
#include <vector>

const double ln2 = log(2);

vector<float> log2_float() {
    int lim = 1 << 24;
    vector<float> table(lim);
    for (int i = 0; i < lim; i++) {
        table[i] = float(log(i) / ln2) - 150;
    }
    return table;
}
const vector<float> LOG2_FLOAT = log2_float();

float fast_log2(float f) {
    uint32_t bits = *reinterpret_cast<uint32_t*>(&f);
    int e = (bits >> 23) & 0xff;
    int m = bits & 0x7fffff;
    return (e == 0 ? LOG2_FLOAT[m << 1] : e + LOG2_FLOAT[m | 0x00800000]);
}

我不会向您展示完整的代码，其余的无关紧要。我已经针对 log2 对我的实现进行了基准测试，发现我的实现快了一个数量级，但我的代码仅适用于浮点数：

fast_log2(3.1415927f) = 1.651497 : 2.222061 nanoseconds
log2(3.1415927f) = 1.651496 : 33.1502 nanoseconds

我实际上并没有使用 pi 进行基准测试，我使用了存储在向量中的随机值，但迭代次数是相同的，并且它们使用来自同一向量的值。

我想知道，如何对 double 执行此操作并获得 10 个正确的小数位？我知道我不能使用查找表，这将需要比我更多的内存。放弃查找表并仍然使用相同的方法将需要调用日志，这会再次减慢代码速度。

我计算出float有24个小数位，所以最大正确的小数位数是7.224719895935548，对于每个正确的小数位数，需要1 / log₁₀(2) = 3.321928094887362位，才能得到8个正确的小数位数，26.575424759098897位是必需的，所以我需要一个包含 2²⁷ 值的表，这将需要 2³⁰ 字节或 1GiB，而我只有 16GiB RAM...

我可以使用什么算法来实现这一目标？

---

原来我的计算是错误的，我只需要 1GiB RAM。当然我实际上并没有创建那么大的表，我只是重新计算了一下。我不是数学家，也没有花太多功夫去计算。

无论如何，我正在做其他事情，几分钟前我尝试了一下并创建了一个双打查找表。为了获得正确的小数点后 7 位，我使用了相同数量的元素，浮点向量需要 32MiB，双精度向量需要 64MiB。

我只是更改了数字，无论出于何种原因它编译成功但我未能得到结果，程序在运行该函数时崩溃了：

vector<double> log2_double() {
    int lim = 1 << 24;
    vector<double> table(lim);
    for (uint64_t i = 0; i < lim; i++) {
        table[i] = log(i << 29) / ln2 - 1075;
    }
    return table;
}
const vector<double> LOG2_DOUBLE = log2_double();

double fast_log2_double(double d) {
    uint64_t bits = *reinterpret_cast<uint64_t*>(&d);
    uint64_t e = (bits >> 52) & 0xfff;
    uint64_t m = bits & 0xfffffffffffff;
    return (e == 0 ? LOG2_DOUBLE[m >> 28] : e + LOG2_DOUBLE[m | 0x20000000000000]);
}

上面的代码有什么问题？

我的目标是获得 10 个正确的十进制数字，使用相同的方法，我需要一个包含 2³⁴ 双精度数的查找表，这需要 128GiB RAM，这太离谱了。

我想要一种快速逼近方法，保证至少 9 位准确的小数位，有时 10 位或更多。

---

我在调试模式下编译了程序，当程序运行时，我得到了类似

IndexingError

的信息，我很惊讶，因为我读到C++更有可能给出一些随机数。

所以我再次查看了链接的答案，以及关于 DPFP 的维基百科文章并再次计算了值，我修复了代码：

vector<double> log2_double() {
    int lim = 1 << 24;
    vector<double> table(lim);
    for (uint64_t i = 0; i < lim; i++) {
        table[i] = log(i << 29) / ln2 - 1075;
    }
    return table;
}
const vector<double> LOG2_DOUBLE = log2_double();

double fast_log2_double(double d) {
    uint64_t bits = *reinterpret_cast<uint64_t*>(&d);
    uint64_t e = (bits >> 52) & 0x7ff;
    uint64_t m = bits & 0xfffffffffffff;
    return (e == 0 ? LOG2_DOUBLE[m >> 28] : e + LOG2_DOUBLE[(m | 0x10000000000000) >> 29]);
}

这次成功了，但正如预期的那样，我只得到了 6 个正确的小数位。我想要至少 9 位正确的小数位，最好是 10 位。我实际上希望使用 a^x = e^{x ln a} 进行小数求幂。因为浮点数的编码方式，我想直接从二进制文件中提取 log₂ a 部分，并使用 ln a = log₂ a * ln 2.

现在我正在考虑在区间 [1, 2] 上使用 e^x 的 Remez 近似，因为二进制编码值中尾数是 1，所以我可以在 Remez 近似和位移位中使用它并做一些基本的操作数学到指数，然后连接两个二进制字符串并转换为 float/double。

我正在努力。

Answer 1

您可以应用与您为 float 展示的方法类似的方法来构造双精度浮点整数的有效 log2 近似值。另一方面，在没有查找表的情况下获得 10 位精确的十进制数字的精度可能很困难，并且可能需要更复杂的近似方法。

多项式近似是一种常用策略。在给定范围内，可以使用多项式来近似 log2 函数。例如，您可以利用泰勒级数展开：

log2(x) = (x - 1) - (1/2)(x - 1)^2 + (1/3)(x - 1)^3 - (1/4)(x - 1)^4 + ...

要实现多项式近似，请将级数截断为给定的项数。近似值越好，提供的项就越多。下面是如何在 C++ 中实现双精度的示例：

#include <iostream>
#include <cmath>

double fast_log2(double x) {
    if (x <= 0.0) {
        // Handle special case for x <= 0
        return -std::numeric_limits<double>::infinity();
    }

    // Ensure x is in the range (0.5, 1]
    if (x < 0.5) {
        x = 1.0 / x;
        return -fast_log2(x);
    }

    // Calculate log2 using a polynomial approximation
    double y = (x - 1.0) / (x + 1.0);
    double result = 0.0;
    double term = y;
    for (int i = 1; i <= 100; i += 2) {
        result += term / i;
        term *= y * y;
    }

    return 2.0 * result;
}

int main() {
    double x = 3.1415927;
    double result = fast_log2(x);
    std::cout << "log2(" << x << ") = " << result << std::endl;
    return 0;
}

对于双精度输入，此方法应该提供足够准确的 log2 近似值，而不需要查找表。近似的精度由级数中的项数决定（在本例中，最多 100）。

根据您的独特需求，您可以改变短语的数量，以在准确性和性能之间取得折衷。更多术语会提高精度，但会减慢过程。