6502.org Forum  Projects  Code  Documents  Tools  Forum
It is currently Thu Mar 28, 2024 3:50 pm

All times are UTC




Post new topic Reply to topic  [ 12 posts ] 
Author Message
PostPosted: Mon Apr 08, 2019 2:24 am 
Offline
User avatar

Joined: Thu Mar 11, 2004 7:42 am
Posts: 362
EhBASIC calculates X^Y using the formula EXP(Y*LOG(X)). While this is mathematically equivalent for nonzero X, in a practical calculation there are some issues.
Code:
PRINT 41*41*41
PRINT 41^3

The first statement outputs the correct value 68921, but the second statement outputs 68921.1. The reason is simple. In a floating point point calculation there are a finite number of values that can be represented exactly. e (2.718...) is an irrational number, so LOG(X), when calculated as a floating point value, is an approximation, Y*LOG(X) is an approximation, and EXP(Y*LOG(X)) is an approximation. Furthermore, most of the terms of the power series used to calculate LOG and EXP are approximations.

Because EhBASIC has a rounding (extra precision) byte for FAC1, the fact that a calculation is really an approximation and not an exact value is usually not an issue, but there are cases, such as the example above, where there is an issue.

One way to address this issue is to handle the special case where Y is a (small enough) integer. The routine below specially handles the case where ABS(Y) < 256 (i.e. the exponent will fit in a single byte). You can extend it to handle larger integers if you wish, but small integers (say Y < 32) are one of the more common cases, and are likely to be where the issue is more visible. So it makes sense to cater to common case.

Another reason to handle integer Y specially, is that particularly for small integers, it will be faster. X^3 is more quickly calculated as X*X*X and Z=X^9 can be calculated as Z=X*X:Z=Z*Z:Z=Z*Z*X which is only 4 multiplications.

In fact, that is basically the algorithm: square the result, and multiply by X where appropriate. Note that the case where Y=0 gets handled before we get to the integer Y special case.

The first thing to do is add the label MUL_Adatal after LAB_2B6E, as follows (note that the only change is to add a new label):
Code:
LAB_2B6E
        STA     Cptrl           ; save count pointer low byte
        STY     Cptrh           ; save count pointer high byte
        JSR     LAB_276E        ; pack FAC1 into Adatal
        LDA     #<Adatal        ; set pointer low byte (Y already $00)
        JSR     LAB_25FB        ; do convert AY, FCA1*(AY)
        JSR     LAB_2B88        ; go do series evaluation
MUL_Adatal
        LDA     #<Adatal        ; pointer to original # low byte
        LDY     #>Adatal        ; pointer to original # high byte
        JMP     LAB_25FB        ; do convert AY, FCA1*(AY) and return

Then update LAB_POWER as follows:
Code:
;
; this part is new
;
POWER_1
        .byte   $FF,$7F,$3F,$1F,$0F,$07,$03,$01,$00
POWER_2
        JSR     LAB_279B        ; copy FAC2 to FAC1
        JSR     LAB_276E        ; pack FAC1 into Adatal
        ASL     func_l+1        ; shift sign into carry, shift exponent of ^ into upper bits
        PHP                     ; save sign
        DEC     func_l          ; decrement number of bits to shift
        BNE     POWER_4         ; always branch
POWER_3
        JSR     LAB_27AB        ; returns with X = $00
        STX     FAC_sc          ; clear sign (sqaures are positive)
        LDA     FAC1_e          ; load FAC1 exponent for LAB_MULTIPLY call
        JSR     LAB_MULTIPLY
        ASL     func_l+1        ; shift exponent of ^
        BCC     POWER_4
        JSR     MUL_Adatal
POWER_4
        DEC     func_l          ; decrement number of bits to shift
        BMI     POWER_3         ; loop until done
        PLP                     ; get sign of exponent
        BCS     POWER_5         ; branch if negative
        RTS
POWER_5
        LDA     #<LAB_259C      ; address of floating point constant 1.0
        LDY     #>LAB_259C
        JMP     LAB_26CA        ; FAC = 1 / FAC

;
; this part is the same as before
;

; perform power function

LAB_POWER
        BEQ     LAB_EXP         ; go do  EXP()
        LDA     FAC2_e          ; get FAC2 exponent
        BNE     LAB_2ABF        ; branch if FAC2<>0
        JMP     LAB_24F3        ; clear FAC1 exponent and sign and return
