;--------------------------------------
; Approximation of pi by checking points
; in a grid if they are lying within
; a quarter of a circle with the radius
; of the grid's length and calculating
; the ratio of points in the circle and
; all points of the grid, known as pi. 
;--------------------------------------
        !to "test.prg", cbm
        !convtab pet
        !sl "test.labels"

        ; Radius and length of the grid.
        RADIUS = 3000
        ; 32 bit fixed point value with two bit fraction.
        HIT_ZONE = (RADIUS * RADIUS) << 2

        ptr1 = $9e
        ptr2 = $a5
        x = $a7         ; 32 bit fixed point with two bit fractional part.
        hit_count = $b0 ; 24 or 32 bit counter.
        i = $f7         ; 16 bit fixed point with one bit fractional part.
        r16a = $f9

        fac1_print = $aabc

        * = $0801

basic_header
        !word .next_line
        !word 2014      ; BASIC line number.
        !pet "ti$", $b2, 34, "000000", 34, ":" ; Set clock to 00:00:00.
        !byte $9e       ; SYS token.
        !byte "0" + (start / 1000), "0" + ((start / 100) % 10)
        !byte "0" + ((start / 10) % 10), "0" + (start % 10)
        !pet ":", $99, "ti$" ; PRINT clock value.
        !byte 0         ; BASIC line end marker.
.next_line
        !word 0

!zone
start:
        jsr calculate_squares
        
        jsr count_hits_on_one_side_of_the_diagonal

        asl hit_count   ; hit_count *= 2
        rol hit_count+1
        rol hit_count+2
        
        jsr count_hits_on_the_diagonal

        ; 
        ; fac1 = 4 * hit_count / int(HIT_ZONE)
        ; 
        asl hit_count   ; hit_count *= 4
        rol hit_count+1
        rol hit_count+2
        rol hit_count+3
        asl hit_count
        rol hit_count+1
        rol hit_count+2
        rol hit_count+3

        lda hit_count+3 ; fac1 = hit_count
        sta $62
        lda hit_count+2
        sta $63
        lda hit_count+1
        sta $64
        lda hit_count
        sta $65
        sec
        lda #0
        ldx #$a0
        jsr $bc4f

        ldx #<.temp_flpt ; .temp_flpt = fac1
        ldy #>.temp_flpt
        jsr $bbd4

                        ; fac1 = int(HIT_ZONE)
        lda #(HIT_ZONE >> 2) AND $ff
        ldx #((HIT_ZONE >> 2) >> 8) AND $ff
        ldy #((HIT_ZONE >> 2) >> 16) AND $ff
        jsr $af87
        sec
        lda #0
        ldx #$a0
        jsr $bc4f

        lda #<.temp_flpt ; fac1 = .temp_flpt / fac1
        ldy #>.temp_flpt
        jsr $bb0f

        jmp fac1_print

.temp_flpt
        !fill 5

;--------------------------------------
!align 255, 0, 0
!zone
calculate_squares
        lda #<values    ; ptr1 = &values
        sta ptr1
        lda #>values
        sta ptr1+1
        ldy #1          ; i = 0.5
        sty i
        dey
        sty i+1

.loop
        lda i           ; r16a = i
        sta r16a
        lda i+1
        sta r16a+1

        lda #0          ; x = i * r16a
        sta x+2
        sta x+3
        ldx #16
-       lsr r16a+1
        ror r16a
        bcc +
        lda x+2
        clc
        adc i
        sta x+2
        lda x+3
        adc i+1
+       ror
        sta x+3
        ror x+2
        ror x+1
        ror x
        dex
        bne -

        ldy #0          ; *ptr1 = x
        lda x
        sta (ptr1),y
        iny
        lda x+1
        sta (ptr1),y
        iny
        lda x+2
        sta (ptr1),y
        iny
        lda x+3
        sta (ptr1),y

        clc             ; i += 1.0
        lda i
        adc #2
        sta i
        bcc +
        inc i+1
+
        clc             ; ptr1 += 4
        lda ptr1
        adc #4
        sta ptr1
        bcc +
        inc ptr1+1
