6502.org Forum  Projects  Code  Documents  Tools  Forum
It is currently Fri Apr 19, 2024 10:59 pm

All times are UTC




Post new topic Reply to topic  [ 15 posts ] 
Author Message
 Post subject: Sine calculation
PostPosted: Mon Dec 04, 2006 7:38 pm 
Offline

Joined: Wed Nov 24, 2004 6:24 pm
Posts: 26
Location: Adelaide
Here's my attempt at making a sine function. If you have any better ones or details on technique please post!

It takes two bytes in zero page for the 16-bit angle (0-360), called THE here. Returns with sine in the accumulator, based on the table. This table only goes from 0-31. You can use a spreadsheet to calculate it for something else.

Code:
SINE    LDA THE+1
        CMP #1
        BEQ SINEH
        LDA THE
        CMP #180
        BCS SIN180
        CMP #90
        BCC SIN0
SIN90   LDA #180
        CLC
        SBC THE
SIN0    TAX
        LDA SINES,X
        RTS
SIN180  LDA THE
        SBC #176
SIN180A TAX
        LDA SINES,X
        EOR #$FF
        CLC
        ADC #01
        RTS
SINEH   LDA THE
        CMP #14
        BCC SIN270
        LDA #105
        CLC
        SBC THE
        JMP SIN180A
SIN270  CLC
        SBC #179
        JMP SIN180A

SINES   .BYT $00 $00 $01 $01 $02 $02 $03 $03 $04 $05
        .BYT $05 $06 $06 $07 $07 $08 $08 $09 $09 $0A
        .BYT $0A $0B $0B $0C $0D $0D $0E $0E $0F $0F
        .BYT $0F $10 $10 $11 $11 $12 $12 $13 $13 $14
        .BYT $14 $14 $15 $15 $16 $16 $17 $17 $17 $18
        .BYT $18 $18 $19 $19 $19 $1A $1A $1A $1B $1B
        .BYT $1B $1B $1C $1C $1C $1D $1D $1D $1D $1D
        .BYT $1E $1E $1E $1E $1E $1E $1F $1F $1F $1F
        .BYT $1F $1F $1F $1F $1F $1F $1F $1F $1F $1F

_________________
Check out my 8080 project: http://kaput.homeunix.org


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Mon Dec 04, 2006 8:33 pm 
Offline
User avatar

Joined: Fri Aug 30, 2002 1:09 am
Posts: 8423
Location: Southern California
Quote:
If you have any better ones or details on technique please post!

What's best depends on what you're after. There are many ways to do it, with the trade-offs being between speed, code size and complexity, precision & accuracy, and whether you already have things like divide and multiply routines in place. Without taking the time to proofread and test your code, I'd have to say it's almost at the extreme of shortness and speed, at the expense of both precision and accuracy.

Getting good precision and accuracy does not require floating-poing however. With 16-bit cells and double-precision (ie, 32-bit) intermediate results, you can represent the full circle in 65,536 increments, which works out nicely with the carry because for example 359°+2°=361° which is the same angle as 1°. The angle can be considered to be represented either signed or unsigned, since for example 359° is the same angle as -1°. Each increment in the 65,536-part circle represents about 1/182 of a degree, or about 20 arc-seconds. If you scale the output also so that ±1 is represented by ±7FFF (in hex), then each increment there represents .00003. Then using the first four terms of the infinite series for sine and cosine in the range of 0-45° but with the coefficients tweaked a little to make up for the fact that you're using only four terms, your output will usually be correct to 16 bits (actually 15 plus sign), and never be off by more than one lsb. Then, with the sin & cos functions from 0 to 45°, you can get the rest of the circle.

If you care to go to this level and either have a higher-level language (even if integer only) or use the multiply and divide routines from the source code repository, I can give you the rest. I started the topic on fixed-point and scaled integer math at viewtopic.php?t=716 but have not gone very far with it yet. [Edit: I took it in a slightly different direction elsewhere. See my article on scaled-integer math and large look-up tables for hyperfast, accurate 16-Bit scaled-integer math, including trig & log functions. You can probably implement them even if your computer is already built up, the address space is full, and your I/O is almost all taken. See how. There in no interpolating required, and the results will be accurate to all 16 bits. In some cases it can be nearly a thousand times as fast as actually calculating the results.]

