1module edsaux;
2
3% Miscellaneous support functions
4
5% Author: David Hartley
6
7% Redistribution and use in source and binary forms, with or without
8% modification, are permitted provided that the following conditions are met:
9%
10%    * Redistributions of source code must retain the relevant copyright
11%      notice, this list of conditions and the following disclaimer.
12%    * Redistributions in binary form must reproduce the above copyright
13%      notice, this list of conditions and the following disclaimer in the
14%      documentation and/or other materials provided with the distribution.
15%
16% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
18% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
19% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR
20% CONTRIBUTORS
21% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27% POSSIBILITY OF SUCH DAMAGE.
28%
29
30
31fluid '(!*edsverbose !*edsdebug);
32
33
34% Operations on eds
35
36% Storing these kernel lists seems to have <1% effect on speed.
37
38symbolic procedure indkrns s;
39   % s:eds -> indkrns:list of lpow pf
40   % independent kernels - leading basis forms in eds_ind s
41%   geteds(s,'indkrns) or
42%   puteds(s,'indkrns,
43      foreach f in eds_ind s join
44            if tc(f := trterm f) = (1 ./ 1) then {tpow f}; %);
45
46
47symbolic procedure depkrns s;
48   % s:eds -> depkrns:list of lpow pf
49   % dependent kernels - a basis for edscob s/eds_ind s
50   % OK?
51%   geteds(s,'depkrns) or
52%   puteds(s,'depkrns,
53      foreach f in  eds_sys s join
54            if lc f = (1 ./ 1) and degreepf f = 1 then {lpow f}; %);
55
56
57symbolic procedure prlkrns s;
58   % s:eds -> prlkrns:list of lpow pf
59   % principal kernels - a basis for the non-linear part of s
60%   geteds(s,'prlkrns) or
61%   puteds(s,'prlkrns,
62      setdiff(edscob s,append(indkrns s,depkrns s)); %);
63
64
65symbolic procedure remkrns s;
66   % s:eds -> remkrns:nil
67   foreach p in {'indkrns,'depkrns,'prlkrns} do rempropeds(s,p);
68
69
70% Operations on sys
71
72
73symbolic procedure degreepart(s,d);
74   % s:sys, d:int -> degreepart:sys
75   foreach f in s join if degreepf f = d then {f};
76
77
78symbolic procedure scalarpart s;
79   % s:sys -> scalarpart:sys
80   degreepart(s,0);
81
82
83symbolic procedure pfaffpart s;
84   % s:sys -> pfaffpart:sys
85   degreepart(s,1);
86
87
88symbolic procedure nonpfaffpart s;
89   % s:sys -> nonpfaffpart:sys
90   foreach f in s join if degreepf f neq 1 then {f};
91
92
93symbolic procedure solvedpart s;
94   % s:sys -> solvedpart:sys
95   foreach f in s join if f and lc f = (1 ./ 1) then {f};
96
97
98symbolic procedure lpows s;
99   % s:list of pf -> lpows:list of lpow pf
100   % returns the leading kernels in s
101   foreach f in s collect lpow f;
102
103
104symbolic procedure singleterms s;
105   % s:list of pf -> singleterms:bool
106   % true if each form in s is a single term
107   null s or null red car s and singleterms cdr s;
108
109
110symbolic procedure !*a2sys u;
111   % u:prefix list -> !*a2sys:list of pf
112   if rlistp u then
113      foreach f in cdr indexexpandeval{u} collect xpartitop f
114   else typerr(u,'list);
115
116
117symbolic procedure !*sys2a u;
118   % u:list of pf -> !*sys2a:prefix list
119   makelist foreach f in u collect !*pf2a f;
120
121
122symbolic procedure !*sys2a1(u,v);
123   % u:list of pf -> !*sys2a:prefix list
124   % !*sys2a with choice of !*sq or true prefix
125   makelist foreach f in u collect !*pf2a1(f,v);
126
127
128% Operations on pf
129
130
131symbolic procedure xpows f;
132   % f:pf -> xpows:list of lpow pf
133   if f then lpow f . xpows red f;
134
135
136symbolic procedure xcoeffs f;
137   % f:pf -> xcoeffs:list of sq
138   if f then lc f . xcoeffs red f;
139
140
141symbolic procedure degreepf f;
142   % f:pf -> degreepf:int
143   % assumes f homogeneous
144   % could replace with xdegree from XIDEAL2
145   if null f then 0
146   else (if null x then 0
147   else if fixp x then x
148   else rerror('eds,130,"Non-integral degree not allowed in EDS"))
149   where x = deg!*form lpow f;
150
151
152symbolic procedure xreorder f;
153   % f:pf -> xreorder:pf
154   if f then
155      addpf(multpfsq(xpartitop lpow f,reordsq lc f),
156            xreorder red f);
157
158
159symbolic procedure xreorder!* f;
160   % f:pf -> xreorder!*:pf
161   % Like xreorder when it is known that only the order between terms
162   % has changed (and not within terms).
163   if f and red f then
164      addpf(lt f .+ nil,xreorder!* red f)
165   else f;
166
167
168symbolic procedure xrepartit f;
169   % f:pf -> xrepartit:pf
170   if f then
171      addpf(wedgepf(xpartitsq subs2 resimp lc f,xpartitop lpow f),
172            xrepartit red f);
173
174
175symbolic procedure xrepartit!* f;
176   % f:pf -> xrepartit!*:pf
177   % Like xrepartit when xvars!* hasn't changed.
178   if f then
179      addpf(multpfsq(xpartitop lpow f,subs2 resimp lc f),
180            xrepartit!* red f);
181
182
183symbolic procedure trterm f;
184   % f:pf -> trterm:lt pf
185   % the trailing term in f
186   if null red f then lt f else trterm red f;
187
188
189symbolic procedure linearpart(f,p);
190   % f:pf, p:list of 1-form kernel -> linearpart:pf
191   % result is the part of f of degree 1 in p
192   if null f then nil
193   else if length intersection(wedgefax lpow f,p) = 1 then
194      lt f .+ linearpart(red f,p)
195   else linearpart(red f,p);
196
197
198symbolic procedure inhomogeneouspart(f,p);
199   % f:pf, p:list of 1-form kernel -> inhomogeneouspart:pf
200   % result is the part of f of degree 0 in p
201   if null f then nil
202   else if length intersection(wedgefax lpow f,p) = 0 then
203      lt f .+ inhomogeneouspart(red f,p)
204   else inhomogeneouspart(red f,p);
205
206
207symbolic procedure xcoeff(f,c);
208   % f:pf, c:pffax -> xcoeff:pf
209   if null f then nil
210   else
211      begin scalar q,s;
212      q := xdiv(c,wedgefax lpow f);
213      if null q then return xcoeff(red f,c);
214      q := mknwedge q;
215      if append(q,c) = lpow f then s := 1 % an easy case
216      else s := numr lc wedgepf(!*k2pf q,!*k2pf mknwedge c);
217      return q .* (if s = 1 then lc f else negsq lc f)
218               .+ xcoeff(red f,c);
219      end;
220
221
222symbolic procedure xvarspf f;
223   % f:pf -> xvarspf:list of kernel
224   if null f then nil
225   else
226      union(foreach k in
227               append(wedgefax lpow f,
228               append(kernels numr lc f,
229                      kernels denr lc f)) join if xvarp k then {k},
230            xvarspf red f);
231
232
233symbolic procedure kernelspf f;
234   % f:pf -> kernelspf:list of kernel
235   % breaks up wedge products
236   if null f then nil
237   else
238      union(append(wedgefax lpow f,
239            append(kernels numr lc f,
240                          kernels denr lc f)),
241            kernelspf red f);
242
243
244symbolic procedure mkform!*(u,p);
245   % u:prefix, p:prefix (usually int) -> mkform!*:prefix
246   % putform with u returned, and covariant flag removed
247   begin
248   putform(u,p);
249   return u;
250   end;
251
252
253% Operations on lists
254
255
256symbolic procedure purge u;
257   % u:list -> purge:list
258   % remove repeated elements from u, leaving last occurence only
259   if null u then nil
260   else if car u member cdr u then purge cdr u
261   else car u . purge cdr u;
262
263
264symbolic procedure rightunion(u,v);
265   % u,v:list -> rightunion:list
266   % Like union, but appends v to right. Ordering of u and v preserved
267   % (last occurence of each element in v used).
268   append(u,foreach x on v join
269      if not(car x member u) and not(car x member cdr x) then {car x});
270
271
272symbolic procedure sublisq(u,v);
273   % u:a-list, v:any -> sublisq:any
274   % like sublis, but leaves structure untouched where possible
275   if null u or null v then v
276   else
277      begin scalar x,y;
278      if (x := atsoc(v,u)) then return cdr x;
279      if atom v then return v;
280      y := sublisq(u,car v) . sublisq(u,cdr v);
281      return if y = v then v else y;
282      end;
283
284
285symbolic procedure ordcomb(u,n);
286   % u:list, n:int -> ordcomb:list of list
287   % List of all combinations of n distinct elements from u.
288   % Order of u is preserved: ordcomb(u,1) = mapcar(u,'list)
289   % which is not true for comb, from which this is copied.
290   begin scalar v; integer m;
291   if n=0 then return {{}}
292   else if (m:=length u-n)<0 then return {}
293   else for i := 1:m do
294    <<v := nconc!*(v,mapcons(ordcomb(cdr u,n-1),car u));
295      u := cdr u>>;
296   return nconc!*(v,{u})
297   end;
298
299
300symbolic procedure updkordl u;
301   % u:list of kernel -> updkordl:list of kernel
302   % list version of updkorder
303   % kernels in u will have highest precedence, order of
304   % other kernels unchanged
305   begin scalar v,w;
306   v := kord!*;
307   w := append(u,setdiff(v,u));
308   if v=w then return v;
309   kord!* := w;
310   alglist!* := nil . nil;        % Since kernel order has changed.
311   return v
312   end;
313
314
315symbolic procedure allequal u;
316   % u:list of any -> allequal:bool
317   % t if all elements of u are equal
318   if length u < 2 then t
319   else car u = cadr u and allequal cdr u;
320
321
322symbolic procedure allexact u;
323   % u:list of kernel -> allexact:bool
324   % t if all elements of u are exact pforms
325   if null u then t
326   else exact car u and allexact cdr u;
327
328
329symbolic procedure coords u;
330   % u:list of kernel -> coords:list of kernel
331   % kernels in u are 1-forms, returns coordinates involved
332   % THIS SHOULD GO
333   foreach k in u join if exact k then {cadr k};
334
335
336symbolic procedure zippf(pfl,sql);
337   % pfl:list of pf, sql:list of sq
338   % -> zippf:pf
339   % Multiply elements of pfl by corresponding elements of sql, and add
340   % results together. If trailing nulls on pfl or sql can be omitted.
341   if null pfl or null sql then nil
342   else if numr car sql and car pfl then
343      addpf(multpfsq(car pfl,car sql),zippf(cdr pfl,cdr sql))
344   else
345      zippf(cdr pfl,cdr sql);
346
347
348% EDS tracing and debugging
349
350
351symbolic procedure errdhh u;
352% Special error call for errors that shouldn't happen
353rerror(eds,999,"Internal EDS error -- please contact David Hartley
354*****" . if atom u then {u} else u);
355
356
357symbolic procedure edsverbose(msg,v,type);
358   % msg:atom or list, v:various, type:id|nil
359   % -> edsverbose:nil
360   % type gives the type of v, one of
361   %    nil  - v is empty
362   %    'sf  - v is a list of sf
363   %    'sq  - v is a list of sq
364   %    'sys - v is a list of pf
365   %    'cob - v is a list of kernel
366   %    'map - v is a map (list of kernel.prefix)
367   %    'xform - v is an xform (list of kernel.pf)
368   %    'prefix- v is prefix
369   if !*edsverbose or !*edsdebug then
370   begin
371   if atom msg then msg := {msg};
372   lpri msg; terpri();
373   if null type then nil
374   else if type = 'sf then
375      foreach f in v do edsmathprint prepf f
376   else if type = 'sq then
377      foreach f in v do edsmathprint prepsq f
378   else if type = 'sys then
379      foreach f in v do edsmathprint preppf f
380   else if type = 'cob then
381      edsmathprint makelist v
382   else if type = 'map then
383      foreach p in v do edsmathprint {'equal,car p,cdr p}
384   else if type = 'rmap then
385   << foreach p in car v do edsmathprint {'equal,car p,cdr p};
386      foreach p in cadr v do edsmathprint {'neq,p,0} >>
387   else if type = 'xform then
388      foreach p in v do edsmathprint {'equal,car p,preppf cdr p}
389   else if type = 'prefix then
390      edsmathprint v
391   else errdhh{"Unrecognised type",type,"in  edsverbose"};
392   end;
393
394
395symbolic procedure edsdebug(msg,v,type);
396   % msg:string, v:various, type:id|nil
397   % -> edsdebug:nil
398   % Like edsverbose, just for debugging.
399   if !*edsdebug then edsverbose(msg,v,type);
400
401
402symbolic procedure edsmathprint f;
403   % f:prefix -> nil
404   % Similar to mathprint, except going in at writepri,
405   % so TRI package picks it up too.
406   <<writepri(mkquote f,'only); terpri!* t; >>;
407
408endmodule;
409
410end;
411