Representation of the formula P(n) = (n!)(6^n) in C I'm trying to represent the following mathematical expression in C:

$$P(n)=(n!)(6^n)$$

The first program should compute the answer to expression when n = 20. I have attempted to create the program in C and it produces the correct answer (when compared with wolfram alpha's answer). The answer is approximately \$10^{34}\$.

Is there any way that I could improve the code?:

#include &lt;math.h&gt;
typedef unsigned int uint;

void main()
{
    uint n=20;
    double F=1,pwr,P;

    pwr=pow(6,n);

    while(n&gt;0){
        F=F*n;
        n=n-1;
    }
    P=F*pwr;
}
Because the result is on the order of \$10^{34}\$, as you have correctly noted, even a double will not be able to accurately represent all digits of such a large result.  We know this because it would take 
$$\frac{34}{\log_{10}2} \approx 113 \text{ bits}$$
to represent such a result, and a double isn't that large on any architecture you're likely to be using.

This then requires that we resort to multiple precision math.  Just as with multiple-digit multiplication you learned when you were very young, multiple-precision multiplication is the same principle: multiply digit by digit, adding the carry from the previous result.  In this sample code, rather than having a "digit" be 10, I've chosen 1000 which results in fewer math operations.

// calculate P(n) = (n!)(6^n)
#include &lt;stdio.h&gt;

typedef unsigned int uint;

const uint OVERFLOWVALUE = 1000;

void multiply(uint *a, unsigned adigits, uint b)
{
    uint carry = 0;
    for ( ; adigits; --adigits, ++a) {
        *a = *a * b + carry;
        if (*a &gt;= OVERFLOWVALUE) {
            carry = *a / OVERFLOWVALUE;
            *a = *a % OVERFLOWVALUE;
        } else {
            carry = 0;
        }
    }
}

void print(uint *P, unsigned digits) 
{
    // print results
    for (unsigned i = digits; i; --i)
        printf("%3.3u ", P[i-1]);
    putchar('\n');
}

int main()
{
    const unsigned digits = 12;
    uint P[digits]; 
    // initialize P
    for (unsigned i = 0; i &lt; digits; ++i)
        P[i] = 0;

    uint n=20;

    printf("{n!*n^6, n=%u} = ", n);
    for (P[0] = 1; n; --n) {
        multiply(P, digits, n*6);
    }
    print(P, digits);
}


The result on this machine:

{n!*n^6, n=20} = 008 895 075 211 041 185 783 708 532 080 640 000 


This agrees precisely with the Wolfram Alpha result.

More mathematical explanation

The multiply routine multiplies a multiple-precision result pointed to by a with adigits values by single value b.  It does not check for mathematical overflow in its calculations, so we need to be careful about the range.  

For calculation purpose, suppose that a uint type has 32 bits.  If we assume that a is well-formed before the calculation, how big can b be without risking numerical overflow?  In this case, the largest number that can be in any "digit" is 1000, so overflow will happen whenever \$b * 1000 &gt; 2^{32}\$.  This implies that to avoid overflow, $$b &lt; \frac{2^{32}}{1000}$$ $$b &lt;= 4294967$$.  With this maximum possible value, the carry value would be 4294, which is greater than 1000, so to be pedantic, the code should actually check each digit before each multiplication to normalize the value. 

Without this precaution, however, the calculation is still free from overflow as long as \$(999+6n)6n &lt; 2^{32}\$.  This is a simple quadratic equation, with the largest positive root being just over 10839, so as long as \$n&lt;10839\$ this code will work correctly as is.