Calculating n factorial for large n ~ 100, but result is incorrect

Question

I am trying to solve the problem using an array where I put the digits of n!

This is what I tried

int main()
{
    int a[1000];
    int n = 1;
    a[0] = 1;
    int limit = 100;
    for (int i = 2; i <= limit; i++) {
        int carry = 0;
        for (int j = 0; j < n; j++) {
            int prod = a[j] * i + carry;
            a[j] = prod % 10;
            carry = prod / 10;
            if (carry) {
                a[j + 1] = carry;
                n = n + 1;
            }
        }
    }
    for (int i = (n - 1); i >= 0; i--)
        printf("%d", a[i]);

    return 0;
}

It works fine for 2! and 3! but for larger value of limit it just prints many random numbers.

For example, the output for 4! is 504.
The expected output is 24.

Answer 1

The problem is in the way you handle carry: you should store the excess digits after finishing the loop over the existing digits. There can be more than 1 extra digit to store so you need a loop.

Also note that you do not need int for the type of the a array, char or unsigned char is sufficient.

Here is a modified version that can handle numbers up to 3248 (less than 10001 digits).

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

int main(int argc, char *argv[]) {
    int limit = (argc > 1) ? atoi(argv[1]) : 100;
    unsigned char a[10000];
    int len = sizeof a / sizeof a[0];
    int n = 1;
    a[0] = 1;
    for (int i = 2; i <= limit; i++) {
        int carry = 0;
        // compute the digits using high school arithmetics
        for (int j = 0; j < n; j++) {
            carry += a[j] * i;
            a[j] = carry % 10;
            carry = carry / 10;
        }
        // store the extra digits
        while (carry > 0) {
            if (n >= len) {
                printf("too many digits\n");
                return 1;
            }
            a[n] = carry % 10;
            carry = carry / 10;
            n = n + 1;
        }
    }
    for (int i = n; i-- > 0;) {
        printf("%d", a[i]);
    }
    printf("\n");
    return 0;
}

Output:

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

For larger factorials, you can use the Stirling formula to approximate the number of decimal digits:

The number of decimal digits comes out as

number_of_digits(n!) ≈ log10(2π)/2 - n.log10(e) + (n+0.5).log10(n)

Here is a modified version, computing 9 digits at a time:

#include <inttypes.h>
#include <math.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>

// approximating the number of digits:
// log10 n! = log10(2*pi)/2 - n*log10(e) + (n+0.5)*log10(n)

int main(int argc, char *argv[]) {
    int limit = (argc > 1) ? atoi(argv[1]) : 100;
    size_t len = (size_t)((limit * 0.434 + (limit + 0.5) * log((double)limit) + 17) / 9);
    uint32_t *a = malloc(len * sizeof(*a));
    if (a == NULL) {
        fprintf(stderr, "cannot allocate %zu bytes\n", len * sizeof(*a));
        return 1;
    }
    size_t n = 1;
    a[0] = 1;
    uint64_t factor = 1;
    for (int i = 2;; i++) {
        if (i > limit || factor * i > UINT32_MAX) {
            uint64_t carry = 0;
            // compute the digits using high school arithmetics
            for (size_t j = 0; j < n; j++) {
                carry += (uint64_t)a[j] * (uint32_t)factor;
                a[j] = carry % 1000000000;
                carry = carry / 1000000000;
            }
            // store the extra digits
            while (carry > 0) {
                if (n >= len) {
                    fprintf(stderr, "too many digits %zu at %d\n", n * 9, i);
                    return 1;
                }
                a[n] = carry % 1000000000;
                carry = carry / 1000000000;
                n = n + 1;
            }
            factor = 1;
            if (i > limit)
                break;
        }
        factor *= i;
    }
    if (n > 1) {
        uint64_t num = a[n - 1] * 1000000000ULL + a[n - 2];
        int exp = (n - 1) * 9;
        // should perform proper rounding
        while (num >= 10000000000) {
            num /= 10;
            exp += 1;
        }
        printf("%d.%09"PRIu64"e%d\n", (int)(num / 1000000000), num % 1000000000, exp);
    } else {
        printf("%"PRIu32"\n", a[n - 1]);
    }
#if 0
    printf("%"PRIu32, a[n - 1]);
    if (n > 1) {
        for (size_t i = n - 1; i-- > 0;) {
            printf("%09"PRIu32, a[i]);
        }
    }
    putchar('\n');
#endif
    free(a);
    return 0;
}

It can compute all digits of large factorials (as can be verified on the Factorial wikipedia page), but its time complexity could be improved using more advanced multiplications techniques, and binary computation leaving base 10 conversion for printing only.

This version is still quite simple and computes the 5565709 digits of 1000000! in 16mns on my M2 macbook air.