How to output integers more than 16 digits correctly in C?

Question

The code here only gives correct output till the factorial of 21 and after that its correct upto 16 digits from the left then remaining digits are just given as zero. I tried changing the type of variable c from double to long double but it just gives errors or doesn't print the factorials.

#include <stdio.h>

FILE *fp;
long double facto(int i);
int main() {
    int n, i;
    double c;
    printf("enter no. to find factorial till:\n");
    scanf("%d", &n);
    fp = fopen("output_of_factorial.txt", "w");
    fputs("Number  |\t Factorial\n\n", fp);

    for (i = 1; i <= n; i++) {
        c = facto(i);
        fprintf(fp, "%d\t|\t %.0Lf\n", i, c);
    }
    fclose(fp);
    return 0;
}

long double facto(int x) {
    if (x == 1)
        return 1;
    else
        return x * facto(x - 1);
}

Answer 1

Tye double only has 53 bits of precision, long double probably has 80 bits on your platform. Using floating point arithmetics will give you an approximate result that is correct for the most significant digits. Using integers give you the exact result, but only if it is smaller than the range of the type.

You can use the type long long that is at least 64-bit wide, thus 19 digits, or for one more bit, type unsigned long long that allows for integers twice as large:

LLONG_MAX = 9223372036854775807  //  > 9.22e19
ULLONG_MAX = 18446744073709551615  // > 1.84e20

Here is a modified version of the code:

#include <stdio.h>

unsigned long long facto(int i);

int main(void) {
    int n, i;
    unsigned long long c;
    FILE *fp;

    printf("enter no. to find factorial till: ");
    if (scanf("%d", &n) == 1) {
        fp = fopen("output_of_factorial.txt", "w");
        if (fp != NULL) {
            fputs("Number |            Factorial\n\n", fp);
            for (i = 1; i <= n; i++) {
                c = facto(i);
                fprintf(fp, "%6d | %20llu\n", i, c);
            }
            fclose(fp);
        }
    }
    return 0;
}

unsigned long long facto(int x) {
    if (x <= 1)
        return 1;
    else
        return x * facto(x - 1);
}

It works all the way to 20:

Number |            Factorial

     1 |                    1
     2 |                    2
     3 |                    6
     4 |                   24
     5 |                  120
     6 |                  720
     7 |                 5040
     8 |                40320
     9 |               362880
    10 |              3628800
    11 |             39916800
    12 |            479001600
    13 |           6227020800
    14 |          87178291200
    15 |        1307674368000
    16 |       20922789888000
    17 |      355687428096000
    18 |     6402373705728000
    19 |   121645100408832000
    20 |  2432902008176640000

But fails for 21 and above because of arithmetic overflow.

To go further, you could use 128-bit integers if they are available on your platform (uint128_t, __uint128 or __uint128_t) but you would need to write your own conversion function to output the decimal representation.

A better approach would be to use multi-precision (aka bignum) packages that can handle extremely large numbers, typically only bound by available memory.

Answer 2

For illustration purposes, here is a simplistic implementation of bignum I just wrote with a printable internal as an ASCII string:

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

typedef struct bignum {
    const char *number;  // pointer to the printable representation
    char *array;         // array of size+1 byte for the digits
    size_t size;         // maximum number of significant digits
    size_t offset;       // offset if the first significant digit
    char buf[4 * sizeof(size_t)];
} bignum;

void bignum_free(bignum *bp) {
    if (bp && bp->array != bp->buf) {
        free(bp->array);
        bp->array = bp->buf;
    }
}

size_t bignum_realloc(bignum *bp) {
    size_t extra = bp->size / 2 + bp->size / 8;  // use pseudo Fibonacci sequence
    char *array, *source = bp->array;
    if (bp->array == bp->buf) {
        array = malloc(bp->size + extra);
    } else {
        source = array = realloc(bp->array, bp->size + extra);
    }
    if (array == NULL)
        return 0;
    memmove(array + extra, source, bp->size + 1);
    bp->array = array;
    bp->offset += extra;
    bp->size += extra;
    return extra;
}

