/* r2x2.h --- 2x2 matrices and operations on them */ /* Last edited on 2024-12-05 10:27:51 by stolfi */ #ifndef r2x2_H #define r2x2_H #include #include #include typedef struct r2x2_t { double c[2][2]; } r2x2_t; /* STANDARD OPERATIONS */ void r2x2_zero(r2x2_t *M); /* Stores in {M} the null matrix. */ void r2x2_ident(r2x2_t *M); /* Stores in {M} the identity matrix. */ void r2x2_transp(r2x2_t *A, r2x2_t *M); /* Sets {M} to the transpose {A^t} of matrix {A}. */ void r2x2_get_row(r2x2_t *A, uint32_t i, r2_t *x); void r2x2_set_row(r2x2_t *A, uint32_t i, r2_t *x); /* These two procedures copy row {i} of matrix {A} to and from vector {x}, respectively. */ void r2x2_get_col(r2x2_t *A, uint32_t j, r2_t *x); void r2x2_set_col(r2x2_t *A, uint32_t j, r2_t *x); /* These two procedures copy column {j} of matrix {A} to and from vector {x}, respectively. */ void r2x2_map_row(r2_t *x, r2x2_t *A, r2_t *r); /* Sets {r} to the product of row vector {x} by matrix {A}. */ void r2x2_map_col(r2x2_t *A, r2_t *x, r2_t *r); /* Sets {r} to the product of matrix {A} by column vector {x}. */ void r2x2_add(r2x2_t *A, r2x2_t *B, r2x2_t *M); /* Sets {M} to the matrix sum {A + B}. */ void r2x2_sub(r2x2_t *A, r2x2_t *B, r2x2_t *M); /* Sets {M} to the matrix difference {A - B}. */ void r2x2_neg(r2x2_t *A, r2x2_t *M); /* Sets {M} to the negated matrix {-A}. */ void r2x2_scale(double s, r2x2_t *A, r2x2_t *M); /* Sets {M} to the product of scalar {s} and matrix {A}. */ void r2x2_mix(double s, r2x2_t *A, double t, r2x2_t *B, r2x2_t *M); /* Sets {M} to the linear combination {s*A + t*B}. */ void r2x2_mul(r2x2_t *A, r2x2_t *B, r2x2_t *M); /* Sets {M} to the product of matrices {A} and {B}. */ void r2x2_mul_tr(r2x2_t *A, r2x2_t *B, r2x2_t *M); /* Computes the matrix product {M = A * B^t}. (In other words, sets {M[i,j]} to the dot product of row {i} of {A} and row {j} of {B}.) */ void r2x2_tr_mul(r2x2_t *A, r2x2_t *B, r2x2_t *M); /* Computes the matrix product {M = A^t * B}. (In other words, sets {M[i,j]} to the dot product of column {i} of {A} and column {j} of {B}.) */ double r2x2_det(r2x2_t *A); /* Returns the determinant of matrix {A}. */ void r2x2_adj(r2x2_t *A, r2x2_t *M); /* Sets {M} to the adjoint of matrix {A}. */ void r2x2_inv(r2x2_t *A, r2x2_t *M); /* Sets {M} to the inverse of matrix {A}. */ void r2x2_sqrt(r2x2_t *A, r2x2_t *M); /* Sets {M} to the square root of {A}. Bombs if the square root is non-real or does not exist. */ double r2x2_norm_sqr(r2x2_t *A); /* Squared Frobenius norm of {A}, i.e. sum of squares of elements. */ double r2x2_norm(r2x2_t *A); /* Frobenius norm of {A}, i.e. square root of sum of squares of elements. */ double r2x2_normalize(r2x2_t *A); /* Divides the matrix {A} by its {r2x2_norm}. If that norm is zero, the matrix is filled with {NAN}. Returns the norm. */ double r2x2_mod_norm_sqr(r2x2_t *A); /* Returns the square of the Frobenius norm of {A-I}. */ double r2x2_L_inf_norm(r2x2_t *A); /* Infinity norm of {A}, i.e. max absolute element value. */ double r2x2_L_inf_normalize(r2x2_t *A); /* Divides the matrix {A} by its {r2x2_L_inf_norm}. If that norm is zero, the matrix is filled with {NAN}. Returns the norm. */ void r2x2_diff_sqr(r2x2_t *A, r2x2_t *B, r2x2_t *R, double *dabs2P, double *drel2P); /* Returns in {*dabs2P} and {*drel2P} (if not NULL) the total squared discrepancy between the matrices {*A} and {*B}, respectively absolute and relative to the matrix {*R}. Namely { dabs2 = SUM { (Ae[s] - Be[s])^2 } } { drel2 = SUM { ((Ae[s] - Be[s])/Re[s])^2 } } where {Ae[s],Be[s],Re[s]} are all corresponding elements of {*A,*B,*C}. excluding those where {Re[s]} is zero. */ bool_t r2x2_is_unif_scaling(r2x2_t *M, double s); /* TRUE iff {M} is a diagonal matrix with all diagonal elements equal to {s}. */ void r2x2_from_rows(r2_t *a, r2_t *b, r2x2_t *M); /* Sets {M} to the matrix whose rows are the vectors {a,b}. */ void r2x2_from_cols(r2_t *a, r2_t *b, r2x2_t *M); /* Sets {M} to the matrix whose columns are the vectors {a,b}. */ /* I/O */ void r2x2_print(FILE *f, r2x2_t *A); /* Prints matrix {A} to file {f}, with default format. */ void r2x2_gen_print ( FILE *f, r2x2_t *A, char *fmt, char *olp, char *osep, char *orp, /* Outer delimiters. */ char *ilp, char *isep, char *irp /* Inner delimiters. */ ); /* Prints matrix {A} to file {f}, using {fmt} for each element. The matrix is bounded by {olp} and {orp}, and rows are separated by {osep}. Each row is bounded by {ilp} and {irp}, and elements are separated by {isep}. Defaults are provided for any of these strings which are NULL. */ /* OPERATIONS SPECIFIC TO 2x2 MATRICES */ void r2x2_rot90(r2x2_t *M); /* Stores in {M} the matrix {[[00,+1],[-1,00]]}, so that {r2x2_map_row(x,M,x)} rotates {x} by 90 degrees counterclockwise. */ void r2x2_rot_and_scale(r2_t *p, r2x2_t *M); /* Stores in {M} the similarity (rotation and uniform scaling) matrix that keeps the origin fixed and maps {(1,0)} to {p}. */ void r2x2_sym_eigen(r2x2_t *A, r2_t *e, r2x2_t *M); /* Computes the eigenvalues {e.c[0],e.c[1]) of the symmetric matrix {A}. The eigenvalues will be sorted in decreasing order of SIGNED value. Actually, the procedure computes the eigenvalues of the matrix {(A+A')/2}, where {A'} is the transpose of {A}. That's always symmetric, and equal to {A} if {A} is symmetric. If {M} is not NULL, the procedure also stores into the rows of {M} the corresponding unit-length eigenvectors. In that case, {M} will be a rotation matrix (i.e. {M'*M = M*M'= IDENT()}, {det(M)=+1}) such that {A = M'*DIAG(e)*M} (or {M*A*M' = DIAG(e)}), where {DIAG(e)} is the diagonal matrix with the elements of {e} along the diagonal. */ void r2x2_from_point_pairs(r2_vec_t *p1, r2_t *bar1, r2_vec_t *p2, r2_t *bar2, r2x2_t *M); /* Returns in {*M} the linear map that takes the points {p1} as close as possible to the points {p2}, in the sense of minimizing the mean squared error. If {bar1} is not NULL, the point {*bar1} is subtracted from every point in {p1}. Ditto for {bar2} and {p2}. The lists {p1 and {p2} must have the same length, and {p1} (after subtracting {bar1}, if not null) must contain at least two linearly independent points; otherwise the result may contain infinite or {NAN} elements. */ #endif