; Linear regression calculator in N8A for n808 VM
; Enter the pairs number by number, end with 0,0 pair
; The program will then output A and B parameters of A + Bx
; linear function and then the correlation coefficient r
; Created by Luxferre in 2025, released into public domain

; variable/constant definitions
#1 xi    ; data x component
#2 yi    ; data y component
#3 xs    ; xi squared
#4 ys    ; yi squared
#5 xy    ; xy
#6 buf   ; buffer
#7 r1    ; r-coefficient buffer 1
#8 r2    ; r-coefficient buffer 2
#10 Sxi  ; x sum
#11 Syi  ; y sum
#12 Sxx  ; x squared sum
#13 Syy  ; y squared sum
#14 Sxy  ; xy sum
#15 n    ; data element counter
#16 A    ; coefficient A
#17 B    ; coefficient B
#18 RC   ; coefficient R   
; zero out all sums
dca 0 @Sxi
dca 0 @Syi
dca 0 @Sxx
dca 0 @Syy
dca 0 @Sxy
dca 0 @n
; data input loop
:lp inp @xi @yi ; loop start, input xi and yi pair
dva @xi @xs     ; prepare x
mul @xs @xs     ; square x
dva @yi @ys     ; prepare y
mul @ys @ys     ; square y
dva @yi @xy     ; prepare y
mul @xi @xy     ; save xy
add @xi @Sxi    ; update x sum
add @yi @Syi    ; update y sum
add @xs @Sxx    ; update x squared sum
add @ys @Syy    ; update y squared sum
add @xy @Sxy    ; update xy sum
inc @n          ; increment element count
add @ys @xs     ; add y-squared to x-squared
jgt @xs :lp     ; loop back if the square sum is over zero
dec @n          ; decrement last n to omit the (0,0) input
; processing and output part
dva @n @buf     ; n => buffer
mul @Sxy @buf   ; n * Sxy => buffer
dva @Syi @r1    ; Syi => r-buffer 1
mul @Sxi @r1    ; Sxi * Syi => r-buffer 1
sub @buf @r1    ; n * Sxy - Sxi * Syi => r-buffer 1 (to be stored)
dva @n @buf     ; n => buffer
mul @Sxx @buf   ; n * Sxx => buffer
dva @Sxi @r2    ; Sxi => r-buffer 2
mul @r2 @r2     ; Sxi squared => r-buffer 2
sub @buf @r2    ; n * Sxx - Sx^2 => r-buffer 2 (to be stored)
dva @r2 @B      ; prepare coefficient B
div @r1 @B      ; store coefficient B
dva @Sxi @buf   ; copy Sxi to buffer
mul @B @buf     ; B * Sxi => buffer
sub @Syi @buf   ; Syi - B * Sxi => buffer
dva @n @A       ; prepare coefficient A
div @buf @A     ; calculate coefficient A
dva @Syi @buf   ; Syi => buffer
mul @Syi @buf   ; Syi squared => buffer
dva @n @xi      ; reuse xi for the second buffer
mul @Syy @xi    ; n * Syy => second buffer
sub @xi @buf    ; n * Syy - Sy^2 => buffer
dva @r2 @xi     ; r-buffer 2 to xi
mul @buf @xi    ; buffer * xi => xi
sqr @xi @RC     ; sqrt(xi) => prepare RC
div @r1 @RC     ; calculate correlation coefficient
out @A @RC      ; output all three resulting numbers