_________________
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  
 Post subject:
PostPosted: Mon Dec 04, 2006 10:05 pm 
Offline

Joined: Fri Aug 30, 2002 2:05 pm
Posts: 347
Location: UK
Here's the code I use ..

Code:
;----------------------------------------------------------------
;
; SIN(A) COS(A) routines. a full circle is represented by $00 to
; $00 in 256 1.40625 degree steps. returned value is signed 15
; bit with X being the high byte and ranges over +/-0.99997

;----------------------------------------------------------------
;
; get COS(A) in AX

cos_A
      CLC                     ; clear carry for add
      ADC   #$40              ; add 1/4 rotation

;----------------------------------------------------------------
;
; get SIN(A) in AX. enter with flags reflecting the contents of A

sin_A
      BPL   sin_cos           ; just get SIN/COS and return if +ve

      AND   #$7F              ; else make +ve
      JSR   sin_cos           ; get SIN/COS
                              ; now do twos complement
      EOR   #$FF              ; toggle low byte
      CLC                     ; clear carry for add
      ADC   #$01              ; add to low byte
      PHA                     ; save low byte
      TXA                     ; copy high byte
      EOR   #$FF              ; toggle high byte
      ADC   #$00              ; add carry from low byte
      TAX                     ; copy back to X
      PLA                     ; restore low byte
      RTS

;----------------------------------------------------------------
;
; get AX from SIN/COS table

sin_cos
      CMP   #$41              ; compare with max+1
      BCC   quadrant          ; branch if less

      EOR   #$7F              ; wrap $41 to $7F ..
      ADC   #$00              ; .. to $3F to $00
quadrant
      ASL                     ; * 2 bytes per value
      TAX                     ; copy to index
      LDA   sintab,X          ; get SIN/COS table value low byte
      PHA                     ; save it
      LDA   sintab+1,X        ; get SIN/COS table value high byte
      TAX                     ; copy to X
      PLA                     ; restore low byte
      RTS

;----------------------------------------------------------------
;
; SIN/COS table. returns values between $0000 and $7FFF

sintab
      .word $0000,$0324,$0647,$096A,$0C8B,$0FAB,$12C8,$15E2
      .word $18F8,$1C0B,$1F19,$2223,$2528,$2826,$2B1F,$2E11
      .word $30FB,$33DE,$36BE,$398C,$3C56,$3F17,$41CE,$447A
      .word $471C,$49B4,$4C3F,$4EBF,$5133,$539B,$55F5,$5842
      .word $5A82,$5CB4,$5ED7,$60EC,$62F2,$64EB,$66CF,$68A6
      .word $6A6D,$6C24,$6DC4,$6F5F,$70E2,$7255,$73B5,$7504
      .word $7641,$776C,$7884,$798A,$7A7D,$7B5D,$7C2A,$7CE3
      .word $7D8A,$7E1D,$7E9D,$7F09,$7F62,$7FA7,$7FD8,$7FF6
      .word $7FFF

;----------------------------------------------------------------

Lee.


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Tue Dec 05, 2006 2:40 am 
Offline

Joined: Wed Nov 24, 2004 6:24 pm
Posts: 26
Location: Adelaide
GARTHWILSON wrote:
Without taking the time to proofread and test your code, I'd have to say it's almost at the extreme of shortness and speed, at the expense of both precision and accuracy.


Indeed. What I was trying to do was something that would return sine multiplied by a constant, in this case 31. So 45 degrees would return 0.5*31, 90 would be 1*31 etc.

_________________
Check out my 8080 project: http://kaput.homeunix.org


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Tue Dec 05, 2006 3:13 am 
Offline

Joined: Fri Aug 30, 2002 2:05 pm
Posts: 347
Location: UK
You can do the same with the code i use. For a range of 0 to 31 just do ..

Code:
      JSR   sin_A             ; or cos_A
      TXA                     ; copy high byte, discard low byte
      CMP   #$80              ; setup for ASR
      ROR                     ; A / 2
      CMP   #$80              ; setup for ASR again
      ROR                     ; A / 4, leaves LSbit in carry

This will give you +/- $00 to $1F with an additional bit in the carry if you want or need it.

