Project: Nixie tube RPN calculator, part 2

If I go to my HP48 calculator (or the Android or iOS versions thereof) and enter `1 3 ÷`, the answer you'll get is `.333333333333`, 12 digits. Now, if I enter `3 x` I get `.999999999999`. And this is correct, since the calculator makes the assumption that all digits past the displayed answer are `0`.

If, however, I do this with 32-bit floating point calculations, we see that the closest floating point representation of `1/3` is `0.3333333432674407958984375` (exactly). How should that be displayed? If I display it as `0.3333333` then I'm lying, because the closest representation of that is actually `0.333333313465118408203125`. But if I display it as `0.33333334` then I'm inaccurate because there is no way that is `1/3`.

Two things come to mind.

So really, we should be using multiple-precision decimal arithmetic, but even that has problems. To solve the `1 3 ÷ 3 x` problem, we could add one guard digit and round, which means `.999999999999(9)` will be displayed as `1`. Even so, there will always be errors, sometimes fatal ones. There is a whole field around this issue.

But this is a calculator. We can provide an answer and assume that the human knows what context their calculation has, so even if they use the HP48 and see `.999999999999`, they'll know when this is `1` and when it is not `1`.

If, however, we want to use scientific notation, then no matter how many digits of precision we use, we run into all the floating point problems again because N digits of precision plus an exponent is the definition of floating point.

Maybe we shouldn't let the point float? What if we want to be able to represent all numbers between `1x10``2` (`100`) and `1x10``-2` (`0.01`) to an accuracy of `1x10``-2`? Then we need 4 digits. If our 4-digit representation is `ABCD`, then the actual number represented is `ABCDx10``-2`, always. If `ABCD = 0001` the number is `1x10``-2` (`0.01`). and if `ABCD = 9999` then the number is `99.99`, or `1x10``2` - `1x10``-2` (close enough, there's always going to be that upper bound missing). The exponent never changes, which means the point does not float. Add a sign bit and you've got arithmetic.

So where do we find a multiple-precision library? The Gnu Multiple Precision library (GMP) holds out some hope. Let's play with it for a bit.

Let's try `1+2`:

```#include <gmp.h>

int main() {
mpz_t x, y, z;

mpz_init(x);
mpz_init(y);
mpz_init(z);

mpz_set_ui(x, 1);
mpz_set_ui(y, 2);

gmp_printf("%Zd\n", z);
}
```

This will simply print `3`. So far, so good. But when we use `mpz_tdiv_q` instead of `mpz_add`, the answer is zero because 1/2 with the result rounded to an integer towards zero is zero.

This is where our assumed exponent comes in. When we divide `Xx10``-2` by `Yx10``-2` we get `[X/Y]x10``0` because mantissas divide while exponents subtract. I use `[X/Y]` to indicate the integer result of `X/Y` rounded towards zero. But we want an answer that is a multiple of `10``-2`. The only way to do this is add more precision to the numerator: `Xx10``2` `x10``-2`. Now we get the correct answer, `[X` `x10``2``/Y]` `x10``-2`.

```#include <gmp.h>

int main() {
mpz_t x, y, z;

mpz_init(x);
mpz_init(y);
mpz_init(z);

mpz_set_ui(x, 10000);
mpz_set_ui(y, 200);
mpz_tdiv_q(z, x, y);

mpz_t intPart, fracPart;
mpz_t hundred;

mpz_init(intPart);
mpz_init(fracPart);
mpz_init_set_ui(hundred, 100);

mpz_tdiv_qr(intPart, fracPart, z, hundred);

gmp_printf("%Zd x 10^-2\n", z);
gmp_printf("%Zd.%02Zd\n", intPart, fracPart);
}
```

The output of this program is:

`50 x 10^-2`
`0.50`

And that is correct. For multiplication, we will end up with twice as many digits as we need, so we just round off.

So much for the four functions. What about exponentiation and roots? Trigonometric functions? For that we'll need the incredible magical CORDIC algorithm, which I'll talk about in part 3.