Community
In my last project I had to create a simple library to perform matrix and fractions calculus in autolisp (work in progress).
It provides matrix and vectors arithmetic, LU and QR decompositions, linear solver, least squares solver, inversion and Moore-Penrose pseudoinvers, adjoint, fractions arithmetic,
conversions between decimal and fraction representation of real numbers, continuous fractions, determinant and a lots of helper functions.
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).
Introduction to continued fractions
Fractions calculations:
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 😎 (18 125) (72 125))
Conversions from fraction to decimal representation of a rational number.
(decimal '(1 2)) →0.5
(to_decimal '((1 2) (1 😎 (18 125) (72 125))) ->(0.500000 0.125000 0.144000 0.576000)
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. In some future reincarnation I'll try to figure out basic mathematical operations
for continuous fractions – it may bring me a Nobel – LOL.
(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 😎
(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 of 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 fulfils 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 😎 (-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))
Miljenko Hatlak
Did your question get successfully answered? Then click on the ACCEPT SOLUTION button.
![Autodesk Design & Make Report](https://forums.autodesk.com/t5/image/serverpage/image-id/1350519iCB55845EBFC43016/image-size/large?v=v2&px=999" "Autodesk Design & Make Report")