You can also increase the resolution to 9 bits, or 0.703125 degree steps, by interpolating the b0 = 1 values from adjacent b0 = 0 values.

Lee.


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Wed Dec 06, 2006 10:51 pm 
Offline

Joined: Tue Nov 18, 2003 8:41 pm
Posts: 250
leeeeee wrote:
Here's the code I use ..

(snipped some)

Quote:
Code:
sin_A
      BPL   sin_cos           ; just get SIN/COS and return if +ve

      AND   #$7F              ; else make +ve
      JSR   sin_cos           ; get SIN/COS
                              ; now do twos complement
      EOR   #$FF              ; toggle low byte
      CLC                     ; clear carry for add
      ADC   #$01              ; add to low byte
      PHA                     ; save low byte
      TXA                     ; copy high byte
      EOR   #$FF              ; toggle high byte
      ADC   #$00              ; add carry from low byte
      TAX                     ; copy back to X
      PLA                     ; restore low byte
      RTS

;----------------------------------------------------------------
;
; get AX from SIN/COS table

sin_cos
      CMP   #$41              ; compare with max+1
      BCC   quadrant          ; branch if less

      EOR   #$7F              ; wrap $41 to $7F ..
      ADC   #$00              ; .. to $3F to $00
quadrant
      ASL                     ; * 2 bytes per value
      TAX                     ; copy to index
      LDA   sintab,X          ; get SIN/COS table value low byte
      PHA                     ; save it
      LDA   sintab+1,X        ; get SIN/COS table value high byte
      TAX                     ; copy to X
      PLA                     ; restore low byte
      RTS

Lee.