bignum *bignum_normalize(bignum *bp) {
    while (bp->offset < bp->size && bp->array[bp->offset] == '0') {
        bp->offset++;
    }
    bp->number = bp->array + bp->offset;
    if (bp->offset >= bp->size) {
        bp->number = (bp->offset == bp->size) ? "0" : "overflow";
    }
    return  bp;
}

bignum *bignum_set_overflow(bignum *bp) {
    bp->offset = bp->size + 1;
    return bignum_normalize(bp);
}

bignum *bignum_init(bignum *bp) {
    size_t size = sizeof(bp->buf) - 1;
    bp->array = bp->buf;
    bp->offset = bp->size = size;
    bp->array[size] = '\0';
    return bignum_normalize(bp);
}

bignum *bignum_set_string(bignum *bp, const char *value) {
    bp->offset = bp->size;  // default value is 0
    if (value) {
        while (*value == '0') {
            value++;
        }
        size_t len = strlen(value);
        while (len > bp->size && bignum_realloc(bp))
            continue;
        if (len > bp->size) {
            bp->offset = bp->size + 1;  // overflow
        } else {
            bp->offset = bp->size - len;
            memcpy(bp->array + bp->offset, value, len);
        }
    }
    return bignum_normalize(bp);
}

bignum *bignum_init_string(bignum *bp, const char *value) {
    return bignum_set_string(bignum_init(bp), value);
}

bignum *bignum_add_uint(bignum *bp, unsigned int num) {
    unsigned long long carry = num;

    if (bp->offset > bp->size || num == 0) {
        return bp;
    }
    for (size_t i = bp->size;;) {
        if (carry == 0) {
            if (i < bp->offset) {
                bp->offset = i;
            }
            break;
        }
        if (i == 0) {
            bp->offset = 0;
            i += bignum_realloc(bp);
            if (i == 0) {
                bp->offset = bp->size + 1;  // overflow
                break;
            }
        }
        i--;
        if (i >= bp->offset) {
            carry += bp->array[i] - '0';
        }
        bp->array[i] = '0' + carry % 10;
        carry /= 10;
    }
    return bignum_normalize(bp);
}

bignum *bignum_set_value(bignum *bp, unsigned int value) {
    bp->offset = bp->size;
    return bignum_add_uint(bp, value);
}

bignum *bignum_init_value(bignum *bp, unsigned int value) {
    return bignum_add_uint(bignum_init(bp), value);
}

bignum *bignum_mul_uint(bignum *bp, unsigned int num) {
    unsigned long long carry = 0;

    if (bp->offset > bp->size || num == 1) {
        return bp;
    }
    if (num == 0) {
        bp->offset = bp->size;  /* special case 0 for speed */
    } else {
        for (size_t i = bp->size;;) {
            if (i == 0) {
                bp->offset = 0;
                if (carry == 0) {
                    break;
                }
                i += bignum_realloc(bp);
                if (i == 0) {
                    bp->offset = bp->size + 1;  // overflow
                    break;
                }
            }
            i--;
            if (i >= bp->offset) {
                carry += (unsigned long long)num * (bp->array[i] - '0');
            } else {
                if (carry == 0) {
                    bp->offset = i + 1;
                    break;
                }
            }
            bp->array[i] = '0' + carry % 10;
            carry /= 10;
        }
    }
    return bignum_normalize(bp);
}

int main(int argc, char *argv[]) {
    unsigned int n0 = 1, n1 = 100, max_digits = -1;
    bignum c;

    if (argc > 1) {
        n0 = n1 = strtol(argv[1], NULL, 0);
        if (argc > 2) {
            n1 = strtol(argv[2], NULL, 0);
            if (argc > 3) {
                max_digits = strtol(argv[3], NULL, 0);
            }
        }
        if (n1 < n0) {
            max_digits = n1;
            n1 = n0;
        }
    }

    bignum_init_value(&c, 1);
    printf("%6s | %s\n", "Number", "Factorial");
    for (unsigned int i = 1; i <= n1; i++) {
        bignum_mul_uint(&c, i);
        if (i >= n0) {
            if (c.size - c.offset > max_digits) {
                printf("%6u | %.1s.%.*se%zu\n",
                       i, c.number, max_digits - 1, c.number + 1,
                       c.size - c.offset - 1);
            } else {
                printf("%6u | %s\n", i, c.number);
            }
        }
    }
    bignum_free(&c);
    return 0;
}

