OKlibrary  0.2.1.6
Basics.mac
Go to the documentation of this file.
00001 /* Oliver Kullmann, 24.1.2008 (Swansea) */
00002 /* Copyright 2008, 2009, 2010, 2011 Oliver Kullmann
00003 This file is part of the OKlibrary. OKlibrary is free software; you can redistribute
00004 it and/or modify it under the terms of the GNU General Public License as published by
00005 the Free Software Foundation and included in this library; either version 3 of the
00006 License, or any later version. */
00007 
00022 oklib_include("OKlib/ComputerAlgebra/Hypergraphs/Lisp/SetSystems.mac")$
00023 oklib_include("OKlib/ComputerAlgebra/DataStructures/Lisp/HashMaps.mac")$
00024 oklib_include("OKlib/ComputerAlgebra/DataStructures/Lisp/Lists.mac")$
00025 
00026 
00027 /* *****************
00028    * Basic notions *
00029    *****************
00030 */
00031 
00032 /*
00033    A "combinatorial matrix" is a triple M = [R,C,f], where R is the set of
00034    row indices, C the set of column indices, and f is a map with domain
00035    R x C. It is a matrix "over H" if all values of f are elements of H.
00036    The abbreviation is "com".
00037 
00038    A "square combinatorial matrix" is a pair [I, f], where f is a map with
00039    domain I x I. Again it is a matrix "over H" if all values of f are elements
00040    of H. The abbreviation is "scom".
00041 
00042    The ordered forms, called "ocom" and "oscom", use lists instead of sets
00043    for the rows and colums.
00044 
00045    "Standardised" versions have row- and column-sets of the form {1,...,n}
00046    for a natural number n >= 0. Abbreviations are "stdcom", "stdscom" and
00047    "ostdcom" resp. "ostdscom".
00048 */
00049 
00050 /* Checking whether fam is a family with index set XxY over Z: */
00051 binfam_p(fam,X,Y,Z) := block(
00052  [e : errcatch(
00053    subsetp(map(lambda([P],apply(fam,P)),cartesian_product(X,Y)),Z))],
00054  not emptyp(e) and e[1])$
00055 
00056 
00057 /* Tests whether M is a combinatorial matrix over H: */
00058 com_p(M,H) := listp(M) and is(length(M) = 3) and setp(M[1]) and setp(M[2]) and
00059  binfam_p(M[3],M[1],M[2],H)$
00060 ocom_p(M,H) := listp(M) and is(length(M) = 3) and listp(M[1]) and listp(M[2])
00061  and is(length(setify(M[1]))=length(M[1]))
00062  and is(length(setify(M[2]))=length(M[2]))
00063  and binfam_p(M[3],setify(M[1]),setify(M[2]),H)$
00064 /* Tests whether M is a square combinatorial matrix over H: */
00065 scom_p(M,H) := listp(M) and is(length(M) = 2) and setp(M[1]) and 
00066   com_p([M[1],M[1],M[2]],H)$
00067 
00068 /* Two types of combinatorial matrices which "look like" square matrices: */
00069 /* A "square-sized combinatorial matrix": */
00070 sqscom_p(M,H) := com_p(M,H) and is(length(M[1]) = length(M[2]))$
00071 /* An "observed square combinatorial matrix": */
00072 osqcom_p(M,H) := com_p(M,H) and is(M[1] = M[2])$
00073 
00074 
00075 /* *********************
00076    * Checking equality *
00077    *********************
00078 */
00079 
00080 /* Equality test for combinatorial matrices: */
00081 /* RENAME: com_equalp */
00082 comequalp(A,B) := is(A[1] = B[1]) and is(A[2] = B[2]) and
00083  block([break : false],
00084   for a in A[1] unless break do for b in A[2] unless break do
00085     if A[3](a,b) # B[3](a,b) then break : true,
00086   return(not break))$
00087 com_equalp(A,B) := comequalp(A,B)$
00088 
00089 /* Equality test for square combinatorial matrices: */
00090 /* RENAME: scom_equalp */
00091 scomequalp(A,B) := is(A[1] = B[1]) and
00092  block([break : false],
00093   for a in A[1] unless break do for b in A[1] unless break do
00094     if A[2](a,b) # B[2](a,b) then break : true,
00095   return(not break))$
00096 scom_equalp(A,B) := scomequalp(A,B)$
00097 
00098 
00099 /* ***************
00100    * Conversions *
00101    ***************
00102 */
00103 
00104 /* Conversions between general and square matrices: */
00105 scom2com(M) := [M[1],M[1],M[2]]$
00106 oscom2ocom(M) := [M[1],M[1],M[2]]$
00107 com2scom(M) := [M[1],M[3]]$
00108 ocom2oscom(M) := [M[1],M[3]]$
00109 /* com2scom, ocom2oscom can only be applied sensibly if M[1] = M[2]. */
00110 
00111 /* Conversions between ordered and unordered matrices: */
00112 ocom2com(M) := [setify(M[1]),setify(M[2]),M[3]]$
00113 oscom2scom(M) := [setify(M[1]),M[2]]$
00114 com2ocom(M) := [listify(M[1]),listify(M[2]),M[3]]$
00115 scom2oscom(M) := [listify(M[1]),M[2]]$
00116 
00117 /* Conversions between combinatorial matrices and Maxima matrices: */
00118 
00119 /* Transforms a combinatorial matrix into a Maxima matrix, using the
00120    implicit order on rows and columns (note that the values are kept): */
00121 com2m(M) := block([rows : listify(M[1]), cols : listify(M[2]), f : M[3]],
00122   apply(matrix, create_list(create_list(f(i,j),j,cols),i,rows)))$
00123 /* Now using the explicit order on rows and columns: */
00124 ocom2m(M) := block([rows : M[1], cols : M[2], f : M[3]],
00125   apply(matrix, create_list(create_list(f(i,j),j,cols),i,rows)))$
00126 /* For square matrices: */
00127 scom2m(M) := com2m(scom2com(M))$
00128 oscom2m(M) := ocom2m(oscom2ocom(M))$
00129 /* For special situations, see lmcom2m, lmscom2m below (w.r.t efficiency). */
00130 
00131 /* Transforms a Maxima matrix into a (standardised) combinatorial matrix
00132    (attention: this constitutes a reference to the original matrix). */
00133 m2com(M) := block([m,n],
00134   [m,n] : matrix_size(M),
00135   [setn(m), setn(n),
00136     buildq([M], lambda([i,j], 'M[i,j]))])$
00137 m2ocom(M) := block([m,n],
00138   [m,n] : matrix_size(M),
00139   [create_list(i,i,1,m), create_list(i,i,1,n),
00140     buildq([M], lambda([i,j], 'M[i,j]))])$
00141 /* RENAME: m2scom */
00142 sm2scom(M) := [setn(matrix_size(M)[1]), buildq([M], lambda([i,j], 'M[i,j]))]$
00143 m2scom(M) := sm2scom(M)$
00144 m2oscom(M) := [create_list(i,i,1,matrix_size(M)[1]), buildq([M], lambda([i,j], 'M[i,j]))]$
00145 /* For combinatorial matrices created in this way, re-extract the
00146    original Maxima matrix:
00147 */
00148 lmcom2m(M) := ev(part(M[3],2,0),eval);
00149 lmscom2m(M) := ev(part(M[2],2,0),eval);
00150 
00151 /* For convenience: Input as for the Maxima function "matrix", but returns
00152    a com resp. scom:
00153 */
00154 com_l([L]) := m2com(apply(matrix,L))$
00155 scom_l([L]) := m2scom(apply(matrix,L))$
00156 
00157 /* Transforms a Maxima matrix M with a given list R of row-names and a
00158    given list C of column-names into an ordered cominatorial matrix.
00159    (Attention: Again, this constitutes a reference to the original matrix.)
00160 */
00161 mrc2ocom(M,R,C) := block(
00162  [hr : osm2hm(l2osm_inv(R)), hc : osm2hm(l2osm_inv(C))],
00163   [R, C,
00164    buildq([M,hr,hc], lambda([a,b], 'M[ev_hm(hr,a),ev_hm(hc,b)]))])$
00165 
00166 /* Access to the matrix-references again by lmcom2m. */
00167 
00168 
00169 /* ***********************
00170    * Basic constructions *
00171    ***********************
00172 */
00173 
00174 /* The empty combinatorial matrix; since no arguments are given, this means
00175    that rowset and columnset are empty; more generally, an "empty com" has no
00176    entries, that is, one of rowset or columnset (at least) is empty: */
00177 emptycom : [{},{},lambda([i,j],false)]$
00178 emptycom_r(R) := [R,{},lambda([i,j],false)]$
00179 emptycom_c(C) := [{},C,lambda([i,j],false)]$
00180 
00181 /* The empty square combinatorial matrix: */
00182 emptyscom : [{},lambda([i,j],false)]$
00183 
00184 /* The zero combinatorial matrix for given index sets: */
00185 zerocom(I,J) := [I,J,lambda([i,j],0)]$
00186 /* The zero square combinatorial matrix for a given index set: */
00187 zeroscom(I) := [I,lambda([i,j],0)]$
00188 
00189 /* The constant combinatorial matrix for given index sets and value: */
00190 /* RENAME: constcom */
00191 constantcom(I,J,a) := [I,J,buildq([a],lambda([i,j],a))]$
00192 constcom(I,J,a) := constantcom(I,J,a)$
00193 conststdcom(n,m,a) := constcom(setn(n),setn(m),a)$
00194 /* The constant square combinatorial matrix for given index set and value: */
00195 /* RENAME: constscom */
00196 constantscom(I,a) := [I,buildq([a],lambda([i,j],a))]$
00197 constscom(I,a) := constantscom(I,a)$
00198 conststdscom(n,a) := constscom(setn(n),a)$
00199 
00200 /* The identity matrix: */
00201 idscom(I) := [I, lambda([i,j],if i=j then 1 else 0)]$
00202 /* The diagonal matrix with constant diagonal value: */
00203 cdiagscom(I,a) := [I,buildq([a],lambda([i,j],if i=j then a else 0))]$
00204 
00205 /* Creating all Maxima matrices of size m*n with elements from
00206    list L (in lexicographical order): */
00207 /* Prerequisite: if m=0 then n=0. */
00208 all_m(m,n,L) := if n=0 then [apply(matrix,create_list([],i,1,m))] else
00209 map(
00210  lambda([L],apply(matrix,partition_elements(L,n))),
00211  cartesian_product_l(create_list(L,i,1,m*n)))$
00212 
00213 /* The corrected Maxima function "genmatrix" for the main case 
00214    (can now handle also genmatrix(f,m,0)): */
00215 genmatrix_m(_fun,m,n) := if n=0 then apply(matrix,create_list([],i,1,m))
00216  else genmatrix(_fun,m,n)$
00217 genmatrix_sm(_fun,n) := if n=0 then matrix()
00218  else genmatrix(_fun,n)$
00219 
00220 
00221 /* *******************
00222    * Random matrices *
00223    *******************
00224 */
00225 
00226 /* A random Maxima matrix of size m x n, with entries random(x)
00227    (as usual, if m=0 then n=0): */
00228 random_m(m,n,x) := genmatrix_m(lambda([i,j],random(x)),m,n)$
00229 /* A random Maxima matrix of size m x n, with random entries
00230    from the integer interval [a,b]: */
00231 random_intiv_m(m,n,a,b) := block([s : b - a + 1],
00232   genmatrix_m(lambda([i,j],random(s)+a),m,n))$
00233 /* A random Maxima matrix of size m x n, with random entries
00234    from the real interval [a,b) (using floating points): */
00235 random_fpiv_m(m,n,a,b) := block([s : float(b - a)],
00236   genmatrix_m(lambda([i,j],random(s)+a),m,n))$
00237 
00238 
00239 /* A random symmetric Maxima matrix: */
00240 random_sm(n,x) := block(
00241  [m : genmatrix_sm(lambda([i,j], 
00242         if i > j then random(x) elseif i=j then random(x)/2 else 0),n)],
00243   m + transpose(m))$
00244 
00245 /* A random permutation of a Maxima matrix, via row and column permutation
00246    (thus the rows and columns themselfes stay intact, they only change
00247    their places and get internally permuted): */
00248 random_permutation_m(M) := block([s : matrix_size(M)],
00249   trans_l2m(random_permutation(create_list(i,i,1,s[1]))) .
00250   M .
00251   trans_l2m(random_permutation(create_list(i,i,1,s[2]))))$
00252 /* A random permutation of the rows of a Maxima matrix: */
00253 random_rpermutation_m(M) := block([s : matrix_size(M)],
00254   trans_l2m(random_permutation(create_list(i,i,1,s[1]))) . M)$
00255 /* A random permutation of the columns of a Maxima matrix: */
00256 random_cpermutation_m(M) := block([s : matrix_size(M)],
00257   M . trans_l2m(random_permutation(create_list(i,i,1,s[2]))))$
00258 
00259 /* A random self-polar matrix : */
00260 random_spolar_m(n,x) := random_permutation_m(random_sm(n,x))$
00261 
00262 /* Takes list of field elements, a set of lists of elements over the field,
00263    a function implementing the field addition, a set of field elements,
00264    a function implementing the field multiplication and returns the linear
00265    hull of adjoin(L,H).
00266 
00267    H is assumed to be a linear hull itself, i.e., all lists in H  must be of
00268    the same size and for any list in H, H contains all linearly dependent
00269    lists of that list.
00270 */
00271 extend_linear_hull(V,H,F,add_f, mul_f) := block([lin_dep_V],
00272   lin_dep_V : create_list(map(lambda([ve], mul_f(ve,e)),V), e, listify(F)),
00273   return(disjoin(create_list(0,i,1,length(V)),
00274       setify(append(lin_dep_V,
00275           create_list(map(add_f,Vl,Vh), Vl, lin_dep_V, Vh, listify(H)))))))$
00276 /* Computes the linear hull of a given set of vectors over the field
00277    defined by the field elements F along with additiona and multiplication
00278    functions. */
00279 linear_hull(L,F,add_f, mul_f) :=
00280   lreduce(lambda([H,V], extend_linear_hull(V,H,F,add_f,mul_f)),
00281     listify(L), {})$
00282 
00283 /* Generates a random invertible matrix with r rows and c columns over
00284    the field with element set F and addition and multiplication functions
00285    add_f and mul_f. */
00286 random_inv_m(r,c,F,add_f,mul_f) := block([avail_rows, rows : [], hull : {}],
00287   avail_rows : listify(disjoin(create_list(0,i,1,c),
00288     apply(cartesian_product, create_list(F,i,1,c)))),
00289   for rp : 1 thru r do block([row],
00290     row : random_element(avail_rows),
00291     rows : cons(row,rows),
00292     hull : extend_linear_hull(row,hull,F,add_f,mul_f),
00293     avail_rows : remove_elements(hull, avail_rows)),
00294   return(apply(matrix, rows)))$
00295   
00296 
00297 /* ****************************
00298    * Rows, columns, diagonals *
00299    ****************************
00300 */
00301 
00302 /* Transforms a combinatorial matrix into a list, row-wise: */
00303 com2l_r(M) := map(lambda([p],apply(M[3],p)), cartesian_product_l(map(listify,rest(M,-1))))$
00304 /* Similarly, but into a list of lists: */
00305 com2ll_r(M) := create_list(create_list(M[3](i,j), j,listify(M[2])), i,listify(M[1]))$
00306 scom2ll_r(M) := com2ll_r(scom2com(M))$
00307 /* Transforms a combinatorial matrix into a list, column-wise: */
00308 com2l_c(M) := com2l_r(trans_com(M))$
00309 /* Similarly, but into a list of lists: */
00310 com2ll_c(M) := create_list(create_list(M[3](i,j), i,listify(M[1])), j,listify(M[2]))$
00311 scom2ll_c(M) := com2ll_c(scom2com(M))$
00312 
00313 /* Remarks: For a Maxima matrix A, by A[i] we obtain the list of the elements
00314    of row i; accordingly, by transpose(A)[i] we obtain column i as list.
00315    On the other hand, by row(A,i) resp. col(A,i) we obtain the rows resp.
00316    columns as (Maxima) matrices.
00317 
00318    For convenience these functions for Maxima-matrices:
00319 */
00320 /* Extracting row i from Maxima-matrix M as list: */
00321 row_m(M,i) := M[i]$
00322 /* Extracting column i from Maxima-matrix M as list: */
00323 column_m(M,i) := transpose(col(M,i))[1]$
00324 
00325 
00326 /* Transforms a combinatorial matrix into an ordered multiset of rows as
00327    ordered multi-sets:
00328 */
00329 com2omsoms_r(M) := list_distribution(map(list_distribution,com2ll_r(M)))$
00330 /* The same, for the columns: */
00331 com2omsoms_c(M) := list_distribution(map(list_distribution,com2ll_c(M)))$
00332 
00333 
00334 /* The main diagonal of a square combinatorial matrix: */
00335 maindiag_scom(M) := create_list(M[2](i,i), i,listify(M[1]))$
00336 /* As an ordered multiset: */
00337 maindiagoms_scom(M) := list_distribution(maindiag_scom(M))$
00338 
00339 
00340 /* *******************
00341    * Transformations *
00342    *******************
00343 */
00344 
00345 /* Transforms a symmetric combinatorial matrix into a graph: */
00346 /* RENAME: scom2g */
00347 com2g(M) := [M[1], subset(powerset(M[1],2), lambda([S],
00348   block([L : listify(S)], is(M[2](L[1],L[2]) # 0))))]$
00349 scom2g(M) := com2g(M)$
00350 /* Transforms a symmetric combinatorial matrix into a graph with loops: */
00351 scom2gl(M) := [M[1], subset(union(powerset(M[1],1),powerset(M[1],2)), 
00352  lambda([S], block([L : listify(S)], 
00353    if length(L)=1 then is(M[2](L[1],L[1]) # 0) else
00354      is(M[2](L[1],L[2]) # 0))))]$
00355 
00356 /* Transforms a square combinatorial matrix into a digraph: */
00357 scom2dg(M) := [M[1], subset(cartesian_product(M[1],M[1]), lambda([L],
00358   is(L[1]#L[2] and M[2](L[1],L[2]) # 0)))]$
00359 /* Transforms a symmetric combinatorial matrix into a digraph with loops: */
00360 scom2dgl(M) := [M[1], subset(cartesian_product(M[1],M[1]), 
00361  lambda([L], is(M[2](L[1],L[2]) # 0)))]$
00362 
00363 
00364 /* Transforms a combinatorial matrix into a general hypergraph: */
00365 /* RENAME: com2ghg */
00366 com2ghyp(M) := [M[1], M[2], buildq([M], 
00367  lambda([j], subset(M[1], lambda([i], is(M[3](i,j) # 0)))))]$
00368 com2ghg(M) := com2ghyp(M)$
00369 
00370 /* Converts a transformation in list-form into a square combinatorial
00371    matrix: */
00372 trans_l2scom(L) := block([A : l2array(L)],[setn(length(L)),
00373   buildq([A],lambda([i,j],if A[j]=i then 1 else 0))]);
00374 trans_l2m(L) := scom2m(trans_l2scom(L))$
00375 
00376 
00377 /* ********************
00378    * Basic operations *
00379    ********************
00380 */
00381 
00382 /* Recall that for a Maxima matrix via matrix_size the pair
00383    [number_rows,number_columns]
00384    is computed.
00385 */
00386 
00387 /* The dimensions of a combinatorial matrix: */
00388 dim_com(M) := [length(M[1]), length(M[2])]$
00389 /* The order of a square combinatorial matrix: */
00390 order_scom(M) := length(M[1])$
00391 
00392 /* Transposition of a combinatorial matrix: */
00393 /* RENAME: trans_com */
00394 trans(M) := [M[2], M[1], buildq([M], lambda([i,j], M[3](j,i)))]$
00395 trans_com(M) := [M[2], M[1], buildq([M], lambda([i,j], M[3](j,i)))]$
00396 /* Transposition of a square combinatorial matrix: */
00397 /* RENAME: trans_scom */
00398 strans(M) := [M[1], buildq([M], lambda([i,j], M[2](j,i)))]$
00399 trans_scom(M) := [M[1], buildq([M], lambda([i,j], M[2](j,i)))]$
00400 
00401 
00402 /* ************************
00403    * Algebraic operations *
00404    ************************
00405 */
00406 
00407 /* The addition of two square combinatorial matrices
00408    Precondition: A[1] = B[1]. */
00409 add_scom(A,B) := [A[1], buildq([A,B],
00410   lambda([i,j], A[2](i,j) + B[2](i,j)))]$
00411 /* The difference of two square combinatorial matrices
00412    Precondition: A[1] = B[1]. */
00413 diff_scom(A,B) := [A[1], buildq([A,B],
00414   lambda([i,j], A[2](i,j) - B[2](i,j)))]$
00415 
00416 /* Multiplying a square combinatorial matrix with a scalar: */
00417 scprod_scom(a,A) := [A[1], buildq([a,A], lambda([i,j], a * A[2](i,j)))]$
00418 
00419 /* The product of two combinatorial matrices.
00420    Precondition: A[2] = B[1]. */
00421 prod_com(A,B) := [A[1], B[2], buildq([A,B],
00422   lambda([i,j], 
00423     sum_l(map(lambda([k],A[3](i,k)),listify(A[2])) * map(lambda([k],B[3](k,j)),listify(A[2])))))]$
00424 
00425 
00426 /* ********************
00427    * Entry statistics *
00428    ********************
00429 */
00430 
00431 /* The minimum over all entries of a non-empty combinatorial matrix: */
00432 min_com(M) := lmin(com2l_r(M))$
00433 /* The maximum over all entries of a non-empty combinatorial matrix: */
00434 max_com(M) := lmax(com2l_r(M))$
00435 
00436 
00437 /* The sum over all entries of a combinatorial matrix */
00438 sum_com(M) := sum_l(com2l_r(M))$
00439 /* The sum over all entries of a square combinatorial matrix */
00440 sum_scom(M) := sum_com(scom2com(M))$
00441 
00442 /* The sum of matrix entries at positions given by L for Maxima matrix A: */
00443 sum_m(A,L) := sum_l(map(lambda([p], A[p[1],p[2]]), L))$
00444 
00445 
00446 /* The entrywise absolute values of a combinatorial matrix */
00447 abs_com(M) := [M[1], M[2], buildq([M],lambda([i,j], abs(M[3](i,j))))]$
00448 /* The entrywise absolute values of a square combinatorial matrix */
00449 abs_scom(M) := [M[1], buildq([M],lambda([i,j], abs(M[2](i,j))))]$
00450 
00451 
00452 /* The distribution of row-sums, that is a list of pairs [possible
00453     row-sum, number of rows], sorted by ascending row-sum. */
00454 rowsums_list_com(M) := block([C : listify(M[2])],
00455  list_distribution(
00456   create_list(gsum_l(lambda([j],M[3](i,j)),C), i,listify(M[1]))))$
00457 
00458 /* The distribution of column-sums, that is a list of pairs [possible
00459     column-sum, number of rows], sorted by ascending column-sum. */
00460 columnsums_list_com(M) := block([R : listify(M[1])],
00461  list_distribution(
00462   create_list(gsum_l(lambda([i],M[3](i,j)),listify(R)), j,listify(M[2]))))$
00463 
00464 /* The distribution of values, as a list of pairs [possible value, number
00465    of entries], sorted by ascending values: */
00466 com_distribution(M) := list_distribution(com2l_r(M))$
00467 
00468 
00469 /* ***************
00470    * Basic tests *
00471    ***************
00472 */
00473 
00474 /* Checking for an "empty" combinatorial matrix (no entries): */
00475 empty_com_p(M) := is(length(M[1]) = 0 or length(M[2]) = 0)$
00476 
00477 every_com_p(_pred,M) := every_s(_pred, com2l_r(M))$
00478 every_scom_p(_pred,M) := every_com_p(_pred,scom2com(M))$
00479 some_com_p(_pred,M) := some_s(_pred, com2l_r(M))$
00480 some_scom_p(_pred,M) := some_com_p(_pred, scom2com(M))$
00481 
00482 /* Every diagonal value is constant = s: */
00483 every_diagc_scom_p(M,s) := 
00484   every_s(lambda([i], is(M[2](i,i)=s)), M[1])$
00485 /* More generally, every diagonal entry fulfils _predx: */
00486 every_diag_scom_p(_predx,M) := 
00487   every_s(lambda([i], _predx(M[2](i,i))), M[1])$
00488 /* Every value outside the diagonal is constant = s: */
00489 every_ndiagc_scom_p(M,s) := block([L : listify(M[1])],
00490   every_s(lambda([p], is(p[1]=p[2] or apply(M[2],p)=s)), cartesian_product_l([L,L])))$
00491 /* More generally, every non-diagonal entry fulfils _predx: */
00492 every_ndiag_scom_p(_predx,M) := block([L : listify(M[1])],
00493   every_s(lambda([p], is(p[1]=p[2] or _predx(apply(M[2],p)))), cartesian_product_l([L,L])))$
00494 /* Some non-diagonal entry fulfils _predx: */
00495 some_ndiag_scom_p(_predx,M) := block([L : listify(M[1])],
00496   some_s(lambda([p], is(p[1]#p[2] and _predx(apply(M[2],p)))), cartesian_product_l([L,L])))$
00497 
00498 /* Externally square combinatorial matrices: */
00499 extscom_com_p(M) := is(length(M[1]) = length(M[2]))$
00500 /* Observedly square combinatorial matrices: */
00501 obsscom_com_p(M) := is(M[1] = M[2])$ 
00502 
00503 symmetric_scom_p(M) := scom_equalp(M,trans_scom(M))$
00504 symmetric_m_p(M) := is(M = transpose(M))$
00505 
00506 /* Whether Maxima matrix A has entries all different at the positions given
00507    by L: */
00508 alldifferent_m_p(A,L) := listnorep_p(map(lambda([p], A[p[1],p[2]]), L))$
00509 
00510 
00511 /* ***************
00512    * Polynomials *
00513    ***************
00514 */
00515 
00516 oklib_plain_include("nchrpl")$
00517 
00518 /* The characteristic polynomial, in the standard sense, where the leading
00519    coefficient is always +1.
00520 */
00521 /* It seems that for all our cases setting ratmx to true yields a *much*
00522    faster computation: */
00523 charpoly_m(M) := block([n : matrix_size(M)[1]],
00524   if n=0 then 1 
00525   elseif n < 30 then
00526     block([ratmx : true], (-1)^(matrix_size(M)[1]) * charpoly(M,x))
00527   else ncharpoly(M,x))$
00528 charpoly_scom(M) := charpoly_m(scom2m(M))$
00529