Edit: There's now a Follow-up question concerning more recent code.From time to time, you see a question over at Stack Overflow about calculating factorials, typically asking why the result soon gets completely wrong when increasing the input number. The typical answers explain how types are limited, how you could approximate a factorial, how you could use a big number library like gmp, and so on. And you will stumble over a well-meant advice not to roll your own bignum implementation.Challenge accepted ;) So the following code is just for training, and, admittedly, for fun, and although I'm quite experienced, it was a learning experience as well. I'm putting the result for review because I might have missed some bugs, some unnecessary inefficiencies, or even might have broken portability in a subtle way, and if so, it would be great if someone could spot this.When I started implementing it, I had the following goals in mind:Handle any size of an unsigned integer, no compile-time limitsDo this in portable codeProvide a reasonable intuitive interface on an opaque typeTry to achieve good performanceFunctions for converting to and from a string are needed and because the naive implementations (using a lot of multiplications and divisions by 10) proved to be horrible bottlenecks, I replaced them with implementations based on the Double dabble algorithm.hugeint.h#ifndef HUGEINT_H#define HUGEINT_H#include <stdarg.h>#include <stddef.h>typedef struct hugeint hugeint;hugeint *hugeint_create(void);hugeint *hugeint_clone(const hugeint *self);hugeint *hugeint_fromUint(unsigned int val);hugeint *hugeint_parse(const char *str);hugeint *hugeint_ladd_cutoverflow(size_t n, const hugeint *const *summands,        unsigned int *residue);hugeint *hugeint_ladd(size_t n, const hugeint *const *summands);hugeint *hugeint_vadd(size_t n, va_list ap);hugeint *hugeint_add(size_t n, ...);hugeint *hugeint_2comp(const hugeint *self);hugeint *hugeint_sub(const hugeint *self, const hugeint *diff);hugeint *hugeint_shiftleft(const hugeint *hi, const hugeint *positions);hugeint *hugeint_shiftright(const hugeint *hi, const hugeint *positions);hugeint *hugeint_mult(const hugeint *hi, const hugeint *factor);hugeint *hugeint_div(const hugeint *hi, const hugeint *divisor,        hugeint **mod);int hugeint_isZero(const hugeint *self);int hugeint_compare(const hugeint *self, const hugeint *other);int hugeint_compareUint(const hugeint *self, unsigned int other);void hugeint_increment(hugeint **self);void hugeint_decrement(hugeint **self);char *hugeint_toString(const hugeint *self);#endifhugeint.c#include <limits.h>#include <stdlib.h>#include <string.h>#include hugeint.h#define HUGEINT_ELEMENT_BITS (CHAR_BIT * sizeof(unsigned int))// start each new hugeint with 256 bits:#define HUGEINT_INITIAL_ELEMENTS (256 / HUGEINT_ELEMENT_BITS)struct hugeint{    size_t n;    unsigned int e[];};static void *xmalloc(size_t size){    void *m = malloc(size);    if (!m) exit(1);    return m;}static void *xrealloc(void *m, size_t size){    void *m2 = realloc(m, size);    if (!m2) exit(1);    return m2;}static size_t copyNum(char **out, const char *str){    const char *p = str;    const char *start;    size_t length = 0;    while (*p && (*p == ' ' || *p == '\t' || *p == '0')) ++p;    if (*p < '0' || *p > '9') return 0;    start = p;    while (*p >= '0' && *p <= '9')    {        ++p;        ++length;    }    *out = xmalloc(length + 1);    (*out)[length] = 0;    memcpy(*out, start, length);    return length;}static hugeint *hugeint_expand(hugeint *self){    hugeint *expanded = xrealloc(self,            sizeof(hugeint) + 2 * self->n * sizeof(unsigned int));    memset(&(expanded->e[expanded->n]), 0, expanded->n * sizeof(unsigned int));    expanded->n *= 2;    return expanded;}hugeint *hugeint_create(void){    hugeint *self = xmalloc(sizeof(hugeint)            + HUGEINT_INITIAL_ELEMENTS * sizeof(unsigned int));    memset(self, 0, sizeof(hugeint)            + HUGEINT_INITIAL_ELEMENTS * sizeof(unsigned int));    self->n = HUGEINT_INITIAL_ELEMENTS;    return self;}hugeint *hugeint_clone(const hugeint *self){    hugeint *clone = xmalloc(sizeof(hugeint) + self->n * sizeof(unsigned int));    memcpy(clone, self, sizeof(hugeint) + self->n * sizeof(unsigned int));    return clone;}hugeint *hugeint_fromUint(unsigned int val){    hugeint *self = hugeint_create();    self->e[0] = val;    return self;}hugeint *hugeint_parse(const char *str){    char *buf;    hugeint *result = hugeint_create();    size_t bcdsize = copyNum(&buf, str);    if (!bcdsize) return result;    size_t scanstart = 0;    size_t n = 0;    size_t i;    unsigned int mask = 1;    for (i = 0; i < bcdsize; ++i) buf[i] -= '0';    while (scanstart < bcdsize)    {        if (buf[bcdsize - 1] & 1) result->e[n] |= mask;        mask <<= 1;        if (!mask)        {            mask = 1;            if (++n == result->n) result = hugeint_expand(result);        }        for (i = bcdsize - 1; i > scanstart; --i)        {            buf[i] >>= 1;            if (buf[i-1] & 1) buf[i] |= 8;        }        buf[scanstart] >>= 1;        while (scanstart < bcdsize && !buf[scanstart]) ++scanstart;        for (i = scanstart; i < bcdsize; ++i)        {            if (buf[i] > 7) buf[i] -= 3;        }    }    free(buf);    return result;}hugeint *hugeint_ladd_cutoverflow(size_t n, const hugeint *const *summands,        unsigned int *residue){    hugeint *result = hugeint_create();    size_t i;    for (i = 0; i < n; ++i)    {        while (summands[i]->n > result->n)        {            result = hugeint_expand(result);        }    }    unsigned int res = 0;    for (i = 0; i < result->n; ++i)    {        for (unsigned int bit = 1; bit; bit <<= 1)        {            for (size_t j = 0; j < n; ++j)            {                if (i >= summands[j]->n) continue;                if (summands[j]->e[i] & bit) ++res;            }            if (res & 1) result->e[i] |= bit;            res >>= 1;        }    }    if (residue) *residue = res;    return result;}hugeint *hugeint_ladd(size_t n, const hugeint *const *summands){    unsigned int res;    hugeint *result = hugeint_ladd_cutoverflow(n, summands, &res);    if (res)    {        size_t i = result->n;        result = hugeint_expand(result);        result->e[i] = res;    }    return result;}hugeint *hugeint_vadd(size_t n, va_list ap){    const hugeint **summands = xmalloc(n * sizeof(hugeint *));    for (size_t i = 0; i < n; ++i)    {        summands[i] = va_arg(ap, hugeint *);    }    hugeint *result = hugeint_ladd(n, summands);    free(summands);    return result;}hugeint *hugeint_add(size_t n, ...){    va_list ap;    va_start(ap, n);    hugeint *result = hugeint_vadd(n, ap);    va_end(ap);    return result;}hugeint *hugeint_2comp(const hugeint *self){    hugeint *tmp = hugeint_clone(self);    if (hugeint_isZero(tmp)) return tmp;    hugeint *one = hugeint_fromUint(1);    for (size_t i = 0; i < tmp->n; ++i)    {        tmp->e[i] = ~(tmp->e[i]);    }    hugeint *result = hugeint_add(2, tmp, one);    free(tmp);    free(one);    return result;}hugeint *hugeint_sub(const hugeint *self, const hugeint *diff){    if (diff->n > self->n) return 0;    int freediff = 0;    if (diff->n < self->n)    {        freediff = 1;        diff = hugeint_clone(diff);        while (diff->n < self->n) diff = hugeint_expand((hugeint *)diff);    }    hugeint *tmp = hugeint_2comp(diff);    if (freediff) free((hugeint *)diff);    unsigned int res;    hugeint *result = hugeint_ladd_cutoverflow(2,            (const hugeint *const []){self, tmp}, &res);    free(tmp);    if (res > 1)    {        free(result);        return 0;    }    return result;}hugeint *hugeint_shiftleft(const hugeint *hi, const hugeint *positions){    hugeint *result = hugeint_clone(hi);    hugeint *count;    if (positions)    {        if (hugeint_isZero(positions)) return result;        count = hugeint_clone(positions);    }    else    {        count = 0;    }    unsigned int highbit = 1U << (HUGEINT_ELEMENT_BITS - 1);    do    {        if (result->e[result->n - 1] & highbit)        {            result = hugeint_expand(result);        }        int incelement = 0;        for (size_t i = 0; i < result->n; ++i)        {            int overflow = !!(result->e[i] & highbit);            result->e[i] <<= 1;            if (incelement) ++result->e[i];            incelement = overflow;        }        if (count) hugeint_decrement(&count);    } while (count && !hugeint_isZero(count));    free(count);    return result;}hugeint *hugeint_shiftright(const hugeint *hi, const hugeint *positions){    hugeint *result = hugeint_clone(hi);    hugeint *count;    if (positions)    {        if (hugeint_isZero(positions)) return result;        count = hugeint_clone(positions);    }    else    {        count = 0;    }    unsigned int highbit = 1U << (HUGEINT_ELEMENT_BITS - 1);    do    {        int incelement = 0;        size_t i = result->n;        while (i > 0)        {            --i;            int overflow = result->e[i] & 1;            result->e[i] >>= 1;            if (incelement) result->e[i] += highbit;            incelement = overflow;        }        if (count) hugeint_decrement(&count);    } while (count && !hugeint_isZero(count));    free(count);    return result;}hugeint *hugeint_mult(const hugeint *hi, const hugeint *factor){    hugeint **summands=xmalloc(factor->n * HUGEINT_ELEMENT_BITS * sizeof(hugeint *));    size_t n = 0;    hugeint *bitnum = hugeint_create();    for (size_t i = 0; i < factor->n; ++i)    {        for (unsigned int bit = 1; bit; bit <<= 1)        {            if (factor->e[i] & bit)            {                hugeint *summand = hugeint_shiftleft(hi, bitnum);                summands[n++] = summand;            }            hugeint_increment(&bitnum);        }    }    free(bitnum);    hugeint *result = hugeint_ladd(n, (const hugeint **)summands);    for (size_t j = 0; j < n; ++j) free(summands[j]);    free(summands);    return result;}hugeint *hugeint_div(const hugeint *hi, const hugeint *divisor,        hugeint **mod){    if (hugeint_isZero(divisor)) return 0;    hugeint *scaled_divisor = hugeint_clone(divisor);    hugeint *remain = hugeint_clone(hi);    hugeint *result = hugeint_create();    hugeint *multiple = hugeint_fromUint(1);    while (hugeint_compare(scaled_divisor, hi) < 0)    {        hugeint *tmp = hugeint_shiftleft(scaled_divisor, 0);        free(scaled_divisor);        scaled_divisor = tmp;        tmp = hugeint_shiftleft(multiple, 0);        free(multiple);        multiple = tmp;    }    do    {        if (hugeint_compare(remain, scaled_divisor) >= 0)        {            hugeint *tmp = hugeint_sub(remain, scaled_divisor);            free(remain);            remain = tmp;            tmp = hugeint_add(2, result, multiple);            free(result);            result = tmp;        }        hugeint *tmp = hugeint_shiftright(scaled_divisor, 0);        free(scaled_divisor);        scaled_divisor = tmp;        tmp = hugeint_shiftright(multiple, 0);        free(multiple);        multiple = tmp;    } while (!hugeint_isZero(multiple));    if (mod) *mod = remain;    else free(remain);    free(multiple);    free(scaled_divisor);    return result;}int hugeint_isZero(const hugeint *self){    for (size_t i = 0; i < self->n; ++i)    {        if (self->e[i]) return 0;    }    return 1;}int hugeint_compare(const hugeint *self, const hugeint *other){    size_t n;    if (self->n > other->n)    {        for (size_t i = other->n; i < self->n; ++i)        {            if (self->e[i]) return 1;        }        n = other->n;    }    else if (self->n < other->n)    {        for (size_t i = self->n; i < other->n; ++i)        {            if (other->e[i]) return -1;        }        n = self->n;    }    else n = self->n;    while (n > 0)    {        --n;        if (self->e[n] > other->e[n]) return 1;        if (self->e[n] < other->e[n]) return -1;    }    return 0;}int hugeint_compareUint(const hugeint *self, unsigned int other){    for (size_t i = self->n - 1; i > 0; --i)    {        if (self->e[i]) return 1;    }    if (self->e[0] > other) return 1;    if (self->e[0] < other) return -1;    return 0;}void hugeint_increment(hugeint **self){    int carry = 0;    for (size_t i = 0; i < (*self)->n; ++i)    {        carry = !++(*self)->e[i];        if (!carry) break;    }    if (carry)    {        size_t n = (*self)->n;        *self = hugeint_expand(*self);        if (*self) (*self)->e[n] = 1;    }}void hugeint_decrement(hugeint **self){    if (hugeint_isZero(*self)) return;    for (size_t i = 0; i < (*self)->n; ++i)    {        if ((*self)->e[i]--) break;    }}char *hugeint_toString(const hugeint *self){    if (hugeint_isZero(self))    {        char *zero = malloc(2);        zero[0] = '0';        zero[1] = 0;        return zero;    }    size_t nbits = HUGEINT_ELEMENT_BITS * self->n;    size_t bcdsize = nbits/3;    size_t scanstart = bcdsize - 2;    char *buf = xmalloc(bcdsize + 1);    memset(buf, 0, bcdsize + 1);    size_t i, j;    i = self->n;    while(i--)    {        unsigned int mask = 1U << (HUGEINT_ELEMENT_BITS - 1);        while (mask)        {            int bit = !!(self->e[i] & mask);            for (j = scanstart; j < bcdsize; ++j)            {                if (buf[j] > 4) buf[j] += 3;            }            if (buf[scanstart] > 7) scanstart -= 1;            for (j = scanstart; j < bcdsize - 1; ++j)            {                buf[j] <<= 1;                buf[j] &= 0xf;                buf[j] |= (buf[j+1] > 7);            }            buf[bcdsize-1] <<= 1;            buf[bcdsize-1] &= 0xf;            buf[bcdsize-1] |= bit;            mask >>= 1;        }    }    for (i = 0; i < bcdsize - 1; ++i)    {        if (buf[i]) break;    }    bcdsize -= i;    memmove(buf, buf + i, bcdsize + 1);    for (i = 0; i < bcdsize; ++i) buf[i] += '0';    buf = xrealloc(buf, bcdsize + 1);    return buf;}Example 1: divide.c#include <stdio.h>#include <stdlib.h>#include hugeint.hint main(int argc, char **argv){    if (argc != 3)    {        fprintf(stderr, Usage: %s [dividend] [divisor]\n, argv[0]);        return 1;    }    hugeint *dividend = hugeint_parse(argv[1]);    hugeint *divisor = hugeint_parse(argv[2]);    hugeint *remain;    hugeint *result = hugeint_div(dividend, divisor, &remain);    free(divisor);    free(dividend);    if (!result)    {        fputs(Error: division unsuccessful.\n, stderr);        return 1;    }    char *resultstr = hugeint_toString(result);    free(result);    char *remainstr = hugeint_toString(remain);    free(remain);    puts(resultstr);    free(resultstr);    fputs(remainder: , stdout);    puts(remainstr);    free(remainstr);    return 0;}Example 2: factorial.c#include <stdio.h>#include <stdlib.h>#include hugeint.hhugeint *factorial(hugeint *self){    hugeint *result = hugeint_fromUint(1);    hugeint *factor = hugeint_clone(self);    while (hugeint_compareUint(factor, 1) > 0)    {        hugeint *tmp = hugeint_mult(result, factor);        free(result);        result = tmp;        hugeint_decrement(&factor);    }    free(factor);    return result;}int main(int argc, char **argv){    if (argc != 2)    {        fprintf(stderr, Usage: %s [number]\n, argv[0]);        return 1;    }    hugeint *number = hugeint_parse(argv[1]);    hugeint *result = factorial(number);    free(number);    char *factstr = hugeint_toString(result);    free(result);    puts(factstr);    free(factstr);    return 0;}The code as presented in the question can be browsed on github.The current HEAD contains improvements mostly based on suggestions/criticism from the answers.
Arbitrary large unsigned integers
c;integer
Looks like fairly clean code, but I see plenty of places it could be improved. Letting my eyes skim over it and stop on the weird bits...hugeint *hugeint_ladd_cutoverflow(size_t n, const hugeint *const *summands,        unsigned int *residue);hugeint *hugeint_ladd(size_t n, const hugeint *const *summands);hugeint *hugeint_vadd(size_t n, va_list ap);hugeint *hugeint_add(size_t n, ...);This is weird. You have a sane-looking sub function, and a sane-looking mult function hugeint *hugeint_sub(const hugeint *self, const hugeint *diff);hugeint *hugeint_mult(const hugeint *hi, const hugeint *factor); yet somehow when it comes to the simplest possible operation, add, you have all these weird helper functions and no way to just add two numbers? I'm looking forhugeint *hugeint_add(const hugeint *a, const hugeint *b);Having to write that function signature made me notice that your sub function's parameters are named self and diff, which makes no sense to me. self is a name usually reserved for the this pointer in OO style, so from the signature of this function I'd be expecting the equivalent ofauto hugeint_sub(auto self, auto diff) {    self -= diff;    return self;}But that's not what the function does  nor what I'd expected it to do before looking at the parameter names. So basically, your parameter names are misleading. The English names for them would be minuend and subtrahend, but personally I'd go with a and b. I think everyone knows what sub(a, b) means. :)So let's look at that weird N-argument add function. Does it work?hugeint *hugeint_ladd(size_t n, const hugeint *const *summands){    unsigned int res;    hugeint *result = hugeint_ladd_cutoverflow(n, summands, &res);    if (res)    {        size_t i = result->n;        result = hugeint_expand(result);        result->e[i] = res;    }    return result;}We can tell instantly that it can't possibly work, because it's got this thing called residue that's basically a carry flag for the addition. The reason we know it can't be right is that this carry flag doesn't have enough bits to store an arbitrarily large number of carries!  Try adding together 4294967297 copies of the hugeint 4294967295 and see what happens. One bit of carry is only enough for 2 addends, and 32 bits of carry is only enough for 4294967296 addends. What you need is a hugeint's worth of carry; or else just get rid of this whole carry business and just expand the result as you go; or else get rid of the notion that it's a good idea to add arbitrarily many addends and just make a hugeint_add(a, b) that does the obvious thing.Speaking of trouble with carries, let's look at the shift-left and shift-right functions. They're interesting because you're using hugeint for the shift count itself, which is unlike any programming language or instruction set I'm aware of. The shift count can only sanely get as high as the log of the left-hand operand, which is to say, it had better fit into size_t. But okay, let's look at right-shift and find the place where you do If the right-hand operand is bigger than 4 billion, just set the result to 0 and return...    if (count) hugeint_decrement(&count);} while (count && !hugeint_isZero(count));...Oh dear. You should definitely give this one a rewrite. At minimum, it should know that bit-shifting by a multiple of 8 is equivalent to a single memmove, and bit-shifting by any other number can be reduced to a memmove plus 17 loop iterations.Speaking of bit-widths, why are you using unsigned int as your unit data type instead of uint64_t or __uint128_t? I guarantee that uint64_t would be faster for most of what you're doing.The other confusing thing about that right-shift code is:if (positions) {    if (hugeint_isZero(positions)) return result;    count = hugeint_clone(positions);} else {    count = 0;}Here we're assigning 0 (an int) to count (a hugeint). This isn't going to do what we wanted! ...Or rather, it is, but only via a very bizarre path. Here 0 means NULL, which because of the if (count) elsewhere in the function, means do only 1 iteration. So this 0 semantically means 1!The ultimate effect of all this spaghetti code is to make it so that if the user accidentally passes in a null pointer hugeint *result = hugeint_shiftright(a, NULL); then his operand will be shifted by 1. If this behavior is desirable, it should be provided as a separate function:hugeint *result = hugeint_shiftright_by_one(a);Providing two separate functions for the two separate functionalities is not only good for the reader's sanity, it's also going to be more efficient, because you won't have all those ifs confusing the CPU's branch predictor.Getting back to add, I see this hugeint_expand function that you're calling all over the place. It doubles the length of the underlying array, which is way overkill. What you should be doing is computing the correct length for the array and then allocating exactly that much: that is, notfor (i = 0; i < n; ++i) {    while (summands[i]->n > result->n) {        result = hugeint_expand(result);    }}but ratherint result_n = result->n;for (size_t i = 0; i < n; ++i) {    result_n = max(result_n, summands[i]->n);}result = hugeint_expand(result, result_n);That is, hugeint_expand needs to take a parameter that tells it the size of the array you want. Compare this API to std::vector::reserve in C++; it should look very very similar.There's definitely more to critique, but this should give you some ideas, anyway.