LowRes Graph C64 - integer math

; 10 SYS2061

*=$0801

        BYTE    $0B, $08, $0A, $00, $9E, $32, $30, $36, $31, $00, $00, $00


        jsr CLRSCR      ; {clrscr}

        ; reset I (I=0)
        lda #0
        sta varI
        sta varI+1

ForLoop
        ; compute P=cos(I/2)
        lda varI+1
        lsr
        sta angle+1
        lda varI
        ror
        sta angle
        ldy #1          ; cos
        jsr SinCos      ; get Cos(I/2) (value by 65535)

        sty varP
        sta varP+1
        lda #0
        rol
        sta sgntrig1    ; save sign of the trig result

        ; compute Y=24*P*Sin(I)+25
        lda varI
        sta angle
        lda varI+1
        sta angle+1
        ldy #0
        jsr SinCos      ; get Sin(I) (value by 65535)

        sty M
        sta M+1
        lda #0
        rol
        eor sgntrig1
        sta sgntrig2    ; flag with the combined signs of the trig results

        lda varP
        ldy varP+1
        sta N
        sty N+1
        jsr mult16

        ldy #3
_lp24   lda div24,y
        sta divisor,y
        dey
        bpl _lp24
        jsr div32       ; division by 65535*65535/24

        lda sgntrig2    ; check if the result has to be subtracted
        beq _sum25
        clc
        lda dividend
        eor #$ff
        ;adc #1
        sta dividend

_sum25  lda #0
        sta varL+1
        clc
        lda #25
        adc dividend
        lsr
        sta varL

        ; compute B=(Y&1)    and    40*int(Y/2)
        rol
        and #1
        sta varB

        lda varL        ; (<= 25)
        asl             ; int(Y/2) * 2  (<= 50)
        asl             ; int(Y/2) * 4  (<= 100)
        adc varL        ; int(Y/2) * 5  (<= 125)
        asl
        ;rol varL+1     ; int(Y/2) * 10 (<= 250)
        asl
        rol varL+1      ; int(Y/2) * 20
        asl
        sta varL
        rol varL+1      ; varL = int(Y/2) * 40


        ; compute X=39*P*Cos(I)+40
        lda varI
        ldy varI+1
        sta angle
        sty angle+1
        ldy #1
        jsr SinCos      ; get Cos(I) (value by 65535)

        sty M
        sta M+1
        lda #0
        rol
        eor sgntrig1
        sta sgntrig2    ; flag with the combined signs of the trig results

        lda varP
        ldy varP+1
        sta N
        sty N+1
        jsr mult16

        ldy #3
_lp39   lda div39,y
        sta divisor,y
        dey
        bpl _lp39
        jsr div32       ; division by 65535*65535/39

        lda sgntrig2    ; check if the result has to be subtracted
        beq _sum40
        clc
        lda dividend
        eor #$ff
        ;adc #1
        sta dividend

_sum40  clc
        lda #40
        adc dividend
        sta dividend

        ; compute 2^((X & 1) + 2*(Y & 1))   and    Scr + int(X/2) + 40*int(Y/2)
        lsr
        clc
        adc varL
        sta varL
        lda #4
        adc varL+1
        sta varL+1

        lsr dividend
        lda varB
        rol
        tay
        lda _pwr2,y
        sta varB        ; varB = 2 ^ (X & 1 + (2 * (Y & 1)))


        ; Draw new point on the screen
        ldy #0
        ldx #15
        lda (varL),y    ; get the byte from screen mem
_loop   cmp _plotab,x   ; seek the value in the table
        beq _found
        dex
        bpl _loop
_found  txa
        ora varB        ; add the offset for the new value
        tax
        lda _plotab,x   ; load the new value
        sta (varL),y    ; on the screen

        ; increment counter and check the loop condition
        inc varI
        bne _sk
        inc varI+1

_sk     lda varI+1      ; loops until I=$2d0 (4*Pi rad)
        cmp #2
        bne _lp
        lda varI
        cmp #$d0
        beq _end

_lp     jmp ForLoop

_end    rts


; memory references
Scr=$400
Fac1=$61
Fac2=$69

; System references
CLRSCR=$e544

; variables
varB=$2                                 ; offset in the _plotab
varL=$fd                                ; Screen Location to plot
varI    byte 00, 00                     ; counter
varP    byte 00, 00

; constants
_plotab byte  32,126,124,226,123, 97,255,236
        byte 108,127,225,251, 98,252,254,160
