6 ms·

Optimizing a bignum library for fun

145 points by azhenley 2y ago

styczen 2y ago
not + -
minimize 2y ago
For maximum efficiency, you should work in binary instead of base 10. Handling carries becomes more straightforward with the right primitives, for example __builtin_addc with GCC: https://gcc.gnu.org/onlinedocs/gcc/Integer-Overflow-Builtins.html https://gcc.gnu.org/onlinedocs/gcc/Integer-Overflow-Builtins... You can also implement it in C if you want a more portable solution: https://github.com/983/bigint/blob/ee0834c65a27d18fa628e6c526a0b83b24db90f9/bigint.c#L27 https://github.com/983/bigint/blob/ee0834c65a27d18fa628e6c52... If you scroll around, you can also find my implementations for multiplication and such.
- rurban 2y ago
  For maximum efficiency you'd avoid addc like hell, because it blocks internal precomputation, and use guarantees like he did, avoiding overflow better. Better use int128 types though. I'd just just stack objects for small typical sizes. And esp. do formal proofs than the random test joke.
  - minimize 2y ago
    > For maximum efficiency you'd avoid addc like hell, because it blocks internal precomputation, and use guarantees like he did, avoiding overflow better. Whether that is the case depends on the CPU. I hit memory bandwidth limits before adc becomes a problem. But I concede that you have a valid point and that there are definitely cases where leaving sufficient space for carries is the better strategy. Anyway, we can probably both agree that "% 1000000000" is not the best way to do it.
  - koverstreet 2y ago
    addc makes all your adds serialize on the flags register, which is really painful. more modern approach is to reserve some bits out of every word for carries, but that drastically complicates things.
    - JonChesterfield 2y ago
      I remembered and ultimately found a source for a workaround for the serialising on flags problem, intel paper at https://web.archive.org/web/20150131061304/http://www.intel.com/content/dam/www/public/us/en/documents/white-papers/ia-large-integer-arithmetic-paper.pdf https://web.archive.org/web/20150131061304/http://www.intel.... amounts to new instructions with better behaviour for ILP
    - haberman 2y ago
      Is this really true? I would intuitively expect that register renaming would apply to eflags too, so that reads from flags don't truly need to be serialized despite nominally writing a bunch of things to the same register. EDIT: this paper (linked in another comment) seems to indicate that this is possible: > An out-of-order machine can look ahead and process the accumulation pass in parallel with the partial sum pass using a renamed eFlags register. https://web.archive.org/web/20150131061304/http://www.intel.com/content/dam/www/public/us/en/documents/white-papers/ia-large-integer-arithmetic-paper.pdf https://web.archive.org/web/20150131061304/http://www.intel....
      - phire 2y ago
        EFLAGS is actually put in the same renaming register as the result, so you get renaming for free. The renaming registers in the Intel Pentium Pro and Pentium II are actually over 80 bits wide. They need to hold a full 80bit float, or 64bit MMX result. The Pentium III extends this to 128bit wide renaming registers to support SSE. This is despite the fact that the P6 architecture only had 32bit bit integer registers until the Core 2 in 2006. So there is plenty of room to store EFLAGS in the same renaming register as the result. This also means that the branch uops point to the result of the most recent flag modifying instruction. It was only with Sandybridge (and the introduction of AVX) that the P6 switched to a PRF design, with separate registers for floats and integers. Of course, Netburst also had a PRF design.
paldepind2 2y ago
Speaking of bignum libraries, I recently watched a talk with Rob Pike where he mentioned that one thing he regretted about Go was not making the default integer implementation arbitrary precision. Supposedly the performance overhead for normal numbers is very small, and you avoid the weirdness and complicated semantics of fixed precision integers. I found that to be quite fascinating, especially coming from a "low-level guy" like Rob Pike. Ever since I've been wanting a language with that feature and to understand how bignum implementations work.
- MaxBarraclough 2y ago
  For what it's worth, Haskell's Integer type is that way. edit To be clear, in Haskell you 'default' to using its arbitrary-precision type, unlike say Java where they're available in the standard library but not in the core primitive types.
- vient 2y ago
  > I've been wanting a language with that feature Python, which realization CPython is mentioned in the article, has arbitrary precision integers.
- omoikane 2y ago
  > I've been wanting a language with that feature GNU MP's C++ interface makes this mostly transparent, expressions like "x = y * z" just work due to operator overloading. The only time you need to know that the underlying class is "mpz_class" instead of "int" is at input/output.
- pfdietz 2y ago
  Common Lisp has been like this forever. It makes testing the language much easier, as you don't have to check for overflows when constructing test code. A dynamic language like CL has the added requirement that the integers fit with everything else under a single top type, which in practice limits fixnums to a few bits less than the word size.
  - globular-toast 2y ago
