1module dpmat; 2 3% Redistribution and use in source and binary forms, with or without 4% modification, are permitted provided that the following conditions are met: 5% 6% * Redistributions of source code must retain the relevant copyright 7% notice, this list of conditions and the following disclaimer. 8% * Redistributions in binary form must reproduce the above copyright 9% notice, this list of conditions and the following disclaimer in the 10% documentation and/or other materials provided with the distribution. 11% 12% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 13% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 14% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 15% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 16% CONTRIBUTORS 17% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 18% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 19% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 20% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 21% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 22% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 23% POSSIBILITY OF SUCH DAMAGE. 24% 25 26 27COMMENT 28 29 ##################### 30 ### ### 31 ### MATRICES ### 32 ### ### 33 ##################### 34 35This module introduces special dpoly matrices with its own matrix 36syntax. 37 38Informal syntax : 39 40 <matrix> ::= list('DPMAT,#rows,#cols,baslist,column_degrees,gb-tag) 41 42Dpmat's are the central data structure exploited in the modules of 43this package. Each such matrix describes a map 44 45 f : R^rows --> R^cols, 46 47that gives rise for the definition of two modules, 48 49 im f = the submodule of R^cols generated by the rows 50 of the matrix 51 52 and coker f = R^cols/im f. 53 54Conceptually dpmat's are identified with im f. 55 56END COMMENT; 57 58% ------------- Reference operators ---------------- 59 60symbolic procedure dpmat_rows m; cadr m; 61symbolic procedure dpmat_cols m; caddr m; 62symbolic procedure dpmat_list m; cadddr m; 63symbolic procedure dpmat_coldegs m; nth(m,5); 64symbolic procedure dpmat_gbtag m; nth(m,6); 65 66% ------------- Elementary operations -------------- 67 68symbolic procedure dpmat_rowdegrees m; 69% Returns the row degrees of the dpmat m as an assoc. list. 70 (for each x in dpmat_list m join 71 if (bas_nr x > 0) and bas_dpoly x then 72 {(bas_nr x).(mo_getdegree(dp_lmon bas_dpoly x,l))}) 73 where l=dpmat_coldegs m; 74 75symbolic procedure dpmat_make(r,c,bas,degs,gbtag); 76 list('dpmat,r,c,bas,degs,gbtag); 77 78symbolic procedure dpmat_element(r,c,mmat); 79% Returns mmat[r,c]. 80 dp_neworder 81 dp_comp(c, bas_dpoly bas_getelement(r,dpmat_list mmat)); 82 83symbolic procedure dpmat_print m; mathprint dpmat_2a m; 84 85symbolic procedure getleadterms!* m; 86% Returns the dpmat with the leading terms of m. 87 (begin scalar b; 88 b:=for each x in dpmat_list m collect 89 bas_make(bas_nr x,list(car bas_dpoly x)); 90 return dpmat_make(dpmat_rows m,dpmat_cols m,b,cali!=degrees,t); 91 end) where cali!=degrees:=dpmat_coldegs m; 92 93% -------- Symbolic mode file transfer -------------- 94 95symbolic procedure savemat!*(m,name); 96% Save the dpmat m under the name <name>. 97 begin scalar nat,c; 98 if not (stringp name or idp name) then typerr(name,"file name"); 99 if not eqcar(m,'dpmat) then typerr(m,"dpmat"); 100 nat:=!*nat; !*nat:=nil; 101 write"Saving as ",name; 102 out name$ 103 write"algebraic(setring "$ 104 105 % mathprint prints lists without terminator, but matrices with 106 % terminator. 107 108 mathprint ring_2a cali!=basering$ write")$"$ 109 write"algebraic(<<basis :="$ dpmat_print m$ 110 if dpmat_cols m=0 then write"$"$ write">>)$"$ 111 if (c:=dpmat_coldegs m) then 112 << write"algebraic(degrees:="$ 113 mathprint moid_2a for each x in c collect cdr x$ write")$"$ 114 >>; 115 write"end$"$ terpri()$ 116 shut name; terpri(); !*nat:=nat; 117 end; 118 119symbolic procedure initmat!* name; 120% Initialize a dpmat from <name>. 121 if not (stringp name or idp name) then typerr(name,"file name") 122 else begin scalar m,c; integer i; 123 write"Initializing ",name; terpri(); 124 in name$ m:=reval 'basis; cali!=degrees:=nil; 125 if eqcar(m,'list) then 126 << m:=bas_from_a m; m:=dpmat_make(length m,0,m,nil,nil)>> 127 else if eqcar(m,'mat) then 128 << c:=moid_from_a reval 'degrees; i:=0; 129 cali!=degrees:=for each x in c collect <<i:=i+1; i . x >>; 130 m:=dpmat_from_a m; 131 >> 132 else typerr(m,"basis or matrix"); 133 dpmat_print m; 134 return m; 135 end; 136 137% ---- Algebraic mode file transfer --------- 138 139symbolic operator savemat; 140symbolic procedure savemat(m,name); 141 if !*mode='algebraic then savemat!*(dpmat_from_a m,name) 142 else savemat!*(m,name); 143 144symbolic operator initmat; 145symbolic procedure initmat name; 146 if !*mode='algebraic then dpmat_2a initmat!* name 147 else initmat!* name; 148 149% --------------- Arithmetics for dpmat's ---------------------- 150 151symbolic procedure dpmat!=dpsubst(a,b); 152% Substitute in the dpoly a each e_i by b_i from the base list b. 153 begin scalar v; 154 for each x in b do 155 v:=dp_sum(v,dp_prod(dp_comp(bas_nr x,a),bas_dpoly x)); 156 return v; 157 end; 158 159symbolic procedure dpmat_mult(a,b); 160% Returns a * b. 161 if not eqn(dpmat_cols a,dpmat_rows b) then 162 rerror('dpmat,1," matrices don't match for MATMULT") 163 else dpmat_make( dpmat_rows a, dpmat_cols b, 164 for each x in dpmat_list a collect 165 bas_make(bas_nr x, 166 dpmat!=dpsubst(bas_dpoly x,dpmat_list b)), 167 cali!=degrees,nil) 168 where cali!=degrees:=dpmat_coldegs b; 169 170symbolic procedure dpmat_times_dpoly(f,m); 171% Returns f * m for the dpoly f and the dpmat m. 172 dpmat_make(dpmat_rows m,dpmat_cols m, 173 for each x in dpmat_list m collect 174 bas_make1(bas_nr x, dp_prod(f,bas_dpoly x), 175 dp_prod(f,bas_rep x)), 176 cali!=degrees,nil) where cali!=degrees:=dpmat_coldegs m; 177 178symbolic procedure dpmat_neg a; 179% Returns - a. 180 dpmat_make( 181 dpmat_rows a, 182 dpmat_cols a, 183 for each x in dpmat_list a collect 184 bas_make1(bas_nr x,dp_neg bas_dpoly x, dp_neg bas_rep x), 185 cali!=degrees,dpmat_gbtag a) 186 where cali!=degrees:=dpmat_coldegs a; 187 188symbolic procedure dpmat_diff(a,b); 189% Returns a - b. 190 dpmat_sum(a,dpmat_neg b); 191 192symbolic procedure dpmat_sum(a,b); 193% Returns a + b. 194 if not (eqn(dpmat_rows a,dpmat_rows b) 195 and eqn(dpmat_cols a, dpmat_cols b) 196 and equal(dpmat_coldegs a,dpmat_coldegs b)) then 197 rerror('dpmat,2,"matrices don't match for MATSUM") 198 else (begin scalar u,v,w; 199 u:=dpmat_list a; v:=dpmat_list b; 200 w:=for i:=1:dpmat_rows a collect 201 (bas_make1(i,dp_sum(bas_dpoly x,bas_dpoly y), 202 dp_sum(bas_rep x,bas_rep y)) 203 where y= bas_getelement(i,v), 204 x= bas_getelement(i,u)); 205 return dpmat_make(dpmat_rows a,dpmat_cols a,w,cali!=degrees, 206 nil); 207 end) where cali!=degrees:=dpmat_coldegs a; 208 209symbolic procedure dpmat_from_dpoly p; 210 if null p then dpmat_make(0,0,nil,nil,t) 211 else dpmat_make(1,0,list bas_make(1,p),nil,t); 212 213symbolic procedure dpmat_unit(n,degs); 214% Returns the unit dpmat of size n. 215 dpmat_make(n,n, for i:=1:n collect bas_make(i,dp_from_ei i),degs,t); 216 217symbolic procedure dpmat_unitideal!? m; 218 (dpmat_cols m = 0) and null matop_pseudomod(dp_fi 1,m); 219 220symbolic procedure dpmat_transpose m; 221% Returns transposed m with consistent column degrees. 222 if (dpmat_cols m = 0) then dpmat!=transpose ideal2mat!* m 223 else dpmat!=transpose m; 224 225symbolic procedure dpmat!=transpose m; 226 (begin scalar b,p,q; 227 cali!=degrees:= 228 for each x in dpmat_rowdegrees m collect 229 (car x).(mo_neg cdr x); 230 for i:=1:dpmat_cols m do 231 << p:=nil; 232 for j:=1:dpmat_rows m do 233 << q:=dpmat_element(j,i,m); 234 if q then p:=dp_sum(p,dp_times_ei(j,q)) 235 >>; 236 if p then b:=bas_make(i,p) . b; 237 >>; 238 return dpmat_make(dpmat_cols m,dpmat_rows m,reverse b, 239 cali!=degrees,nil); 240 end) where cali!=degrees:=cali!=degrees; 241 242symbolic procedure ideal2mat!* u; 243% Returns u as column vector if dpmat_cols u = 0. 244 if dpmat_cols u neq 0 then 245 rerror('dpmat,4,"IDEAL2MAT only for ideal bases") 246 else dpmat_make(dpmat_rows u,1, 247 for each x in dpmat_list u collect 248 bas_make(bas_nr x,dp_times_ei(1,bas_dpoly x)), 249 nil,dpmat_gbtag u) where cali!=degrees:=nil; 250 251symbolic procedure dpmat_renumber old; 252% Renumber dpmat_list old. 253% Returns (new . change) with new = change * old. 254 if null dpmat_list old then (old . dpmat_unit(dpmat_rows old,nil)) 255 else (begin scalar i,u,v,w; 256 cali!=degrees:=dpmat_rowdegrees old; 257 i:=0; u:=dpmat_list old; 258 while u do 259 <<i:=i+1; v:=bas_newnumber(i,car u) . v; 260 w:=bas_make(i,dp_from_ei bas_nr car u) . w ; u:=cdr u>>; 261 return dpmat_make(i,dpmat_cols old, 262 reverse v,dpmat_coldegs old,dpmat_gbtag old) . 263 dpmat_make(i,dpmat_rows old,reverse w,cali!=degrees,t); 264 end) where cali!=degrees:=cali!=degrees; 265 266symbolic procedure mathomogenize!*(m,var); 267% Returns m with homogenized rows using the var. name var. 268 dpmat_make(dpmat_rows m, dpmat_cols m, 269 bas_homogenize(dpmat_list m,var), cali!=degrees,nil) 270 where cali!=degrees:=dpmat_coldegs m; 271 272symbolic operator mathomogenize; 273symbolic procedure mathomogenize(m,v); 274% Returns the homogenized matrix of m with respect to the variable v. 275 if !*mode='algebraic then 276 dpmat_2a mathomogenize!*(dpmat_from_a reval m,v) 277 else matdehomogenize!*(m,v); 278 279symbolic procedure matdehomogenize!*(m,var); 280% Returns m with var. name var set equal to one. 281 dpmat_make(dpmat_rows m, dpmat_cols m, 282 bas_dehomogenize(dpmat_list m,var), cali!=degrees,nil) 283 where cali!=degrees:=dpmat_coldegs m; 284 285symbolic procedure dpmat_sieve(m,vars,gbtag); 286% Apply bas_sieve to dpmat_list m. The gbtag slot allows to set the 287% gbtag of the result. 288 dpmat_make(length x,dpmat_cols m,x,cali!=degrees,gbtag) 289 where x=bas_sieve(dpmat_list m,vars) 290 where cali!=degrees:=dpmat_coldegs m; 291 292symbolic procedure dpmat_neworder(m,gbtag); 293% Apply bas_neworder to dpmat_list m with current cali!=degrees. 294% The gbtag sets the gbtag part of the result. 295 dpmat_make(dpmat_rows m,dpmat_cols m, 296 bas_neworder dpmat_list m,cali!=degrees,gbtag); 297 298symbolic procedure dpmat_zero!? m; 299% Test whether m is a zero dpmat. 300 bas_zero!? dpmat_list m; 301 302symbolic procedure dpmat_project(m,k); 303% Project the dpmat m onto its first k components. 304 dpmat_make(dpmat_rows m,k, 305 for each x in dpmat_list m collect 306 bas_make(bas_nr x,dp_project(bas_dpoly x,k)), 307 dpmat_coldegs m,nil); 308 309% ---------- Interface to algebraic mode 310 311symbolic procedure dpmat_2a m; 312% Convert the dpmat m to a matrix (c>0) or a polynomial list (c=0) in 313% algebraic (pseudo)prefix form. 314 if dpmat_cols m=0 then bas_2a dpmat_list m 315 else 'mat . 316 if dpmat_rows m=0 then list for j:=1:dpmat_cols m collect 0 317 else for i:=1:dpmat_rows m collect 318 for j:=1:dpmat_cols m collect 319 dp_2a dpmat_element(i,j,m); 320 321symbolic procedure dpmat_from_a m; 322% Convert an algebraic polynomial list or matrix expression into a 323% dpmat with respect to the current setting of cali!=degrees. 324 if eqcar(m,'mat) then 325 begin integer i; scalar u,p; m:=cdr m; 326 for each x in m do 327 << i:=1; p:=nil; 328 for each y in x do 329 << p:=dp_sum(p,dp_times_ei(i,dp_from_a reval y)); i:=i+1 >>; 330 u:=bas_make(0,p).u 331 >>; 332 return dpmat_make(length m,length car m, 333 bas_renumber reversip u, cali!=degrees,nil); 334 end 335 else if eqcar(m,'list) then 336 ((begin scalar x; x:=bas_from_a reval m; 337 return dpmat_make(length x,0,x,nil,nil) 338 end) where cali!=degrees:=nil) 339 else typerr(m,"polynomial list or matrix"); 340 341% ---- Substitution in dpmats -------------- 342 343symbolic procedure dpmat_sub(a,m); 344% a=list of (var . alg. prefix form) to be substituted into the dpmat 345% m. 346 dpmat_from_a subeval1(a,dpmat_2a m) 347 where cali!=degrees:=dpmat_coldegs m; 348 349% ------------- Determinant ------------------------ 350 351symbolic procedure dpmat_det m; 352% Returns the determinant of the dpmat m. 353 if dpmat_rows m neq dpmat_cols m then rederr "non-square matrix" 354 else dp_from_a prepf numr detq matsm dpmat_2a m; 355 356endmodule; % dpmat 357 358end; 359