arbitrary-precision | :100 : Arbitrary precision | Apps library

CPython's int, in Python 3, is represented as a sign-magnitude array of values, where each element in the array represents 15 or 30 bits of the magnitude, for 32 and 64 bit Python builds respectively. This is an implementation detail, but a longstanding one (it was originally 15 all the time, but it was found to be an easy win to double the size and number of used bits per "digit" in the array when working on 64 bit systems). It has optimizations for ints that fit into a single (or sometimes two) such array values (it extracts the raw value from the array and performs a single CPU operation, skipping loops and algorithms that apply to the arbitrary length case), and on 64 bit builds of CPython, that currently means values with magnitude of 30 bits or less are generally optimized specially (with 60 bit magnitudes occasionally having fast paths).

That said, there's rarely a reason to take this into account; CPython interpreter overhead is pretty high, and it's pretty hard to imagine a scenario where manually breaking down a larger operation into smaller ones (incurring more interpreter overhead to do many small operations at the Python layer) would outweigh the much smaller cost of Python doing the array-based operations (even without special fast paths) at the C layer. The exceptions to this rule would all rely on non-Python-int-solutions, using fixed size numpy arrays to vectorize the work, and largely following C rules at that point (since numpy arrays are wrappers around raw C arrays most of the time).

In answer to your question why?, it is by design, and the rationale for this is in PEP 238:

Floor division will be implemented in all the Python numeric types, and will have the semantics of:

a // b == floor(a/b)

except that the result type will be the common type into which a and b are coerced before the operation.

Specifically, if a and b are of the same type, a//b will be of that type too. If the inputs are of different types, they are first coerced to a common type using the same rules used for all other arithmetic operators.

...

For floating point inputs, the result is a float. For example:

3.5//2.0 == 1.0

For complex numbers, // raises an exception, since floor() of a complex number is not allowed.

This PEP dates back to Python 2.2. I've suppressed a paragraph that discusses the now obsolete distinction between int and long.

Dyalog APL includes the operator big in the dfns workspace which allow arbitrary-precision arithmetic:

Just like with Python's decimal.Decimal, you can configure the precision of Julia's BigFloat:

That algorithm is horrible.

The first step should be to determine the sign of the result (from the sign of the numerator and divisor); then find the magnitude of the numerator and divisor (and have short-cuts for the "numerator is 0" and "abs(divisor) is 0 or 1" cases where actually doing a division is avoidable or impossible) so that the code that does the actual division only ever deals with positive numbers.

The second step should be to determine if the divisor is small enough to fit in a single digit (with digits that are in whatever base is the largest that your language/environment supports - e.g. for C with 32-bit integers it might be "base 65536", and for 64-bit 80x86 assembly language (where you can use 128-bit numerators) it might be "base 18446744073709551616"). At this point you branch to one of 2 completely different algorithms.

Small Divisor

This is a relatively trivial "for each digit in numerator { divide digit by divisor to find the digit in the result}" loop (followed by fixing up the sign of the result that you determined at the start).

Large Divisor

For this I'd use binary division. The idea is to shift the numerator and divisor left/right so that divisor becomes as large as possible without being larger than the numerator. Then you subtract divisor from numerator and set the bit in the result that corresponds to how much you shifted left/right; and repeat this, so that (after the initial shifting) it ends up being "shift whatever remains of the numerator left; then compare to divisor, and if numerator is larger than divisor subtract divisor from numerator and set bit in result" in a loop that terminates when there's nothing remaining of the numerator).

Off-Topic Alternative

For most cases where you need arbitrary-precision division; it's better to use rational numbers where each number is stored as three (arbitrary sized) integers - numerator, divisor and exponent (like number = numerator/divisor * (1 << exponent)). In this case you never need to divide - you only multiply by the reciprocal. This makes it better for performance, but also significantly better for precision (e.g. you can calculate (1/3) * 6 and guarantee that there will be no precision loss).

