Fun with CBM arithmetics

Basic and Machine Language

Moderator: Moderators

User avatar
Mike
Herr VC
Posts: 4816
Joined: Wed Dec 01, 2004 1:57 pm
Location: Munich, Germany
Occupation: electrical engineer

Fun with CBM arithmetics

Post by Mike »

Just thought I'd like to share this one:

Code: Select all

1 FORT=-100TO100
2 PRINT3/4+T/1E9
3 NEXT
Look what happens after '.75' is output. :lol:
wimoos
Vic 20 Afficionado
Posts: 345
Joined: Tue Apr 14, 2009 8:15 am
Website: http://wimbasic.webs.com
Location: Netherlands
Occupation: farmer

Re: Fun with CBM arithmetics

Post by wimoos »

I reckon this behaviour has to do with rounding the accumulator.

The remarkable thing is that when you reverse the calculation, and do a little rounding, the oddity is compensated:

Code: Select all

X=3/4+T/1E9
X=X-3/4
X=X*1E9
X=SGN(X)*INT(ABS(X)+.5)
Returns X equal to T.

Regards,

Wim
VICE; selfwritten 65asmgen; tasm; maintainer of WimBasic
User avatar
Jeff-20
Denial Founder
Posts: 5759
Joined: Wed Dec 31, 1969 6:00 pm

Re: Fun with CBM arithmetics

Post by Jeff-20 »

I guess I shouldn't calculate my taxes on the VIC.
High Scores, Links, and Jeff's Basic Games page.
User avatar
Mike
Herr VC
Posts: 4816
Joined: Wed Dec 01, 2004 1:57 pm
Location: Munich, Germany
Occupation: electrical engineer

Re: Fun with CBM arithmetics

Post by Mike »

Jeff-20 wrote:I guess I shouldn't calculate my taxes on the VIC.
For the moment being, yes.

However, I've found strong hints, that only the conversion of CBM floats to decimal output is affected. C64, C16/C116/+4 and even the V4 machines show the same error, but the C128 gets it right, so Commodore got aware of the issue at some time.

I'm investigating further, and want to provide a solution in the next time.
FD22
Vic 20 Hobbyist
Posts: 148
Joined: Mon Feb 15, 2010 12:31 pm

Re: Fun with CBM arithmetics

Post by FD22 »

I'd be inclined to suspect Monte Davidoff, since he wrote the FP maths package for Microsofts' Altair BASIC, and it's that code that Commodore licensed for the majority of their 8-bit machines - with the exception being the C128 which had BASIC 7 largely re-written from scratch (and that, presumably, included the math logic).
groepaz
Vic 20 Scientist
Posts: 1180
Joined: Wed Aug 25, 2010 5:30 pm

Re: Fun with CBM arithmetics

Post by groepaz »

it's probably just a different way of implementing floats - with different kind of errors (every implementation has those by principle)
I guess I shouldn't calculate my taxes on the VIC.
you should never do that on anything that uses floating point arithmetics. financial packages usually use so called "bignum" libraries for that (which work eg in BCD with a fixed resolution)

you should be able to find demonstrations on how these errors show in excel too, for example :)
I'm just a Software Guy who has no Idea how the Hardware works. Don't listen to me.
User avatar
Mike
Herr VC
Posts: 4816
Joined: Wed Dec 01, 2004 1:57 pm
Location: Munich, Germany
Occupation: electrical engineer

Re: Fun with CBM arithmetics

Post by Mike »

groepaz,

this is not just the occasional rounding error with 1, maybe 2 units wrong in the lowest place we're talking about here. It's a *defect* in the float to ASCII conversion which happens under reproducable circumstances and which, when happens, outputs a grossly wrong result.

If you enter PRINT.750000059 (five zeroes) and get back .75000003 as a result, that is just not acceptable.
groepaz
Vic 20 Scientist
Posts: 1180
Joined: Wed Aug 25, 2010 5:30 pm

Re: Fun with CBM arithmetics

Post by groepaz »

well, i dont know according to what standard (if any) the floats in cbm basic have been implemented, so its hard to tell whether this is acceptable or not.

and dont forget that the typical error is 1 or 2 _binary_ units (powers of two!) - and the decimal/float conversion involves multiplications and divisions - errors in the range you are showing there are well expected. (again: you can find similar "strange errors" in excel. or the windows calc.exe. or pretty much everything =P)
I'm just a Software Guy who has no Idea how the Hardware works. Don't listen to me.
User avatar
Mike
Herr VC
Posts: 4816
Joined: Wed Dec 01, 2004 1:57 pm
Location: Munich, Germany
Occupation: electrical engineer

Re: Fun with CBM arithmetics

Post by Mike »