+
        lda ptr1
        cmp #<(RADIUS * 4 + values)
        bne .loop
        lda ptr1+1
        cmp #>(RADIUS * 4 + values)
        bne .loop

        rts

;--------------------------------------
!align 255, 0, 0
!zone
count_hits_on_one_side_of_the_diagonal
        lda #0          ; hit_count = 0
        sta hit_count
        sta hit_count+1
        sta hit_count+2
        sta hit_count+3
        
        lda #<values    ; ptr1 = &values
        sta ptr1
        lda #>values
        sta ptr1+1

.outer_loop
        inc $d020

        lda ptr1
        sta ptr2
        lda ptr1+1
        sta ptr2+1
        jmp .inner_continue

; Make sure the inner loop is within one page.
!align 255, 0, 0

.inner_loop
        clc             ; x = *ptr1 + *ptr2
        ldy #0
        lda (ptr1),y
        adc (ptr2),y
        sta x
        iny
        lda (ptr1),y
        adc (ptr2),y
        sta x+1
        iny
        lda (ptr1),y
        adc (ptr2),y
        sta x+2
        iny
        lda (ptr1),y
        adc (ptr2),y
        ; sta x+3

        ; lda x+3
        cmp #(HIT_ZONE >> 24) AND $ff
        bcc +
        bne .skip
        lda x+2
        cmp #(HIT_ZONE >> 16) AND $ff
        bcc +
        bne .skip
        lda x+1
        cmp #(HIT_ZONE >> 8) AND $ff
        bcc +
        bne .skip
        lda x
        cmp #HIT_ZONE AND $ff
        bcs .skip
+
        inc hit_count   ; ... ++hit_count
        bne +
        inc hit_count+1
        bne +
        inc hit_count+2
+
.skip

.inner_continue
        clc             ; ptr2 += 4
        lda ptr2
        adc #4
        sta ptr2
        bcc +
        inc ptr2+1
+
        lda ptr2
        cmp #<(RADIUS * 4 + values)
        bne .inner_loop
        lda ptr2+1
        cmp #>(RADIUS * 4 + values)
        bne .inner_loop

.outer_continue
        clc             ; ptr1 += 4
        lda ptr1
        adc #4
        sta ptr1
        bcc +
        inc ptr1+1
+
        lda ptr1
        cmp #<(RADIUS * 4 + values)
        bne .to_outer_loop
        lda ptr1+1
        cmp #>(RADIUS * 4 + values)
        beq +
.to_outer_loop
        jmp .outer_loop
+
        rts

;--------------------------------------
!align 255, 0, 0
!zone
count_hits_on_the_diagonal
        lda #<values    ; ptr1 = &values
        sta ptr1
        lda #>values
        sta ptr1+1

.loop
        ldy #0          ; x = *ptr1 * 2
        lda (ptr1),y
        asl
        sta x
        iny
        lda (ptr1),y
        rol
        sta x+1
        iny
        lda (ptr1),y
        rol
        sta x+2
        iny
        lda (ptr1),y
        rol
        ; sta x+3

        ; lda x+3
        cmp #(HIT_ZONE >> 24) AND $ff
        bcc +
        bne .skip
        lda x+2
        cmp #(HIT_ZONE >> 16) AND $ff
        bcc +
        bne .skip
        lda x+1
        cmp #(HIT_ZONE >> 8) AND $ff
        bcc +
        bne .skip
        lda x
        cmp #HIT_ZONE AND $ff
        bcs .skip
+
        inc hit_count   ; ... ++hit_count
        bne +
        inc hit_count+1
        bne +
        inc hit_count+2
+
.skip
        clc             ; ptr1 += 4
        lda ptr1
        adc #4
        sta ptr1
        bcc +
        inc ptr1+1
+
        lda ptr1+1
        cmp #>(RADIUS * 4 + values)
        bcc .loop
        lda ptr1
        cmp #<(RADIUS * 4 + values)
        bcc .loop

        rts

;--------------------------------------
; Variables:

; Make sure no indexed access to a 32 bit value crosses page boundaries.
!align 4, 0, 0  
; 32 bit fixed point values with 2 bits fractional part.
values

; Mind the BASIC ROM.
!if values + RADIUS * 4 > $9fff {
        !error "`values` array does not fit into memory"
}