# Arbitrary large unsigned integers

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 limits
• Do this in portable code
• Provide a reasonable "intuitive" interface on an opaque type
• Try to achieve good performance

Functions 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_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);

#endif


hugeint.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 = 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;

for (i = 0; i < bcdsize; ++i) buf[i] -= '0';

while (scanstart < bcdsize)
{
if (buf[bcdsize - 1] & 1) result->e[n] |= mask;
{
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;
}

{
const hugeint **summands = xmalloc(n * sizeof(hugeint *));

for (size_t i = 0; i < n; ++i)
{
summands[i] = va_arg(ap, hugeint *);
}

free(summands);

return result;
}

{
va_list ap;
va_start(ap, n);
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;
(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;
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 > other) return 1;
if (self->e < 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';
zero = 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);
{
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;
}
}

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.h"

int main(int argc, char **argv)
{
if (argc != 3)
{
fprintf(stderr, "Usage: %s [dividend] [divisor]\n", argv);
return 1;
}
hugeint *dividend = hugeint_parse(argv);
hugeint *divisor = hugeint_parse(argv);
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.h"

hugeint *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);
return 1;
}
hugeint *number = hugeint_parse(argv);
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.

• Note that allocating 0 bytes and receiving a NULL pointer does not indicate out-of-memory. if (!m) exit(1); --> if (!m && size) exit(1); This concern may not be possible with OP's code though. Jun 1 '17 at 3:05
• @chux, yes, I consider an attempt to allocate "nothing" a bug in this code, it's designed so this shouldn't happen. Jun 4 '17 at 9:19

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);


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 for

hugeint *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 of

auto 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 1–7 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, not

for (i = 0; i < n; ++i) {
while (summands[i]->n > result->n) {
result = hugeint_expand(result);
}
}


but rather

int 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.

• self is a name usually reserved for the "this pointer" --> Anything in C specified that supports this? Jun 1 '17 at 3:07
• You spotted an inconsistent change in design. The original idea was to have objects with mutator methods, self and the expand logic are leftovers from this. So thanks for bringing it to my attention! Also good catch about "add", i'll change this as well, and probably make the non-expanding version "private" (static) exclusively for subtracting. Jun 1 '17 at 5:53
• I addressed some of your remarks here, but now I have a question regarding the base type: I started using uintmax_t and later changed it to unsigned int under the assumption that calculations in the natural word size of the platform should perform best. So why would a fixed size (uint64_t) perform better? Isn't this possibly implemented in software in the compiler runtime on some platforms, so the software aritmetics would "cascade"? Jun 1 '17 at 11:08
• You can't "guarantee that uint64_t would be faster" than unsigned int, as implementations are permitted where they are the same type.... Jun 1 '17 at 12:50
• @chux: It's a keyword in Objective-C (and Objective-C++), Smalltalk, maybe some others. For its use in C code, you'd have to google around; I don't have any links on hand. (It's a pretty hard term to google for! Maybe just start with "C object-oriented programming" and see if you get some hits right away.) Jun 1 '17 at 18:08
• Multiplication cries for optimization (e.g. Karatzuba or Toom-Cook)

• Division also can be optimized; given a fast multiplication a most straightforward approach is Newton-Raphson.

• The semantics of hugeint_sub and hugeint_decrement is questionable, or, if I may say so, surprising.

• hugeint_ladd_cutoverflow assumes that the carry fits in an unsigned int. It is not necessarily a case.

• I don't think it is a good idea to interpret a position being NULL as an implicit shift by 1.

• Thanks, your links made me curious enough to invest some more time in this toy project. I now have a multiplication method roughly using karatsuba, which made better shift, add and sub routines necessary. Still it's not a huge performance boost with the example naively calculating factorials, but I guess that's because the naive factorial algorithm very often multiplies a huge number with a small one ... Jun 4 '17 at 9:07