groepaz wrote:and dont forget that the typical error is 1 or 2 _binary_ units (powers of two!)
Yes, that's what I meant in the first place.
well, i dont know according to what standard (if any) the floats in cbm basic have been implemented, so its hard to tell whether this is acceptable or not.
The implementation surely predates IEEE 754 by a few years, but generally, the routines for basic arithmetic in CBM BASIC are just fine. They round correctly to the last digit for all basic arithmetic functions (i.e. addition, subtraction, multiplication and division).
- and the decimal/float conversion involves multiplications and divisions - errors in the range you are showing there are well expected.
No, *they* *are* *not*.

All conversion routines worth their salt don't use the standard float routines for digit extraction, only for scaling the number with the necessary powers of 10 to get an integer number. The digit extraction then is done with integer arithmetic. And then there's no other reason to introduce an error of 29 decimal units in the last place - as shown in the example -, other than a plain defect.
groepaz
Vic 20 Scientist
Posts: 1180
Joined: Wed Aug 25, 2010 5:30 pm

Re: Fun with CBM arithmetics

Post by groepaz »

but generally, the routines for basic arithmetic in CBM BASIC are just fine.
depends on your expectations :) when i made a wrapper for them and used them in C, they behaved surprisingly different to what you'd expect in some corner cases. really impossible to say if the float routines are to blame here - since we dont know what to expect =P
All conversion routines worth their salt don't use the standard float routines for digit extraction, only for scaling the number with the necessary powers of 10 to get an integer number. The digit extraction then is done with integer arithmetic. And then there's no other reason to introduce an error of 29 decimal units in the last place - as shown in the example -, other than a plain defect.
time for disassembling :)
I'm just a Software Guy who has no Idea how the Hardware works. Don't listen to me.
Commodore Explorer
Vic 20 Newbie
Posts: 11
Joined: Thu Jan 31, 2008 12:37 pm

Re: Fun with CBM arithmetics

Post by Commodore Explorer »

It's not just immediately after ".75" is output. Keep reading:

Code: Select all

 .75000027
 .75000027
 .75000028
 .75000028
 .75000029
 .75000029
 .7500006
 .75000061
 .75000062
 .75000063
At some point, a threshold is crossed, and it starts rounding differently. Very strange.

I must confess, I have never investigated floating-point arithmetic very closely, so this should be interesting. (I did try writing my own printf() function in C and the digit extraction for the %f specifier was a pain to implement. Maybe I was doing it incorrectly?)

This reminds me of a joke about a bug in the "Calculator" program that shipped with Windows 3.1.

Q: What is the difference between Windows 3.11 and Windows 3.1?
A: None! Enter 3.11 - 3.1 on Windows calculator. Press = and instead of showing the correct answer, 0.01, it shows 0.00.

At some point (starting with the "Calculator" that shipped with Windows 95, I think?), Microsoft fixed this bug.
User avatar
Kweepa
Vic 20 Scientist
Posts: 1314
Joined: Fri Jan 04, 2008 5:11 pm
Location: Austin, Texas
Occupation: Game maker

Re: Fun with CBM arithmetics

Post by Kweepa »

I did a little investigation.
It seems to be the function that converts FAC1 to a string.
PRINT .750000059 has FAC1 E = 80 M = C0 00 00 FD on entering $DDDD which is correct.
If 0.5 <= FAC1 < 1.0, FAC1 is multiplied by 1000000000 using the routine at $DA28.
At this point, the value in FAC1 is 750000030 - so the multiply seems to lose a little accuracy.

You can see the same thing by typing PRINT 1000000000*.750000059-750000000, which returns 29.5.

Could be something to do with the FAC rounding byte...
User avatar
Mike
Herr VC
Posts: 4816
Joined: Wed Dec 01, 2004 1:57 pm
Location: Munich, Germany
Occupation: electrical engineer

Re: Fun with CBM arithmetics

Post by Mike »

Seems like we're zeroing in on the bug. I had hopes it would 'only' be the conversion routine which was broken, but it's (sadly) most probably the multiplication at fault: :(

Code: Select all

1 A=.75+59E-9
2 PRINTA*1E9-7.5E8,1E9*A-7.5E8
3 PRINT(A+0)*(1E9+0)-7.5E8,(1E9+0)*(A+0)-7.5E8
4 PRINT(A*1)*(1E9*1)-7.5E8,(1E9*1)*(A*1)-7.5E8
The left column gives the correct result of 59, while the right column errorneously prints 29.5. In lines 3 and 4 the '+0' and '*1' are two possible ways to force an operand out of the FAC into memory and back again, thus the rounding byte is not directly involved. Nonetheless, the multiplication is not commutative (which I already found out about in an earlier thread), I hadn't expected the current issue though.

Other numbers with both the middle bytes 0 in the mantissa seem to be affected the same way.
User avatar
Mike
Herr VC
Posts: 4816
Joined: Wed Dec 01, 2004 1:57 pm
Location: Munich, Germany
Occupation: electrical engineer

Re: Fun with CBM arithmetics

Post by Mike »

I've found the bug.

The circumstances that trigger the bug are a bit complicated, but I'll try to explain.

During multiplication, this routine is called 4 times for each byte of the mantissa of one of the factors. The current mantissa byte is contained in A, and the Z flag has been set accordingly:

Code: Select all

.DA59  D0 03     BNE $DA5E
.DA5B  4C 83 D9  JMP $D983
.DA5E  4A        LSR A
.DA5F  09 80     ORA #$80
.DA61  A8        TAY
.DA62  90 19     BCC $DA7D
.DA64  18        CLC
.DA65  A5 29     LDA $29
.DA67  65 6D     ADC $6D
.DA69  85 29     STA $29
.DA6B  A5 28     LDA $28
.DA6D  65 6C     ADC $6C
.DA6F  85 28     STA $28
.DA71  A5 27     LDA $27
.DA73  65 6B     ADC $6B
.DA75  85 27     STA $27
.DA77  A5 26     LDA $26
.DA79  65 6A     ADC $6A
.DA7B  85 26     STA $26
.DA7D  66 26     ROR $26
.DA7F  66 27     ROR $27
.DA81  66 28     ROR $28
.DA83  66 29     ROR $29
.DA85  66 70     ROR $70
.DA87  98        TYA
.DA88  4A        LSR A
.DA89  D0 D6     BNE $DA61
.DA8B  60        RTS
Now, the first two instructions, BNE $DA5E and JMP $D983 are supposed to execute a shortcut in case the mantissa byte to multiply with is 0. It is sensible to optimize for this, especially because small integer constants do contain zeroes in their lower significant parts of the mantissa.

The routine would still work without that optimization, but the branch at $DA62 would always be executed (for 8 times in total), skipping the additions to $26 .. $29 - and all the routine then does is painstakingly move the bytes $26 .. $29 over to $27 .. $29, $70, one bit at a time. Of course this can be done much faster with just four load and store instrúctions. For this another routine at $D983 is 'reused', which normally 'normalizes' the mantissa:

Code: Select all

.D983  A2 25     LDX #$25
.D985  B4 04     LDY $04,X
.D987  84 70     STY $70
.D989  B4 03     LDY $03,X
.D98B  94 04     STY $04,X
.D98D  B4 02     LDY $02,X
.D98F  94 03     STY $03,X
.D991  B4 01     LDY $01,X
.D993  94 02     STY $02,X
.D995  A4 68     LDY $68
.D997  94 01     STY $01,X
.D999  69 08     ADC #$08
.D99B  30 E8     BMI $D985
.D99D  F0 E6     BEQ $D985
.D99F  E9 08     SBC #$08
.D9A1  A8        TAY
.D9A2  A5 70     LDA $70
.D9A4  B0 14     BCS $D9BA ---+ normally,
.D9A6  16 01     ASL $01,X    | this branch
.D9A8  90 02     BCC $D9AC    | is supposed
.D9AA  F6 01     INC $01,X    | to happen.
.D9AC  76 01     ROR $01,X    |
.D9AE  76 01     ROR $01,X    |
.D9B0  76 02     ROR $02,X    |
.D9B2  76 03     ROR $03,X    |
.D9B4  76 04     ROR $04,X    |
.D9B6  6A        ROR A        |
.D9B7  C8        INY          |
.D9B8  D0 EC     BNE $D9A6    |
.D9BA  18        CLC       <--+
.D9BB  60        RTS
The instructions from $D983 to $D993 do as supposed, however already beginning at $D995, LDY $68 with STY $01,X is a little bit dubious: Without the optimization, a 0 would end up at $26 as the result of being rotated to the right 8 times. Whether there is a 0 at $68 at all times the routine is called with this purpose cannot be guaranteed. But anyway, that is not what triggers the bug in the first place.

The routine at $DA59 is first called with a non-0 byte - for the examples in the earlier posts this is ultimately the result of adding something around 1..59 divided by 10^9 to the constant 0.75. The routine exits with C=1, which will become important now!

The next mantissa byte is zero, so now the shortcut at $D983 is called. With A=0 and C=1 on entry, the two instructions ADC #$08 and SBC #$08 result in A=0 and C=1 again. At $D9A4, the instruction BCS $D9BA skips the second half of the routine, which is a good thing, however it also executes a CLC at $D9BA, and from there everything goes downhill.

The third mantissa byte is *also* zero, so now the shortcut gets called a second time. This time, however C is 0 (with A again being 0), which results in C=0 and A=255 after SBC #$08. The instructions from $D9A6 .. $D9B8 are now executed, which shift the whole mantissa at least one bit to the right!

For the remaining mantissa byte, the check for a shortcut is skipped. Whatever is finally added to the resulting mantissa, the earlier parts have been inadvertently divided by 2 before, which is exactly what can be seen as false result.
User avatar
Kweepa
Vic 20 Scientist
Posts: 1314
Joined: Fri Jan 04, 2008 5:11 pm
Location: Austin, Texas
Occupation: Game maker

Re: Fun with CBM arithmetics

Post by Kweepa »

Nice find!
The same bug is of course in the c64...
Post Reply