4
\$\begingroup\$

This is a followup question to Arbitrary large unsigned integers. The feedback was very helpful and made me redesign large portions of the code.

I changed the struct to distinguish between the number of meaningful elements (n) and the number of elements there's space allocated for (s):

struct hugeint
{
    size_t s;
    size_t n;
    unsigned int e[];
};

And the helper function to change the size of an instance became a little more complex:

static hugeint *hugeint_scale(hugeint *self, size_t newSize)
{
    if (newSize == self->n) return self;
    if (newSize > self->s)
    {
        size_t s = self->s;
        while (newSize > s) s *= 2;
        self = xrealloc(self,
                sizeof(hugeint) + s * sizeof(unsigned int));
        self->s = s;
    }
    if (newSize > self->n)
    {
        memset(&(self->e[self->n]), 0,
                (self->s - self->n) * sizeof(unsigned int));
    }
    else if (newSize < self->n)
    {
        memset(&(self->e[newSize]), 0,
                (self->n - newSize) * sizeof(unsigned int));
    }
    self->n = newSize;
    return self;
}

Now I have the following function for an "in place addition":

void hugeint_addToSelf(hugeint **self, const hugeint *other)
{
    if ((*self)->n < other->n)
    {
        *self = hugeint_scale(*self, other->n);
    }

    unsigned int carry = 0;
    size_t i;
    for (i = 0; i < other->n; ++i)
    {
        unsigned int v = (*self)->e[i] + other->e[i];
        unsigned int nextCarry = v < other->e[i];
        (*self)->e[i] = v + carry;
        if (!carry || (*self)->e[i]) carry = nextCarry;
    }
    for (; i < (*self)->n; ++i)
    {
        (*self)->e[i] = (*self)->e[i] + carry;
        carry = carry && !(*self)->e[i];
    }
    if (carry)
    {
        *self = hugeint_scale(*self, (*self)->n + 1);
        (*self)->e[(*self)->n - 1] = 1;
    }
}

and a very similar function for "in place subtraction":

void hugeint_subFromSelf(hugeint **self, const hugeint *other)
{
    if (hugeint_isZero(other)) return;
    int comp = hugeint_compare(*self, other);
    if (comp < 0)
    {
        free(*self);
        *self = 0;
        return;
    }
    if (comp == 0)
    {
        *self = hugeint_scale(*self, 1);
        (*self)->e[0] = 0;
        return;
    }

    unsigned int carry = 1;
    for (size_t i = 0; i < (*self)->n; ++i)
    {
        unsigned int v = (*self)->e[i] + (i < other->n ? ~(other->e[i]) : ~0U);
        unsigned int nextCarry = v < (*self)->e[i];
        (*self)->e[i] = v + carry;
        if (!carry || (*self)->e[i]) carry = nextCarry;
    }
    hugeint_autoscale(self);
}

Basically, I have two questions left (and more ideas to improve the code are of course welcome as well!):

  • Is there some better way to handle carry? As I aim for a pure C solution, I can't take advantage of a carry flag, but this "soft carry" kind of looks over-complicated?
  • As pointed out in the comments on my earlier question, using unsigned int I do all calculations in 32bit on x86_64, therefore missing on better performance by using 64bit arithmetic instructions of the processor. Can somebody suggest me a good way to pick the "ideal" base int type, so I use the widest type with direct CPU support on any platform?

For reference, the whole code:

hugeint.h:

#ifndef HUGEINT_H
#define HUGEINT_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_parseHex(const char *str);

hugeint *hugeint_add(const hugeint *a, const hugeint *b);
hugeint *hugeint_sub(const hugeint *minuend, const hugeint *subtrahend);
hugeint *hugeint_mult(const hugeint *a, const hugeint *b);
hugeint *hugeint_div(const hugeint *dividend, const hugeint *divisor,
        hugeint **remainder);

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);
void hugeint_addToSelf(hugeint **self, const hugeint *other);
void hugeint_subFromSelf(hugeint **self, const hugeint *other);
void hugeint_shiftLeft(hugeint **self, size_t positions);
void hugeint_shiftRight(hugeint **self, size_t positions);

