; Linear regression calculator in MU8A for mu808 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

set 0 0 10     ; x-sum to loc 10
set 0 0 11     ; y-sum to loc 11
set 0 0 12     ; x-square-sum to loc 12
set 0 0 13     ; y-square-sum to loc 13
set 0 0 14     ; xy-sum to loc 14
set 0 0 15     ; n to loc 15
set 1 0 50     ; constant 1 to loc 50
:lp inp 1 2 0  ; loop start; input xi and yi into loc 1 and loc 2
cpy 1 1 3      ; copy x into loc 3
fma 0 3 3      ; save x-squared into loc 3
cpy 1 2 4      ; copy y into loc 4
fma 0 4 4      ; save y-squared into loc y
cpy 1 2 5      ; copy y into loc 5
fma 0 1 5      ; save xy into loc 5
fma 1 50 10    ; update x-sum (Sx) => loc 10
fma 2 50 11    ; update y-sum (Sy) => loc 11
fma 3 50 12    ; update x-square-sum (Sxx) => loc 12
fma 4 50 13    ; update y-square-sum (Syy) => loc 13
fma 5 50 14    ; update xy-sum (Sxy) => loc 14
fma 50 50 15   ; increment n at loc 15
cpy 1 3 6      ; copy x-squared into loc 6
fma 4 50 6     ; add x-squared and y-squared into loc 6
jmp 1 6 :lp    ; loop back to the input if the square sum is over zero
sub 15 50 15   ; decrement last n at loc 15 to omit the (0,0) input
cpy 1 15 1     ; copy n to loc 1
fma 0 14 1     ; n * Sxy => loc 1
cpy 1 11 5     ; copy Sy to loc 5
fma 0 10 5     ; Sx * Sy => loc 5
sub 1 5 6      ; n * Sxy - Sx * Sy => loc 6 (to be stored for r calculation)
cpy 1 15 3     ; copy n to loc 3
fma 0 12 3     ; n * Sxx => loc 3
cpy 1 10 4     ; copy Sx to loc 4
fma 0 10 4     ; Sx squared => loc 4
sub 3 4 3      ; n * Sxx - Sx^2 => loc 3 (to be stored for r calculation)
div 6 3 21     ; coefficient B => loc 21
cpy 1 10 1     ; copy Sx to loc 1
fma 0 21 1     ; B * Sx => loc 1
sub 11 1 1     ; Sy - B * Sx => loc 1
div 1 15 20    ; coefficient A => loc 20
cpy 1 11 2     ; copy Sy to loc 2
fma 0 11 2     ; (Sy) ^ 2 => loc 2
cpy 1 15 5     ; n to loc 5
fma 0 13 5     ; n * Syy => loc 5
sub 5 2 2      ; n * Syy - Sy^2 => loc 2
cpy 1 3 1      ; loc 3 => loc 1
fma 0 2 1      ; loc 3 * loc 2 => loc 1
sqr 0 1 1      ; sqrt(loc 1) => loc 1
div 6 1 22     ; correlation coefficient r => loc 22
out 20 22 0    ; output all three resulting numbers