How can I understand numerical precision of values in Microsoft BASIC (on the Dragon 32)?

Question

I implemented a differential equation solver on my Dragon 32 (which uses Microsoft Extended Color BASIC) but the results I get quickly diverge from those I expect, which are the same as those I get when I write the solver in Python. I suspect that the divergence is due to numerical precision, but I don't know how to model the precision of numbers stored in this BASIC (meanwhile, modelling imprecision in IEEE854 floating point numbers is well-documented).

I don't even know whether non-integer numbers are stored using a floating point representation: all I know (which I worked out by PRINT MEM before and after DIMing some of my data, to ensure the problem would fit in memory) is that five bytes are used per number. None of the books I own on the Dragon 32 talk about numerical precision at all, most of them deal only with integer maths (including the Dragon manual, which only has a little diversion into the trig functions at the end).

Is the memory layout of numbers in Microsoft Extended BASIC documented anywhere? What limits on precision are there in use?

Answer 1

Microsoft Extended BASIC, as used by the Dragon, uses a 40 bit (5 byte) float format(*1):

 Field          Size (Bits)  
 Exponent Sign      1 
 Exponent           7
 Mantissa Sign      1
 Mantissa          31
                   --
                   40

Numbers are always normalized with the top mantissa bit removed.

As an oddity, they are stored in a 6 byte structure with the 6th byte unused.

I suspect that the divergence is due to numerical precision,

Most likely, as Microsofts 40 bit format differs from IEEE-754 single (32 Bit) as well as from double (64 bit). Since Python uses double IEEE format, its results will differ in precision and these differences may add up.

In addition, if you use any trigonometric function or square-root, results may wander off even more due to a different approximation.

---

Some unordered thoughts:

One solution would be to use single precision (with numpy) with your Python solver, as Jean-François Fabre sugsests. But it's still a different precision, so depending on the task the results may still vary a lot.

This could be levelled by 'amputating' the BASIC numbers down to single precision as well by clearing the lowest mantissa byte after each computational assignment. Doing this in BASIC alone will be rather cumbersome (*2), but a quite short machine routine called with USR will do the trick.

100 A=1+1 : A=USR0(A)

USR will lookup the variable and move it to the Floating Point Accumulator at $004F and the assignment stores the result in A again. Of course any other variable can be uses as well (*3). Since BASIC is doing most of the Stuff, the machine routine is essentially just two instructions:

    B7 53 00   STAA   $4F    ; Clear the lowest mantissa byte (*4)
    39         RTS           ; Return to BASIC

Of course this has to be stored somewhere. For routines that short, the 6 unused (*5) bytes at $011A work fine, so the add

10 POKE &H011A,&HB7: POKE  &H011B,....
20 DEFUSR0=&H011A

to initialise the USR function.

(CAVEAT: this is a quick hack from my faint memory while peeking at the memory map).

Of course, this still doesn't solve possible different approximations of build in function. Only precision should now be aligned.

(And sorry for being carried away, it was too tempting :))

---

*1 - It's the so called Extended Format (sic) added with the 8 KiB extended Basic.

*2 - Or not, as the Dragon knows VATPTR ... so forget all that follows and do simply A=1+1 : POKE VARPTR(A)+4,0. With preparing the addresses ahead of time (AA=VARPTR(A)+4 ... POKE AA,0) this might be as fast as the machine language routine :(

*3 - BTW, to speed up calculations it's helpful to define the most used variables first.

*4 - Confused ... well, this is using a neat side effect. BASIC passes the type of the parameter in A - 1 for a string and 0 for a float - so we receive a cleared A, ready to be used :))

*5 - Only with a Dragon 32.

Answer 2

As far as I can see, any Microsoft BASIC at that time should have used Microsoft Binary Format for floating-point numbers. Since you found out that numbers are occupying five bytes, that should mean that your implementation is using the 40-bit format.

This format is similar enough to IEEE-754 binary formats in that it has a sign bit, a base-2, biased exponent, and a mantissa with an implicit 1-bit to the left (but there is no support for denormals). The 40-bit Microsoft format has an 8-bit exponent and a 31-bit mantissa, as compared to single-precision IEEE-754's 8 exponent + 23 mantissa, or double-precision IEEE-754's 11 exponent + 52 mantissa, so its precision (as long as you don't touch on the range of denormals) should lie between the two.

Answer 3

You may want to look at the PARANOIA floating point test suite, which tests for quite a few characteristics of the floating point format and implementation (range, bits of precision, rounding, guard bits, etc.).

It is said that it was originally written in BASIC, but I could not find the BASIC version (EDIT: hat tip to @scruss, here it is). The C variant, as follows from the comments, is a rewrite of the Pascal version, which is a rewrite of the original BASIC version.

The coding style, though, is very easy to follow to re-convert back to BASIC. For example, the numerical radix of the f. p. representation is calculated as follows (note that literal constants are replaced with variables to inhibit optimization, but it should not be a problem in BASIC, you could write 1.0 instead of One, etc.):

printf ( "Searching for Radix and Precision.\n" );
W = One;

do {
    W = W + W;
    Y = W + One;
    Z = Y - W;
    Y = Z - One;
} while (MinusOne + FABS(Y) < Zero);

/*
  Now W is just big enough that |((W+1)-W)-1| >= 1.
*/
Precision = Zero;
Y = One;
do {
    Radix = W + Y;
    Y = Y + Y;
    Radix = Radix - W;
} while ( Radix == Zero);

if ( Radix < Two ) {
    Radix = One;
}

printf ( "Radix = %f\n", Radix );

In the pre-IEEE854 formats, radix was not necessarily 2. For example, it was 16 in the IBM floating point.

Other specific tests are equally short and lucid.