char *hugeint_toString(const hugeint *self);
char *hugeint_toHexString(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))
#define HUGEINT_INITIAL_ELEMENTS (256 / HUGEINT_ELEMENT_BITS)

struct hugeint
{
    size_t s;
    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_scale(hugeint *self, size_t newSize)
{
    if (newSize == self->n) return self;
    if (newSize > self->s)
    {
        size_t s = self->s;
        while (newSize > s) s *= 2;
        self = xrealloc(self,
                sizeof(hugeint) + s * sizeof(unsigned int));
        self->s = s;
    }
    if (newSize > self->n)
    {
        memset(&(self->e[self->n]), 0,
                (self->s - self->n) * sizeof(unsigned int));
    }
    else if (newSize < self->n)
    {
        memset(&(self->e[newSize]), 0,
                (self->n - newSize) * sizeof(unsigned int));
    }
    self->n = newSize;
    return self;
}

static void hugeint_autoscale(hugeint **self)
{
    (*self)->n = (*self)->s;
    while ((*self)->n > 1 && !(*self)->e[(*self)->n-1]) --(*self)->n;
}

static hugeint *hugeint_createSized(size_t size)
{
    size_t s = size;
    if (s < HUGEINT_INITIAL_ELEMENTS) s = HUGEINT_INITIAL_ELEMENTS;
    hugeint *self = xmalloc(sizeof(hugeint) + s * sizeof(unsigned int));
    memset(self, 0, sizeof(hugeint) + s * sizeof(unsigned int));
    self->s = s;
    self->n = size;
    return self;
}

hugeint *hugeint_create(void)
{
    return hugeint_createSized(1);
}

hugeint *hugeint_clone(const hugeint *self)
{
    hugeint *clone = xmalloc(sizeof(hugeint) + self->s * sizeof(unsigned int));
    memcpy(clone, self, sizeof(hugeint) + self->s * 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_scale(result, result->n + 1);
        }
        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;
}

static unsigned char hexNibble(const char *str)
{
    unsigned char hex = *str;
    if (hex >= 'a' && hex <='f')
    {
        hex -= ('a' - 10);
    }
    else if (hex >= 'A' && hex <= 'F')
    {
        hex -= ('A' - 10);
    }
    else if (hex >= '0' && hex <= '9')
    {
        hex -= '0';
    }
    else
    {
        hex = 0;
    }
    return hex;
}

hugeint *hugeint_parseHex(const char *str)
{
    size_t len = strlen(str);
    size_t bits = len * 4;
    size_t n = HUGEINT_ELEMENT_BITS / bits;
    size_t leading = HUGEINT_ELEMENT_BITS % bits;
    size_t i = n;
    if (leading) ++n;
    hugeint *result = hugeint_createSized(n);
    if (leading)
    {
        unsigned int shift = leading;
        while (shift)
        {
            shift -= 4;
            unsigned char nibble = hexNibble(str++);
            result->e[i] |= ((unsigned int)nibble << shift);
        }
    }
    while (i)
    {
        --i;
        unsigned int shift = HUGEINT_ELEMENT_BITS;
        while (shift)
        {
            shift -= 4;
            unsigned char nibble = hexNibble(str++);
            result->e[i] |= ((unsigned int)nibble << shift);
        }
    }
    return result;
}

hugeint *hugeint_add(const hugeint *a, const hugeint *b)
{
    if (a->n < b->n)
    {
        const hugeint *tmp = a;
        a = b;
        b = tmp;
    }

    hugeint *result = hugeint_clone(a);
    hugeint_addToSelf(&result, b);
    return result;
}

hugeint *hugeint_sub(const hugeint *minuend, const hugeint *subtrahend)
{
    hugeint *result = hugeint_clone(minuend);
    hugeint_subFromSelf(&result, subtrahend);
    return result;
}

static hugeint *multUint(unsigned int a, unsigned int b)
{
    if (!a || !b) return hugeint_create();
    if (a < b)
    {
        unsigned int tmp = a;
        a = b;
        b = tmp;
    }
    unsigned int halfMask = (1 << (HUGEINT_ELEMENT_BITS / 2)) - 1;
    unsigned int maskA = halfMask;
    unsigned int maskB = halfMask;
    unsigned int direct = 0;
    while (maskB)
    {
        if (a > maskA && b > maskB) break;
        if (a <= maskA && b <= maskB)
        {
            direct = 1;
            break;
        }
        maskA = (maskA << 1) | 1;
        maskB >>= 1;
    }
    if (direct)
    {
        return hugeint_fromUint(a * b);
    }
    unsigned int ah = a >> (HUGEINT_ELEMENT_BITS / 2);
    unsigned int al = a & halfMask;
    unsigned int bh = b >> (HUGEINT_ELEMENT_BITS / 2);
    unsigned int bl = b & halfMask;

    hugeint *p1 = multUint(ah, bh);
    hugeint *p2 = multUint(al, bl);
    unsigned int p3f1u = ah + al;
    unsigned int p3f2u = bh + bl;
    hugeint *p3;
    if (p3f1u < ah || p3f2u < bh)
    {
        hugeint *p3f1;
        hugeint *p3f2;
        if (p3f1u < ah)
        {
            p3f1 = hugeint_createSized(2);
            p3f1->e[0] = p3f1u;
            p3f1->e[1] = 1;
        }
        else p3f1 = hugeint_fromUint(p3f1u);
        if (p3f2u < bh)
        {
            p3f2 = hugeint_createSized(2);
            p3f2->e[0] = p3f2u;
            p3f2->e[1] = 1;
        }
        else p3f2 = hugeint_fromUint(p3f2u);
        p3 = hugeint_mult(p3f1, p3f2);
        free(p3f1);
        free(p3f2);
    }
    else
    {
        p3 = multUint(p3f1u, p3f2u);
    }
    hugeint_subFromSelf(&p3, p2);
    hugeint_subFromSelf(&p3, p1);
    hugeint_shiftLeft(&p3, HUGEINT_ELEMENT_BITS - (HUGEINT_ELEMENT_BITS / 2));
    hugeint_shiftLeft(&p1, HUGEINT_ELEMENT_BITS);
    hugeint_addToSelf(&p3, p1);
    hugeint_addToSelf(&p3, p2);
    free(p2);
    free(p1);
    return p3;
}

hugeint *hugeint_mult(const hugeint *a, const hugeint *b)
{
    if (hugeint_isZero(a) || hugeint_isZero(b)) return hugeint_create();
    if (b->n > a->n)
    {
        const hugeint *tmp = a;
        a = b;
        b = tmp;
    }
    if (a->n == 1) return multUint(a->e[0], b->e[0]);

    size_t nh = a->n / 2;
    size_t nl = a->n - nh;
    size_t bnl = nl;
    if (b->n < bnl) bnl = b->n;

    hugeint *ah = hugeint_createSized(nh);
    memcpy(&(ah->e), &(a->e[nl]), nh * sizeof(unsigned int));
    hugeint *al = hugeint_createSized(nl);
    memcpy(&(al->e), &(a->e), nl * sizeof(unsigned int));
    hugeint *bh;
    if (b->n > nl)
    {
        bh = hugeint_createSized(b->n - nl);
        memcpy(&(bh->e), &(b->e[nl]), (b->n - nl) * sizeof(unsigned int));
    }
    else
    {
        bh = hugeint_create();
    }
    hugeint *bl = hugeint_createSized(bnl);
    memcpy(&(bl->e), &(b->e), bnl * sizeof(unsigned int));

    hugeint_autoscale(&ah);
    hugeint_autoscale(&al);
    hugeint_autoscale(&bh);
    hugeint_autoscale(&bl);

    hugeint *p1 = hugeint_mult(ah, bh);
    hugeint *p2 = hugeint_mult(al, bl);
    hugeint_addToSelf(&ah, al);
    hugeint_addToSelf(&bh, bl);
    free(al);
    free(bl);
    hugeint *p3 = hugeint_mult(ah, bh);
    free(ah);
    free(bh);
    hugeint_subFromSelf(&p3, p2);
    hugeint_subFromSelf(&p3, p1);
    hugeint_shiftLeft(&p3, nl * HUGEINT_ELEMENT_BITS);
    hugeint_shiftLeft(&p1, 2 * nl * HUGEINT_ELEMENT_BITS);
    hugeint_addToSelf(&p3, p1);
    hugeint_addToSelf(&p3, p2);
    free(p2);
    free(p1);
    return p3;
}

hugeint *hugeint_div(const hugeint *dividend, const hugeint *divisor,
        hugeint **remainder)
{
    if (hugeint_isZero(divisor)) return 0;

    hugeint *scaledDivisor = hugeint_clone(divisor);
    hugeint *remain = hugeint_clone(dividend);
    hugeint *result = hugeint_create();
    hugeint *multiple = hugeint_fromUint(1);

    while (hugeint_compare(scaledDivisor, dividend) < 0)
    {
        hugeint_shiftLeft(&scaledDivisor, 1);
        hugeint_shiftLeft(&multiple, 1);
    }

    do
    {
        if (hugeint_compare(remain, scaledDivisor) >= 0)
        {
            hugeint_subFromSelf(&remain, scaledDivisor);
            hugeint_addToSelf(&result, multiple);
        }
        hugeint_shiftRight(&scaledDivisor, 1);
        hugeint_shiftRight(&multiple, 1);
    } while (!hugeint_isZero(multiple));

    if (remainder) *remainder = remain;
    else free(remain);
    free(multiple);
    free(scaledDivisor);
    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_scale(*self, n + 1);
        (*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;
    }
    hugeint_autoscale(self);
}

void hugeint_addToSelf(hugeint **self, const hugeint *other)
{
    if ((*self)->n < other->n)
    {
        *self = hugeint_scale(*self, other->n);
    }

    unsigned int carry = 0;
    size_t i;
    for (i = 0; i < other->n; ++i)
    {
        unsigned int v = (*self)->e[i] + other->e[i];
        unsigned int nextCarry = v < other->e[i];
        (*self)->e[i] = v + carry;
        if (!carry || (*self)->e[i]) carry = nextCarry;
    }
    for (; i < (*self)->n; ++i)
    {
        (*self)->e[i] = (*self)->e[i] + carry;
        carry = carry && !(*self)->e[i];
    }
    if (carry)
    {
        *self = hugeint_scale(*self, (*self)->n + 1);
        (*self)->e[(*self)->n - 1] = 1;
    }
}

void hugeint_subFromSelf(hugeint **self, const hugeint *other)
{
    if (hugeint_isZero(other)) return;
    int comp = hugeint_compare(*self, other);
    if (comp < 0)
    {
        free(*self);
        *self = 0;
        return;
    }
    if (comp == 0)
    {
        *self = hugeint_scale(*self, 1);
        (*self)->e[0] = 0;
        return;
    }

    unsigned int carry = 1;
    for (size_t i = 0; i < (*self)->n; ++i)
    {
        unsigned int v = (*self)->e[i] + (i < other->n ? ~(other->e[i]) : ~0U);
        unsigned int nextCarry = v < (*self)->e[i];
        (*self)->e[i] = v + carry;
        if (!carry || (*self)->e[i]) carry = nextCarry;
    }
    hugeint_autoscale(self);
}

void hugeint_shiftLeft(hugeint **self, size_t positions)
{
    if (!positions) return;
    if (hugeint_isZero(*self)) return;
    size_t shiftElements = positions / HUGEINT_ELEMENT_BITS;
    unsigned int shiftBits = positions % HUGEINT_ELEMENT_BITS;
    size_t oldSize = (*self)->n;
    unsigned int topBits = 0;
    if (shiftBits) topBits = (*self)->e[oldSize - 1]
            >> (HUGEINT_ELEMENT_BITS - shiftBits);
    size_t newSize = oldSize + shiftElements + !!topBits;

    if (newSize > oldSize) *self = hugeint_scale(*self, newSize);

    if (shiftElements)
    {
        memmove(&((*self)->e[shiftElements]), &((*self)->e[0]),
                oldSize * sizeof(unsigned int));
        memset(&((*self)->e[0]), 0, shiftElements * sizeof(unsigned int));
    }

    if (shiftBits)
    {
        unsigned int overflowBits = 0;
        for (size_t i = shiftElements; i < newSize; ++i)
        {
            unsigned int nextOverflow = (*self)->e[i]
                    >> (HUGEINT_ELEMENT_BITS - shiftBits);
            (*self)->e[i] <<= shiftBits;
            (*self)->e[i] |= overflowBits;
            overflowBits = nextOverflow;
        }
    }
    hugeint_autoscale(self);
}

void hugeint_shiftRight(hugeint **self, size_t positions)
{
    if (!positions) return;
    if (hugeint_isZero(*self)) return;
    size_t shiftElements = positions / HUGEINT_ELEMENT_BITS;
    if (shiftElements >= (*self)->n)
    {
        *self = hugeint_scale(*self, 1);
        (*self)->e[0] = 0;
        return;
    }

    unsigned int shiftBits = positions % HUGEINT_ELEMENT_BITS;

    if (shiftElements)
    {
        memmove(&((*self)->e[0]), &((*self)->e[shiftElements]),
                ((*self)->n - shiftElements) * sizeof(unsigned int));
    }

    if (shiftBits)
    {
        unsigned int overflowBits = 0;
        size_t i = (*self)->n - shiftElements;
        while (i)
        {
            --i;
            unsigned int nextOverflow = (*self)->e[i]
                    << (HUGEINT_ELEMENT_BITS - shiftBits);
            (*self)->e[i] >>= shiftBits;
            (*self)->e[i] |= overflowBits;
            overflowBits = nextOverflow;
        }
    }
    hugeint_autoscale(self);
}

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 = 1 << (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;
}

char *hugeint_toHexString(const hugeint *self)
{
    size_t len = self->n * HUGEINT_ELEMENT_BITS / 4;

    unsigned int mask = 0xfU << (HUGEINT_ELEMENT_BITS - 4);
    while (mask && !(self->e[self->n-1] & mask))
    {
        --len;
        mask >>= 4;
    }
    if (!len)
    {
        char *result = malloc(2);
        result[0] = '0';
        result[1] = 0;
        return result;
    }

    char *result = malloc(len + 1);
    result[len] = 0;

    size_t i = 0;
    while (len)
    {
        unsigned int v = self->e[i];
        for (unsigned int j = 0; j < HUGEINT_ELEMENT_BITS / 4; ++j)
        {
            unsigned char nibble = v & 0xf;
            result[--len] = nibble > 9 ? nibble - 10 + 'a' : nibble + '0';
            if (!len) break;
            v >>= 4;
        }
        ++i;
    }
    return result;
}

The code is also available on Github, in the revision corresponding to this question and in the current revision, maybe based on answers and comments.

\$\endgroup\$
  • \$\begingroup\$ Why not typedef the "base" type and make it so that you can easily profile which type works best on your current hardware? \$\endgroup\$ – mvds Jun 7 '17 at 1:08
  • \$\begingroup\$ It could be interesting to see if you can use SIMD instructions and/or registers, if performance is a concern. From my little experience in this field, long ago, I remember this: you can either hand-craft C code with SIMD intrinsics, or provide gentle nudges (in code) and flags to the compiler so that it starts using them. Such a nudge would be to make any array a multiple of 8 elements, and make that explicit in your code, e.g. so that a loop can process elements in groups of 8. (i.e. for ( int i=0;i<count&~7;i++ ) { } signals to the compiler that the loop always runs a multiple of 8 times) \$\endgroup\$ – mvds Jun 7 '17 at 1:20
  • \$\begingroup\$ @mvds: the typedefis already there in the newest revision :) Thanks for the further idea, this indeed sounds interesting! \$\endgroup\$ – Felix Palmen Jun 7 '17 at 7:28