Given you need precise decimal places, all of the libraries you mentioned are either out of bounds or will require setting absurd levels of precision to guarantee your results are rounded accurately.

The problem is they're all based on binary floating point, so when you specify a precision, e.g. with mpf_set_default_prec in GMP/MPIR (or the equivalent call in MPFR, a fork of GMP that provides more guarantees on exact, portable precision), you're saying how many bits will be used for the significand (roughly, how many bits of integer precision exist before it is scaled by the exponent). While there is a rough correspondence to decimal precision (so using four bits of significand for every one decimal digit of precision needed is usually good enough for most purposes), the lack of decimal precision means explicit output rounding is always needed to avoid false precision, and there is always the risk of inconsistent rounding: a number that logically ends in 5 and should be subjected to rounding rules like round half-even, round half-up or round half-down might end up with false precision that makes it appear to not end exactly with 5, so the rounding rules aren't applied.

What you want is a decimal math library, like Java's BigDecimal, Python's decimal, etc. The C/C++ standard library for this is libmpdec, the core of the mpdecimal package (on Python 3.3+, this is the backend for the decimal module, so you can experiment with its behavior on Python before committing to using it in your C/C++ project). Since it's decimal based, you have guaranteed precision levels; if you set the precision to 4 and mode to ROUND_HALF_EVEN, then parsing "XX.X45" will consistently produce "XX.X4", while "XX.X75" will consistently produce "XX.X8" (the final 5 rounds the next digit to even numbers). You can also "quantize" to round to a specific precision after the decimal point (so your precision might be 100 to allow for 100 digits combined on left and right hand sides of the decimal point, but you can quantize relative to "0.0000" to force rounding of the right hand of the decimal to exactly four digits).

Point is, if you want to perform exact decimal rounding of strings representing numbers in the reals in C/C++, take a look at libmpdec.

For your particular use case, a Python version of the code (using the moneyfmt recipe, because formatting numbers like "1.1E-26" to avoid scientific notation is a pain) would look roughly like:

You cannot detect signed int overflow. You have to write your code to avoid it.

Signed int overflow is Undefined Behaviour and if it is present in your program, the program is invalid and the compiler is not required to generate any specific behaviour.

As we can see in the reference source, BigInteger in .NET uses a fairly slow multiplication algorithm, the usual quadratic time algorithm using 32x32->64 multiplies. But it is written with low overhead: iterative, few allocations, and no calls to non-inlinable ASM procedures. Partial products are added into the result immediately rather than materialized separately.

The non-inlinable ASM procedure can be replaced with the _umul128 intrinsic. The manual carry calculations (both the conditional +1 and determining the output carry) can be replaced by the _addcarry_u64 intrinsic.

Fancier algorithms such as Karatsuba multiplication and Toom-Cook multiplication can be effective, but not when the recursion is done all the way down to the single limb level - that is far past the point where the overhead outweighs the saved elementary multiplications. As a concrete example, this implementation of Java's BigInteger switches to Karatsuba for 80 limbs (2560 bits because they use 32 bit limbs), and to 3-way Toom-Cook for 240 limbs. Given that threshold of 80, with only 64 limbs I would not expect too much gain there anyway, if any.

First, remember the sign (+ or −) and change the number to its absolute value.

Normalize the significand to be in [1, 2) and calculate the exponent:

• If the number is greater than two, divide it by two until it is in [1, 2). (That is the interval from 1 to 2 that includes 1 but not 2.) The number of times you divided will be the exponent.

• If the number is less than one, multiply it by two until it is in [1, 2). The exponent will be the negative of the number of times you multiplied.

• If the number is already in [1, 2), the exponent is zero. Convert the number in [1, 2) that you ended up with to binary.

Then round the number to fit in the significand of the floating-point format:

• Round the number to as many bits as the significand allows. If this pushes the number up to two, divide it by two and add one to the exponent.

Now you have the sign, the exponent, and the significand.