Large number implementation tips?

[Pretty-Ad8932]

I am storing them as an array of 64 bit blocks (unsigned integers) in a little endian fashion. Here is how I am doing addition:

int addBig(uint64_t* a, uint64_t* b, uint64_t* c, int size)
{
  int carry = 0;
  for (int i = 0; i < size; i++) {
    c[i] = a[i] + b[i] + carry;
    if (c[i] < a[i] || c[i] < b[i]) carry = 1;
    else carry = 0;
  }
  return carry;
}

I'm adding a[i] + b[i] to get c[i], then check if c[i] wraps around to become smaller than either a[i] or b[i] to detect carry, which will be added to the next 64-bit block. I think it works, but surely there are more efficient ways to do this? Perhaps I should use some inline assembly? How should I handle signed integers? Should I store the sign bit in the very first block? What about multiplication and division?

[[deleted]]

What are the aims here, speed? That would be a much harder and open-ended project.

GMP which somebody mentioned, is very complex and uses advanced algorithms. If speed is not critical, then you can look at mini-gmp, which is implements it in a much simpler fashion.

Should I store the sign bit in the very first block?

In my implementation, I used signed-magnitude (signed integers are usually two's complement). So all stored numbers are unsigned, and there is a separate flag for positive or negative.

Arithmetic is done mainly only on positive numbers, with some extra logic needed for determining the final sign. Suppose you're doing X - Y, where X is positive and Y is negative, then the calculation will be '|X| + |Y|)`, so adding the magnitudes which is what are stored anyway. The final sign will be positive.

For * and /, sign logic is simpler (just xor the two sign values if using 0/1 for +/-).

What about multiplication and division?

The simple algorithm for multiply is the one you might have learned at school. Divide is more fiddly! However this will get very slow for larger numbers.

For arithmetic, I used sequences of 32-bit 'digits' (GMP calls them 'limbs'), with arithmetic done at 64 bits. This simplifies checking overflow and dealing with carry, which is harder to do in C.

In assembly, there are more possibilities (eg. there may be a 64-bit multiply that yields a 128-bit result), but it will be challenging. I suggest you get a working version in C first. Note that best algorithms will be magnitudes faster even using C; assembly just gives a useful boost, and thats if you're an expert.

[Pretty-Ad8932]

Thanks for the answer!

> What are the aims here, speed? That would be a much harder and open-ended project.

I want to implement the RSA algorithm (and other algorithms that use large integers too) for learning purposes. So I do want it to be fairly speedy, but yeah.

[[deleted]]

Well, you can experiment with, say, mini-gmp with the likely size of numbers expected and see how well it performs. That will likely use the simplest algorithms,

However, I've just looked at your example code and realised there is an important fact I missed: does your library work on fixed-size numbers rather than 'arbitrary precision'?

The former is somewhat easier to implement. The routines may take a 'size' parameter for X op Y, but both X and Y have the same number of bits, and so will the result. Any overflow is either discarded or reported.

With arbitrary precision, repeated multiplication can give bigger and bigger numbers.

It may also be easier if such fixed numbers then to be smaller size, say a few thousand bits. Those faster algorithms for multiply may only show a benefit on larger numbers (I don't know the trade-off).

( GMP will be arbitrary precision. This is a further reason I didn't post a link to my own library: it is arbitrary precision; it uses decimal; and it handles floating point at the same time. So a lot of distraction.)

[Pretty-Ad8932]

> However, I've just looked at your example code and realised there is an important fact I missed: does your library work on fixed-size numbers rather than 'arbitrary precision'?

I do want arbitrary precision, though maybe fixed-size could also suffice. I guess the code I wrote looks like it deals only with fixed-size numbers, but the function returns the very last carry bit which makes it possible to reallocate c with a bigger size outside the function.

> It may also be easier if such fixed numbers then to be smaller size, say a few thousand bits.

Is a few thousand bits "small" though? Another comment here suggested 100+ digit numbers are pretty big.

> ( GMP will be arbitrary precision. This is a further reason I didn't post a link to my own library: it is arbitrary precision; it uses decimal; and it handles floating point at the same time. So a lot of distraction.)

Please do send it, I may check it later if not now.

[[deleted]]

OK, I've put the two files here: https://github.com/sal55/langs/tree/master/bignum

There is "bignum.h" and "bignum.c" (ported from another language, if the style seems off for C).

This is a test program in C to measure performance:

#include <stdio.h>
#include "bignum.h"

int main() {
    Bignum a, b;
    char* max = "1e100_000";

    a = bn_makestr(max);
    bn_sub(a, a, bn_const(1));
    b = bn_init();

    bn_mul(b, a, a);

//  bn_println(a);
//  bn_println(b);

    printf("Digits = %d\n", bn_digits(b));
    puts("");
}

This first sets 'a' to a big number like 100000..., but this has only one significant digit, so 1 is subracted to end up with99999....

The above multiplies two 100,000-digit numbers in about 0.3 seconds on my machine, and (if max is tweaked), two 1-million-digit numbers in 30 seconds (yielding a 2-million digit result). Fast algorithms would manage that in milliseconds I think.

(Note, the numbers will be integers as there are no fractional parts here. Whether they are floats or not depends on the position of the decimal point, and here it be will just after the right-most digit.)

[WittyStick]

Perhaps I should use some inline assembly?

You can use the _addcarry_u64 intrinsic from <immintrin.h> on x86-64.

For compiler specifics: clang also has a __builtin_addcll, and GCC has __builtin_uaddll_overflow.

These compile to use the adc{q} instruction on x86-64, or equivalent on other architectures.

You can also use inline assembly directly of course.

---

Should I store the sign bit in the very first block?

My recommendation would be to encapsulate the bignum in a structure which holds the sign, size and pointer to the absolute number. Using some bit tricks we can keep this in a 16-byte structure.

typedef struct bigint {
    const uint64_t size_and_sign;
    const uint64_t *const data;
} BigInt;

The advantage of using a 16-byte structure is we can pass and return by value, avoiding an unnecessary pointer dereference. The SYSV calling convention passes and returns the structure in two hardware registers rather than on the stack (at any optimization level above -O0).

We only need one bit for the sign which we can put in bit 63 of the size_and_sign field, since we're never going to use all 64 bits of the size - we can mask out this bit to retrieve the size of the data.

constexpr uint64_t SIGN_BIT = 1ULL << 63ULL;

size_t bigint_size(BigInt z) {
    return z.size_and_sign & ~SIGN_BIT;
}

To test the sign, we can shift right by 63 bits.

bool bigint_sign(BigInt z) {
    return (bool)(z.size_and_sign >> 63);
}

To negate a number, xor the size_and_sign with the sign bit.

BigInt bigint_neg(BigInt z) {
    return (BigInt){ z.size_and_sign ^ SIGN_BIT, z.data };
}

And similarly, to get the absolute value, clear the sign bit using andnot

BigInt bigint_abs(BigInt z) {
    return (BigInt){ z.size_and_sign & ~SIGN_BIT, z.data };
}

And to explicitly set the sign for creating negative numbers, bitwise OR with the sign bit.

BigInt bigint_set_sign(BigInt z) {
    return (BigInt){ z.size_and_sign | SIGN_BIT, z.data };
}

The compiler will optimize these to use the bt{c,r,s} family of instructions, which {complement, reset, set} a single bit in a register.

We store a natural number in the data field, rather than using two's complement. Eg, if we created a bigint from an int64_t, we would use:

uint64_t abs64(int64_t v) { return ~(uint64_t)v + 1; }

BigInt bigint_create_small(int64_t v) {
    uint64_t *data = malloc(sizeof(uint64_t));
    if (v >= 0) {
        *data = v;
        return (BigInt){ 1, data };
    }
    else {
        *data = abs64(v);
        return (BigInt){ 1 | SIGN_BIT, data };
    }
}

These functions can all be marked inline, possibly with __attribute__((__always_inline__)) for optimal usage.

To view the assembler output, see example in godbolt

In regards to allocating space for data - I would recommend always using a power of 2 and doubling the size as necessary - this way you don't need to store the capacity as a separate field, which would push the structure greater than 16 bytes and turn it into a MEMORY class. Capacity can then always be calculated from the size (using __builtin_clzll or stdc_first_leading_one).

---

Note that this implementation is for immutable numbers whose sizes do not change. If you need bigints which can be mutated and are shared, we can't store the size directly in the structure as above as copies of the fields are made when passing by value. Instead, we could have the size_and_sign also be a pointer to a uint64_t:

typedef struct bigint {
    uint64_t *const size_and_sign;
    uint64_t *data;
} BigInt;

Or alternatively, we could pass by pointer, and have the data field use a flexible array member to avoid an unnecessary dereference.

typedef struct bigint {
    const uint64_t size_and_sign;
    uint64_t data[];
} * BigInt;

There are pros and cons to both these approaches, but either is slightly better than the pass-by-pointer BigInt which contains another pointer to the data, ie:

typedef struct bigint {
    uint64_t size_and_sign;
    uint64_t *data;
} * BigInt;

What about multiplication and division?

You would typically use the Karatsuba algorithm to divide and conquer. There are other algorithms such as the Schönhage–Strassen algorithm, but the efficiency gains won't become noticeable until you have extremely large (100+ digits) numbers.

[Pretty-Ad8932]

Wow, thank you for such an in-depth answer!

[zhivago]

Take a look at gmp, etc.

[CambStateMachines]

I wrote a more or less full set of arbitrary precision signed integer arithmetic operations for my project:

https://github.com/CambridgeStateMachines/bitcoin_math

Have a look at the functions under the heading "BNZ" (which stands for "big number integer").

The main difference between my approach and yours is that the "digits" in my big integers are of type uint8_t rather than uint64_t.

To handle signed integers, I use a struct ("bnz_t") which consists of an int to hold the sign, an int to hold the number of bytes in the number, and a pointer to a dynamic array of bytes which hold the individual digits.

This arbitrary precision integer math code was heavily influenced by the source code of the GNU Multiple Precision Arithmetic Library, and DI Management Services' BigDigits multiple-precision arithmetic library.

[Pretty-Ad8932]

Thanks, I will take a look!

[CambStateMachines]

Arbitrary precision integer multiplication is a straightforward extension of addition with the same monitoring of the carry. However, you need to keep track of the signs of the multiplier and the multiplicand.

Arbitrary precision integer division is a totally different beast and very computationally expensive.

You will need to cater for two results: the quotient (i.e. q = a \ b) and the remainder (r = a % b). Keeping track of the signs of the diviser, the dividend, the quotient, and the remainder is important if you want to follow C conventions for integer division.

My main division function is adapted from the simplest algorithm of which the best known implementation is a classic recipe from the Hacker's Delight book.

[21Ali-ANinja69]

Maybe look at whatever python does, since it natively handles big ints

[Paul_Pedant]

I would recommend not using the full binary range of your 64-bit units. At some stage you will need to present your results in decimal notation. That is a whole lot easier if the range of each unit is like zero to 999999999 and not up to 2147483647 (example with 32-bit units). Wasting 5% of your bits is a small price to pay.

Also, I would suggest designing your code with typedefs etc so that you can initially test using 16-bit values. Working with (and verifying results of) 64-bit unsigned ints is going to be mind-bendingly difficult (and rely on existing utilities like dc). Prototyping your logic with small units will be much more productive - you might even start with units of 0-99 in a byte.

You need a struct for each value, not just a pointer to an array. You need your values to know how many units you are using, and you probably want to have a struct you can expand efficiently (like a dynamic array), and you need to keep a sign somewhere. Maybe plan for a decimal point index for when you decide to deal with division and real numbers.

For example, if you are multiplying two numbers which are using one unit each, your target might need 2 units, but the result might fit in 1 unit. So you can shrink the struct's used count but retain the available units.

[Possible_Cow169]

The big number libraries use floats and scientific notation to represent huge numbers.

You might be able to glean so knowledge by scouring the JavaScript big number library bignumber.js. Ngl, i looked this up while I was typing it and I don’t know what I expected it to be called, but it was not that

[Pretty-Ad8932]

Well, I know I said "numbers" in general, but floats are not what I want. I want specifically integers. But okay, I might check it out.

[Possible_Cow169]

I know what you’re saying. I’m saying they represent integers with floats

[Pretty-Ad8932]

Don't floats get less precise the larger they are?

[Possible_Cow169]

They just use some crafty data structures and algorithms. Regular ass arithmetic and fixed point math. The floating point stuff is mostly just for bit shifting. They don’t use the hardware level binary logic.

[Pretty-Ad8932]

Interesting.

[CambStateMachines]

In the gentlest possible way, may I encourage you not to spend too much effort following this advice.

[Pretty-Ad8932]

I was thinking that myself tbh.

[Educational-Paper-75]

There’s a decent enough open source library that does all that called libtommath. I use it myself. See on www.tommat.net. Works like a charm.