1
\$\begingroup\$

Recommend separate allocation

Currently, you are using a trick where your data structure has variable length:

struct hugeint
{
    size_t s;
    size_t n;
    unsigned int e[];
};

I would recommend changing it to a fixed size structure with a separately allocated array:

struct hugeint
{
    size_t s;
    size_t n;
    unsigned int *e;
};

The reason I recommend this is because with your current data structure, you can only ever have one useful reference (i.e. pointer) to any hugeint, because modifying a hugeint requires reallocating the entire data structure. Any extra pointers to a hugeint become invalid the moment the hugeint is modified. For example:

hugeint *x = hugeint_parse(str);
hugeint *y = x;

hugeint_increment(&x);
// Now y is invalid, because x may have been reallocated

You might think that you would never need to have two pointers to the same hugeint, but I could imagine cases where it might be handy to be able to do so.

An extra benefit of this change is that it would remove one level of pointer indirection from many of your functions, leaving them more readable. For example, you would write self->e[i] instead of (*self)->e[i].

Answers to questions

Is there a better way to handle carries?

GMP uses something called "nails" to handle carries efficiently. Essentially, you would use only 31 out of each 32 bits for data, and reserve the upper bit for carries. That way, you can just add two 31 bit numbers and read the carry out of the top bit.

How do I pick the best type?