LAB_2ABF
        LDX     #<func_l        ; set destination pointer low byte
        LDY     #>func_l        ; set destination pointer high byte
        JSR     LAB_2778        ; pack FAC1 into (XY)

;
; this part is new
;

; $81 %10000000 $00 $00 = 1
; $82 %1x000000 $00 $00 = 2 to 3
; $83 %1xx00000 $00 $00 = 4 to 7
; $84 %1xxx0000 $00 $00 = 8 to 15
; $85 %1xxxx000 $00 $00 = 16 to 31
; $86 %1xxxxx00 $00 $00 = 32 to 63
; $87 %1xxxxxx0 $00 $00 = 64 to 127
; $88 %1xxxxxxx $00 $00 = 128 to 255
;
; Note that the high bit of FAC1_1 is 1 for a normalized number, so FAC1_e =
; $80 is rejected by masking off no bits (i.e. ANDing with $FF) of FAC1_1
;
        LDX     FAC1_e          ; check if $80 <= FAC1 exponent < $89
        BPL     POWER_6
        CPX     #$89
        BCS     POWER_6
        LDA     POWER_1-$80,X   ; check if lower bits of FAC1_1 are 0
        AND     FAC1_1
        ORA     FAC1_2          ; check if lower 2 bytes of FAC1 are 0
        ORA     FAC1_3
        BEQ     POWER_2         ; branch if ABS(FAC) < 256
POWER_6

;
; this part is the same as before
;
        LDA     FAC2_s          ; get FAC2 sign (b7)
        BPL     LAB_2AD9        ; branch if FAC2>0
                                ; else FAC2 is -ve and can only be raised to an
                                ; integer power which gives an x +j0 result
        JSR     LAB_INT         ; perform INT
        LDA     #<func_l        ; set source pointer low byte
        LDY     #>func_l        ; set source pointer high byte
        JSR     LAB_27F8        ; compare FAC1 with (AY)
        BNE     LAB_2AD9        ; branch if FAC1 <> (AY) to allow Function Call error
                                ; this will leave FAC1 -ve and cause a Function Call
                                ; error when LOG() is called
        TYA                     ; clear sign b7
        LDY     Temp3           ; save mantissa 3 from INT() function as sign in Y
                                ; for possible later negation, b0
LAB_2AD9
        JSR     LAB_279D        ; save FAC1 sign and copy ABS(FAC2) to FAC1
        TYA                     ; copy sign back ..
        PHA                     ; .. and save it
        JSR     LAB_LOG         ; do LOG(n)
        LDA     #<garb_l        ; set pointer low byte
        LDY     #>garb_l        ; set pointer high byte
        JSR     LAB_25FB        ; do convert AY, FCA1*(AY) (square the value)
        JSR     LAB_EXP         ; go do EXP(n)
        PLA                     ; pull sign from stack
        LSR                     ; b0 is to be tested, shift to Cb
        BCC     LAB_2AF9        ; if no bit then exit

; fall thru to LAB_GTHAN (same as before)


Top
 Profile  
Reply with quote  
PostPosted: Mon Apr 08, 2019 7:10 am 
Offline
User avatar

Joined: Thu Dec 11, 2008 1:28 pm
Posts: 10760
Location: England
Nicely done! Thanks for sharing.


Top
 Profile  
Reply with quote  
PostPosted: Mon Apr 08, 2019 8:43 am 
Offline
User avatar

Joined: Wed Feb 14, 2018 2:33 pm
Posts: 1392
Location: Scotland
Interesting stuff - It's easy to get complacent and "trust the computer" - part of my course a uni a long time back was dealing with floating point, rounding and accuracy too...

I also remember schools having a bit of a calculator accuracy test back then too - probably not really relevant, but it as along the lines of: Start with 49 (degrees), then sin, cos, tan, atan, acos, asin then subtract 49...


Relevant to me as I recently started to look into some FP routines for my 6502 project, (I need some 32-bit ones, ideally IEEE 754)

But decided to try this on what I'm currently using: BBC Basic 4:

Code:
>P. 41*41*41
     68921
>P. 41^3
     68921
>P. EXP (3 * LN (41))
     68921


It uses a 5-byte format and may well have better rounding, precision, etc. and it looks OK.

What I was (still am) after was some IEEE 754 32-bit routines and so-far have not found any other than the 128-bit ones here by Marco Granati which I'm sure can be adapted however my plan B was to use the ATmega host processor to my system as a math co-processor which might be a bit cheaty, but seems to work (and is obviously somewhat faster). I just checked it's accuracy with:

