Math128
[Numeric: Rational Number Handling with Rounding Error Control]

---

Detailed Description

Quick-n-dirty 128-bit integer math lib. Things seem to mostly work, and have been tested, but not comprehensively tested.

Data Structures
struct	QofInt128
Functions
gboolean	equal128 (QofInt128 a, QofInt128 b)
gint	cmp128 (QofInt128 a, QofInt128 b)
QofInt128	shift128 (QofInt128 x)
QofInt128	shiftleft128 (QofInt128 x)
QofInt128	inc128 (QofInt128 a)
QofInt128	add128 (QofInt128 a, QofInt128 b)
QofInt128	mult128 (gint64 a, gint64 b)
QofInt128	div128 (QofInt128 n, gint64 d)
gint64	rem128 (QofInt128 n, gint64 d)
guint64	gcf64 (guint64 num, guint64 denom)
QofInt128	lcm128 (guint64 a, guint64 b)

---

Function Documentation

QofInt128 add128	(	QofInt128	a,
		QofInt128	b
	)			[inline]

Add a pair of 128-bit numbers, returning a 128-bit number

Definition at line 308 of file qofmath128.c.

{
    QofInt128 sum;
    if (a.isneg == b.isneg)
    {
        sum.isneg = a.isneg;
        sum.hi = a.hi + b.hi;
        sum.lo = a.lo + b.lo;
        if ((sum.lo < a.lo) || (sum.lo < b.lo))
        {
            sum.hi++;
        }
        sum.isbig = sum.hi || (sum.lo >> 63);
        return sum;
    }
    if ((b.hi > a.hi) || ((b.hi == a.hi) && (b.lo > a.lo)))
    {
        QofInt128 tmp = a;
        a = b;
        b = tmp;
    }

    sum.isneg = a.isneg;
    sum.hi = a.hi - b.hi;
    sum.lo = a.lo - b.lo;

    if (sum.lo > a.lo)
    {
        sum.hi--;
    }

    sum.isbig = sum.hi || (sum.lo >> 63);
    return sum;
}

gint cmp128	(	QofInt128	a,
		QofInt128	b
	)			[inline]

Return returns 1 if a>b, -1 if b>a, 0 if a == b

Definition at line 246 of file qofmath128.c.

{
    if ((0 == a.isneg) && b.isneg)
        return 1;
    if (a.isneg && (0 == b.isneg))
        return -1;
    if (0 == a.isneg)
    {
        if (a.hi > b.hi)
            return 1;
        if (a.hi < b.hi)
            return -1;
        if (a.lo > b.lo)
            return 1;
        if (a.lo < b.lo)
            return -1;
        return 0;
    }

    if (a.hi > b.hi)
        return -1;
    if (a.hi < b.hi)
        return 1;
    if (a.lo > b.lo)
        return -1;
    if (a.lo < b.lo)
        return 1;
    return 0;
}

QofInt128 div128	(	QofInt128	n,
		gint64	d
	)			[inline]

Divide a signed 128-bit number by a signed 64-bit, returning a signed 128-bit number.

Definition at line 180 of file qofmath128.c.

{
    QofInt128 quotient;
    int i;
    gint64 remainder = 0;

    quotient = n;
    if (0 > d)
    {
        d = -d;
        quotient.isneg = !quotient.isneg;
    }

    /* Use grade-school long division algorithm */
    for (i = 0; i < 128; i++)
    {
        guint64 sbit = HIBIT & quotient.hi;
        remainder <<= 1;
        if (sbit)
            remainder |= 1;
        quotient = shiftleft128 (quotient);
        if (remainder >= d)
        {
            remainder -= d;
            quotient.lo |= 1;
        }
    }

    /* compute the carry situation */
    quotient.isbig = (quotient.hi || (quotient.lo >> 63));

    return quotient;
}

gboolean equal128	(	QofInt128	a,
		QofInt128	b
	)			[inline]

Return true of two numbers are equal

Definition at line 233 of file qofmath128.c.

{
    if (a.lo != b.lo)
        return 0;
    if (a.hi != b.hi)
        return 0;
    if (a.isneg != b.isneg)
        return 0;
    return 1;
}

guint64 gcf64	(	guint64	num,
		guint64	denom
	)			[inline]

Return the greatest common factor of two 64-bit numbers

Definition at line 278 of file qofmath128.c.

