FAST get 10 to the power n(10^n) in C/C++

Question

I want to calculate 10(yes only 10) to the power n[0..308] fast. i came up with some methods.

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;
}

time: O(logn), could it be faster? ( // f1() can do a little bit optimization but still O(logn))

2)

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

time: O(1), very good. But space: 309 * 8 bytes = 2472 bytes.. whoops it's too huge...

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 combines f1 and f2 to avoid multiplication such as 1e128*1e128, i wished it faster, but.. actually f3 is slower than f2.. because of ++i i guess..

well i almost gave up before i typed these codes,

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

and compiled it to .s by g++

LCPI0_0:
    .quad   0x4059000000000000              ## double 100

HOW does the compiler know 1e+2 is 0x4059000000000000?

i mean all i want to get is a double values 1e+n. but when a compiler compiles "double d = 1e+2", it know d should be 0x4059000000000000. can i just use some method to directly return something like 1e+n. OR can i do some things beyond C/C++ to get my value?

thanks a lot. Please point out if there is something wrong or unclear.

Answer 1

This is fast pow10-function in C++20:

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;
};

Answer 2

HOW does the compiler know 1e+2 is 0x4059000000000000?

Because 1e+2 is a literal (so a compile-time constant) and the compiler know the target architecture. Thus it can directly store the constant in the target program. Note that it can compute constants like pow(10.0, n) if it can deduce the value of n at compile time and optimizations are enabled. Because, n is likely a variable, 10e+n needs to be computed at runtime (eg. using pow(10.0, n) like you did for example). If you know that the input n is always a compile-time constant, then you can use constexpr (or templates).

f3 combines f1 and f2 to avoid multiplication such as 1e128*1e128, i wished it faster, but.. actually f3 is slower than f2.. because of ++i i guess..

No, it is a bit more complex. f3 is not very fast because of the loop-carried dependency that prevent your processor to execute the loop efficiently. Indeed, modern processors are superscalar and they can execute many loop iterations in parallel as long as there is no dependency between the instructions. In this case, the bottleneck come from the (fully sequential) modification of res in each iteration combined with a high latency of the associated floating-point instruction.

Note that the predictability of the conditionals and loops also plays a significant role.

I want to calculate 10(yes only 10) to the power n[0..308] fast

f2 can be very fast if the table is already stored in the L1 cache (and so can fit inside). Almost all mainstream modern processors have at least 16 KiB of cache nowadays. f2 can be slower than f3 if the table is stored in RAM and needs to be retrieved to to the cache (causing a very slow cache miss due to the high latency of the RAM).

You can speed up the computation of f3 by a large margin using a manual reduction that can be easily unrolled by the compiler. Here is an example of code:

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];
}

The compiler Clang produces relatively good instructions for f4 resulting in a much faster execution. Unfortunately, some compilers like GCC generate a quite bad code for f4: they use slow conditional jumps. Still, this reduction-based approach is far from being cheap.

As proposed by @MarcGlisse in the comments, you can do 2 lookups on very-small tables so to get a very fast implementation:

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];
}