worstspotgain 2y ago
I've had to implement BigNum on an embedded system that had very little RAM and where the initial optimized version of ModExp took many minutes to complete. After much hair-pulling, the final version took 40 seconds. First, you should work with unsigned numbers, and use a power of 2 as your word size. The fastest choice for a word size is very operation- and CPU-dependent. A key trick is to lay out your bignums in the same order as the endianity of each word in memory, e.g. least-significant word first on little-endian systems. This will allow you to choose your word size dynamically for each operation: in memory, a number with M words of N bits each is identical to a number with M / 2 words of N * 2 bits each. For multiplication, identify the CPU instruction with the widest result, then use half that size as your word size. Each step through the arrays generates a word result in the low half and a word carry in the top half. The carry gets added to the result of the next step, possibly overflowing. For addition, use the widest result as your word size. This can also overflow. How you deal with overflows is very CPU-dependent. You can use adc/addc as someone else mentioned, which will be faster on embedded and may be faster on fatter chips. Alternatively, you can halve the word size and use the top half as the carry. If addc is not available, you can test for overflows as follows: uint32_t a = ..., b = ...; uint32_t res = a + b; uint32_t carry = res < a; On overflow, res must necessarily be less than both a and b, so no need to check b. If SIMD instructions are available, that will almost always be the fastest choice by far. While it doesn't change the above guidelines in principle, there are often e.g. optimized overflow mechanisms.
- less_less 2y ago
  """ uint32_t a = ..., b = ...; uint32_t res = a + b; uint32_t carry = res < a; """ But note well: this trick doesn't work when there is also a carry-in. That is, uint32_t res = a + b + carry; carry = res < a; won't work when b = 0xFFFFFFFF and carry = 1, so you have to check for carries twice, or set `carry = res < a || (res == a && carry);` or similar. So if you are implementing bignums on an architecture that doesn't have ADC (such as RISC-V), or on a vector architecture where the performance would be bad, then you often need an alternative strategy such as using less than 32/64 bits per limb (whether that's half-size, or 30 bits or 51 bits or whatever).
nj5rq 2y ago
Fascinating article, I have always been wondering how these big number libraries worked. As a side question, does anyone the program that the author used when making that "addition" and "multiplication" performance graph? Thanks.
- nerdponx 2y ago
  I don't know for sure, but you can do that kind of thing pretty easily with Matplotlib in Python. Or in R base graphics, with more effort to get it looking pretty.
- jaeckel 2y ago
  If you're interested in even more details you could also have a look at Tom's book about bignum arithmetics. C.f. third paragraph of https://github.com/libtom/libtommath#summary https://github.com/libtom/libtommath#summary The library evolved since then, so the algorithm documentation only exists in code nowadays. FTR I'm involved in the project as maintainer.
kazinator 2y ago
LOL, I've been there: https://www.kylheku.com/cgit/txr/commit/mpi?id=98dedd310b1d5d876b0fbb0ebd6c4df9bd7b2d88 https://www.kylheku.com/cgit/txr/commit/mpi?id=98dedd310b1d5... (Patch was originally from 2011; it was bugfixed once, and then in 2015 converted to git commit: https://www.kylheku.com/cgit/txr/commit/mpi-patches/faster-square-root?id=ef47dfe4fcb7c1be369ae83221386b9da6474a1e https://www.kylheku.com/cgit/txr/commit/mpi-patches/faster-s... ) Another one: faster highest-bit search: https://www.kylheku.com/cgit/txr/tree/mpi-patches/bit-search-optimizations?id=124e7dd6977a0853d7a8399921e31fd1ccde2dcb https://www.kylheku.com/cgit/txr/tree/mpi-patches/bit-search... That does use GCC built-ins today: https://www.kylheku.com/cgit/txr/commit/?id=15b7c542dc44899e8db7addfcc2f1c1c4a188b49 https://www.kylheku.com/cgit/txr/commit/?id=15b7c542dc44899e...
haberman 2y ago
Coincidentally I was just writing a bignum library from scratch two weeks ago. A few interesting things I learned: 1. Knuth volume 2 has an extensive discussion of the problem space. I've only gotten a chance to skim it so far, but it looks interesting and approachable. 2. I need to support bitwise operations, which operate on a two's complement representation, so I figured it would be simpler to use two's complement internally, despite seeing that most (all?) bignum libraries use signed magnitude. I'm starting to regret this: two's complement introduces a lot of complexity. The most fun thing about working on a bignum library is that it makes the algorithms you learned in grade school for add/subtract/multiply/divide relevant again. The basic ("classical") algorithms on bignum are basically exactly the same thing you learned in grade school, except on a much larger base than base-10.
- vintermann 2y ago
  Knuth has a "hugenum" scheme, TCALC, which is based on trees rather than byte arrays. There are many ways to do it, depending on which type of numbers you need most - for TCALC, it's numbers close to towers of exponents. Of course for every big number encoding that can beat an array of bytes for some numbers, it also has to lose to byte arrays for some numbers... But it's possible to bound this inefficiency, so that it's never more than twice as bad as byte arrays, but potentially exponentially more efficient for the favoured numbers.
  - lysp 2y ago
    Twice as bad being a timeout and 2nd call?
    - vintermann 2y ago
      I don't quite get what you mean. Computation time is also bounded with the size of the tree.
