How can I calculate 2^n for large n?

The problem is not to compute 2 to a high power, but to convert this number to a decimal representation:

  • Let's represent large numbers with arrays of unsigned 32-bit integers.
  • Computing 2n is as easy as setting a single bit.
  • Converting to binary can be performed by repeatedly dividing this number by 1000000000, producing 9 digits at a time.

Here is a simple but fast implementation:

#include <stdint.h>
#include <stdio.h>

void print_2_pow_n(int n) {
    int i, j, blen = n / 32 + 1, dlen = n / 29 + 1;
    uint32_t bin[blen], dec[dlen];
    uint64_t num;

    for (i = 0; i < blen; i++)
        bin[i] = 0;
    bin[n / 32] = (uint32_t)1 << (n % 32);

    for (j = 0; blen > 0; ) {
        for (num = 0, i = blen; i-- > 0;) {
            num = (num << 32) | bin[i];
            bin[i] = num / 1000000000;
            num = num % 1000000000;
        }
        dec[j++] = (uint32_t)num;
        while (blen > 0 && bin[blen - 1] == 0)
            blen--;
    }
    printf("2^%d = %u", n, dec[--j]);
    while (j-- > 0)
        printf("%09u", dec[j]);
    printf("\n");
}

int main() {
    int i;
    for (i = 0; i <= 100; i += 5)
        print_2_pow_n(i);
    print_2_pow_n(1000);
    print_2_pow_n(10000);
    print_2_pow_n(100000);
    return 0;
}

Output:

2^0 = 1
2^5 = 32
2^10 = 1024
2^15 = 32768
2^20 = 1048576
2^25 = 33554432
2^30 = 1073741824
2^35 = 34359738368
2^40 = 1099511627776
2^45 = 35184372088832
2^50 = 1125899906842624
2^55 = 36028797018963968
2^60 = 1152921504606846976
2^65 = 36893488147419103232
2^70 = 1180591620717411303424
2^75 = 37778931862957161709568
2^80 = 1208925819614629174706176
2^85 = 38685626227668133590597632
2^90 = 1237940039285380274899124224
2^95 = 39614081257132168796771975168
2^100 = 1267650600228229401496703205376
2^1000 = 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376
2^10000 = 1995063116880758384883742<...>91511681774304792596709376
2^100000 = 9990020930143845079440327<...>97025155304734389883109376

2100000 has 30103 digits, which is exactly floor(100000 * log10(2)). It executes in 33 milliseconds on my old laptop.


Simply make a bit array and set the nth-bit. Then divide by 10 as if the bit array were a little-endian number and print the remainders in reverse to get the base-10 representation of your nth-power of two.

This quick program below does it and it's giving me the same results as bc, so I guess it works. The printing routine could use some tuning.

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>

uint_least32_t div32(size_t N, uint_least32_t Z[/*N*/], uint_least32_t X[/*N*/], uint_least32_t Y)
{
    uint_least64_t carry; size_t i;
    for(carry=0, i = N-1; i!=-1; i--)
        carry = (carry << 32) + X[i], Z[i] = carry/Y, carry %= Y;
    return carry;
}

void pr10(uint_least32_t *X, size_t N)
{
    /*very quick and dirty; based on recursion*/
    uint_least32_t rem=0;
    if(!X[N?N-1:0]) return;
    rem = div32(N,X,X,10);
    while(N && !X[N-1]) N--;
    pr10(X,N);
    putchar(rem+'0');
}
int main(int C, char **V)
{
    uint_least32_t exp = atoi(V[1]);
    size_t nrcells = exp/32+1;
    uint_least32_t *pow  = calloc(sizeof(uint_least32_t),nrcells);
    if(!pow) return perror(0),1;
    else pow[exp/32] = UINT32_C(1)<<(exp%32);
    pr10(pow,nrcells);

}

Example run:

$ ./a.out 100
1267650600228229401496703205376

Step 1: Decide how you're going to represent bignums

There exist libraries out there already for this. The GNU Multiple Precision Integer library is a commonly used option. (But according to your edit, that isn't an option. You might still glance at some of them to see how they do things, but it isn't necessary.)

If you want to roll your own, I do not recommend storing the decimal digits. If you do that, you'll need to convert to and from a binary representation every time you want to do arithmetic on the components. Better to have something like a linked list of uint32_ts, along with a sign bit. You can convert from/to decimal when you want to read and write, but do your math in binary.

Step 2: Implement exponentiation

I'll be assuming the linked list bignum implementation here; you can adapt the algorithms as needed.

If you're just calculating a power of 2, it's easy. It's a 1 followed by N 0s, so if each block stores M bits and you want to represent 2^N, then just have floor(N/M) blocks of all 0s, and store 1 << (N % M) in the most significant block.

If you want to be able to do exponentiation with arbitrary bases in an efficient manner, you should use exponentiation by squaring. The idea behind this is that if you want to compute 3^20, you don't multiply 3 * 3 * 3 * ... * 3. Rather, you compute 3^2 = 3 * 3. Then 3^4 = 3^2 * 3^2. 3^8 = 3^4 * 3^4. 3^16 = 3^8 * 3^8. And you store each of these intermediate results as you go. Then once you reach the point where squaring it again would result in a larger number than the one you want, you stop squaring and assemble the final result from the pieces you have. In this case, 3^20 = 3^16 * 3^4.

This approach computes the final result in 5 steps instead of 20, and since the time is logarithmic in terms of the exponent, the speed gain gets more pronounced the larger the exponent. Even computing 3^100000 only takes 21 multiplications.

There isn't a clever approach to the multiplication that I know of; you can probably just do something along the lines of the basic long multiplication algorithm you learned in elementary school, but at the level of blocks: the reason we used uint32_ts earlier instead of uint64_t`s is so we can cast the operands to the larger type and multiply them without risk of losing the carry bits to overflow.

Convert from binary to decimal for printing

First, find the largest multiple of 10 less than your number.
I leave doing this efficiently as an exercise for the reader, but you can probably manage it by doing exponentiation by squaring to find an upper bound, then subtracting various stored intermediate values to get down to the actual value faster than you would by dividing by 10 repeatedly.

Or you can just find the number by repeatedly multiplying by 10; the rest of this is going to be linear no matter how the first part is handled.

But however you get it, you have a q such that q = k * 10, 10 * q > n, q <= n, you can just cycle through one decimal digit at a time:

for (; q; q /= 10) {
   int digit = n / q; //truncated down to floor(n/q)
   printf("%d", digit);
   n -= digit * q;
}

It's possible that there's a more efficient method in the literature somewhere, but I'm not familiar with one offhand. But it's not a major deal so long as we only need to do the inefficient part when writing output; that's slow no matter the algorithm. By which I mean, it might take a millisecond or two to print all 100,000 digits. That doesn't matter when we're displaying the number for human consumption, but if we had to wait a millisecond as part of a calculation in a loop somewhere, it'd add up and become terribly inefficient. That's why we never store numbers in a decimal representation: by representing it as binary internally, we do the inefficient parts once on input and once on output, but everything in between is fast.

Tags:

C