How can I calculate 2^n for large n?

Question

I am trying to write a program that takes a number, n, as input, and outputs the result of 2 to the power of n. The problem is, n can be very large (up to 100,000). Essentially, I am trying to calculate pow(2, n); for very large numbers.

I think the way to do this is to store the digits in an array, since there is no built-in numeric type that can hold values that are this large.

The digits are in decimal format (base 10).

I am using C, not C++, so I cannot use STL vectors and other C++ containers. I also cannot use external libraries, like GMP. I need to implement the algorithm manually in pure C.

Answer 1

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 2ⁿ 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

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

Answer 2

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

Answer 3

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.

Answer 4

This is a rather naive and inefficient solution. As asked for, the numbers are represented in an array of decimal digits. We compute the exponential 2ⁿ by repeated addition of the number 2 with itself: start with e := 2 and repeat e := e + e n times.

To come up with an upper bound for the length of our digits array we use the following approach:

• The representation of number x in base b has ⌈log_b(x)⌉ digits.
• Comparing the number of digits between the binary and decimal representation of any number x we find that they are only off by a constant factor if we neglect the rounding (⌈⌉).
  log₂(x) / log₁₀(x) = 1 / log₁₀(2) = 3.3219... > 3
• 2ⁿ has log₂(2ⁿ) = n binary digits.
• Therefore 2ⁿ has around n / 3 decimal digits. Because of the rounding issues we add +1 to that.

void print(int digits[], int length) {
    for (int i = length - 1; i >= 0; --i)
        printf("%d", digits[i]);
    printf("\n");
}
void times2(int digits[], int length) {
    int carry = 0;
    for (int i = 0; i < length; ++i) {
        int d = 2 * digits[i] + carry;
        digits[i] = d % 10;
        carry = d / 10;
    }
}
int lengthOfPow2(int exponent) {
    return exponent / 3 + 1;
}
// works only for epxonents > 0
void pow2(int digits[], int length, int exponent) {
    memset(digits, 0, sizeof(int) * length);
    digits[0] = 2;
    for (int i = 1; i < exponent; ++i)
        times2(digits, length);
}
int main() {
    int n = 100000;
    int length = lengthOfPow2(n);
    int digits[length];
    pow2(digits, length, n);
    print(digits, length);
    return 0;
}

On unix-like systems you can check correctness for a fixed n using

diff \
  <(compiledProgram | sed 's/^0*//' | tr -d '\n') \
  <(bc <<< '2^100000' | tr -d '\n\\')

As already pointed out, this solution is not very efficient. Compiled with clang -O2 computing 2^100'000 took 8 seconds on an Intel i5-4570 (3,2GHz).

The next step to speed this up would be to repeatedly cube your number instead of repeatedly multiplying by 2. Even with a naive implementation of the cube step this should be faster than the implementation presented in this answer.

If you need to be even more efficient you can implement the cube step using something like Karatsuba's algorithm or even fast fourier transformation (FFT). With the cubing approach and FFT you can compute 2ⁿ in around O(n·log(n)) (there may be an additional log(log(n)) factor due to rounding issues in FFT).

Answer 5

I wasn't able to find a solution in logarithmic complexity(exponentiation by squaring) but I did manage to code a naive implementation with time complexity O(noOfDigits*pow), noOfDigits in 2^n will be n*log10(2)+1;

I only checked the answer with the first few digits of https://www.mathsisfun.com/calculator-precision.html and it seems to be correct.

#include <stdio.h>
#include <math.h>
//MAX is no of digits in 2^1000000
#define MAX 30103
int a[MAX];
int n;
void ipow(int base, int exp,int maxdigits)
{
    a[0]=1;
    for (;exp>0;exp--){
            int b=0;
            for(int i=0;i<maxdigits;i++){
                a[i]*=base;
                a[i]+=b;
                b=a[i]/10;
                a[i]%=10;
            }
    }
}
int main()
{
    int base=2;
    int pow=100000;
    n=log10(2)*pow+1;
    printf("Digits=%d\n",n);
    ipow(base,pow,n);
    for(int i=n-1;i>=0;i--){
        printf("%d",a[i]);
    }
    return 0;
}

I also wrote code for exponentiation by squaring but with an unoptimized multiplication function. Which seems to be faster than above implementation.

#define MAX 30103
int a[MAX];
int b[MAX];
int z[MAX];
//stores product in x[]; mul of large arrays implemented in n^2 complexity
//n and m are no of digits in x[] and y[]
//returns no of digits in product
int mul(int x[],int y[],int n,int m){
    for(int i=0;i<n+m;i++)
        z[i]=0;
    for(int j=0;j<m;j++){
        int c=0;
        for(int i=0;i<n+m;i++){
            z[i+j]+=x[i]*y[j];
            z[i+j]+=c;
            c=z[i+j]/10;
            z[i+j]%=10;
        }
    }
    for(int i=0;i<n+m;i++){
            x[i]=z[i];
    }
    if(x[n+m-1]==0)
        return n+m-1;
    return n+m;
}
//stores answer in x[]
int ipow(int base, int exp)
{
    int n=1,m=0;
    for(int i=0;base>0;i++){
        b[i]=base%10;
        base/=10;
        m++;
    }
    a[0]=1;
    for (;;)
    {
        if (exp & 1)
            n=mul(a,b,n,m);
        exp >>= 1;
        if (!exp)
            break;
        m=mul(b,b,m,m);
    }
}
int main()
{
    int base=2;
    int pow=100000;
    n=log10(2)*pow+1;
    printf("Digits=%d\n",n);
    ipow(base,pow);
    printf("\n");
    for(int i=n-1;i>=0;i--){
        printf("%d",a[i]);
    }
    return 0;
}

Answer 6

Since the original problem statement does not specify the output base, here is a joke implementation:

#include <stdio.h>

void print_2_pow_n(int n) {
    printf("2^%d = 0x%d%.*d\n", n, 1 << (n % 4), n / 4, 0);
}

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

Output:

2^0 = 0x1
2^1 = 0x2
2^2 = 0x4
2^3 = 0x8
2^4 = 0x10
2^5 = 0x20
2^6 = 0x40
2^7 = 0x80
2^8 = 0x100
2^9 = 0x200
2^10 = 0x400
2^11 = 0x800
2^12 = 0x1000
2^13 = 0x2000
2^14 = 0x4000
2^15 = 0x8000
2^100 = 0x10000000000000000000000000
2^100000 = 0x10...<0 repeated 24998 times>...0