OKlibrary  0.2.1.6
QuadraticForms.mac
Go to the documentation of this file.
00001 /* Oliver Kullmann, 25.1.2008 (Swansea) */
00002 /* Copyright 2008 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/CombinatorialMatrices/Lisp/Basics.mac")$
00023 
00024 
00025 /* ***********************************************
00026    * Indices of positivity, negativity, nullity *
00027    ***********************************************
00028 */
00029 
00030 /* Auxiliary function: Find an element in a symmetric matrix which is largest 
00031    with respect to its absolute value, where diagonal elements are preferred.
00032    Return the index-pair. Assumes the matrix is non-empty. */
00033 pivot_d(M) := block([A : abs(M), m : matrix_size(M)[1], l, il : 1, jl : 1, a],
00034   l : A[il,jl],
00035   for i : 2 thru m do (a : A[i,i], if a > l then (l : a, il : i, jl : i)),
00036   for i : 1 thru m-1 do for j : i+1 thru m do
00037     (a : A[i,j], if a > l then (l : a, il : i, jl : j)),
00038   return([il,jl]))$
00039 
00040 /* For a symmetric real matrix M compute [index of positivity, index of
00041   negativity, index of nullity]. */
00042 /* M is a Maxima matrix. */
00043 pnn_indices(M) := block([m : matrix_size(M)[1], pi, pj, s, res],
00044   if m = 0 then return([0,0,0]),
00045   [pi,pj] : pivot_d(M),
00046   if [pi,pj] # [1,1] then (
00047     M : rowswap(M,1,pi), M : columnswap(M,1,pi),
00048     if pi # pj then (
00049       M : columnop(M,1,pj,-1), M : rowop(M,1,pj,-1)
00050     )),
00051   s : sign(M[1,1]),
00052   if s # zero then (
00053     for i : 2 thru m do block([a : M[1,i] / M[1,1]],
00054       M : rowop(M,i,1,a), M : columnop(M,i,1,a))),
00055   M : minor(M,1,1),
00056   res : pnn_indices(M),
00057   if s = zero then return(res + [0,0,1])
00058   elseif s = pos then return(res + [1,0,0])
00059   else return(res + [0,1,0]))$
00060 
00061 /* The hermitian rank of a symmetric real matrix */
00062 hermitian_rank(M) := block([pnn : pnn_indices(M)],
00063   max(pnn[1],pnn[2]))$
00064 
00065 /* The hermitian defect of a symmetric real matrix */
00066 hermitian_def(M) := matrix_size(M)[1] - hermitian_rank(M)$
00067 
00068 /* Computing the hermitian rank via (numerical!) eigenvalues. */
00069 hermitian_rank_eig(M) := if emptyp(M) then 0 else
00070 block(
00071  [eigvals : eigens_by_jacobi(M,bigfloatfield)[1], iplus : 0, iminus : 0],
00072   for e in eigvals do block([s, eps : 10^(1-fpprec)],
00073     s : if abs(e) <= eps then zero else sign(e),
00074     if s = pos then iplus : iplus + 1
00075     elseif s = neg then iminus : iminus + 1),
00076   return(max(iplus,iminus)))$
00077 
00078 /* Computing the hermitian rank via the characteristic polynomial. */
00079 /* This method seems to be fastest currently. */
00080 hermitian_rank_charpoly(M) := if emptyp(M) then 0 else
00081 block(
00082  [char_poly, highpow, ipluss : 0, iminuss : 0, s1, x],
00083   char_poly : charpoly_m(M),
00084   highpow : hipow(char_poly,x),
00085   for i : highpow step -1 thru 0 do block(
00086    [s : sign(coeff(char_poly,x,i))],
00087     if s # zero then s1 : s),
00088   for j : 0 thru highpow do block([s2 : sign(coeff(char_poly,x,j))],
00089     if s1 # s2 and s2 # zero then (s1 : s2, ipluss : ipluss + 1)),
00090   iminuss : matrix_size(M)[1] - ipluss - lopow(char_poly,x),
00091   return(max(ipluss,iminuss)))$
00092 
00093