tgot 2y ago
I think that your description is almost excellent, but that you're fundamentally misleading in describing what you are doing as a "30-bit" digit. It's a 10^9 digit mathematically, occupying 30 bits of storage. You do briefly mention that it's 10^9, but repeatedly say 30-bits.
- samatman 2y ago
  It isn't misleading at all. A hexadecimal digit has 4 bits of entropy. You can guess it correctly by chance one in sixteen times. Calling that a four-bit digit is correct. Same with a 30 bit digit, all that changes is the magnitude. The "storage bits" aren't a detail, they're an essential property. Stored in balanced ternary it would still have 30 bits of entropy.
  - derdi 2y ago
    I think you misunderstood the parent. If you speak of a number in the range 0-15 stored in 4 bits, we can all agree that "4 bit digit" is the appropriate term for it. But what about a number also stored in 4 bits but restricted to the range 0-13? It's a digit that fits in 4 bits, but calling it a "4 bit digit" without further qualification would omit relevant information.
    - samatman 2y ago
      That hadn't occurred to me, I had thought (having read the Fine Article, but not closely, nor really more than glanced at the code) that it was using exactly 30 bits and reserving two for, IDK, carry detection, and that 10^9 was just approximating 1073741824. We do that every time we call 1024 a KB and not a KiB, so there's precedent. Ok. Fair. It's a ~28.9 bit digit taking up 32 bits of storage. It turns out I was mislead!
parentheses 2y ago
I am so surprised that there's no exploration of Karatsuba's algorithm. That's what makes the Python implementation perform. I actually came here hoping to find discussion on Karatsuba. https://en.m.wikipedia.org/wiki/Karatsuba_algorithm https://en.m.wikipedia.org/wiki/Karatsuba_algorithm
- jmts 2y ago
  Please see the section titled "Faster multiplication"
