在 C/C++ 中快速求 10 的 n(10^n) 次方

Question

我想快速计算 10（是的，只有 10）的 n[0..308] 次方。我想出了一些方法。

1)

double f(int n) {
  return pow(10.0, n);
}
double f1(int n) {
  double a = 10.0;
  double res = 1.0;
  while(n) {
    if(n&1) res *= a;
    a *= a;
    n >>= 1;
  }
  return res;
}

时间：O(logn)，还能更快吗？ ( // f1() 可以做一点优化，但仍然是 O(logn))

2)

double f2(int n) {
  static const double e[] = { 1e+0, 1e+1, 1e+2, ..., 1e+308 };
  return e[n];
}

时间：O(1)，非常好。 但是空间：309 * 8 字节 = 2472 字节..哎哟，它太大了...

3)

double f3(int n){
    static const double e[] = {
        1e+1, 1e+2, 1e+4, 1e+8, 1e+16, 1e+32, 1e+64, 1e+128, 1e+256
    };
    double res = 1.0;
    for(int i = 0; n; ++i){
        if(n & 1){
            res *= e[i];
        }
        n >>= 1;
    }
    return res;
}

f3 结合了 f1 和 f2 以避免乘法，例如 1e128*1e128，我希望它更快，但是..实际上 f3 比 f2 慢..因为我猜是 ++i..

好吧，在我输入这些代码之前我几乎放弃了，

int main(){
    double d = 1e+2;
    return 0;
}

并通过g++将其编译为.s

LCPI0_0:
    .quad   0x4059000000000000              ## double 100

编译器如何知道 1e+2 是 0x4059000000000000？

我的意思是我想要得到的只是一个双精度值 1e+n。但是当编译器编译“double d = 1e+2”时，它知道d应该是0x4059000000000000。我可以使用某种方法直接返回 1e+n 这样的东西吗？或者我可以做一些超越 C/C++ 的事情来获得我的价值吗？

非常感谢。如有错误或不清楚的地方请指出。

Answer 1

编译器如何知道 1e+2 是 0x4059000000000000？

因为

1e+2

是一个文字（因此是一个编译时常量）并且编译器知道目标体系结构。这样就可以将常量直接存储到目标程序中。请注意，如果它可以在编译时推导出

pow(10.0, n)

的值并且启用优化，则它可以计算像

n

这样的常量。因为，

n

可能是一个变量，所以需要在运行时计算10e+n（例如，像您一样使用

pow(10.0, n)

）。如果您知道输入

n

始终是编译时常量，那么您可以使用

constexpr

（或模板）。

f3 结合了 f1 和 f2 以避免乘法，例如 1e128*1e128，我希望它更快，但是..实际上 f3 比 f2 慢..因为我猜是 ++i..

不，它有点复杂。

f3

速度不是很快，因为循环携带的依赖性会阻止处理器有效地执行循环。事实上，现代处理器是超标量，只要指令之间不存在依赖性，它们就可以并行执行许多循环迭代。在这种情况下，瓶颈来自每次迭代中对

res

的（完全顺序）修改以及相关浮点指令的高延迟。

请注意，条件和循环的可预测性也发挥着重要作用。

我想快速计算 10（是的，只有 10）的 n[0..308] 次方

如果表已经存储在 L1 缓存中（因此可以放入其中），

f2

会非常快。现在几乎所有主流现代处理器都至少有 16 KiB 的缓存。如果表存储在 RAM 中并且需要检索到缓存，

f2

可能会比

f3

慢（由于 RAM 的高延迟，导致非常慢的 缓存未命中）。

您可以使用手动

reduction

大幅加快 f3 的计算速度，编译器可以轻松地展开。这是代码示例：

double f4(uint32_t n) {
    static const double e[] = {
        1e+1, 1e+2, 1e+4, 1e+8, 1e+16, 1e+32, 1e+64, 1e+128, 1e+256
    };

    double p[9];

    for(int i=0 ; i<9 ; ++i)
        p[i] = (n & (1 << i)) ? e[i] : 1.0;

    return ((p[0] * p[4]) * (p[1] * p[5])) * ((p[2] * p[6]) * (p[3] * p[7])) * p[8];
}

编译器 Clang 为

f4

生成相对较好的指令，从而导致执行速度更快。不幸的是，一些编译器（例如 GCC）会为

f4

生成非常糟糕的代码：它们使用缓慢的条件跳转。尽管如此，这种基于减少的方法并不便宜。

正如 @MarcGlisse 在评论中所建议的，您可以在非常小的表上进行 2 次查找，以便获得非常快速的实现：

double f6(uint32_t n) {
    static const double eLow[] = { 1e+0, 1e+1, 1e+2, 1e+3, 1e+4, 1e+5, 1e+6, 1e+7, 1e+8, 1e+9, 1e+10, 1e+11, 1e+12, 1e+13, 1e+14, 1e+15, 1e+16, 1e+17, 1e+18, 1e+19, 1e+20, 1e+21, 1e+22, 1e+23, 1e+24, 1e+25, 1e+26, 1e+27, 1e+28, 1e+29, 1e+30, 1e+31 };
    static const double eHigh[] = { 1e+0, 1e+32, 1e+64, 1e+96, 1e+128, 1e+160, 1e+192, 1e+224, 1e+256, 1e+288 };
    return eLow[n & 0x1F] * eHigh[n >> 5];
}

Answer 2

这是 C++20 中的快速 pow10 函数：

double pow10( int64_t exp )
{
    constexpr uint64_t EXP_MASK = 0x7FFull << 52;
    // table for binary exponentation with 10 ^ (2 ^ N)
    static array<double, 64> tenPows;
    // table initialized ?
    if( static atomic_bool once( false ); !once.load( memory_order_acquire ) )
    {
        static mutex mtx;
        // lock mtx and check once-flag again
        lock_guard lock( mtx );
        if( !once.load( memory_order_relaxed ) )
        {
            // strongly no: calculate table
            for( double p10x2xN = 10.0; double &pow : tenPows )
                pow = p10x2xN,
                p10x2xN *= p10x2xN;
            // set initialized flag with release semantics
            once.store( true, memory_order_release );
        }
    }
    // begin with 1.0 since x ^ 0 = 1
    double result = 1.0;
    // unsigned exponent
    uint64_t uExp = exp >= 0 ? exp : -exp;
    // iterator to highest exponent
    auto itExp = tenPows.rbegin();
    // as long as there are bits set in uExp
    for( size_t gap ; uExp; uExp <<= gap, uExp <<= 1 )
    {
        // exponent bits gap
        gap = countl_zero( uExp );
        // get next exponent
        itExp += gap;
        // multiply result by next pow10 exponent
        result *= *itExp++;
        // overlow / underflow ?
        if( (bit_cast<uint64_t>( result ) & EXP_MASK) == EXP_MASK ) [[unlikely]]
            // yes: result wouldn't change furhter; stop
            return exp >= 0 ? result : 0.0;
    }
    return exp >= 0 ? result : 1.0 / result;
};