_pwr2   byte   1, 2, 4, 8

div24 byte $55, $95, $aa, $0a           ; 65535*65535/24 = 178951509
div39 byte $e6, $5b, $90, $06           ; 65535*65535/39 = 110124005


P    =fac1
M    =fac1+4
N    =fac1+6

mult16
        lda #$00
        sta P+2
        sta P+3
        ldx #$10
Lp      lsr N+1
        ror N
        bcc Sk
        lda P+2
        clc
        adc M
        sta P+2
        lda P+3
        adc M+1
Sk      ror
        sta P+3
        ror P+2
        ror P+1
        ror P
        dex
        bne Lp
        rts

div32   lda #0          ;preset remainder to 0
        sta remainder
        sta remainder+1
        sta remainder+2
        sta remainder+3
        ldx #32         ;repeat for each bit: ...

divloop asl dividend    ;dividend lb & hb*2, msb -> Carry
        rol dividend+1
        rol dividend+2
        rol dividend+3
        rol remainder   ;remainder lb & hb * 2 + msb from carry
        rol remainder+1
        rol remainder+2
        rol remainder+3
        lda remainder
        sec
        sbc divisor     ;substract divisor to see if it fits in
        tay             ;lb result -> Y, for we may need it later
        lda remainder+1
        sbc divisor+1
        sta pztemp
        lda remainder+2
        sbc divisor+2
        sta pztemp2
        lda remainder+3
        sbc divisor+3
        bcc skip        ;if carry=0 then divisor didn't fit in yet

        sta remainder+3 ;else save substraction result as new remainder,
        lda pztemp2
        sta remainder+2
        lda pztemp
        sta remainder+1
        sty remainder
        inc dividend    ;and INCrement result cause divisor fit in 1 times

skip    dex
        bne divloop
        rts

sgntrig1 =fac2+5
sgntrig2 =fac2+4
dividend =fac1
divisor  =fac1+4
remainder=fac2
pztemp   =$fb
pztemp2  =$fc
angle    =$fb

; *** SinCos ************
; call with Y=0 for Sin
;           Y=1 for Cos
;
; Result in A/Y (hi/lo)
;           sign in C
; ***********************
SinCos
; check angle≤90
ChkAng lda angle+1
       bne ChQuad
       lda angle
       cmp #91
       bcc FrstQ

; change quadrant
ChQuad iny
       sec
       lda angle
       sbc #90
       sta angle
       bcs ChkAng
       dec angle+1
       bcc ChQuad
FrstQ  tax
       tya
       lsr
       bcc SinX

; will get 90-α
       tay
       txa
       eor #$FF
       ;sec
       adc #90
       tax
       tya

; res: A/Y  sign: C
SinX   lsr
       lda SinHi,X
       ldy SinLo,X
       rts
SinHi   ; hibytes of Sin(angle)×65535
  byte 0, 4, 8, 13, 17, 22, 26, 31, 35, 40, 44, 48, 53
  byte 57, 61, 66, 70, 74, 79, 83, 87, 91, 95, 100, 104
  byte 108, 112, 116, 120, 124, 127, 131, 135, 139, 143
  byte 146, 150, 154, 157, 161, 164, 167, 171, 174, 177
  byte 181, 184, 187, 190, 193, 196, 198, 201, 204
  byte 207, 209, 212, 214, 217, 219, 221, 223, 226
  byte 228, 230, 232, 233, 235, 237, 238, 240, 242
  byte 243, 244, 246, 247, 248, 249, 250, 251, 252
  byte 252, 253, 254, 254, 255, 255, 255, 255, 255, 255
SinLo   ; lobytes of Sin(angle)×65535
  byte 0, 120, 239, 102, 219, 80, 194, 51, 161, 12
  byte 116, 217, 57, 150, 238, 66, 144, 217, 27, 88
  byte 142, 190, 230, 7, 31, 48, 57, 56, 47, 28, 255
  byte 217, 168, 109, 39, 213, 121, 16, 155, 27, 141
  byte 243, 75, 151, 212, 4, 38, 57, 62, 52, 27, 242
  byte 186, 115, 27, 179, 59, 178, 25, 110, 179, 230, 8
  byte 24, 22, 3, 221, 165, 91, 254, 143, 13
  byte 119, 207, 20, 70, 100, 111, 103, 75, 27
  byte 216, 129, 23, 152, 6, 95, 165, 215, 245, 255