globular-toast 2y ago
Funny, I did this myself 10 years ago. Shit, it's been that long..? For my undergrad project I wrote a computer algebra system to do symbolic integration. The supervisor was a hardcore, old school C guy, so naturally I was going to just use C and no libraries. He told me I'd need bignums first, so I got to work (this is because many algorithms like polynomial GCD create massive numbers as they go, even if the inputs and final outputs are generally very small). I just couldn't figure out how to do better than the largest power of 10 per digit at the time. Working with non base 10 arithmetic was a mind fuck for me at the time. So I did it with digits holding 10^9 and the classical algorithms from Knuth. Division is the hardest! At some point I discovered the GNU multiple precision library (GMP) and made my program work with that instead of mine. I was shocked at how much faster GMP was! I finished my project with my own code, but I knew I had to come back to do it better. The breakthrough came when I got a copy of Hacker's Delight. It has stuff like how to detect overflow after it's happened (in C). Something twigged and then I just understood how to fill each word completely rather than use a power of 10. I don't know what confused me before. But, of course, the real way to do it is to use assembly. You can't get close to high performance in C alone. In assembly you get the overflow bit. It's actually easier in a way! So you write tiny platform specific bits for the digits and build on that in C. My add and subtract were then as fast as GMP. I lost interest when it came to implement faster multiplication algorithms. Code in case anyone is interested: https://github.com/georgek/bignums https://github.com/georgek/bignums
chubot 2y ago
Very useful post! Also it's cool to see how many people in this thread have worked on this problem -- lots of new info here I haven't seen I wonder if anyone is interested in implementing a big numbers in Oils? It's a Unix shell with TWO complete implementations - the "executable spec" in Python, and an automatic translation to pure C++ (which is 30x-50x faster) We currently use 64-bit integers in C++, but big nums are a better semantic. Some trivia about bad shell semantics here: Integers - Don't do whatever Python or C++ does - https://www.oilshell.org/blog/2024/03/release-0.21.0.html#integers-dont-do-whatever-python-or-c-does https://www.oilshell.org/blog/2024/03/release-0.21.0.html#in... This is a very self-contained project: the interface is defined by 200 lines of Python: https://github.com/oilshell/oil/blob/master/mycpp/mops.py https://github.com/oilshell/oil/blob/master/mycpp/mops.py and the trivial 64-bit overflowing implementation is also about 200 lines: https://github.com/oilshell/oil/blob/master/mycpp/gc_mops.h https://github.com/oilshell/oil/blob/master/mycpp/gc_mops.h https://github.com/oilshell/oil/blob/master/mycpp/gc_mops.cc https://github.com/oilshell/oil/blob/master/mycpp/gc_mops.cc (We have a fast Ninja-based build system, so you can probably iterate on this in 100 milliseconds or less -- it should be fun for the right person) --- I think the main reason it is specific to Oils is that the bigger number should become GC objects. Details on our GC here: Pictures of a Working Garbage Collector - https://www.oilshell.org/blog/2023/01/garbage-collector.html https://www.oilshell.org/blog/2023/01/garbage-collector.html It's been very solid for the last 18 months, basically because it's well tested by ASAN and #ifdef testing modes. The main thing I'd be concerned with is how to TEST that big number operations are correct. I think there are probably some interesting strategies, which I'd love to discuss. You're of course welcome to adapt existing open source code, including code you've already written -- I probably even prefer that, i.e. something that has had some real world testing. We want all the operations in that file, and it should be integrated with our GC. --- We've had ~6 contributors funded by grants from https://nlnet.nl https://nlnet.nl for the past couple years, so you can even be paid (there's a chance it depends on the country you live in, but generally a wide range of situations is OK). Contact me at andy at oilshell.org or https://oilshell.zulipchat.com/ https://oilshell.zulipchat.com/ if interested!
- JonChesterfield 2y ago
  It doesn't look especially specific to me - slightly hack the structure in https://github.com/libtom/libtommath https://github.com/libtom/libtommath to make sense to your garbage collector and you're pretty much done. That somewhat punts on the testing question but I doubt add/mul/bitops from that library are buggy.
amelius 2y ago
The title is a little bit dishonest, because they are optimizing their own bignum library.
Kimberrly 2y ago
[dead]
booleandilemma 2y ago
Except he's not doing it for fun at all. He's doing it for clout and publicity on HN.
- JonChesterfield 2y ago
  Looks like play to me. The earlier post doesn't have a HN reference and there's not much correlation with their publications list. People legitimately do write things like this as entertainment.
guyomes 2y ago
For a small survey of practical efficient methods for bignum arithmetic operations, the Algorithms section of the documentation of GMP [1] is excellent. [1]: https://gmplib.org/manual/Algorithms https://gmplib.org/manual/Algorithms
lifthrasiir 2y ago
Modulo is surprisingly expensive even when you combine it with a quotient. It is almost always better to use binary "limbs", in this case 31 or 32 bits wide, because decimal parsing and printing should be much rarer than individual operations in general.
- colonwqbang 2y ago
  This is a bit misleading. Quotient and modulo with constant right-hand side is readily optimised into faster operations, like multiplication and bit shifts. "All" optimising compilers do this and don't emit a true divide instruction in such cases. See for instance: https://godbolt.org/z/vxnTPvY6M https://godbolt.org/z/vxnTPvY6M Quotient and modulo with a variable RHS is expensive, but this isn't necessary for the algorithm showed in the article. Although I agree that binary wins over decimal anyway. We are using binary computers after all, for the last 50 years or so.
  - Someone 2y ago
    > Quotient and modulo with a variable RHS is expensive For dividing ‘big’ bignums by ‘small’ divisors, I would use/tweak libdivide (https://libdivide.com/ https://libdivide.com/), and spend the time to optimize the division. For dividing ‘big’ bignums by ‘big’ divisors, that may not be worth it. (Determining the correct values of the various ‘big’ and ‘small’ values left as an exercise for the reader. Those will depend on the hardware at hand)
  - worstspotgain 2y ago
    The reason a division by 1000000 can be efficiently turned into a multiplication is because you're targeting a system that has a 32x32->64 multiplication instruction, which has plenty of room to spare. But if you have that, you would be better off using e.g. 2^32 as your modulo.
    - colonwqbang 2y ago
      I agree that 2^32 is better, but the reason is not that division is too expensive because OPs algorithm after standard optimisation does not need any division. Clearly choosing a base type that is fast on the machine, is key. But I feel even a 64/32 divide with constant operand will be optimised into something reasonably fast on most systems. If you have a counterexample it would be interesting to hear. A system without fast divide likely doesn't have any floats either. So all arithmetic would be fixed point and 32x32->64 multiplication would be a pretty common operation.

145 points by azhenley 2y ago

not + -

For maximum efficiency, you should work in binary instead of base 10. Handling carries becomes more straightforward with the right primitives, for example __builtin_addc with GCC: https://gcc.gnu.org/onlinedocs/gcc/Integer-Overflow-Builtins.html https://gcc.gnu.org/onlinedocs/gcc/Integer-Overflow-Builtins... You can also implement it in C if you want a more portable solution: https://github.com/983/bigint/blob/ee0834c65a27d18fa628e6c526a0b83b24db90f9/bigint.c#L27 https://github.com/983/bigint/blob/ee0834c65a27d18fa628e6c52... If you scroll around, you can also find my implementations for multiplication and such.

rurban 2y ago
For maximum efficiency you'd avoid addc like hell, because it blocks internal precomputation, and use guarantees like he did, avoiding overflow better. Better use int128 types though. I'd just just stack objects for small typical sizes. And esp. do formal proofs than the random test joke.
- minimize 2y ago
  > For maximum efficiency you'd avoid addc like hell, because it blocks internal precomputation, and use guarantees like he did, avoiding overflow better. Whether that is the case depends on the CPU. I hit memory bandwidth limits before adc becomes a problem. But I concede that you have a valid point and that there are definitely cases where leaving sufficient space for carries is the better strategy. Anyway, we can probably both agree that "% 1000000000" is not the best way to do it.
- koverstreet 2y ago
  addc makes all your adds serialize on the flags register, which is really painful. more modern approach is to reserve some bits out of every word for carries, but that drastically complicates things.
  - JonChesterfield 2y ago
    I remembered and ultimately found a source for a workaround for the serialising on flags problem, intel paper at https://web.archive.org/web/20150131061304/http://www.intel.com/content/dam/www/public/us/en/documents/white-papers/ia-large-integer-arithmetic-paper.pdf https://web.archive.org/web/20150131061304/http://www.intel.... amounts to new instructions with better behaviour for ILP
  - haberman 2y ago
    Is this really true? I would intuitively expect that register renaming would apply to eflags too, so that reads from flags don't truly need to be serialized despite nominally writing a bunch of things to the same register. EDIT: this paper (linked in another comment) seems to indicate that this is possible: > An out-of-order machine can look ahead and process the accumulation pass in parallel with the partial sum pass using a renamed eFlags register. https://web.archive.org/web/20150131061304/http://www.intel.com/content/dam/www/public/us/en/documents/white-papers/ia-large-integer-arithmetic-paper.pdf https://web.archive.org/web/20150131061304/http://www.intel....
    - phire 2y ago
      EFLAGS is actually put in the same renaming register as the result, so you get renaming for free. The renaming registers in the Intel Pentium Pro and Pentium II are actually over 80 bits wide. They need to hold a full 80bit float, or 64bit MMX result. The Pentium III extends this to 128bit wide renaming registers to support SSE. This is despite the fact that the P6 architecture only had 32bit bit integer registers until the Core 2 in 2006. So there is plenty of room to store EFLAGS in the same renaming register as the result. This also means that the branch uops point to the result of the most recent flag modifying instruction. It was only with Sandybridge (and the introduction of AVX) that the P6 switched to a PRF design, with separate registers for floats and integers. Of course, Netburst also had a PRF design.

paldepind2 2y ago

Speaking of bignum libraries, I recently watched a talk with Rob Pike where he mentioned that one thing he regretted about Go was not making the default integer implementation arbitrary precision. Supposedly the performance overhead for normal numbers is very small, and you avoid the weirdness and complicated semantics of fixed precision integers. I found that to be quite fascinating, especially coming from a "low-level guy" like Rob Pike. Ever since I've been wanting a language with that feature and to understand how bignum implementations work.

MaxBarraclough 2y ago
For what it's worth, Haskell's Integer type is that way. edit To be clear, in Haskell you 'default' to using its arbitrary-precision type, unlike say Java where they're available in the standard library but not in the core primitive types.
vient 2y ago
> I've been wanting a language with that feature Python, which realization CPython is mentioned in the article, has arbitrary precision integers.
omoikane 2y ago
> I've been wanting a language with that feature GNU MP's C++ interface makes this mostly transparent, expressions like "x = y * z" just work due to operator overloading. The only time you need to know that the underlying class is "mpz_class" instead of "int" is at input/output.
pfdietz 2y ago
Common Lisp has been like this forever. It makes testing the language much easier, as you don't have to check for overflows when constructing test code. A dynamic language like CL has the added requirement that the integers fit with everything else under a single top type, which in practice limits fixnums to a few bits less than the word size.
- globular-toast 2y ago

worstspotgain 2y ago

I've had to implement BigNum on an embedded system that had very little RAM and where the initial optimized version of ModExp took many minutes to complete. After much hair-pulling, the final version took 40 seconds. First, you should work with unsigned numbers, and use a power of 2 as your word size. The fastest choice for a word size is very operation- and CPU-dependent. A key trick is to lay out your bignums in the same order as the endianity of each word in memory, e.g. least-significant word first on little-endian systems. This will allow you to choose your word size dynamically for each operation: in memory, a number with M words of N bits each is identical to a number with M / 2 words of N * 2 bits each. For multiplication, identify the CPU instruction with the widest result, then use half that size as your word size. Each step through the arrays generates a word result in the low half and a word carry in the top half. The carry gets added to the result of the next step, possibly overflowing. For addition, use the widest result as your word size. This can also overflow. How you deal with overflows is very CPU-dependent. You can use adc/addc as someone else mentioned, which will be faster on embedded and may be faster on fatter chips. Alternatively, you can halve the word size and use the top half as the carry. If addc is not available, you can test for overflows as follows: uint32_t a = ..., b = ...; uint32_t res = a + b; uint32_t carry = res < a; On overflow, res must necessarily be less than both a and b, so no need to check b. If SIMD instructions are available, that will almost always be the fastest choice by far. While it doesn't change the above guidelines in principle, there are often e.g. optimized overflow mechanisms.

less_less 2y ago
""" uint32_t a = ..., b = ...; uint32_t res = a + b; uint32_t carry = res < a; """ But note well: this trick doesn't work when there is also a carry-in. That is, uint32_t res = a + b + carry; carry = res < a; won't work when b = 0xFFFFFFFF and carry = 1, so you have to check for carries twice, or set `carry = res < a || (res == a && carry);` or similar. So if you are implementing bignums on an architecture that doesn't have ADC (such as RISC-V), or on a vector architecture where the performance would be bad, then you often need an alternative strategy such as using less than 32/64 bits per limb (whether that's half-size, or 30 bits or 51 bits or whatever).

nj5rq 2y ago

Fascinating article, I have always been wondering how these big number libraries worked. As a side question, does anyone the program that the author used when making that "addition" and "multiplication" performance graph? Thanks.

nerdponx 2y ago
I don't know for sure, but you can do that kind of thing pretty easily with Matplotlib in Python. Or in R base graphics, with more effort to get it looking pretty.
jaeckel 2y ago
If you're interested in even more details you could also have a look at Tom's book about bignum arithmetics. C.f. third paragraph of https://github.com/libtom/libtommath#summary https://github.com/libtom/libtommath#summary The library evolved since then, so the algorithm documentation only exists in code nowadays. FTR I'm involved in the project as maintainer.

kazinator 2y ago

LOL, I've been there: https://www.kylheku.com/cgit/txr/commit/mpi?id=98dedd310b1d5d876b0fbb0ebd6c4df9bd7b2d88 https://www.kylheku.com/cgit/txr/commit/mpi?id=98dedd310b1d5... (Patch was originally from 2011; it was bugfixed once, and then in 2015 converted to git commit: https://www.kylheku.com/cgit/txr/commit/mpi-patches/faster-square-root?id=ef47dfe4fcb7c1be369ae83221386b9da6474a1e https://www.kylheku.com/cgit/txr/commit/mpi-patches/faster-s... ) Another one: faster highest-bit search: https://www.kylheku.com/cgit/txr/tree/mpi-patches/bit-search-optimizations?id=124e7dd6977a0853d7a8399921e31fd1ccde2dcb https://www.kylheku.com/cgit/txr/tree/mpi-patches/bit-search... That does use GCC built-ins today: https://www.kylheku.com/cgit/txr/commit/?id=15b7c542dc44899e8db7addfcc2f1c1c4a188b49 https://www.kylheku.com/cgit/txr/commit/?id=15b7c542dc44899e...

haberman 2y ago

Coincidentally I was just writing a bignum library from scratch two weeks ago. A few interesting things I learned: 1. Knuth volume 2 has an extensive discussion of the problem space. I've only gotten a chance to skim it so far, but it looks interesting and approachable. 2. I need to support bitwise operations, which operate on a two's complement representation, so I figured it would be simpler to use two's complement internally, despite seeing that most (all?) bignum libraries use signed magnitude. I'm starting to regret this: two's complement introduces a lot of complexity. The most fun thing about working on a bignum library is that it makes the algorithms you learned in grade school for add/subtract/multiply/divide relevant again. The basic ("classical") algorithms on bignum are basically exactly the same thing you learned in grade school, except on a much larger base than base-10.

vintermann 2y ago
Knuth has a "hugenum" scheme, TCALC, which is based on trees rather than byte arrays. There are many ways to do it, depending on which type of numbers you need most - for TCALC, it's numbers close to towers of exponents. Of course for every big number encoding that can beat an array of bytes for some numbers, it also has to lose to byte arrays for some numbers... But it's possible to bound this inefficiency, so that it's never more than twice as bad as byte arrays, but potentially exponentially more efficient for the favoured numbers.
- lysp 2y ago
  Twice as bad being a timeout and 2nd call?
  - vintermann 2y ago
    I don't quite get what you mean. Computation time is also bounded with the size of the tree.

tgot 2y ago

I think that your description is almost excellent, but that you're fundamentally misleading in describing what you are doing as a "30-bit" digit. It's a 10^9 digit mathematically, occupying 30 bits of storage. You do briefly mention that it's 10^9, but repeatedly say 30-bits.

samatman 2y ago
It isn't misleading at all. A hexadecimal digit has 4 bits of entropy. You can guess it correctly by chance one in sixteen times. Calling that a four-bit digit is correct. Same with a 30 bit digit, all that changes is the magnitude. The "storage bits" aren't a detail, they're an essential property. Stored in balanced ternary it would still have 30 bits of entropy.
- derdi 2y ago
  I think you misunderstood the parent. If you speak of a number in the range 0-15 stored in 4 bits, we can all agree that "4 bit digit" is the appropriate term for it. But what about a number also stored in 4 bits but restricted to the range 0-13? It's a digit that fits in 4 bits, but calling it a "4 bit digit" without further qualification would omit relevant information.
  - samatman 2y ago
    That hadn't occurred to me, I had thought (having read the Fine Article, but not closely, nor really more than glanced at the code) that it was using exactly 30 bits and reserving two for, IDK, carry detection, and that 10^9 was just approximating 1073741824. We do that every time we call 1024 a KB and not a KiB, so there's precedent. Ok. Fair. It's a ~28.9 bit digit taking up 32 bits of storage. It turns out I was mislead!

parentheses 2y ago

I am so surprised that there's no exploration of Karatsuba's algorithm. That's what makes the Python implementation perform. I actually came here hoping to find discussion on Karatsuba. https://en.m.wikipedia.org/wiki/Karatsuba_algorithm https://en.m.wikipedia.org/wiki/Karatsuba_algorithm

jmts 2y ago
Please see the section titled "Faster multiplication"

globular-toast 2y ago

Funny, I did this myself 10 years ago. Shit, it's been that long..? For my undergrad project I wrote a computer algebra system to do symbolic integration. The supervisor was a hardcore, old school C guy, so naturally I was going to just use C and no libraries. He told me I'd need bignums first, so I got to work (this is because many algorithms like polynomial GCD create massive numbers as they go, even if the inputs and final outputs are generally very small). I just couldn't figure out how to do better than the largest power of 10 per digit at the time. Working with non base 10 arithmetic was a mind fuck for me at the time. So I did it with digits holding 10^9 and the classical algorithms from Knuth. Division is the hardest! At some point I discovered the GNU multiple precision library (GMP) and made my program work with that instead of mine. I was shocked at how much faster GMP was! I finished my project with my own code, but I knew I had to come back to do it better. The breakthrough came when I got a copy of Hacker's Delight. It has stuff like how to detect overflow after it's happened (in C). Something twigged and then I just understood how to fill each word completely rather than use a power of 10. I don't know what confused me before. But, of course, the real way to do it is to use assembly. You can't get close to high performance in C alone. In assembly you get the overflow bit. It's actually easier in a way! So you write tiny platform specific bits for the digits and build on that in C. My add and subtract were then as fast as GMP. I lost interest when it came to implement faster multiplication algorithms. Code in case anyone is interested: https://github.com/georgek/bignums https://github.com/georgek/bignums

chubot 2y ago