Code:
double x = 41.0 ;
printf ("%g\n", x * x * x) ;
printf ("%g\n", pow (41.0, 3.0)) ;


and it printed the expected 68921 twice. (doubles are really 32-bit floats in the ATmega)

Back to the MS BASICs though - it does make you wonder just how many programs were written (and might still be in-use!) to deal with math issues that may not be quite as accurate as it ought to be. (not counting the Pentium FDIV bug, of-course!)

Doing the old school calculator test on BBC Basic yields some rounding though:

Code:
>X=RAD(49) : REM No DEG or RAD mode in BBC Basic - all in RADs, so built-in functions
>P.X
0.855211333
>Y=TAN(COS(SIN(X)))
>Z=ASN(ACS(ATN(Y)))
>Q=DEG(Z)
>P.Q
        49
>P.Q-49
1.49011612E-8



Cheers,

-Gordon

_________________
--
Gordon Henderson.
See my Ruby 6502 and 65816 SBC projects here: https://projects.drogon.net/ruby/


Top
 Profile  
Reply with quote  
PostPosted: Mon Apr 08, 2019 9:05 am 
Offline
User avatar

Joined: Fri Aug 30, 2002 1:09 am
Posts: 8412
Location: Southern California
Quote:
I also remember schools having a bit of a calculator accuracy test back then too - probably not really relevant, but it as along the lines of: Start with 49 (degrees), then sin, cos, tan, atan, acos, asin then subtract 49...

Here's a repeat of a post of mine on a calculator forum a couple of years ago in response to a complaint that the HP-41 calculator did not have good accuracy because alog(log(60)) showed 59.999999950:

    This is something slide-rule users will have a better understanding of. It's not really any bug in the machine, just a limitation of the number of digits it uses.

    Consider: ln(60) is 4.09434456222, according to my HP-71B which uses more digits than the 41 does. Using only 10 digits, e^4.094344562 is 59.9999999867, and e^4.94344563 (ie, incrementing the least-significant digit by 1) is 60.0000000467, both quite a ways off from 60 in the 10th digit. IOW, there is no 10-digit number that will give you 60.00000000 when you take the natural antilog.

    You find this more in for example cos(1°), only one digit, where you get .99984769516; but using only four digits of it to convert back to the angle, .9998 gets you over 14% high, at 1.1459..., and .9999 gets you .81029..., almost 19% low. There is no 4-digit number which, when taking its ACOS, will result in 1.0° (two digits), let alone 1.000° (four digits).

    In my electronics engineering work, I've never needed even six digits. It's easy for math enthusiasts to get carried away with precision if they don't have a feel for what it means in the physical world. Although a few functions will make unacceptable error levels accumulate (for example, for a bank making amortization calculations, ten digits is nowhere near enough, even if the interest rate is expressed with only four digits), most calculators have far, far more precision than most real-life situations have any use for. Audio quality better than the finest cassettes can be had with only three digits to the samples. (999 is close to 1023 which is the maximum that 10 bits can represent; and the digital will have less distortion, and with an adequate sampling rate, better frequency response.) A lot of digital signal processing is even done in 16-bit scaled-integer arithmetic.

    There is no limit to what can be done with numbers, but their relevance to what they measure, control, etc. is another matter. It is difficult or impossible to measure or control pressure, weight, fluid flow, temperature, light intensity, position/displacement, field strength, etc. to precisions anywhere near what common calculators handle. When I need A/D or D/A converters for experiments or development on the workbench, 8-bit is normally enough with the appropriate scaling. Although I'm glad to have some guard digits, I usually use my HP-41 calculator in ENG3 (engineering notation with three digits after the first significant one) display mode.

_________________
http://WilsonMinesCo.com/ lots of 6502 resources
The "second front page" is http://wilsonminesco.com/links.html .
What's an additional VIA among friends, anyhow?


Top
 Profile  
Reply with quote  
PostPosted: Mon Apr 08, 2019 10:11 am 
Offline

Joined: Mon Jun 24, 2013 8:18 am
Posts: 83
Location: Italy
I agree with GARTH, but with an additional comment.

In some applications losing accurance in computation can give very bad pratical result; in GPS the cancellation error due to the computing differences of numbers about fractions of uS, can give a wrong position about hundreds meters (for no say about the error caused neglecting relativity effects). Another example was the loss of significance in patriot missile (incident, 1991), due to a cumulative relative error.