{
    guint64 t;

    t = num % denom;
    num = denom;
    denom = t;

    /* The strategy is to use Euclid's algorithm */
    while (0 != denom)
    {
        t = num % denom;
        num = denom;
        denom = t;
    }
    /* num now holds the GCD (Greatest Common Divisor) */
    return num;
}

QofInt128 inc128

(

QofInt128

a

)

[inline]

Increment by one

increment a 128-bit number by one

Definition at line 153 of file qofmath128.c.

{
    if (0 == a.isneg)
    {
        a.lo++;
        if (0 == a.lo)
        {
            a.hi++;
        }
    }
    else
    {
        if (0 == a.lo)
        {
            a.hi--;
        }
        a.lo--;
    }

    a.isbig = (a.hi != 0) || (a.lo >> 63);
    return a;
}

QofInt128 lcm128	(	guint64	a,
		guint64	b
	)			[inline]

Return the least common multiple of two 64-bit numbers.

Definition at line 299 of file qofmath128.c.

{
    guint64 gcf = gcf64 (a, b);
    b /= gcf;
    return mult128 (a, b);
}

QofInt128 mult128	(	gint64	a,
		gint64	b
	)			[inline]

Multiply a pair of signed 64-bit numbers, returning a signed 128-bit number.

Definition at line 41 of file qofmath128.c.

{
    QofInt128 prod;
    guint64 a0, a1;
    guint64 b0, b1;
    guint64 d, d0, d1;
    guint64 e, e0, e1;
    guint64 f, f0, f1;
    guint64 g, g0, g1;
    guint64 sum, carry, roll, pmax;

    prod.isneg = 0;
    if (0 > a)
    {
        prod.isneg = !prod.isneg;
        a = -a;
    }

    if (0 > b)
    {
        prod.isneg = !prod.isneg;
        b = -b;
    }

    a1 = a >> 32;
    a0 = a - (a1 << 32);

    b1 = b >> 32;
    b0 = b - (b1 << 32);

    d = a0 * b0;
    d1 = d >> 32;
    d0 = d - (d1 << 32);

    e = a0 * b1;
    e1 = e >> 32;
    e0 = e - (e1 << 32);

    f = a1 * b0;
    f1 = f >> 32;
    f0 = f - (f1 << 32);

    g = a1 * b1;
    g1 = g >> 32;
    g0 = g - (g1 << 32);

    sum = d1 + e0 + f0;
    carry = 0;
    /* Can't say 1<<32 cause cpp will goof it up; 1ULL<<32 might work */
    roll = 1 << 30;
    roll <<= 2;

    pmax = roll - 1;
    while (pmax < sum)
    {
        sum -= roll;
        carry++;
    }

    prod.lo = d0 + (sum << 32);
    prod.hi = carry + e1 + f1 + g0 + (g1 << 32);
    // prod.isbig = (prod.hi || (sum >> 31));
    prod.isbig = prod.hi || (prod.lo >> 63);

    return prod;
}

gint64 rem128	(	QofInt128	n,
		gint64	d
	)			[inline]

Return the remainder of a signed 128-bit number modulo a signed 64-bit. That is, return nd in 128-bit math. I beleive that ths algo is overflow-free, but should be audited some more ...

Definition at line 220 of file qofmath128.c.

{
    QofInt128 quotient = div128 (n, d);

    QofInt128 mu = mult128 (quotient.lo, d);

    gint64 nn = 0x7fffffffffffffffULL & n.lo;
    gint64 rr = 0x7fffffffffffffffULL & mu.lo;
    return nn - rr;
}

QofInt128 shift128

(

QofInt128

x

)

[inline]

Shift right by one bit (i.e. divide by two)

Definition at line 110 of file qofmath128.c.

{
    guint64 sbit = x.hi & 0x1;
    x.hi >>= 1;
    x.lo >>= 1;
    x.isbig = 0;
    if (sbit)
    {
        x.lo |= HIBIT;
        x.isbig = 1;
        return x;
    }
    if (x.hi)
    {
        x.isbig = 1;
    }
    return x;
}

QofInt128 shiftleft128

(

QofInt128

x

)

[inline]

Shift left by one bit (i.e. multiply by two)

Shift leftt by one bit (i.e. multiply by two)

Definition at line 131 of file qofmath128.c.

{
    guint64 sbit;
    sbit = x.lo & HIBIT;
    x.hi <<= 1;
    x.lo <<= 1;
    x.isbig = 0;
    if (sbit)
    {
        x.hi |= 1;
        x.isbig = 1;
        return x;
    }
    if (x.hi)
    {
        x.isbig = 1;
    }
    return x;
}

---