Very useful post! Also it's cool to see how many people in this thread have worked on this problem -- lots of new info here I haven't seen I wonder if anyone is interested in implementing a big numbers in Oils? It's a Unix shell with TWO complete implementations - the "executable spec" in Python, and an automatic translation to pure C++ (which is 30x-50x faster) We currently use 64-bit integers in C++, but big nums are a better semantic. Some trivia about bad shell semantics here: Integers - Don't do whatever Python or C++ does - https://www.oilshell.org/blog/2024/03/release-0.21.0.html#integers-dont-do-whatever-python-or-c-does https://www.oilshell.org/blog/2024/03/release-0.21.0.html#in... This is a very self-contained project: the interface is defined by 200 lines of Python: https://github.com/oilshell/oil/blob/master/mycpp/mops.py https://github.com/oilshell/oil/blob/master/mycpp/mops.py and the trivial 64-bit overflowing implementation is also about 200 lines: https://github.com/oilshell/oil/blob/master/mycpp/gc_mops.h https://github.com/oilshell/oil/blob/master/mycpp/gc_mops.h https://github.com/oilshell/oil/blob/master/mycpp/gc_mops.cc https://github.com/oilshell/oil/blob/master/mycpp/gc_mops.cc (We have a fast Ninja-based build system, so you can probably iterate on this in 100 milliseconds or less -- it should be fun for the right person) --- I think the main reason it is specific to Oils is that the bigger number should become GC objects. Details on our GC here: Pictures of a Working Garbage Collector - https://www.oilshell.org/blog/2023/01/garbage-collector.html https://www.oilshell.org/blog/2023/01/garbage-collector.html It's been very solid for the last 18 months, basically because it's well tested by ASAN and #ifdef testing modes. The main thing I'd be concerned with is how to TEST that big number operations are correct. I think there are probably some interesting strategies, which I'd love to discuss. You're of course welcome to adapt existing open source code, including code you've already written -- I probably even prefer that, i.e. something that has had some real world testing. We want all the operations in that file, and it should be integrated with our GC. --- We've had ~6 contributors funded by grants from https://nlnet.nl https://nlnet.nl for the past couple years, so you can even be paid (there's a chance it depends on the country you live in, but generally a wide range of situations is OK). Contact me at andy at oilshell.org or https://oilshell.zulipchat.com/ https://oilshell.zulipchat.com/ if interested!

JonChesterfield 2y ago
It doesn't look especially specific to me - slightly hack the structure in https://github.com/libtom/libtommath https://github.com/libtom/libtommath to make sense to your garbage collector and you're pretty much done. That somewhat punts on the testing question but I doubt add/mul/bitops from that library are buggy.

amelius 2y ago

The title is a little bit dishonest, because they are optimizing their own bignum library.

Kimberrly 2y ago

[dead]

booleandilemma 2y ago

Except he's not doing it for fun at all. He's doing it for clout and publicity on HN.

JonChesterfield 2y ago
Looks like play to me. The earlier post doesn't have a HN reference and there's not much correlation with their publications list. People legitimately do write things like this as entertainment.

guyomes 2y ago

For a small survey of practical efficient methods for bignum arithmetic operations, the Algorithms section of the documentation of GMP [1] is excellent. [1]: https://gmplib.org/manual/Algorithms https://gmplib.org/manual/Algorithms

lifthrasiir 2y ago

Modulo is surprisingly expensive even when you combine it with a quotient. It is almost always better to use binary "limbs", in this case 31 or 32 bits wide, because decimal parsing and printing should be much rarer than individual operations in general.

colonwqbang 2y ago
This is a bit misleading. Quotient and modulo with constant right-hand side is readily optimised into faster operations, like multiplication and bit shifts. "All" optimising compilers do this and don't emit a true divide instruction in such cases. See for instance: https://godbolt.org/z/vxnTPvY6M https://godbolt.org/z/vxnTPvY6M Quotient and modulo with a variable RHS is expensive, but this isn't necessary for the algorithm showed in the article. Although I agree that binary wins over decimal anyway. We are using binary computers after all, for the last 50 years or so.
- Someone 2y ago
  > Quotient and modulo with a variable RHS is expensive For dividing ‘big’ bignums by ‘small’ divisors, I would use/tweak libdivide (https://libdivide.com/ https://libdivide.com/), and spend the time to optimize the division. For dividing ‘big’ bignums by ‘big’ divisors, that may not be worth it. (Determining the correct values of the various ‘big’ and ‘small’ values left as an exercise for the reader. Those will depend on the hardware at hand)
- worstspotgain 2y ago
  The reason a division by 1000000 can be efficiently turned into a multiplication is because you're targeting a system that has a 32x32->64 multiplication instruction, which has plenty of room to spare. But if you have that, you would be better off using e.g. 2^32 as your modulo.
  - colonwqbang 2y ago
    I agree that 2^32 is better, but the reason is not that division is too expensive because OPs algorithm after standard optimisation does not need any division. Clearly choosing a base type that is fast on the machine, is key. But I feel even a 64/32 divide with constant operand will be optimised into something reasonably fast on most systems. If you have a counterexample it would be interesting to hear. A system without fast divide likely doesn't have any floats either. So all arithmetic would be fixed point and 32x32->64 multiplication would be a pretty common operation.