One way would be to always use uint64_t. If the hardware is 64-bit then you are fine, and if not you can assume that the compiler will generate efficient code to handle the 64-bit math (for example using add-with-carry instructions when adding 32-bit amounts).

Alternatively, you could also use uintptr_t, which is the size of a pointer. This typically should match the "register size" of your target. There may be cases where it doesn't, though.

\$\endgroup\$
  • \$\begingroup\$ This is indeed a good argument, thanks. My intention was to minimize dynamic allocations, but maybe I could just make the struct itself public... \$\endgroup\$ – Felix Palmen Jun 6 '17 at 9:11
1
\$\begingroup\$

Consider using "half-size" unsigned for the element size. This avoids the need for complex code to handle carries and overflow.

The following is certainly the work-horse loop for addition.

for (i = 0; i < other->n; ++i) {
    unsigned int v = (*self)->e[i] + other->e[i];
    unsigned int nextCarry = v < other->e[i];
    (*self)->e[i] = v + carry;
    if (!carry || (*self)->e[i]) carry = nextCarry;
}

If e[i] was half the width of unsigned, say as unsigned_half, then

unsigned carry = 0;
for (i = 0; i < other->n; ++i) {
    carry += (*self)->e[i] + other->e[i];
    (*self)->e[i] = (unsigned_half) carry;
    carry %= UNSIGNED_HALF_MAX + 1; 
}

I submit that there is a good chance that although the loop needs to iterate twice as often, the overall time will be less that the original code. Profiling will help determine that.

Further, for multiplying this will certainly be fastest to use half-sized than OP's unposted multiply routine.

// something like
unsigned_half c[a->n + b->n] = {0};
for (i = 0; i < a->n; ++i) {
    unsigned carry = 0;
    for (j = 0; j < b->n; ++j) {
      carry += 1u*a[i]*b[j] + c[i+j];
      c[i+j] = (unsigned_half) carry;
      carry %= UNSIGNED_HALF_MAX + 1; 
    }
    c[i+j] = (unsigned_half) carry;
}
\$\endgroup\$
  • \$\begingroup\$ Thanks for the answer. I'll probably do some experimentation concerning carry. The multiply routine is there (see the box with the whole code, hugeint_mult()), and it could be adapted to only use half the bits of the underlying type -- I'll check. \$\endgroup\$ – Felix Palmen Jun 7 '17 at 7:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.