epicartificials blog

Arbitrary precision numbers

20.04.2020

For this week's homework we got to implement numbers with arbitrary precision. In programming we store numbers in binary format. This means also the decimal numbers are stored in binary. This means we need to store all the digits in a finite amount of disk space (usually 32-bits or 64-bits). The result of that is that we store numbers such as $\frac{1}{3}$ as $0.33333333333$. But when we add the number $\frac{1}{3}$ three times we expect to get $1$ but thanks to the "rounding" we get $0.99999999$. Another problem we need to solve is that programming languages use limited data types (in C++ int is stored in 4 bytes, that means a range from $<-2147483648, 2147483647>$). We have to be able to store much larger numbers.

Our school assignment:

This assignment will be about numbers: your goal will be to
implement a class, ‹number›, which represents a real number with an
arbitrary number of digits (also known as an ‘arbitrary precision’
number). The class should have operators for addition, subtraction,
multiplication, division and comparisons, as well as unary minus.

Additionally, add ‹power( int )› and ‹sqrt( int )› as ‹const›
methods. The argument to ‹power› is the exponent. The argument of
‹sqrt› gives the number of «decimal» digits that should be exact in
the result (see also below).

The constructor of ‹number› takes a single integer argument and
constructs an instance of ‹number› with the argument as its value
(this constructor should allow implicit conversions from ‹int›). A
default-constructed ‹number› should be 0.

As an example, all of the following should be valid code:

    number a( 10 ), b( 25 );
    number c = a + b;
    number d = c / a;
    number e = d * a;
    assert( e == a + b );
    assert( e != a );
    assert( c > a );
    assert( a < b );
    assert( c.power( 2 ) > c );
    c = number( 10 ).power( -5 );
    assert( c > c.power( 2 ) );

The decimal digits supplied to ‹sqrt› should be interpreted as
follows:

    number s = number( 145 ).sqrt( 3 ); /* 3 fractional digits */
    /* the exact result rounded to 3 fractional places is 12.042 */
    number lower = number( 120415 ) * number( 10 ).power( -4 );
    number upper = number( 120425 ) * number( 10 ).power( -4 );
    assert( s > lower );
    assert( s < upper );

Or in other words, if your result is (single) rounded to the given
number of decimal places, it must agree in all those decimal places
with a correctly rounded exact result.

Division by zero and square roots of negative numbers are undefined.

PS: A reference solution is about 250 lines of not very dense code.

At first this looks quite simple. Implement $+, -, *, /, =, \neq$ operations with power and square rooot. Turns out this was more difficult then expected becuase again, we need a good time complexity. $\mathcal{O}(n^2)$ wasn't enough.

So how do we store the numbers?

Since we know that the division is not precise enough we will store our number as fractions $\frac{p}{q}$. This means both $p$ and $q$ will be whole numbers.

The whole number can't be an already existing datatype because any type would be too small. But as we all know we can represent a number as its digits multiplied by the base $1234 = 110^3 + 210^2 + 310^1 + 410^0$. We can emulate this behaviour byt storing the digits in a list.

But now that we store the digits in a list, what about negative numbers? We can store the sign separately in our number and let our whole number always be positive.

Now we need to understand how to store the numbers we can look on how the operations $+, -, *, /, =, \neq$ really work to implement them. For the complexities and algorithms that can be used you can look at this table

($+$) whole plus

For plus we only need to implement the school method which has $\mathcal{O}(n)$ complexity. Allign both numbers under each other justified right. Add the digits up and if the digit is greater than base, then subtract base and carry 1 to the left.

($-$) whole minus

Minus, same as plus, has $\mathcal{O}(n)$ complexity. Allign both numbers under each other justified right. subtract the digits up and if the digit is smaller than 0, then add base and borrow 1 from the left.

($*$) whole times

School method of multiplication is unfortunately slow. Its time complexity is $\mathcal{O}(n^2)$. But that is good enough on small numbers (by small lets say $<$ 20 digits long). Can we do any better?

Yes. One faster algorithm is called Karatsuba algorithm which has $\mathcal{O}(n^{log_2(3)}) \approx \mathcal{O}(n^{1.58})$. Karatsuba uses divide-and-conquer method to split the number into halves and split into three resursive calls and then combines the result together.

There are methods faster than this, but for this assignment $\mathcal{O}(n^{1.58})$ was enough.

($/$ and $\%$) whole division

For this the school method was satisfactory as well. This division is also very small $O(n^2)$, but we can use it less often in our code to slow other operations.

When implementing the division operator we can return the final reminder to obtain also the operator for (%) modulo.

Thats for whole numbers, but what about the rest?

As I said earlier, we will store our number as a $\frac{p}{q}$ where $p$ and $q$ are whole numbers. Now we know how to operate on those whole numbers but what now?

Since we store our number as a fraction we will implement our operations as fraction operations.

($+$) fractional addition

As we all know from school we add fractions as follows: $\frac{p}{q} + \frac{r}{s} = \frac{p * s + q * r}{q * s}$. Although, $q * s$ isn't the smallest common denominator. To calculate the smallest common denominator we use $lcm(q, s)$. Using this we get the result as $\frac{p * \frac{lcm(q, s)}{q} + \frac{lcm(q,s)}{s} * r}{lcm(q, s)}$.

($-$) fractional addition

Same as with addition we get $\frac{p}{q} - \frac{r}{s} = \frac{p * s - q * r}{q * s} \Rightarrow \frac{p * \frac{lcm(q, s)}{q} - \frac{lcm(q,s)}{s} * r}{lcm(q, s)}$.

($*$) fraction subtraction

Fraction multiplication is quite simple $\frac{p}{q} * \frac{r}{s} = \frac{p * r}{q * s}$. This, without simplification might lead to some big numbers $\frac{1}{12345} * \frac{12345}{1} = \frac{12345}{12345}$ which can be simplified to just $\frac{1}{1}$. But how do we acieve this?

To simplify a fraction we calculate the $gcd(p, q)$ and divide both numerator and denumerator by the $gcd$ value.

($/$) fraction division

Since division is not precise, we will multiply instead. To multiply fractions we multiply by the inverse of the fraction so $\frac{p}{q} / \frac{r}{s} = \frac{p}{q} * \frac{s}{r} = \frac{p * s}{q * r}$.

But wait, can I add a negative number?

Since we store the sign separately in the number we can get operations such as $a + (-b)$. We can solve these situations by flipping the sign and calling $a - b$ instead.

We also need to be careful about subtraction. Rememner, that we use whole numbers. To account for this we swap the operands and flip the sign $a - b, a<b$ becomes $-(b - a)$.

Now that we got the mathematical operations figured out, its time to look at the power and square root.

$n^{exp}$

One prerequisite is that it needs to be faster than $\mathcal{O}(exp)$.

To achieve faster results we istead of multiplying $n$ $exp$ times we square $n$ until the exponents match. You can find the code on gekksforgeeks.org. This method has $\mathcal{O}(\log(exp))$ complexity which is much better.

$sqrt(n, prec)$

Calculating precise square root is probably impossible, so we need to calculate square root to $prec$ decimal places. I found two methods of calculating the square root.

First one was the Babylonian method. Someone calculated here that this method has $\mathcal{O}(\log(\log(n/m))$ where $n$ is input and $m$ is the starting seed. If you read more about this method you will find that depending on the starting value you are able to find the square root faster or slower. Since I didn't choose the right starting value the calculation got stuck on computing numbers with precision $> 11$.

The second option was a binary search algorithm. This has a time complexity of $\mathcal{O}(\sqrt{n})$ and has worked much better and reliable for me.

Conclusion

We have learned how we can implement an arbitrary large number, how addition, subtraction, multiplication and division work on both whole and fracional numbers, how to implement square root and nth-power.

I hope that you have learned something new and I wish you a nice day.

#math #C++ #algorithm