; 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
:lp inp 1 2 0  ; loop start; input xi and yi into loc 1 and loc 2
dva 1 3        ; copy x into loc 3
mul 3 3        ; save x-squared into loc 3
dva 2 4        ; copy y into loc 4
mul 4 4        ; save y-squared into loc y
dva 2 5        ; copy y into loc 5
mul 1 5        ; save xy into loc 5
fma 1 127 10   ; update x-sum (Sx) => loc 10
fma 2 127 11   ; update y-sum (Sy) => loc 11
fma 3 127 12   ; update x-square-sum (Sxx) => loc 12
fma 4 127 13   ; update y-square-sum (Syy) => loc 13
fma 5 127 14   ; update xy-sum (Sxy) => loc 14
inc 15         ; increment n at loc 15
dva 3 6        ; copy x-squared into loc 6
fma 4 127 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 127 15  ; decrement last n at loc 15 to omit the (0,0) input
dva 15 1       ; copy n to loc 1
mul 14 1       ; n * Sxy => loc 1
dva 11 5       ; copy Sy to loc 5
mul 10 5       ; Sx * Sy => loc 5
sub 1 5 6      ; n * Sxy - Sx * Sy => loc 6 (to be stored for r calculation)
dva 15 3       ; copy n to loc 3
mul 12 3       ; n * Sxx => loc 3
dva 10 4       ; copy Sx to loc 4
mul 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
dva 10 1       ; copy Sx to loc 1
mul 21 1       ; B * Sx => loc 1
sub 11 1 1     ; Sy - B * Sx => loc 1
div 1 15 20    ; coefficient A => loc 20
dva 11 2       ; copy Sy to loc 2
mul 11 2       ; (Sy) ^ 2 => loc 2
dva 15 5       ; n to loc 5
mul 13 5       ; n * Syy => loc 5
sub 5 2 2      ; n * Syy - Sy^2 => loc 2
dva 3 1        ; loc 3 => loc 1
mul 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