You can pass arguments on the command line:

chqrlie@macbook ~/dev/stackoverflow > ./factomp 1 60

Number | Factorial
     1 | 1
     2 | 2
     3 | 6
     4 | 24
     5 | 120
     6 | 720
     7 | 5040
     8 | 40320
     9 | 362880
    10 | 3628800
    11 | 39916800
    12 | 479001600
    13 | 6227020800
    14 | 87178291200
    15 | 1307674368000
    16 | 20922789888000
    17 | 355687428096000
    18 | 6402373705728000
    19 | 121645100408832000
    20 | 2432902008176640000
    21 | 51090942171709440000
    22 | 1124000727777607680000
    23 | 25852016738884976640000
    24 | 620448401733239439360000
    25 | 15511210043330985984000000
    26 | 403291461126605635584000000
    27 | 10888869450418352160768000000
    28 | 304888344611713860501504000000
    29 | 8841761993739701954543616000000
    30 | 265252859812191058636308480000000
    31 | 8222838654177922817725562880000000
    32 | 263130836933693530167218012160000000
    33 | 8683317618811886495518194401280000000
    34 | 295232799039604140847618609643520000000
    35 | 10333147966386144929666651337523200000000
    36 | 371993326789901217467999448150835200000000
    37 | 13763753091226345046315979581580902400000000
    38 | 523022617466601111760007224100074291200000000
    39 | 20397882081197443358640281739902897356800000000
    40 | 815915283247897734345611269596115894272000000000
    41 | 33452526613163807108170062053440751665152000000000
    42 | 1405006117752879898543142606244511569936384000000000
    43 | 60415263063373835637355132068513997507264512000000000
    44 | 2658271574788448768043625811014615890319638528000000000
    45 | 119622220865480194561963161495657715064383733760000000000
    46 | 5502622159812088949850305428800254892961651752960000000000
    47 | 258623241511168180642964355153611979969197632389120000000000
    48 | 12413915592536072670862289047373375038521486354677760000000000
    49 | 608281864034267560872252163321295376887552831379210240000000000
    50 | 30414093201713378043612608166064768844377641568960512000000000000
    51 | 1551118753287382280224243016469303211063259720016986112000000000000
    52 | 80658175170943878571660636856403766975289505440883277824000000000000
    53 | 4274883284060025564298013753389399649690343788366813724672000000000000
    54 | 230843697339241380472092742683027581083278564571807941132288000000000000
    55 | 12696403353658275925965100847566516959580321051449436762275840000000000000
    56 | 710998587804863451854045647463724949736497978881168458687447040000000000000
    57 | 40526919504877216755680601905432322134980384796226602145184481280000000000000
    58 | 2350561331282878571829474910515074683828862318181142924420699914240000000000000
    59 | 138683118545689835737939019720389406345902876772687432540821294940160000000000000
    60 | 8320987112741390144276341183223364380754172606361245952449277696409600000000000000

The program is only limited by available memory, but it takes 54 seconds to compute factorial(100000), a number with 456574 decimal digits. You can pass a limit to the number of digits:

chqrlie@macbook ~/dev/stackoverflow > ./factomp 100000 60

Number | Factorial
100000 | 2.82422940796034787429342157802453551847749492609122485057891e456573

Answer 3

C normally uses the machines word size for integer storage. With factorials you end up running out of bits for storage quite quickly. At some point you will need to store the result in an array and account for carries. Assembly might be some help (since you might find multiplication with carries) or alternatively you can use the top half of your integer as the and write a multiplication algorithm similar to thecway you were taught to multiply in school.

If all this sounds tedious, you can use an arbitrary precision library.

Just so you know long long is 64 bits on many platforms and gcc supports int_128 on some platforms. There however will give limited milage with factorials.