1module edssolve;
2
3% Specialised solvers for EDS
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
31COMMENT. The EDS solve routines are interfaces to the general REDUCE
32solver, presenting the results in a form more useful in EDS. The most
33primitive is edssolve, which basically just turns off arbvars and
34reorganises the answer, and the most sophisticated is edssolvegraded,
35which tries to solve the principal part first if the equations are
36semilinear. This is useful since the principal part is often less
37complex than the rest. Solutions are returned as rmaps, giving
38information about both the solution and its applicability.
39
40endcomment;
41
42
43fluid '(!*edsverbose !*edsdebug !*arbvars !*varopt !*groebopt
44        !*solveinconsistent depl!* cfrmcrd!* cfrmrsx!* xvars!*
45        !*edssloppy !*edsdisjoint);
46
47
48symbolic procedure edssolvegraded(xl,vl,rsx);
49   % xl:list of sq, vl:list of list of kernel, rsx:rsx
50   % -> edssolvegraded:list of (t.rmap)|(nil.list of sq)
51   (begin scalar fl,z;
52   z := foreach l in vl join append(l,{});
53   if vl and semilinearsql(xl,car vl,if !*edssloppy then car vl else z)
54     then return edssolvesemi(xl,vl,rsx);
55   foreach l on vl do
56      foreach y in car l do foreach ll in cdr l do
57          foreach x in ll do depend1(y,x,t);
58   xl := edssolve(xl,z);
59   fl := {};
60   xl := foreach s in xl join
61      if null car s  then <<fl := s . fl; nil>>
62      else if not member(0,z := pullbackrsx(rsx,cadr s)) then
63          {{cadr s,purgersx append(z,caddr s)}};
64   return append(reverse fl,foreach s in edsdisjoin xl collect t . s);
65   end) where depl!* = depl!* ;
66
67
68symbolic procedure semilinearsql(xl,vl,kl);
69   % xl:list of sq, vl,kl:list of kernel -> semilinearsql:bool
70   null xl
71      or semilinearsq(car xl,vl,kl) and semilinearsql(cdr xl,vl,kl);
72
73
74symbolic procedure semilinearsq(x,vl,kl);
75   % x:sq, vl,kl:list of kernel -> semilinearsq:bool
76   % True if x is linear in vl with coefficients independent of kl.
77   % Assumes that kl includes vl.
78   semilinearpf0(xpartitsq x,vl,kl) where xvars!* = vl;
79
80
81symbolic procedure semilinearpf0(f,vl,kl);
82   % f:pf, vl,kl:list of kernel  -> semilinearpf0:bool
83   % Works when xvars!* allows 0-forms as well - used in solvegraded.
84   % True if x is linear in vl with coefficients independent of kl.
85   % Assumes that kl includes vl.
86   null f or
87      (lpow f = 1 and
88       freeoffl(denr lc f,vl) and   % because vl might occur inside
89                                    % operators
90       freeoffl(numr lc f,vl)
91                 or
92       length wedgefax lpow f = 1 and
93       freeoffl(denr lc f,kl) and   % because vl might occur inside
94                                    % operators
95       freeoffl(numr lc f,kl)) and
96      semilinearpf0(red f,vl,kl);
97
98
99symbolic procedure edssolvesemi(xl,vl,rsx);
100   % xl:list of sq, vl:list of list of kernel, rsx:rsx
101   % -> edssolvesemi:list of (t.rmap)|(nil.list of sq)
102   % xl is known to be semilinear
103   begin scalar sl,nl;
104   xl := edspartsolve(xl,car vl);
105   foreach x in car xl do
106      if not(car x memq cadr xl) then
107         sl := (car x . mk!*sq subs2 subsq(simp!* cdr x,caddr xl)) . sl
108      else if numr
109      << x := addsq(negsq !*k2q car x,simp!* cdr x);
110         x := subs2 subsq(x,caddr xl) >>
111      then nl := x . nl;
112   sl := !*map2rmap reversip sl;
113   if 0 member (rsx := pullbackrsx(rsx,car sl)) then return nil;
114   rsx := purgersx append(rsx,cadr sl); sl := car sl;
115   if null nl then return {t . {sl,rsx}};
116   xl := edssolvegraded(nl,cdr vl,rsx);
117   return foreach s in xl collect
118      if null car s then
119         nil . append(cdr s,foreach x in sl collect
120                                             addsq(negsq !*k2q car x,simp!* cdr x))
121      else
122         t . {append(pullbackmap(sl,cadr s),cadr s),caddr s};
123   end;
124
125
126symbolic procedure edssolve(xl,vl);
127   % xl:list of sq, vl:list of kernel -> edssolve: list of tag.value
128   % where tag.value is one of
129   %        t.rmap          an explicit solution
130   %        nil.list of sq    a (partial) system which couldn't be solved
131   % This is an interface to the REDUCE solve routines, for use by other
132   % eds routines. It switches off arbvars and returns the result in a
133   % form more easily used in eds. If the system is inconsistent, an
134   % empty list {} is returned.
135   begin scalar kl,sl,msg;
136         scalar !*arbvars;  % stop arbcomplex's being generated
137   msg := {"Solving",length xl,"equations in",length vl,"variables"};
138   foreach q in xl do kl := union(kernels numr q,kl);
139   kl := expanddependence kl;
140   vl := intersection(vl,kl);        % must preserve order of vl
141   msg := append(msg,{{length vl,"present"}}); % {a,b} prints as (a b)
142   edsdebug(msg,nil,nil);
143   sl := edssolvesys(foreach q in xl collect numr q,vl);
144   if eqcar(sl,'inconsistent) then
145      sl := {}
146   else if eqcar(sl,'failed) or eqcar(sl,nil) then
147      sl := {nil . xl}
148   else if eqcar(sl,t) then
149      sl := foreach s in cdr sl join
150               if not explicitsolutionp s then
151            {nil . foreach pr in pair(cadr s,car s) collect
152               addsq(negsq !*k2q car pr,cdr pr)}
153               else % need to reorder return value from edssolvesys!
154            {t . !*map2rmap pair(cadr s,for each q in car s
155                                           collect mk!*sq reordsq q)}
156   else errdhh{"Unexpected result from solvesys:",sl};
157   return sl;
158   end;
159
160
161symbolic procedure expanddependence kl;
162   % kl: list of kernel -> expanddependence:list of kernel
163   % expand kl recursively to include all kernels explicitly appearing
164   % as operator arguments. Should we check implicit dependence also?
165   if null kl then nil
166   else if atomf car kl then
167      car kl . expanddependence cdr kl
168   else % must be operator, so add all arguments to list
169      car kl . expanddependence append(cdar kl,cdr kl);
170
171
172symbolic procedure edssolvesys(xl,vl);
173   % xl:list of sf, vl:list of kernel -> edssolvesys: list of tag.value
174   % where tag.value is given in solve.red. This just calls solvesys.
175   solvesys(xl,vl);
176
177
178symbolic procedure explicitsolutionp s;
179   % s:solve solution -> explicitsolutionp:bool
180   not smember('root_of,s) and not smember('one_of,s);
181
182
183symbolic procedure edspartsolve(xl,vl);
184   % x:list of sq,
185   % vl:list of kernel -> edspartsolve:{map,list of kernel,map}
186   % Solves the equations xl for the variables vl by splitting into
187   % homogeneous and inhomogeneous parts.  The results are only
188   % guaranteed for linear equations whose coefficient are independent
189   % of the inhomogeneous parts.  The solution is returned in terms of
190   % some temporary variables representing the inhomogeneous parts.
191   % The temporary variables are given in the original variables by the
192   % third return value.
193   (begin scalar al,il,l;
194%        scalar depl!*;       % preserve dependencies
195   xl := splitlinearequations(xl,vl);
196   xl := foreach p in xl collect
197                  if null numr cdr p then
198               car p
199            else if (l := mk!*sq cdr p member il) then
200               addsq(car p,!*k2q nth(al,1 + length il - length l))
201            else
202            << al := intern gensym() . al;
203               il := mk!*sq cdr p . il;
204               addsq(car p,!*k2q car al) >>;
205   al := reversip al; il := reversip il;
206   foreach y in vl do
207      foreach z in al do depend1(y,z,t);
208   xl := edssolve(xl,append(vl,al));
209   if length xl neq 1 or null car(xl := car xl) then
210      errdhh{"Bad solution in edspartsolve",xl};
211   return {cadr xl,al,pair(al,il)};
212   end) where depl!* = depl!*;     % preserve dependencies
213
214
215symbolic procedure splitlinearequations(v,c);
216   % v:list of sq, c:list of kernel
217   % -> splitlinearequations:{list of sq.sq}
218   % Splits rational expressions in v into homogeneous and
219   % inhomogeneous parts wrt variables in c.
220   begin scalar ok,g,h;
221         scalar !*exp;
222   !*exp := t;
223   ok := setkorder c;
224   v := foreach q in v collect
225        << g := h := nil;
226                 while not domainp f and mvar f memq c do
227                    << g := lt f . g; f := red f >>;
228           foreach u in g do h := u .+ h;
229                 cancel(h ./ d) . cancel(f ./ d) >>
230         where f = reorder numr q,
231               d = reorder denr q;
232   setkorder ok;
233   return foreach p in v collect reordsq car p . reordsq cdr p;
234   end;
235
236
237symbolic procedure edsgradecoords(crd,jet0);
238   % crd,jet0:list of kernel -> edsgradecoords:list of list of kernel
239   % grade coordinates according to jet order, highest jet first
240   if null jet0 then {crd}
241   else
242      begin scalar u;
243      foreach c in crd do
244         begin scalar j0,c0;
245         j0 := jet0;
246         while j0 and null(c0 := splitoffindices(car j0,c)) do
247            j0 := cdr j0;
248         if j0 := assoc(length c0,u) then nconc(j0,{c})
249         else u := (length c0 . {c}) . u;
250         end;
251      return foreach v in sort(u,function greaterpcar) collect cdr v;
252      end;
253
254
255
256COMMENT. The routine solvepfsys tries to bring a system into solved form
257in the current environment specified by cfrmcrd!* and cfrmrsx!*.
258
259endcomment;
260
261symbolic procedure solvepfsys sys;
262   % sys:sys -> solvepfsys:{sys,sys}
263   % Bring sys into solved form as far as possible.
264   solvepfsys1(sys,{});
265
266symbolic procedure solvepfsys1(sys,vars);
267   % sys:sys, vars:list of lpow pf -> solvepfsys1:{sys,sys}
268   % Solve sys for vars. If vars is {} then solve for anything. Kernel
269   % ordering changed so that solved set comes ahead of rest. Ordering
270   % within solved set and others is unchanged. The solved part is
271   % returned sorted in decreasing order of lpow.  NBB. Kernel ordering
272   % changed!!
273   begin scalar ok,sl,nl;
274   % save old kernel ordering
275   ok := updkordl nil;
276   nl := foreach f in sys collect subs2pf!* f;
277   % First try for constant coefficients
278   if nl then
279      begin
280%      nl := solvepfsyseasy(nl,vars,'domainp);
281      nl := solvepfsyseasy(nl,vars,'cfrmconstant);
282      if car nl then
283          edsdebug("Solved with constant coefficients:",car nl,'sys);
284      if cadr nl then
285          edsdebug("Unsolved with constant coefficients:",cadr nl,'sys);
286      sl := append(sl,car nl);
287      nl := cadr nl;
288      end;
289   % If that's not enough, try for nowhere-zero coefficients
290   if nl then
291      begin
292      nl := solvepfsyseasy(nl,vars,'cfrmnowherezero);
293      if car nl then
294          edsdebug("Solved with nonzero coefficients:",car nl,'sys);
295      if cadr nl then
296          edsdebug("Unsolved with nonzero coefficients:",cadr nl,'sys);
297      sl := append(sl,car nl);
298      nl := cadr nl;
299      end;
300   % If that's not enough, try Cramer's rule.
301   if nl then
302      begin
303      nl := solvepfsyshard(nl,vars);
304      if car nl then
305          edsdebug("Solved with Cramer's rule:",car nl,'sys);
306      if cadr nl then
307          edsdebug("Unsolved with Cramer's rule:",cadr nl,'sys);
308      sl := append(sl,car nl);
309      nl := cadr nl;
310      end;
311   % Fix up kernel ordering and back-substitute solved forms
312   setkorder ok;
313   updkordl lpows sl;
314   sl := xautoreduce1 xreordersys sl;
315   nl := xreordersys nl;
316   % Block structure may have messed up order of solved kernels
317   sl := sortsys(sl,ok);
318   updkordl lpows sl;
319   return {sl,nl};
320   end;
321
322
323symbolic procedure solvepfsyseasy(sys,vars,test);
324   % sys:sys, vars:list of lpow pf,
325   % test:function -> solvepfsyseasy:{sys,sys}
326   % Recursively bring sys into weakly solved form as far as possible
327   % just by changing kernel ordering and normalising.  The solved part
328   % is in normalised upper triangular form, and korder is completely
329   % messed up - results must be reordered under lpows solved part.
330   begin scalar sl,nl,kl;
331   kl := purge foreach f in sys join xsolveables(f,test);
332   kl := sort(if vars then intersection(kl,vars) else kl,'termordp);
333   if null kl then return {{},sys};
334   updkordl kl;
335   foreach f in sys do
336      if (f := xreduce(xreorder f,sl)) and apply1(test,numr lc f) then
337         sl := xnormalise f . sl
338      else if f then
339         nl := f . nl;
340   sl := reversip sl; % sl in upper triangular form
341   nl := reversip foreach f in nl join
342      if f := subs2pf!* xreduce(f,sl) then {f};
343   if null sl or null nl then return {sl,nl};
344   nl := solvepfsyseasy(nl,vars,test);
345   return {append(sl,car nl),cadr nl};
346   end;
347
348
349symbolic procedure xsolveables(f,test);
350   % f:pf, test:function -> xsolveables:list of lpow pf
351   % All powers in f whose coefficient numerators satisfy test
352   if null f then nil
353   else if apply1(test,numr lc f) then lpow f . xsolveables(red f,test)
354   else xsolveables(red f,test);
355
356
357symbolic procedure solvepfsyshard(sys,vars);
358   % sys:sys, vars:list of lpow pf -> solvepfsyshard:{sys,sys}
359   % Bring sys into weakly solved form by Cramer's rule.  The solved
360   % part is in normalised upper triangular form, and korder is
361   % completely messed up - results must be reordered under lpows
362   % solved part.  Allowing !*edssloppy to work here for the first time
363   % minimises the number of restrictions added.
364   begin scalar sl,nl,kl,w;
365   if null sys then return {{},{}};
366   if vars then updkordl vars;
367   w := !*k2pf 1;
368   foreach f in sys do
369      if f := subs2pf wedgepf(if vars then xreorder f else f,w) then
370         if null vars or subsetp(wedgefax lpow f,vars) then
371            w := f
372         else
373         << if degreepf f neq 1 then
374               f := xcoeff(f,intersection(wedgefax lpow f,vars));
375            nl := multpfsq(f,invsq lc w) . nl >>; % exact divisions
376   kl := xsolveables(w,'cfrmnowherezero);
377   if null kl and !*edssloppy then kl := xpows w;
378   if vars then while kl and not subsetp(wedgefax car kl,vars) do
379         kl := cdr kl;
380   if null kl or (car kl = 1) then return {{},sys};
381   kl := wedgefax car kl;
382   if !*edssloppy then
383      cfrmrsx!* := xpartitsq(numr lc xcoeff(w,kl) ./ 1) . cfrmrsx!*
384                                     where xvars!* = cfrmcrd!*;
385   sl := foreach k in kl collect xcoeff(w,delete(k,kl));
386   updkordl kl;
387   return {foreach f in sl collect xnormalise xreorder f,
388                 foreach f in nl collect xrepartit!* f};
389   end;
390
391endmodule;
392
393end;
394