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