I rearranged your code for no better reason than it's not the way I would have done it.
(I think mine's prettier :) )
(I didn't actually test it to see if I got it right, so take it with a grain of salt ;) )

Code:

 ASL                ; angle x 2 for indexing 2 consecutive
                    ; bytes and to move the upper 2 bits
                    ; to where we can test for
                    ; which quadrant we're in

 BPL skip           ; even quadrant? if so skip complement

 PHP                ; odd quadrant save carry to test sine sign
 EOR #FE            ; and complement the angle in 64 (x 2)
 CLC                ;
 ADC #02            ; lsb = 0 leave it alone
 PLP                ; restore the carry

skip
 TAX                ; copy angle to x for indexing
 BCC pos_sine       ; do we want a positive sine?       

 LDA #00            ; negative sine so subtract from 0
 SBC sinetab,X      ; to negate
 PHA
 LDA #00
 SBC sinetab+1,X
 TAX
 PLA
 RTS
         
pos_sine
 LDA sinetab,X
 PHA
 LDA sinetab+1,x
 TAX
 PLA
 RTS



now I'm curious if there are unstated reasons for the way you did it
(like it's easier to add interpolation if you have a subroutine that just
returns positive values, or something like that) ?


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Wed Dec 06, 2006 11:46 pm 
Offline

Joined: Fri Aug 30, 2002 2:05 pm
Posts: 347
Location: UK
Quote:
now I'm curious if there are unstated reasons for the way you did it

Like Topsy it just growed. First make a routine that returns the first quadrant, then add code to handle the first two quadrants. Next add the complement code to handle the other two quadrants and finally add the COS() entry point.

The code to interpolate for a 9 bit entry value is another step up.

I keep it like that because its simple to change to suit my needs. If I change the code to get the value from the table I only have one routine to change and debug, you have two.

Lee.


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Mon Jan 29, 2007 8:11 am 
Offline

Joined: Sun Aug 24, 2003 7:17 pm
Posts: 111
The only systematic approach to mathematics is to use floating point. On this WEB page you have http://6502.org/source/floats/wozfp1.txt with floating point routines that can be used for this. But unfortunately SIN is missing there, the only functions supported are LOG and EXP. The approach to compute

J = 32767 * SIN( PI * I/32767) -32767 <= I <= 32767

would normally be

1 - Compute F = 32767 as floating point number using subroutine FLOAT
2 - Compute F1 = 31416 as floating point number using subroutine FLOAT
3 - Compute F2 = 10000 as floating point number using subroutine FLOAT
4 - Compute G = (F1/F2) /F as floating point number using subroutine FDIV twice

With the floating point numbers F and G available the procedure would then be:

1 - Compute F1 = I as floating point number using subroutine FLOAT
2 - Compute F2 = G * F1 using subroutine FMUL
3 - Compute F3 = SIN(F2) as floating point number using the "missing" routine
4 - Compute F4 = F * F3 as floating point number using FMUL
5 - Compute 16 bit signed integer J as J = F4 using subroutine FIX

But this is presently theory as the SIN function is missing. The point is, floating point is there to be used!


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Mon Jan 29, 2007 9:37 am 
Offline
User avatar

Joined: Fri Aug 30, 2002 1:09 am
Posts: 8423
Location: Southern California
Quote:
The only systematic approach to mathematics is to use floating point.

I'm not sure what you mean by "systematic," but as I said in my post above, I'm doing sines in 16-bit scaled-integer (actually 15 bits plus sign), and getting all 16 bits correct in most cases, and never getting an error of more than 1 lsb. I used to be very skeptical about going without floating-point, but now I'm a believer, for most (not all) applications.


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Mon Jan 29, 2007 10:48 am 
Offline

Joined: Wed Oct 22, 2003 4:07 am
Posts: 51
Location: Norway
Why not keep it simple?
Code:
; 8 bit angle in x, 8 bit result in A
lda    sinlut,x

...or if you need more precision....
Code:
; 9 bit angle, sign in carry, remaining 8 bits in x
; 9 bit result, sign in carry, remaining 8 bits in A (as a Sign and Magnitude number)
lda    sinlut,x

:lol:


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Mon Jan 29, 2007 1:09 pm 
Offline

Joined: Sun Aug 24, 2003 7:17 pm
Posts: 111
But is there any reason to try to avoid floating point? Just to make some advertisment for floating point, see the code for the floating point computation proposed in my post!

Code:
SIGN = $00
X2   = $01
M2   = $02
X1   = $05
M1   = $06
E    = $09
Z    = $0D
T    = $11
SEXP = $15
INT  = $19
      *=$0300
;CONVERT I1 TO FLOAT
     LDA I1
     STA M1
     LDA I1+1
     STA M1+1
     JSR FLOAT
;MOVE TO X2,M2 POSITION
     LDX #$04
L1   LDA X1-1,X
     STA X2-1,X
     DEX
     BNE L1
;CONVERT I2 TO FLOAT
     LDA I2
     STA M1
     LDA I2+1
     STA M1+1
     JSR FLOAT
     JSR FDIV
;MOVE TO X2,M2 POSITION
     LDX #$04
L2   LDA X1-1,X
     STA X2-1,X
     DEX
     BNE L2
;CONVERT IFF TO FLOAT
     LDA IFF
     STA M1
     LDA IFF+1
     STA M1+1
     JSR FLOAT
;SAVE FLOATING POINT NUMBER F
     LDX #$04
L3   LDA X1-1,X
     STA F-1,X
     DEX
     BNE L3
     JSR FDIV
;SAVE FLOATING POINT NUMBER G
     LDX #$04
L4   LDA X1-1,X
     STA G-1,X
     DEX
     BNE L4
;   
;COMPUTE PI*I/32767
;
;CONVERT I TO FLOAT
     LDA I
     STA M1
     LDA I+1
     STA M1+1
     JSR FLOAT
;RETRIEVE FLOATING POINT NUMBER G
     LDX #$04
L5   LDA G-1,X
     STA X2-1,X
     DEX
     BNE L5
     JSR FMUL
;    JSR SIN   ;MISSING FUNCTION
;RETRIEVE FLOATING POINT NUMBER F
     LDX #$04
L6   LDA F-1,X
     STA X2-1,X
     DEX
     BNE L6
     JSR FMUL
     JSR FIX
     BRK
IFF   .DB  $7F  ;IFF = 32767
      .DB  $FF
I1    .DB  $7A  ;I1=31416
      .DB  $B8
I2    .DB  $27  ;I2=10000
      .DB  $10
I     .DB  $16  ;I=5825 (32 DEG)
      .DB  $C1
F     *=*+4
F1    *=*+4
F2    *=*+4
G     *=*+4
;     BASIC FLOATING POINT ROUTINES
;
ADD    CLC         ;CLEAR CARRY
       LDX #$02    ;INDEX FOR 3-BYTE ADD
ADD1   LDA M1,X
       ADC M2,X    ;ADD A BYTE OF MANT2 TO MANT1
       STA M1,X
       DEX         ;ADVANCE INDEX TO NEXT MORE SIGNIF.BYTE
       BPL ADD1    ;LOOP UNTIL DONE.
       RTS         ;RETURN
MD1    ASL SIGN    ;CLEAR LSB OF SIGN
       JSR ABSWAP  ;ABS VAL OF MANT1, THEN SWAP MANT2
ABSWAP BIT M1      ;MANT1 NEG?
       BPL ABSWP1  ;NO,SWAP WITH MANT2 AND RETURN
       JSR FCOMPL  ;YES, COMPLIMENT IT.
       INC SIGN    ;INCR SIGN, COMPLEMENTING LSB
ABSWP1 SEC         ;SET CARRY FOR RETURN TO MUL/DIV
;
;     SWAP EXP/MANT1 WITH EXP/MANT2
;
SWAP   LDX #$04    ;INDEX FOR 4-BYTE SWAP.
SWAP1  STY E-1,X
       LDA X1-1,X  ;SWAP A BYTE OF EXP/MANT1 WITH
       LDY X2-1,X  ;EXP/MANT2 AND LEAVEA COPY OF
       STY X1-1,X  ;MANT1 IN E(3BYTES). E+3 USED.
       STA X2-1,X
       DEX         ;ADVANCE INDEX TO NEXT BYTE
       BNE SWAP1   ;LOOP UNTIL DONE.
       RTS
;
;
;
;     CONVERT 16 BIT INTEGER IN M1(HIGH) AND M1+1(LOW) TO F.P.
;     RESULT IN EXP/MANT1.  EXP/MANT2 UNEFFECTED
;
;
FLOAT  LDA #$8E
       STA X1      ;SET EXPN TO 14 DEC
       LDA #$00    ;CLEAR LOW ORDER BYTE
       STA M1+2
       BEQ NORM    ;NORMALIZE RESULT
NORM1  DEC X1      ;DECREMENT EXP1
       ASL M1+2
       ROL M1+1    ;SHIFT MANT1 (3 BYTES) LEFT
       ROL M1
NORM   LDA M1      ;HIGH ORDER MANT1 BYTE
       ASL         ;UPPER TWO BITS UNEQUAL?
       EOR M1
       BMI RTS1    ;YES,RETURN WITH MANT1 NORMALIZED
       LDA X1      ;EXP1 ZERO?
       BNE NORM1   ;NO, CONTINUE NORMALIZING
RTS1   RTS         ;RETURN
;
;
;     EXP/MANT2-EXP/MANT1 RESULT IN EXP/MANT1
;
FSUB   JSR FCOMPL  ;CMPL MANT1 CLEARS CARRY UNLESS ZERO
SWPALG JSR ALGNSW  ;RIGHT SHIFT MANT1 OR SWAP WITH MANT2 ON CARRY
;
;     ADD EXP/MANT1 AND EXP/MANT2 RESULT IN EXP/MANT1
;
FADD   LDA X2
       CMP X1      ;COMPARE EXP1 WITH EXP2
       BNE SWPALG  ;IF UNEQUAL, SWAP ADDENDS OR ALIGN MANTISSAS
       JSR ADD     ;ADD ALIGNED MANTISSAS
ADDEND BVC NORM    ;NO OVERFLOW, NORMALIZE RESULTS
       BVS RTLOG   ;OV: SHIFT MANT1 RIGHT. NOTE CARRY IS CORRECT SIGN
ALGNSW BCC SWAP    ;SWAP IF CARRY CLEAR, ELSE SHIFT RIGHT ARITH.
RTAR   LDA M1      ;SIGN OF MANT1 INTO CARRY FOR
       ASL         ;RIGHT ARITH SHIFT
RTLOG  INC X1      ;INCR EXP1 TO COMPENSATE FOR RT SHIFT
       BEQ OVFL    ;EXP1 OUT OF RANGE.
RTLOG1 LDX #$FA    ;INDEX FOR 6 BYTE RIGHT SHIFT
ROR1   LDA #$80
       BCS ROR2
       ASL
ROR2   LSR E+3,X   ;SIMULATE ROR E+3,X
       ORA E+3,X
       STA E+3,X
       INX         ;NEXT BYTE OF SHIFT
       BNE ROR1    ;LOOP UNTIL DONE
       RTS         ;RETURN
;
;
;     EXP/MANT1 X EXP/MANT2 RESULT IN EXP/MANT1
;
FMUL   JSR MD1     ;ABS. VAL OF MANT1, MANT2
       ADC X1      ;ADD EXP1 TO EXP2 FOR PRODUCT EXPONENT
       JSR MD2     ;CHECK PRODUCT EXP AND PREPARE FOR MUL
       CLC         ;CLEAR CARRY
MUL1   JSR RTLOG1  ;MANT1 AND E RIGHT.(PRODUCT AND MPLIER)
       BCC MUL2    ;IF CARRY CLEAR, SKIP PARTIAL PRODUCT
       JSR ADD     ;ADD MULTIPLICAN TO PRODUCT
MUL2   DEY         ;NEXT MUL ITERATION
       BPL MUL1    ;LOOP UNTIL DONE
MDEND  LSR SIGN    ;TEST SIGN (EVEN/ODD)
NORMX  BCC NORM    ;IF EXEN, NORMALIZE PRODUCT, ELSE COMPLEMENT
FCOMPL SEC         ;SET CARRY FOR SUBTRACT
       LDX #$03    ;INDEX FOR 3 BYTE SUBTRACTION
COMPL1 LDA #$00    ;CLEAR A
       SBC X1,X    ;SUBTRACT BYTE OF EXP1
       STA X1,X    ;RESTORE IT
       DEX         ;NEXT MORE SIGNIFICANT BYTE
       BNE COMPL1  ;LOOP UNTIL DONE
       BEQ ADDEND  ;NORMALIZE (OR SHIFT RIGHT IF OVERFLOW)
;
;
;     EXP/MANT2 / EXP/MANT1 RESULT IN EXP/MANT1
;
FDIV   JSR MD1     ;TAKE ABS VAL OF MANT1, MANT2
       SBC X1      ;SUBTRACT EXP1 FROM EXP2
       JSR MD2     ;SAVE AS QUOTIENT EXP
DIV1   SEC         ;SET CARRY FOR SUBTRACT
       LDX #$02    ;INDEX FOR 3-BYTE INSTRUCTION
DIV2   LDA M2,X
       SBC E,X     ;SUBTRACT A BYTE OF E FROM MANT2
       PHA         ;SAVE ON STACK
       DEX         ;NEXT MORE SIGNIF BYTE
       BPL DIV2    ;LOOP UNTIL DONE
       LDX #$FD    ;INDEX FOR 3-BYTE CONDITIONAL MOVE
DIV3   PLA         ;PULL A BYTE OF DIFFERENCE OFF STACK
       BCC DIV4    ;IF MANT2<E THEN DONT RESTORE MANT2
       STA M2+3,X
DIV4   INX         ;NEXT LESS SIGNIF BYTE
       BNE DIV3    ;LOOP UNTIL DONE
       ROL M1+2
       ROL M1+1    ;ROLL QUOTIENT LEFT, CARRY INTO LSB
       ROL M1
       ASL M2+2
       ROL M2+1    ;SHIFT DIVIDEND LEFT
       ROL M2
       BCS OVFL    ;OVERFLOW IS DUE TO UNNORMALIZED DIVISOR
       DEY         ;NEXT DIVIDE ITERATION
       BNE DIV1    ;LOOP UNTIL DONE 23 ITERATIONS
       BEQ MDEND   ;NORMALIZE QUOTIENT AND CORRECT SIGN
MD2    STX M1+2
       STX M1+1    ;CLR MANT1 (3 BYTES) FOR MUL/DIV
       STX M1
       BCS OVCHK   ;IF EXP CALC SET CARRY, CHECK FOR OVFL
       BMI MD3     ;IF NEG NO UNDERFLOW
       PLA         ;POP ONE
       PLA         ;RETURN LEVEL
       BCC NORMX   ;CLEAR X1 AND RETURN
MD3    EOR #$80    ;COMPLIMENT SIGN BIT OF EXP
       STA X1      ;STORE IT
       LDY #$17    ;COUNT FOR 24 MUL OR 23 DIV ITERATIONS
       RTS         ;RETURN
OVCHK  BPL MD3     ;IF POS EXP THEN NO OVERFLOW
OVFL   BRK
;
;
;     CONVERT EXP/MANT1 TO INTEGER IN M1 (HIGH) AND M1+1(LOW)
;      EXP/MANT2 UNEFFECTED
;
       JSR RTAR    ;SHIFT MANT1 RT AND INCREMENT EXPNT
FIX    LDA X1      ;CHECK EXPONENT
       CMP #$8E    ;IS EXPONENT 14?
       BNE FIX-3   ;NO, SHIFT
RTRN   RTS         ;RETURN


This is it! I see no reason to avoid floating point!

Code:
31416            $8E7AB8
10000            $8D4E20
PI               $816487FC

32767            $8E7FFF
PI / 32767       $726488C4

PI * I / 32767   $7F477C5F (for I=5825)

(PI *I / 32767) * 32767      $8E477BD0

(= PI * 5825 = 18299 as the missing "SIN" is commented out)


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Mon Jan 29, 2007 6:35 pm 
Offline
User avatar

Joined: Fri Aug 30, 2002 1:09 am
Posts: 8423
Location: Southern California
Quote:
But is there any reason to try to avoid floating point?

Just that it is much slower for a processor that does not have floating-point hardware, and takes more memory space. In the first embedded application I wrote, I did write a set of floating-point routines, thinking it was the only way to do what we wanted. Later I was introduced to the world of fixed-point and scaled-integer math and found out how much faster a 6502 could be without giving up much at all. My 5MHz 65c02 workbench computer with 16K of RAM does a 2K complex FFT in about 5 seconds in Forth (not even assembly language), using fixed-point/scaled-integer. When I go to a 16MHz 65816 and optimize the routine, I know I can get it under half a second.

Thankyou for posting your FP material though. I may make use of it someday. What I wrote before was in decimal because of how much human interface there was, and I didn't want to have to be continually converting between decimal and hex.


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Tue Jan 30, 2007 10:14 am 
Offline

Joined: Sun Aug 24, 2003 7:17 pm
Posts: 111
My conclusion is rather that it does not matter that the 6502 only has 8 bit add (ADC and SBC) in hardware, i.e. no floating point operations, as this can be replaced by a few routines. For my work I often do heavy computations using double precision floating point (64 bit) to propagate spacecraft orbits numerically over a long time. For this I am happy to use a high performance workstation rather then a 6502 with floating point in software running at 1 MHz(or 8 MHz)! But for any (control) application were you seriously would consider to use a 6502 I cannot imagine that such a high speed is needed or even useful!

Obviously it is nothing new that this can be done with a 6502, Commodore & Apple made this 30 years ago! But we need the source code for this to be able to easily burn corresponding software into our ROMs.

I tried to "hack" the Commodore floating point code out of my PETs ROM. Commodore uses 4 bits mantissa and a format that (essentially) is the modern IEEE format. But this "hacking" is more difficult then to write the code new from scratch what I also have done. My code for Commodore floating point format is straightforward and correct but tests have shown that the orignal code often is faster, sometimes in the order of 30%. They certainly use some clever "tricks" to save cycles! It would certainly be useful/interesting to have the original source code (with comments!) for Commodore floating point software available as well. With SIN and the other standard functions which all are available in BASIC.


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Tue Jan 30, 2007 10:18 am 
Offline

Joined: Sun Aug 24, 2003 7:17 pm
Posts: 111
Sorry
Quote:
4 bits mantissa
, read 4 byte mantissa


Top
 Profile  
Reply with quote  
 Post subject:
PostPosted: Tue Jan 30, 2007 5:08 pm 
Offline

Joined: Sat Jan 04, 2003 10:03 pm
Posts: 1706
The reason Commodore's routines are 30% faster is because they pack everything into 4 bits. That allows them to work with the whole number in a single CPU instruction. Of course, this would be *blindingly* fast, were it not for the REST of the instructions, which compress and decompress 24 bit mantissas into those 4 bits. ;D Nobody has yet reverse engineered that compression algorithm too. ;D

(Sorry, Mats! I just had to run with it.)


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

All times are UTC


Who is online

Users browsing this forum: No registered users and 4 guests


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: