10

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?

  • Before going down this road, double check that you are not bitten by operator precedence differences. – Thorbjørn Ravn Andersen Aug 10 at 20:18
  • 2
    @ThorbjørnRavnAndersen definitely not in my case, as I hate remembering that stuff so use parentheses liberally ;-) – Graham Lee Aug 10 at 20:29
  • 2
    Microsoft Extended BASIC, as used by the Dragon uses a 5 byte float format, with a 1 byte exponent and 4 byte mantissa. Let me check if I got some definition somewhere. – Raffzahn Aug 10 at 20:56
  • 2
    Quite interesting that when MS Basic was first introduced too my Math dept, the poo-pooed it as a toy until they heard that it could do 40 bit floating point. That was far more accurate than the 32-bit floats (these were 16-bit minis) that they were used to. The room went silent and then everyone was suddenly interested. – cup Aug 12 at 7:24
11

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.

| improve this answer | |
  • 1
    yes, python uses the machine IEEE representation, double for native python, and can use single if numpy is used. – Jean-François Fabre Aug 10 at 21:13
  • 1
    The rounding mode and whether hidden guard bits are used for rounding tie-breaks might matter for the numerical stability as well. – Leo B. Aug 10 at 21:58
  • 1
    Thanks! That's what I was looking for. I'm seeing errors become meaningful fast in numbers that are ~10^-4, and I'm just using basic maths operations. I'm a little surprised how "bad" my results are, when I've seen pretty good Mandelbrot plots in MS BASIC. It's possible - nay probable - I'm holding it wrong. – Graham Lee Aug 10 at 22:04
  • 2
    @GrahamLee just a suggestion for a different cause of the problem: in Python you have two different local variables with the same name in different functions, but in Basic you have one global variable. – alephzero Aug 11 at 14:50
  • 1
    Unsurprisingly the error was that I was doing the wrong thing :). "My Python implementation" was actually "my using a Python library" that had a much more nuanced implementation of the algorithm than I thought. If I implement my rough BASIC algorithm in Python I get equally bad answers, so now I need to implement the higher-order algorithm in BASIC. Thanks, everyone. – Graham Lee Aug 12 at 11:15
5

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.

| improve this answer | |
  • Around 1983-84, Microsoft produced several interpreters that used BCD floating point. These were noticeably slower than MBF interpreters, but tended to produce fewer spurious results. If a BASIC interpreter has a maximum floating point exponent that ends in 99, it's likely a BCD interpreter. The Broucke test may also identify BCD BASIC. The Tandy 100, MSX computers and one of the Macintosh MS interpreters used BCD instead of MBF. – scruss Aug 15 at 16:28
  • @scruss interesting :) – hobbs Aug 15 at 17:26
4

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.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.