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)
Login
Register
FAQ
Search