_________________
http://65xx.unet.bz/ - Hardware & Software 65XX family


Top
 Profile  
Reply with quote  
PostPosted: Mon Apr 08, 2019 6:27 pm 
Offline

Joined: Sat Dec 13, 2003 3:37 pm
Posts: 1004
drogon wrote:
Back to the MS BASICs though - it does make you wonder just how many programs were written (and might still be in-use!) to deal with math issues that may not be quite as accurate as it ought to be. (not counting the Pentium FDIV bug, of-course!)

I certainly don't lose sleep over it.

The simple truth is for the vast majority of applications, FP math is "good enough". Those that care about those edge case issues run in to those edge case issues, and fix things.

If someone is running MS-BASIC on a Mars Orbiter, well, that's their folly, frankly.

We've all heard horror stories about miscalculations having dire outcomes, but I don't know if I've heard any issues on subtle inaccuracies done within an FP implementation that could should have been mitigated in other ways.

We've know for a long time the problems of accumulating rounding error over time with FP (and non-FP) implementations, and there are techniques aside from underlying FP raw accuracy to mitigate those problems.

I know I had to fight such issues in the past during long running financial reports. Today, I think I would have tried different techniques to mitigate the issues I was having then, but, in the end, floating point is simply not the best option for financials if you can avoid it. Adequate, but imperfect. I had a large GL report that would just refuse to balance. Just one. Everything else worked great, but this one, large report was just a stubborn cuss.


Top
 Profile  
Reply with quote  
PostPosted: Mon Apr 15, 2019 12:54 am 
Offline
User avatar

Joined: Thu Mar 11, 2004 7:42 am
Posts: 362
The issue you tend to run into with EXP(Y*LOG(X)) is not so much how much precision is good enough, but instead the case where you're expecting an exact value (e.g. an integer) and you don't get one. Example:
Code:
10 N=300^2
20 S=10000
30 FOR I=S TO N STEP S:PRINT I:NEXT
40 FOR I=20*S TO N STEP -S:PRINT I:NEXT

Because 300^2 turns out to be a little more than the correct value, the increasing loop in line 30 works, but the decreasing loop in line 40 loops one fewer time than expected. Change 300^2 to 300*300 in line 10 and then line 40 behaves as expected.

You don't need a large result to encounter this issue. PRINT 3^6 prints 729.001. An example that you're more likely to run into is 2^15, which is slightly larger than the correct value. PRINT 2^15 = 32768 prints 0 because the two values are not equal.


Top
 Profile  
Reply with quote  
PostPosted: Mon Apr 15, 2019 4:39 pm 
Offline
User avatar

Joined: Wed Feb 14, 2018 2:33 pm
Posts: 1392
Location: Scotland
dclxvi wrote:
The issue you tend to run into with EXP(Y*LOG(X)) is not so much how much precision is good enough, but instead the case where you're expecting an exact value (e.g. an integer) and you don't get one. Example:
Code:
10 N=300^2
20 S=10000
30 FOR I=S TO N STEP S:PRINT I:NEXT
40 FOR I=20*S TO N STEP -S:PRINT I:NEXT

Because 300^2 turns out to be a little more than the correct value, the increasing loop in line 30 works, but the decreasing loop in line 40 loops one fewer time than expected. Change 300^2 to 300*300 in line 10 and then line 40 behaves as expected.

You don't need a large result to encounter this issue. PRINT 3^6 prints 729.001. An example that you're more likely to run into is 2^15, which is slightly larger than the correct value. PRINT 2^15 = 32768 prints 0 because the two values are not equal.


It's easy to over-generalise. One BASIC may not behave correctly, but there are other 8-bit,6502 BASICs that work just fine. This is BBC Basic 4:
Code:
>L.
   10 N=300^2
   20 S=10000
   25 PRINT "Loop1"
   30 FOR I=S TO N STEP S:PRINT I:NEXT
   35 PRINT "Loop2"
   40 FOR I=20*S TO N STEP -S:PRINT I:NEXT
>RUN
Loop1
     10000
     20000
     30000
     40000
     50000
     60000
     70000
     80000
     90000
Loop2
    200000
    190000
    180000
    170000
    160000
    150000
    140000
    130000
    120000
    110000
    100000
     90000
>P.3^6
       729
>P.3^6-729
         0
>


Nothings ever quite the bed of roses though:

