Commodore BASIC and binary floating point precision

Question

I am mildly curious that though the 6502 provides BCD arithmetic which would be useful for implementing decimal floating point, Commodore BASIC uses, like all (?) Micro-Soft BASIC, binary floating point instead.

Are there any easy example that show precision errors in Commodore BASIC, that would not be present if it would be based on decimal FP?

A classic test of the difference is 0.1 + 0.2 = 0.3; this evaluates to false in pretty much every modern language (since almost all of them use the IEEE floating point that is built into modern hardware). There is even a website devoted to this oddity: 0.30000000000000004.com

I tried this on a C64 emulator (which uses the same BASIC as the PET) and to my astonishment, it correctly evaluated to true. So did some other obvious tests like 0.1 * 10 = 1 and 0.1 + 0.9 = 1 but they worked as well.

What test would give a wrong answer on Commodore BASIC? That is, I'm not asking for a way to get it to demonstrate rounding errors per se; that much is trivial. I'm asking for a way to get it to give a wrong answer, not because it lacks infinite precision, but specifically for(simple) cases where decimal arithmetic would give the right answer. Some Commodore BASIC (MS-BASIC) equivalent to the 0.1 + 0.2 = 0.3 test on IEEE 754.

Answer 1

This example reveals a rounding error under Commodore BASIC V2.0:

  A=0.3:B=0.6:IF A+B<>0.9 THEN PRINT A+B-0.9

Running this on a C64 yields a difference of 2.32830644e-10. Other pairs that fail are 0.4+0.5, 0.6+0.1 and 0.8+0.1. Please note that also the order in which the numbers are summed up affects the result. 0.6+0.1-0.7 yields a difference, while 0.1+0.6-0.7 results to 0.

Answer 2

Here is my favourite example for this problem. I often use it to show Excel's mathematical shortcomings, but not surprisingly it works the same in the C64:

10 A = 0.1
20 B = 0.1
30 FOR I = 1 TO 10
40 D = B
50 B = 20 * A - 19 * B
60 PRINT B
70 A = D
80 NEXT I

In every iteration, the algorithm should be doing 20 * 0.1 - 19 * 0.1 = 0.1, but the output on this simulator is

 .0999999999
 .100000002
 .0999999578
 .100000845
 .0999831052
 .100337895
 .0932421037
 .235157925
-2.60315849
 54.1631699

Answer 3

Might I suggest you try 0.11+0.12?

I believe IEEE754 will in fact give the right answer on 0.1+0.2=0.3, using standard single precision. It is, however, not difficult to provoke IEEE754 failures, for instance on 0.11+0.12. The C program below show the raw bin32 representations of the relevant IEEE754 numbers, the program output is:

a  :3dcccccd
b  :3e4ccccd
a+b:3e99999a
c  :3e99999a
IEEE754 copes
a  :3de147ae
b  :3df5c28f
a+b:3e6b851e
c  :3e6b851f
IEEE754 fails

Program:

#include <stdio.h>
#include <stdint.h>
int main( void ) {
    float a = 0.1;
    float b = 0.2;
    float c = 0.3;
    float apb = a+b;
    printf( "a  :%x\n", *(uint32_t *)&a);
    printf( "b  :%x\n", *(uint32_t *)&b);
    printf( "a+b:%x\n", *(uint32_t *)&apb);
    printf( "c  :%x\n", *(uint32_t *)&c);
    if ( a + b == c ) {
            printf( "IEEE754 copes\n" );
    } else {
            printf( "IEEE754 fails\n" );
    }

    a = 0.11;
    b = 0.12;
    c = 0.23;
    apb = a+b;
    printf( "a  :%x\n", *(uint32_t *)&a);
    printf( "b  :%x\n", *(uint32_t *)&b);
    printf( "a+b:%x\n", *(uint32_t *)&apb);
    printf( "c  :%x\n", *(uint32_t *)&c);
    if ( a + b == c ) {
            printf( "IEEE754 copes\n" );
    } else {
            printf( "IEEE754 fails\n" );
    }

    return 0;
}

Answer 4

Although not a direct answer to the question, I believe this is still worth pointing out:

Atari BASIC used BCD for its FP implementation. It was free from the sorts of errors being outlined above. I have found several other BCD implementations, but nothing "mainstream" to this degree.

The flip side is that BCD is theoretically slightly slower, you have to be careful turning the BCD mode on and off, and you get slightly fewer digits per 40-bit store.

On top of that, the Atari implementation was notoriously slow, but that was not specifically because it was BCD, it was just bad.

Answer 5

Probably the simplest pattern to look for is 0.1+2-2-0.1. Depending upon precision, 0.1 is going to be (the .... is some number of repetitions of 0011)

0.00011....00110011
 0.00011....0011010
  0.00011....001101
   0.00011....00110

Adding 2 is going to reduce the available precision at the low end by five bits (more than a whole repetition of the 0011 pattern), forcing some kind of rounding there. After 2 is subtracted, zeroes would be shifted into the low-order bits, and subtracting 0.1 would leave values in those bits.

When using decimal floating point, all calculations would be exact, and thus the final residue would be zero.

Incidentally, with regard to speed, the 6502 has patented circuitry to perform BCD addition and subtraction as fast as binary (interestingly, the CMOS versions of the chip do not have such circuitry, causing BCD addition and subtraction to be slower), but performing most other kinds of BCD math efficiently requires some substantial lookup tables.

Answer 6

Another precision test, previously mentioned in a comment by Tim Locke on this site, was published in Antic magazine Vol 1 No.4. It was submitted by one "R. Broucke" (the late Roger A. Broucke, then at UT Austin):

10 S=0
20 X=0
30 FOR N=1 TO 1000
40 S=S+X*X
50 X=X+0.00123
60 NEXT N
70 PRINT S,X
80 PRINT "CORRECT RESULT: 503.54380215, 1.23"

When run on a binary floating point interpreter (cbmbasic) it produces:

 503.543832            1.23000004

On a decimal floating point interpreter (tibasic):

 503.5438022   1.23 

Typically (but not always) binary floating point interpreters add spurious digits after the 1.23 result. Microsoft produced a few interpreters that used decimal floating point around 1983-4. These include BASIC for the Tandy 100 portable and MSX BASIC. Microsoft BASIC for Macintosh included two interpreters, one binary, one decimal.

Answer 7

A shorter version of this answer is simply:

PRINT (0.3+0.6)-0.9

...which will result in 2.32830644E-10