;------------- POČETAK - PODNAPUTBINE ZA REDONIZOVE, NIZOVE, RAZLOMKE ----------------------------------------------
; Copyright (C) 2018 Miljenko Hatlak (hak_vz)
; https://forums.autodesk.com/t5/user/viewprofilepage/user-id/5530556
; https://forums.autodesk.com/t5/visual-lisp-autolisp-and-general/sharing-my-matrix-and-fractions-library/td-p/8164199
; It provides matrix and vectors arithmetic, LU and QR decompositions, linear and least squares solver, matrix inversion and Moore-Penrose pseudoinverse,
; adjoint, fractions arithmetics, conversions between decimal and fraction representation of real numbers, continuous fractions, determinant
; and a lots of helper functions and more.
; Fractions calculations:
; Fractions calculations are performed using continued fractions. If you round your decimals to fixed number of digits you will always receive exact result. For larger decimal numbers calculate fraction for decimal part only (restrictions of continuous fractions and internal integer representation in autolisp).
; Decimal number is in its fraction form represented as a list i.e. 0.5 (1 2). It may also be represented as a list of two decimal number (2.5 5), a combination of a decimal number and a list (2.5 ( 5 1)), or a combination of a lists – double fractions ((5 2)(5 1))
; Conversions from decimal to fraction representation of a rational number.
; (frac 0.5) → (1 2)
; (frac '(2.5 5)) → (1 2)
; (frac '(2.5 ( 5 1)) → (1 2)
; (frac '((5 2)(5 1))) → (1 2)
; (frac 0) → (0 0)
; (to fractions '(0.5 0.125 0.144 0.576)) → ((1 2) (1 8) (18 125) (72 125))
; Conversions from fraction to decimal representation of a rational number.
; (decimal '(1 2)) → 0.5
; (to_decimal '((1 2) (1 8) (18 125) (72 125))) ->((1.00000 2.00000) (1.00000 8.00000) (18.0000 125.000) (72.0000 125.000))
; Conversing between two representations is performed using continuous fractions so decimal number 0.128 is represented in continuous fraction form as ((1) (0 7 1 4 3)), and -0.125 in a form ((-1) (0 8)). First sublist represents positive or negative number and second sublist contains continuous fractions representation of decimal number.
; (con_frac 0.125) ->((-1) (0 8))
; (con_frac_to_frac '((-1) (1 8 2 5))) → (104 93)
; (decimal '(104 93)) → 1.8828
; Correct result is 1,1182795698924731182795698924731 but autocad shows only five significant digits although if uses double precision that is correct to at least 10-14 . To fetch more significant digits you may use (rtos(decimal '(104 93))2 15) → "1.118279569892473". In technical applications more than five significant digits is rarely needed.
; (con_frac_to_frac(con_frac (sqrt 2))) → (47321 33461)
; (con_frac_to_frac(con_frac pi )) → (80143857 25510582)
; (con_frac_to_frac(con_frac 2.718281828459045235360287)) → (25946 9545)
; Fractions operations f+, f-, f*, f/:
; (f+ 0.125 '(1 2)) → (5 8)
; (decimal (f+ 0.125 '(1 2))) → 0.625
; (f+ '((1 2)(3 4))'((-1 4)(2 3))) → (7 24)
; (f+ '((1 2.5)(3.2 4))'((-1.25 4)(2.25 3))) → (1 12)
; (f- 0.12 0.132) → (-3 250)
; (f/ 0 0) → error: Undefined (/ 0 0)
; (f/ 3 0) → error: Cannot divide by zero
; In fractions representations zero is represented as (0 0) but in calculations operations are restricted regarding operations with zero.
; Full strength of this conversations will come to light later in a text when dealing with matrices, linear solvers …..
; Vector functions:
; (v+ v1 v2) – addition of two vectors
; (v- v1 v2) – subtraction of two vectors
; (v*s v s) – multiply vector with scalar s
; (v/s v s) – divide vector with scalar s
; (unit_vect v) – calculates unit vector of vector v
; (angles_vect_deg v) – returns angles of vector v relative to axes X Y Z in deg
; (angles_vect v) – returns angles of vector v relative to axes X Y Z in radians
; (vect_dot v1 v2) – dot product of vectors v1 i v2
; (mod_vect v) – vector modulus (absolute length)
; (vect_L1 v) (vect_L2 v) (vect_L_inf) – vector norms
; Matrix functions:
; cm - creates named matrix
; (create_unit_matrix n) - create unit matrix (square) with n rows
; (matrix_size m) - returns matrix size
; (is_diagonal_matrix m) - checks if matrix is diagonal. Return list which first element is boolean – if true matrix is diagonal, and pivot of a matrix
; (is_scalar_matrix m) - checks if matrix is scalar. Return list which first element is boolean – if true matrix is scalar, and pivot of a matrix
; (is_unit_matrix m) - checks if matrix is unite matrix (eye)
; (is_null_matrix m) - checks if matrix is null (zero) matrix
; (is_symmetric m) - checks if matrix is symmetric
; (get_matrix_column m i) - fetches i-th column of a matrix m - returns vector
; (get_matrix_column_transp m i) - fetches i-th column of a matrix m - returns matrix
; (get_matrix_row m i) - fetches i-th row of a matrix m
; (swap_rows matrix m n) - swaps rows m and n in matrix
; (swap_cols matrix m n) - swaps columns m and n in matrix
; (mult_nth_row mat n val) - multiply nth row (n) of a matrix mat with a scalar value val
; (set_m_n mat m n val) – sets element in m-th row at nth position im matrix mat to value val
; (get_m_n mat m n) – gets value of a element in m-th row at nth position of a matrix mat
; (set_nth_row mat n row) – sets nth row of a matrix
; (set_nth_column mat n col) – sets nth column of a matrix
; (m+ m1 m2) – matrix addition
; (m- m1 m2) – matrix subtraction
; (m* m1 m2) – matrix multiplication
; (m*v mat vect) – matrix to vector multiplication
; (v*m vect mat) – vector to matrix multiplication
; (m+s m s) – add scalar value s to all elements of a matrix
; (m-s m s) – subtract scalar value s from all elements of a matrix
; (m*s m s) – multiply all elements of a matrix with scalar s
; (m/s m s) – divide all elements of a matrix with scalar s
; (remove_row mat row) – removes nth row of a matrix
; (remove_column mat col) – removes nth column of a matrix
; (submatrix mat row col) – return submatrix (minor) of a matrix excluding i-th row and column
; (trace mat) – returns matrix trace (sum of all elements at main diagonal)
; (main_diagonal mat) – return elements ath matrix main diagonal
; (determ mat) – calculates determinant value of a square matrix using Laplace method
; (pivot mat) - returns pivoted unit matrix, to shuffle matrix rows and prevent zero values at main diagonal by multiplication of a matrix with its pivoted unit matrix
; (lu_dec mat) - returns LU decomposition of a matrix mat with pivoting using Doolitle algorithm, returns vector containing L,U, and P matrix
; (lu mat) – performs LU decomposition of a matrix mat and store values in variables L,U and P
; (linear_solve A b) – solves linear system by using LU decomposition of a matrix A
; (inverse mat) – calculates inverse of a matrix trough LU decomposition, return list containing inverse of a matrix and its pivot
; (transp mat) – returns transposed matrix (code by Dough Wilson)
; (inv mat) – returns non-pivoted inverse of a matrix
; (adjoint mat) – calculates matrix adjoint
; (inverse_adj mat) – calculates inverse of a matrix using its adjoint
; (pinv mat) – calculates Moore – Penrose pseudoinverse (sometimes called general inverse) of matrix mat. If matrix is square then returned matrix is a true inverse, while for non-square matrices it is a matrix that fulfills nearly all prerequisites of a matrix inverse.
; (QR_decomposition mat) – QR decomposition of a matrix using Wedderburn's rank reduction method; returns a list containing matrices Q and R
; (QR a) – QR decomposition of a matrix, results are stored in global variables Q and R
; (qr_least_squares A b) – calculates least squares solution of a system ;minimizes ||Ax-b||^2 aka minimize a sum of squares of affine functions
; Examples:
; vectors:
; (setq v '(1 2 3))
; (setq u (unit_vector v)) → (0.267261 0.534522 0.801784)
; (angles_vect_deg u) → (74.4986 57.6885 36.6992)
; (mod_vect v) → 3.74166
; matrices:
; A - generated using command cm
; A -
; (
; (1.00000 2.00000 3.00000 4.00000)
; (2.00000 0.000000 4.00000 -4.00000)
; (0.000000 1.00000 2.00000 -1.00000)
; (2.00000 1.00000 0.000000 3.00000)
; )
; (determ A) - 48
; (inv A)
; (
; (-0.166667 0.166667 -0.0833333 0.416667)
; (-0.416667 -0.333333 1.29167 0.541667)
; (0.333333 0.166667 -0.333333 -0.333333)
; (0.250000 0.000000 -0.375000 -0.125000)
; )
; (m*m A invA)
; (
; (1.00000 0.000000 0.000000 0.000000)
; (0.000000 1.00000 0.000000 0.000000)
; (0.000000 0.000000 1.00000 0.000000)
; (0.000000 0.000000 0.000000 1.00000)
; )
; (m*m invA A)
; (
; (1.00000 0.000000 0.000000 0.000000)
; (0.000000 1.00000 0.000000 0.000000)
; (0.000000 0.000000 1.00000 0.000000)
; (0.000000 0.000000 0.000000 1.00000)
; )
; (to_fractions invA)
; (
; ((-1 6) (1 6) (-1 12) (5 12))
; ((-5 12) (-1 3) (31 24) (13 24))
; ((1 3) (1 6) (-1 3) (-1 3))
; ((1 4) (0 0) (-3 8) (-1 8))
; )
; (setq b '(30 -2 4 16))
; (m*v A b)
; ((102.000) (12.0000) (-10.0000) (106.000))
; (lu A)
; !L
; (
; (1.00000 0 0 0)
; (0.500000 1.00000 0 0)
; (0.000000 0.500000 1.00000 0)
; (1.00000 0.500000 -3.00000 1.00000)
; )
; !U
; (
; (2.00000 0.000000 4.00000 -4.00000)
; (0 2.00000 1.00000 6.00000)
; (0 0 1.50000 -4.00000)
; (0 0 0 -8.00000)
; )
; !P
; (
; (0 1 0 0)
; (1 0 0 0)
; (0 0 1 0)
; (0 0 0 1)
; )
; (m*m P (m*m L U))
; (
; (1.00000 2.00000 3.00000 4.00000)
; (2.00000 0.000000 4.00000 -4.00000)
; (0.000000 1.00000 2.00000 -1.00000)
; (2.00000 1.00000 0.000000 3.00000)
; )
; (linear_solve A b)
; (1.00000 2.00000 3.00000 4.00000)
; (qr A)
; !Q
; (
; (0.333333 0.757033 0.173422 0.534522)
; (0.666667 -0.432590 0.606977 0.000000)
; (0.000000 0.486664 0.346844 -0.801784)
; (0.666667 0.0540738 -0.693688 -0.267261)
; )
; !R
; (
; (3.00000 1.33333 3.66667 0.666667)
; (0.000000 2.05480 1.51407 4.43405)
; (0.000000 0.000000 3.64186 -4.16213)
; (0.000000 0.000000 0.000000 2.13809)
; )
; (qr_least_squares A b)
; ((1.00000) (2.00000) (3.00000) (4.00000))
; (pinv a)
; (
; (-0.166667 0.166667 -0.0833333 0.416667)
; (-0.416667 -0.333333 1.29167 0.541667)
; (0.333333 0.166667 -0.333333 -0.333333)
; (0.250000 0.000000 -0.375000 -0.125000)
; )
;usage of other helpers is more or less intuitive
; Code starts here
;----------------------------------------------------------
;Raises x to potential that is represented as a fraction
; (^ 3 '(2 5)) -> 1.55185
(defun ^ (x y) (expt x (decimal (frac y))))
(defun mappend (fn lst) ; P. Norvig ??
;"Append the results of calling fn on each element of list.
;Like mapcon, but uses append instead of nconc."
;One thing to notice is that fn must return a list, otherwise, it will go wrong.
; usage: (mappend '(lambda (x) (list x (* x x))) '(1 2 3))
(apply 'append (mapcar fn lst))
)
(defun mklist (x)
"If x is a list return it, otherwise return the list of x"
(if (listp x) x (list x)))
(defun flatten (exp)
;"Get rid of imbedded lists (to one level only)."
(mappend 'mklist exp))
(defun *error* (msg)
;;check for contents of msg
(cond
;;if msg is nil, do nothing
((not msg))
;;if msg is "quit / exit abort" (the content you want to trap from
;;the EXIT function, do nothing
;;NOTE: you may get this if the user cancels the routine as well
((equal (strcase msg t) "quit / exit abort"))
;;if msg is within this list, user may have canceled the function
((member (strcase msg t) '("console break" "function canceled"))
;;inform the user the routine was canceled
(princ "*Cancel*\n")
)
;;something else happened, print the error
((princ (strcat "\n; error: " msg "\n")))
)
;;'silent' exit
(princ)
);;end of *Error*
(defun log10 (n) (/ (log n) (log 10.0)))
; Random number generator
(defun rnd (/ modulus multiplier increment random)
(if (not seed)
(setq seed (getvar "DATE"))
)
(setq modulus 65536
multiplier 25173
increment 13849
seed (rem (+ (* multiplier seed) increment) modulus)
random (/ seed modulus)
)
)
;Generates random matrix
(defun randmatr (n m / row mat)
(repeat n
(repeat n
(setq a (* (rnd) 10))
(if (> a 5) (setq a 1.0)(setq a -1.0))
(setq row (cons (* a (* (rnd) m)) row))
)
(setq mat (cons row mat) row nil)
)
mat
)
;Number testing functions
(defun plusp (num) (cond ((numberp num) (>= num 0.0))))
(defun minusp (num) (cond ((numberp num) (< num 0.0))))
(defun oddp (num)(cond ((numberp num) (/= (rem num 2.0) 0.0))))
(defun evenp (num)(cond ((numberp num) (= (rem num 2.0) 0.0))))
(defun is_string (ent) (eq (type ent) 'str))
;Converts a number or a list containing numbers to float
;Used when integer value is needed in floating point representation
;in conversion to continuous fractions or similar
(defun to_float (e / ret k)
(if (listp e)
(progn
(while e
(setq
k (car e)
e (cdr e)
)
(if(listp k)
(setq k (mapcar 'float k)) ; don't join
(setq k (float k))
)
(setq ret (append ret (list k)))
)
(reverse ret)
)
(setq ret (float e))
)
ret
)
; (fix_numbers '(1 2 (3 5))) -> (1.00000 2.00000 (3.00000 5.00000))
(defun fix_numbers (e / ret k z r fix_num);
(defun fix_num (num / f r)
(setq f (fix num) r (- num f))
(if (zerop r) (setq num (float f)))
num
)
(if (listp e)
(progn
(while e
(setq
k (car e)
e (cdr e)
)
(if(listp k)
(setq k (mapcar 'fix_num k))
(setq k (fix_num k))
)
(setq ret (append ret (list k)))
)
;(reverse ret)
)
(setq ret(fix_num e))
)
ret
)
;Removes n-th element of a list
;Miljenko Hatlak
(defun rem_nth (lst n)
(cond
((and (plusp n) (< n (length lst)))
(cond
((zerop n)
(cdr lst)
)
(t (cons (car lst) (rem_nth (cdr lst) (1- n))))
)
)
(T lst)
)
)
;Sets n-th element of a list to new value
;(set_nth '(1 2 3 4) 2 5) -> (1 2 5 4)
(defun set_nth (lst n value)
(cond
((and (plusp n) (<= n (length lst)))
(cond
((zerop n)
(cons value (cdr lst))
)
(t (cons (car lst) (set_nth (cdr lst) (1- n) value)))
)
)
)
)
;Returns first n elemenst of a list
;(first_n '((1 2) 3 4 5) 2)-> ((1 2) 3)
(defun first_n (vect n / ret)
(cond
((and (plusp n) (<= n (length vect)))
(repeat n
(setq ret (cons (car vect) ret) vect (cdr vect))
)
)
)
(reverse ret)
)
;Returns last n elemenst of a list
;(last_n '((1 2) 3 4 5) 2)-> (4 5)
(defun last_n (vect n)
(reverse(first_n (reverse vect) n))
)
(defun rad_to_deg (rad)(* 180.0 (/ rad pi)))
(defun deg_to_rad (deg)(* pi (/ deg 180.0)))
(defun asin (x)
(cond
((and(> x -1.0)(< x 1.0)) (atan (/ x (sqrt (- 1.0 (* x x))))))
((= x -1.0) (* -1.0 (/ pi 2)))
((= x 1) (/ pi 2))
)
)
(defun acos (x)
(cond
((and(>= x -1.0)(<= x 1.0)) (-(* pi 0.5) (asin x)))
)
)
(defun rad_to_deg (rad)(* 180.0 (/ rad pi)))
(defun deg_to_rad (deg)(* pi (/ deg 180.0)))
;Creates matrix assignet to global variable
(defun c:cm( / rows cols matrix_name)
(setq
matrix_name (getstring "\nMatrix variable name >")
rows (getint "\nNumber of rows >")
cols (getint "\nNumber of cols >")
)
(create_named_matrix_ex matrix_name rows cols )
)
;Returns unnamed matrix element in iterative form
(defun create_matrix (rows cols / i j rows cols matrix row)
(setq i 0 j 0 matrix nil)
(repeat rows
(repeat cols
(setq row (append row (list(getreal (strcat "\n Element " (itoa (+ i 1)) " x " (itoa (+ j 1)) " > ")))))
(setq j (+ j 1))
)
(setq matrix (append matrix (list row)))
(setq j 0 i (+ i 1) row '())
)
matrix
)
;Returns named matrix element in iterative form
(defun create_named_matrix ( / rows cols matrix_name)
(setq matrix_name (getstring "\n Matrix variable name? > "))
(setq rows (getint "\n Number of rows? > "))
(setq cols (getint "\n Number of columns? > "))
(set (read matrix_name)(create_matrix rows cols))
)
;Returns named matrix element in iterative form
(defun create_named_matrix_ex (matrix_name rows cols)
(set (read matrix_name)(create_matrix rows cols))
)
(defun vect_shift_right (vect / a )
(setq a (cdr vect)) (append a (list(car vect)))
)
;Doug Wilson
;Transpose of a list i.e. matrix
(defun transp ( m )(apply 'mapcar (cons 'list m)))
;Creates unit matrix of a size n
(defun create_unit_matrix (n / unit)
(setq vect (cons 1 (create_null_vector (- n 1))))
(repeat n
(setq vect (vect_shift_right vect) unit (cons vect unit ))
)
)
;Returns list containing number of rows and columns of a matrix (rows cols)
(defun matrix_size (m)(list (length m) (length (car m))))
;Tests if list is a proper representation of a matrix
(defun is_matrix (m / elements_per_row elements_per_col equal_cols equal_rows not_nil ret)
(if (and m) ;nil is also list
(if (listp m)
(progn
(setq
elements_per_row (mapcar 'length m)
elements_per_col (mapcar 'length (transp m))
equal_rows (apply '= elements_per_row)
equal_cols (apply '= elements_per_col)
all_numbers T
)
(while m
(setq
all_numbers (and all_numbers (apply 'and (mapcar 'numberp (car m))))
m (cdr m)
)
)
(setq ret(and equal_cols equal_rows all_numbers))
)
(*error* "This isn't a matrix")
)
)
ret
)
;Tests if matrix is diagonal
(defun is_diagonal_matrix (mat / is i k dia );
(cond
((is_square_matrix mat)
(setq p (pivot mat) mat (m*m p mat))
(setq
dia (main_diagonal mat)
i 0)
(while (< i (length mat))
(setq
mat (set_m_n mat i i 0)
k (cons (apply '= (nth i mat)) k)
i (+ i 1))
)
(if (apply 'and k) (setq is T))
))
(list is p)
)
;Tests if matrix is scalar
(defun is_scalar_matrix (mat / is i k dia);
(cond
((is_square_matrix mat)
(setq p (pivot mat) mat (m*m p mat))
(setq
dia (main_diagonal mat)
i 0)
(while (< i (length mat))
(setq
mat (set_m_n mat i i 0)
k (cons (apply '= (nth i mat)) k)
i (+ i 1))
)
(if (and (apply '= dia) (apply 'and k)) (setq is T))
))
(list is p)
)
;Tests if matrix is unite matrix (eye)
(defun is_unit_matrix (mat / is_diag main is_unit ret)
(setq is_diag (is_diagonal_matrix mat))
(if (and (car is_diag))
(setq
mat (m*m (cadr is_diag) mat)
main (main_diagonal mat)
is_unit (and (= (car main) 1)(apply '= main))
ret (list is_unit (pivot mat)))
)
ret
)
;Tests if matrix is null (zero) matrix
(defun is_null_matrix (mat / i dia m k)
(setq i 0 dia (main_diagonal mat))
(while (< i (length mat))
(setq
m (cons m (= (get_m_n mat i i) 0))
k (cons (apply '= (nth i mat)) k)
i (+ i 1))
)
(and (apply 'and dia) (apply 'and k))
)
;Tests if matrix is square
(defun is_square_matrix(m)(cond((is_matrix m) (apply '= (matrix_size m)))))
;symmetric matrix generator (to be implemented)
;Tests if matrix is symetric
(defun is_symmetric (mat / tr ret)
(if (is_square_matrix mat)
(apply 'and (mapcar '= (flatten mat) (flatten(transp mat))))
)
)
;Returns i-th column of a matrix
(defun get_matrix_column (matrix i)(mapcar 'car(transp (list(nth i (transp matrix))))))
;Returns i-th row of a matrix
(defun get_matrix_row (matrix i) (nth i matrix))
;Returns i-th row of a matrix in matrix representation
(defun get_matrix_column_transp (matrix i)(mapcar 'list (get_matrix_column matrix i)))
;Swaps m-th and n-th rows of a matrix mat
;(swap_rows '((1 2)(3 4)) 0 1) ->((3 4)(1 2))
(defun swap_rows (mat m n / nel mel)
(setq
nel (nth n mat)
mel (nth m mat)
mat (set_nth mat m nel)
mat (set_nth mat n mel)
)
)
;Swaps m-th and n-th column of a matrix mat
;(swap_cols '((1 2)(3 4)) 0 1) ->((2 1)(4 3))
(defun swap_cols (mat m n)
(transp (swap_rows (transp mat) m n))
)
;Multiplies n-th row of a matrix with a scalar
;(mult_nth_row '((1 2)(3 4)) 0 3) ->((3 6) (3 4))
(defun mult_nth_row (mat n val)
(setq
nel (v*s (nth n mat) val)
mat (set_nth mat n nel)
)
)
;Adds two rows of a matrix
;(add_rows '((1 2)(3 4)) 1 0) ->((4 6) (3 4))
(defun add_rows (mat m n)
(setq
nel (nth n mat)
mel (nth m mat)
nel (mapcar '+ mel nel)
mel (set_nth mat n nel)
)
)
;Updates non-square matrix to full square by adding zero rows or columns
;(update_to_square '((1 2 3)(4 5 6))) ->((1 2 3) (4 5 6) (0 0 0))
(defun update_to_square (mat / size rows cols n )
(setq
size (matrix_size mat)
rows (car size)
cols (cadr size)
n (- rows cols)
)
(if (< n 0) ;update rows
(repeat (abs n)
(setq mat (append mat (list (create_null_vector cols))))
)
)
(if (> n 0) ;update cols
(progn
(setq mat (transp mat))
(repeat n
(setq mat (append mat (list (create_null_vector rows))))
)
(setq mat (transp mat))
)
)
mat
)
;Sets element of a matrix at m-th row and nth column to new value
;(set_m_n '((1 2)(3 4)) 0 1 10) -> ((1 10) (3 4))
(defun set_m_n (mat m n val / row el)
(setq
row (nth m mat)
row (set_nth row n val)
mat (set_nth mat m row)
)
)
;Get element of a matrix in mth row and nth column
(defun get_m_n (mat m n / row col)
(setq
row (nth m mat)
col (nth n row)
)
)
;Replaces row of a matrix with new one
;(set_nth_row '((1 2)(3 4)) 0 '(5 6)) -> ((5 6) (3 4))
(defun set_nth_row (mat index row)
(if (is_matrix mat)
(set_nth mat index row)
(list row)
)
)
;Replaces column of a matrix with new one
;(set_nth_column '((1 2)(3 4)) 0 '(5 6)) -> (((5 2) (6 4)))
(defun set_nth_column (mat index col)
(if (is_matrix mat)
(setq mat (transp mat)
mat(transp(set_nth mat index col))
)
(setq mat(transp (list col)))
)
)
;Multiplies two matrices that fulfill rules of matrix multiplication
(defun mmult (mat1 mat2 / rows_m1 rows_m2 row cols_m1 cols_m2 col i j aij mrow mult);
(if (and (is_matrix mat1)(is_matrix mat2))
(progn
(setq
mat2 (transp mat2)
rows_m1 (length mat1)
cols_m1 (length (car mat1))
rows_m2 (length mat2)
cols_m2 (length (car mat2))
i 0
j 0
)
(if (= cols_m1 cols_m2)
(repeat rows_m1
(setq row (nth i mat1))
(repeat rows_m2
(setq col (get_matrix_row mat2 j)
aij (apply '+ (mapcar '* row col))
)
(if (< (abs aij) 1e-14) (setq mrow (append mrow (list 0)))(setq mrow (append mrow (list aij))))
(setq j (+ j 1))
)
(setq mult (append mult (list mrow)))
(setq j 0 i(+ i 1) mrow nil)
)
)
)
)
mult
)
;Creates null vector of length n
;(create_null_vector 5) -> (0 0 0 0 0)
(defun create_null_vector (n / r ) (repeat n (setq r (cons 0 r))) r)
;Creates value vector of length n
;(create_value_vector 5 1) -> (1 1 1 1 1)
(defun create_value_vector (n val / r ) (repeat n (setq r (cons val r))) r)
;Creates vector of length n with start value start of length n
;(create_ordered_vector 5 0 2) -> (0 2 4 6 8)
(defun create_ordered_vector (n start step / r i) (setq i 0)(repeat n (setq r (cons (+ start (* i step)) r) i (+ i 1))(reverse r)))
;Creates null matrix
;(create_null_matrix 2 3) -> ((0.000000 0.000000 0.000000) (0.000000 0.000000 0.000000))
(defun create_null_matrix (rows cols / mat r)
(setq r (create_null_vector cols))
(repeat rows (setq mat (cons r mat)))
mat
)
;Mapper for matrix algebra
(defun matrix_matrix_operations (mat1 mat2 operation)
(if (and (is_matrix mat1)(is_matrix mat2))
(if (apply 'and (mapcar '= (matrix_size mat1)(matrix_size mat1)))
(fix_numbers (mapcar '(lambda ( r s ) (mapcar operation r s)) mat1 mat2))
)
)
)
;Matrix algebra function
(defun m+m (m1 m2) (matrix_matrix_operations m1 m2 '+))
(defun m-m (m1 m2) (matrix_matrix_operations m1 m2 '-))
(defun m*m (m1 m2) (fix_numbers (mmult m1 m2)))
;Multiplies matrix with a vector right side
;(m*v '((1 2)(3 4)) '(5 10)) -> ((25.0000) (55.0000))
(defun m*v (mat vect / r i )
(if
(and (eq (length(car mat)) (length vect)))
(progn
(setq i 0)
(repeat (length mat)
(setq r (append r (list(list (apply '+ (mapcar '* vect (nth i mat)))))) i (+ i 1))
)
(fix_numbers r)
)
)
)
;Multiplies vector with a matrix left side
;(v*m '(5 10) '((1 2)(3 4))) -> ((35.0000) (50.0000))
(defun v*m (vect mat / r i )
(if
(and (eq (length vect) (length mat) ))
(progn
(setq i 0)
(repeat (length vect)
(setq r (append r (list(list (apply '+ (mapcar '* vect (get_matrix_column mat i)))))) i (+ i 1))
)
(fix_numbers r)
)
)
)
;Mapper for matrix scalar operations
(defun matrix_scalar_operations (mat skal operation / r sv i row col);
(setq row (length mat) col (length (car mat)) sv (create_value_vector col skal) i 0)
(repeat row
(setq r (append r (list (mapcar operation (nth i mat) sv ))) i (+ i 1))
)
(fix_numbers r)
)
;Matrix to scalar operations
;(m*s '((1 2)(3 4)) 4) ->((4.00000 8.00000) (12.0000 16.0000))
(defun m+s (m s) (matrix_scalar_operations m s '+))
(defun m-s (m s) (matrix_scalar_operations m s '-))
(defun m*s (m s) (matrix_scalar_operations m s '*))
(defun m/s (m s) (matrix_scalar_operations m s '/))
;Vector algebra
(defun v+ (v1 v2)(mapcar '+ v1 v2))
(defun v- (v1 v2)(mapcar '- v1 v2))
(defun v*s (v s)(mapcar '* v (create_value_vector (length v) s)))
(defun v/s (v s)(mapcar '/ v (create_value_vector (length v) s)))
;Return unite vector of a vector v
(defun unit_vect (v)(mapcar '* v (create_value_vector (length v) (/ 1 (mod_vect v)))))
;Returns angles of a vector in degrees relative to horizontal plane
(defun angles_vect_deg (v)(mapcar 'rad_to_deg (mapcar 'acos (unit_vect v))))
;Returns angles of a vector in radians relative to horizontal plane
(defun angles_vect (v)(mapcar 'acos (unit_vect v)))
;Dot (inner) produvt of two vectors
(defun vect_dot (v1 v2)(apply '+ (mapcar '* v1 v2)))
;Vector modulus - absolute length
(defun mod_vect (v)(sqrt(vect_dot v v)))
;Various norms of a vector
(defun vect_L1(v) (apply '+ (mapcar 'abs v)))
(defun vect_L2(v) (mod_vect v)) ;Euclidean norm
(defun vect_L_inf (v) (apply 'max (mapcar 'abs v)))
;Removes nth column of a matrix
;(remove_column '((1 2)(3 4)) 0) ->((2) (4))
(defun remove_column (mat col)
(transp(rem_nth (transp mat) col))
)
;Removes nth row of a matrix
;(remove_row '((1 2)(3 4)) 0) ->((3 4))
(defun remove_row (mat row)
(rem_nth mat row)
)
;Returns submatrix - minor of a matrix
;(submatrix '((1 2 3)(4 5 6)(7 8 9)) 1 1) -> ((1 3) (7 9))
(defun submatrix (mat row col / ret)
(remove_row (remove_column mat col) row)
)
;Returns matrix trace - sum of all elements at main diagonal of a matrix
(defun trace (m / ret)
(if (is_square_matrix m)
(setq ret (apply '+ (main_diagonal m)))
)
ret
)
;Returns main diagonal of a matrix
(defun main_diagonal (mat / ret)
(while (and mat)
(setq
ret (cons (caar mat) ret)
mat (cdr (mapcar 'cdr mat))
)
)
(reverse ret)
)
;Returns second diagonal of a matrix
(defun sub_diagonal (mat) (reverse (main_diagonal (mapcar 'reverse mat))))
;Returns cofactor multiplier vector used in determinante calculations
;(cofact 5) -> (1 -1 1 -1 1)
(defun cofact (n / a cof)
(setq a (list 1 -1))
(repeat n
(setq cof (cons (car a) cof) a (reverse a))
)
(reverse cof)
)
;Calculates matrix determinante using Laplace method
;Author Miljenko Hatlak
(defun determ (mat / rows ret det first_row cof cofact subs i);
(defun det (mat)
(if(= (length mat) 1)
; don't join sets
(setq ret (car mat))
(setq ret
(cons(-(*(caar mat)(cadadr mat))(* (caadr mat)(cadar mat))) ret)
)
)
)
(defun cof (n / a cofs)
(setq a (list 1 -1))
(repeat n
(setq cofs (cons (car a) cofs) a (reverse a))
)
(reverse cofs)
)
(setq rows (length mat) cols (length (car mat)))
(if (is_square_matrix mat)
(if (<= rows 2)
(det mat)
(progn
(setq
first_row (car mat)
cofact (cof (length first_row))
first_row (mapcar '* first_row cofact)
i 0
)
(repeat (length mat)
(setq
subs (cons (submatrix mat 0 i) subs)
i (1+ i)
)
)
(setq subs (reverse subs)
ret (mapcar '* first_row (mapcar 'determ subs))
)
)
)
(*error* "Matrix is not a square!")
)
(apply '+ ret)
)
;Returns index of a first largest element in a vector
(defun vector_max_index (v / mx mi i k)
(setq mx (car v) v (cdr v) mi 0 i 0)
(while v
(setq i (+ i 1))
(setq k (car v))
(if (> k mx) (setq mx k mi i))
(setq v (cdr v))
)
mi
)
;Returns index of a first smallest element in a vector
(defun vector_min_index (v / mn mi i k)
(setq mn (car v) v (cdr v) mi 0 i 0)
(repeat (length v)
(setq i (+ i 1))
(setq k (car v))
(if (< k mn) (setq mn k mi i))
(setq v (cdr v))
)
mi
)
;returns pivoted unit matrix, to shuffle matrix rows and prevent zero values at diagonal places
;Author Miljenko Hatlak
(defun pivot (mat / pivot i mx n col);
(setq
n (length mat)
pivot (create_unit_matrix n)
i 0
eta 1
)
(cond
((and (= n 2) (= (caar mat) 0))
(setq
mat (reverse mat)
pivot (reverse pivot)
)
(if (zerop(caar mat)) (*error* "All elements in first matrix column are zero!! "))
)
)
(while (< i (- n 2)) ;ne diraj
(setq col(mapcar 'car mat)
mx (+ i (vector_max_index col))
)
(if (zerop (nth (- mx i) col))
(setq mx (+ i (vector_max_index (mapcar 'abs col))))
) ; prevent that max elemnt in col is 0
(if (> mx i)
(setq
pivot (swap_rows pivot mx i)
mat (swap_rows mat mx i)
)
)
(setq
mat (mapcar 'cdr (cdr mat))
i(+ i 1)
)
)
pivot
)
;Calculates LU decomposition of a square matrix using Doolitle's algorithm
;with rows pivoting of a initial matrix
;Assigns matrices L U and P to global variables L U and P
;Author: Miljenko Hatlak
(defun lu (mat)(mapcar 'set '(L U p) (lu_dec mat))) ; L U P are global variables
;Calculates LU decomposition of a square matrix using Doolitle's algorithm
; with rows pivoting of a initial matrix
;Returns a list containing matrices L and U with a pivot matrix P
;Author: Miljenko Hatlak
;(mapcar 'set '(L U P) (lu_dec a))
(defun lu_dec (mat / n L U p i j k sum)
(setq
n (length mat)
L (create_null_matrix n n)
U (create_null_matrix n n)
mat (to_float mat)
p (pivot mat)
mat (m*m p mat)
i 0 j 0 k 0
)
(while (< i n)
(setq k i)
(while (< k n)
(setq
sum 0
j 0
)
(while (< j i)
(setq
sum (+ sum (* (get_m_n L i j)(get_m_n U j k)))
j (+ j 1)
)
)
(setq
u (set_m_n U i k (- (get_m_n mat i k) sum))
k (+ k 1)
)
)
(setq k i)
(while (< k n)
(if (= i k)
(setq l (set_m_n L i i 1.0))
(progn
(setq sum 0 j 0)
(while (< j i)
(setq
sum (+ sum (* (get_m_n L k j)(get_m_n U j i)))
j (+ j 1)
)
)
(setq l (set_m_n L k i (/ (- (get_m_n mat k i) sum) (get_m_n U i i))))
)
)
(setq k (+ k 1))
)
(setq i (+ i 1))
)
(list L U p)
)
;Solves determinate system of linear equations
;(linear_solve '((1 2)(3 4)) '(5 6)) -> (-4.00000 4.50000)
(defun linear_solve (A b / n lu L U P b x y i );
(setq
n (length b)
lu (lu_dec A)
L (car lu)
U (cadr lu)
P (caddr lu)
b (to_float(car (transp (m*v p b))))
x (create_null_vector n)
y x
y (set_nth y 0 (car b))
i 1
)
(while (< i n)
(setq
y (set_nth y i (- (nth i b) (vect_dot (nth i L) y)))
i (+ i 1)
)
)
(setq
U (reverse (mapcar 'reverse U))
y (reverse y)
x (set_nth x 0 (/ (car y) (caar U)))
i 1
)
(while (< i n)
(setq
x (set_nth x i (/ (- (nth i y) (vect_dot (nth i U) x))(get_m_n U i i)))
i (+ i 1)
)
)
(fix_numbers (reverse x))
)
;return pivoted inverse matrix trough LU decomposition and partial (row) pivoting -> returns list (A-1,p)
(defun inverse (mat / n p unit i inv eow)
(setq
n (length mat)
p (pivot mat)
mat (m*m p mat)
unit (create_unit_matrix n)
i 0
)
(while unit
(setq
row (car unit)
unit (cdr unit)
inv (append inv (list(linear_solve mat row)))
)
)
(list (transp inv) p)
)
;returns only matrix inverse
(defun inv (mat / inv_m)(setq inv_m (inverse mat))(m*m (car inv_m) (cadr inv_m)))
;Calculates adjoint ( hermitian conjugate) of a matrix
;Autor Miljenko Hatlak
(defun adjoint (mat / size rows cols adj i j cof)
(setq
size (matrix_size mat)
rows (car size)
cols (cadr size)
)
(cond ((is_square_matrix mat)
(setq
adj (create_null_matrix rows cols)
i 0
j 0
)
(while (< i rows)
(while (< j rows)
(setq
cof (expt -1.0 (+ i j))
adj (set_m_n adj i j (* cof (determ (submatrix mat i j))))
j (+ j 1)
)
)
(setq i (+ i 1) j 0)
)
))
(if (and adj) (transp adj) (*error* "Matrix is not square !"))
)
;Calculates inverse of a matrix by using its adjoint
(defun inverse_adj (mat / inv det size )
(setq det (determ mat) size (matrix_size mat))
(cond ((= (car size)(cadr size))
(if (/= det 0.0) (setq inv(m*s (adjoint mat) (/ 1.0 det))))
))
inv
)
;Calculates Moore - Pentose pseudoinverse of a matrix
;If matrix is square then returned matrix is a true inverse, while for non-square matrices it is a matrix
;that fulfils nearly all prerequisites of a matrix inverse.
;https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse#Definition
(defun pinv (mat / tr ret)
(if (is_matrix mat)
(setq
tr (transp mat)
mat (m*m mat tr)
mat (inv mat)
mat (m*m tr mat)
)
(setq mat nil)
)
mat
)
;Wedderburn's rank reduction to find the QR factorization (decomposition) of a matrix A
;Returns a list containing matrices Q and R
(defun QR_decomposition (A / p size rows cols b q r qq rr i B Q_R)
(setq
size (matrix_size A)
rows (car size)
cols (cadr size)
A (update_to_square A)
B A
i 0
)
(repeat (min rows cols)
(setq
qq (get_matrix_column B i)
qq (fix_numbers (v*s qq (/ 1.0 (vect_L2 qq))))
Q (set_nth_column Q i qq)
rr (fix_numbers (first_n (car(transp(v*m qq B))) cols))
R (set_nth_row R i rr)
Q_R (fix_numbers (m*m Q R))
B (fix_numbers (m-m A Q_R))
i (+ i 1)
)
)
(if (< rows cols)(setq q (first_n q rows)))
(list Q R)
)
;Performs QR decomposition of a matrix and assign resulting matrices to global variables Q and R
(defun qr (mat)(mapcar 'set '(q r) (QR_decomposition mat)))
;Least squares solution to the system of equations using QR decomposition
(defun qr_least_squares (A b / size rows cols len QR_dec Q R x);
(setq
size (matrix_size A)
rows (car size)
cols (cadr size)
len (length b)
)
(if (= rows len)
(setq
QR_dec (QR_decomposition A)
Q (car QR_dec)
R (cadr QR_dec)
x (m*v (m*m (pinv R)(transp Q)) b)
)
)
x
)
;Fractions calculations
;Fetches numerator and denumerator of a fraction represented in a form of a lisr (numerator denumerator)
(defun numer (frac) (car frac))
(defun denom (frac) (cadr frac))
;Tests if element is a fraction or multiple fraction in any allowed form
;and performs fraction rationalizations to form (numerator denumerator)
;(frac '(1 2)) -> 0.5
;(frac'((1 2)(3 4))) -> (2 3)
;(frac'((1 2.5)(3.6 4))) -> (4 9)
(defun frac (el / a b c d e f g h ret)
(cond
((numberp el)
(setq ret (con_frac_to_frac(con_frac el)))
)
((and(listp el)(= (length el) 1))
(setq ret (con_frac_to_frac(con_frac (car el))))
)
((and(listp el)(=(length el) 2))
(setq
a (car el)
b (cadr el)
)
(if (numberp a) (setq c (con_frac_to_frac(con_frac a))))
(if (and(listp a)(= (length a) 1))(setq c(con_frac_to_frac(con_frac (car a)))))
(if (and(listp a)(= (length a) 2))
(setq
e(con_frac_to_frac(con_frac (car a)))
f(con_frac_to_frac(con_frac (cadr a)))
g(* (numer e) (denom f)) h (* (denom e) (numer f))
c (list g h)
)
)
(if (numberp b) (setq d (con_frac_to_frac(con_frac b))))
(if (and(listp b)(= (length b) 1))(setq d(con_frac_to_frac(con_frac (car b)))))
(if (and (listp b)(= (length b) 2))
(setq
e(con_frac_to_frac(con_frac (car b)))
f(con_frac_to_frac(con_frac (cadr b)))
g(* (numer e) (denom f)) h (* (denom e) (numer f))
d (list g h)
)
)
(setq e(* (numer c) (denom d)) f (* (denom c) (numer d)))
(setq ret (list e f ))
)
)
(frac_rationalize ret)
)
;Converts a rational number to its continuous fractions form
;(con_frac -3.456) ->((-1) (3 2 5 5 2))
(defun con_frac (num / ret sign recip i) ;
(setq
num (fix_numbers num)
sign (minusp num)
num (abs num)
ret (cons (fix num) ret)
num (- num (fix num))
i 0
)
(if (> num 0)
(setq recip (fix (/ 1 num))))
(cond ((and recip)
(while (and (< (car ret) 500)(< i 12))
(if (not (zerop num))
(setq num (/ 1.0 num) ret (cons (fix num) ret) num (- num (car ret))))
(setq i (+ i 1))
)
(if(> (car ret)500)
(setq ret (cdr ret)))
))
(if (and (= (car ret) 1)(> (length ret) 1))
(setq ret (cdr ret) ret(set_nth ret 0 (+ (car ret)1))))
(setq ret (reverse ret))
(if (and recip (= (length ret) 1))
(setq ret (append ret (list recip))))
(if (and (not recip)(= (length ret) 1))
(setq ret ret))
(if (and sign)
(setq sign -1)
(setq sign 1))
(list (list sign) ret)
)
;Converts continuous fraction to fraction
;(con_frac_to_frac '((-1) (3 2 5 5 2))) -> (-432 125)
;(decimal '(-432 125)) -> -3.45600
(defun con_frac_to_frac (con / sign z ret )
(if (= (length (cadr con)) 1)
(setq ret (list (* (caar con) (caadr con)) 1))
(progn
(setq
sign (car con)
con (reverse(cadr con))
z (list 1 (car con))
con (cdr con)
)
(repeat (- (length con) 1)
(setq
z (f/ 1 (f+ (car con) z))
con (cdr con)
)
)
(setq
ret (f+ (car con) z)
ret (set_nth ret 0 (* (car sign)(car ret)))
)
))
ret
)
;Tests fraction for illegal operations i.e dividing with zero .....
;and rationalizes it to form (numerator denumerator)
(defun frac_rationalize (fr / e f g _gcd);
(setq e (numer fr) f (denom fr))
(cond
((and (= e 0)(= f 0)) (setq fr (list 0 0))) ;undefined as a fraction but defined as a zero representation
((and (= e 0)(/= f 0)) (setq fr (list 0 0)))
((and (/= e 0)(= f 0)) (*error* "Cannot divide by zero!"))
((and (/= e 0)(/= f 0)) (setq fr fr))
)
(cond
((listp fr)
(setq e (numer fr) f (denom fr))
(if (and (zerop e) (zerop f))
(setq fr (list 0 0))
(progn
(setq
_gcd (gcd (abs e)(abs f))
_gcd (list _gcd _gcd)
g (mapcar '/ fr _gcd)
e (numer g) f (denom g)
)
(if (and (minusp e)(minusp f))
(setq e (* -1 e) f(* -1 f)))
(if (and (not(minusp e))(minusp f))
(setq e (* -1 e) f (* -1 f)))
(setq fr (list e f))
)
)
)
)
fr
)
;Fractions algebra
;(f+ '(1 2) '(3 4)) -> (5 4)
;(f+ 0.5 '(3 4)) -> (5 4)
;(f+ '(0.5 (2 5))'((3 4)(5.2 0.6))) -> (139 104)
(defun f+ (x y / a b ret);
(setq
x (frac x)
y (frac y)
a (apply 'and (mapcar 'zerop x))
b(apply 'and(mapcar 'zerop y))
)
(cond
((and a) (setq ret (frac_rationalize y)))
((and b) (setq ret (frac_rationalize x)))
((and (not a)(not b))
(setq ret
(frac_rationalize
(list
(+ (* (numer x) (denom y))(* (numer y) (denom x)))
(* (denom x) (denom y))
)
)
))
)
ret
)
(defun f- (x y / a b ret);
(setq
x (frac x)
y (frac y)
a (apply 'and(mapcar 'zerop x))
b (apply 'and(mapcar 'zerop y))
)
(cond
((and a) (setq ret (frac_rationalize (f* -1 y))))
((and b) (setq ret (frac_rationalize x)))
((and (not a)(not b))
(setq ret
(frac_rationalize
(list
(- (* (numer x) (denom y))(* (numer y) (denom x)))
(* (denom x) (denom y))
)
)
))
)
ret
)
(defun f* ( x y)
(setq x (frac x) y (frac y))
(frac_rationalize
(list
(* (numer x) (numer y))(* (denom x) (denom y))
)
)
)
(defun f/ (x y / a b ret);
(setq
x (frac x)
y (frac y)
a (apply 'and(mapcar 'zerop x))
b (apply 'and(mapcar 'zerop y))
)
(cond
((and a b) (*error* "Undefined (/ 0 0)"))
((and a)(setq ret (list 0 0)))
((and b)(*error* "Cannot divide by zero"))
((and (not a)(not b))
(setq ret
(frac_rationalize
(list
(* (numer x) (denom y))(* (denom x) (numer y))
)
)
))
)
ret
)
;Tests if two fractions are equals
;(f= '(1 2)'(2.5 5)) -> T
(defun f= (x y)
(setq x (frac x) y (frac y))
(= (* (numer x) (denom y))
(* (numer y) (denom x))
)
)
;Converts a list of numbers of complex list of numbers to its fractions form
(defun to_fractions (lst / x z w)
(foreach x lst
(if (listp x) (setq z (mapcar 'frac x)))
(if (numberp x) (setq z (frac x)))
(if (and z)(setq w (cons z w)))
;(setq z nil)
)
(reverse w)
)
;Conversion from fraction to decimal form
;(decimal '(139 104))-> 1.33654
(defun decimal (fr / n d)
(setq fr (frac fr)
n (numer fr)
d (denom fr)
)
(cond
((and (= n 0) (= d 0)) 0.0)
(T (/ (to_float n) d))
)
)
;Converts a list of fractions of complex list of fractions to its decimal form
(defun to_decimal (frac_lst / x z)
(while (and frac_lst)
(setq x (car frac_lst) frac_lst (cdr frac_lst))
(setq z (cons (decimal x) z))
)
(reverse z)
)
;Clean exit at load
(princ)
;------------- KRAJ - PODNAPUTBINE ZA REDONIZOVE, NIZOVE, RAZLOMKE ----------------------------------------------