Code:
>P. EXP (2 * LN (300))
89999.9998
>P. 300^2
     90000
>P. 300^2-90000
         0
>


I wonder what other 6502 BASICs there are, other than the generic (& enhanced) MS BASIC and BBC Basic...

-Gordon

_________________
--
Gordon Henderson.
See my Ruby 6502 and 65816 SBC projects here: https://projects.drogon.net/ruby/


Top
 Profile  
Reply with quote  
PostPosted: Mon Apr 15, 2019 5:17 pm 
Offline
User avatar

Joined: Thu Dec 11, 2008 1:28 pm
Posts: 10760
Location: England
drogon wrote:
I wonder what other 6502 BASICs there are, other than the generic (& enhanced) MS BASIC and BBC Basic...

Quite a few: see here.

It's good if X^Y can square (or cube) a negative number.


Top
 Profile  
Reply with quote  
PostPosted: Mon Apr 22, 2019 2:45 am 
Offline
User avatar

Joined: Thu Mar 11, 2004 7:42 am
Posts: 362
Different implementations are inexact in different places. 7^2 is exact for EhBASIC but not for Applesoft. 2^6 is exact for Applesoft but not for EhBASIC. 2^15 is not exact for both; 3^6 is also not exact for both.

For division, it doesn't take long for most people to realize that A=1/3 is going to be an approximation of 1/3 and not exactly 1/3. Less obvious is that .1 is also an approximation, but once you know the floating point format, that the exponent is 2^N instead of 10^N, it then becomes clear that the mantissa is a fraction whose denominator is 2^N. And at that point you more or less have the whole picture, and it's relatively simple to know when a result will be exact. The inexact cases aren't always convenient, of course, but when they are limited and predictable, it is easier to deal with them when they occur.

For exponentiation, even if you know about EXP(Y*LOG(X)), there are only a few cases that clearly will be exact regardless of implementation, e.g. 0^N or 1^N (LOG(1) will be exactly 0 and EXP(0) will be exactly 1).

Handling (some) integer exponents as a special case not only results in a speed improvement, but also gives you a range of useful cases that clearly will be exact (as long as there are enough bits in the mantissa), i.e. X^2, X^3, 2^Y, or 10^Y are among the more common cases. Of course, you'll still have cases that clearly will be an approximation, e.g. 2^(1/3).


Top
 Profile  
Reply with quote  
PostPosted: Tue Apr 30, 2019 2:27 am 
Offline

Joined: Mon May 21, 2018 8:09 pm
Posts: 1462
A common situation where loss of precision leads to subtle and hard-to-debug (but actually quite annoying) problems is in game physics.

Have you ever noticed that most games limit map size to, at most, a few kilometres across? That's the limit beyond which the precision of position vectors becomes too coarse to ensure reasonable physics behaviour - *if* the game uses a single-precision physics engine, as most do. IEEE-754 single precision can represent 1mm increments at up to 4km from the origin on any axis, and that's the level of margin often programmed in to ensure stability.

Most games that implement larger game worlds (such as flight sims) end up moving the physics origin around with the player, so that sufficient precision is retained to stop things exploding when they're on a 4000km flight. This, if done badly, results in a noticeable "bump" periodically, when the origin is moved but the physics vectors are not entirely consistent between the new and old positions. Compounding this problem is that origin moves often coincide with scenery loading.

A better solution would be to put the origin at the centre of the Earth and implement double-precision physics. This allows 1µm positioning precision out to 8 million km radius, which is vastly more than needed to present objects on or near the Earth's surface (or even in orbit around it). Modern CPUs (and even recent GPUs) have sufficiently fast double-precision implementations that the main penalty is a doubling of memory usage for physics data.


Top
 Profile  
Reply with quote  
PostPosted: Sun Oct 06, 2019 1:00 am 
Offline

Joined: Sun Feb 22, 2004 9:01 pm
Posts: 78
I think I'm probably lucky in doing O-level maths at the same time as investigating floating point coding, so the obviousness of multiplication being mantissa*mantissa and exponent+exponent, and powers being exponent*exponent came simultaneously. Some years ago I was digging my old Sallas, Hille & Anderson out when coding the floating point code in PDP11 BASIC. ;)

_________________
--
JGH - http://mdfs.net


Top
 Profile  
Reply with quote  
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 12 posts ] 

All times are UTC


Who is online

Users browsing this forum: No registered users and 1 guest


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to: