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 1x102 (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 1x102 - 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]x100 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: Xx102 x10-2. Now we get the correct answer, [X x102/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.