1module edsuser;
2
3% Miscellaneous user 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
31
32fluid '(alglist!* subfg!* !*arbvars xvars!*);
33
34global '(indxl!*);
35
36
37put('index_expand,'rtypefn,'quotelist);
38put('index_expand,'listfn,'indexexpandeval0);
39
40symbolic procedure indexexpandeval0(u,v);
41   % u:list of prefix, v:bool -> indexexpandeval0:prefix list
42   % kludge to add v argument to index_expand listfn
43   makelist foreach p in getrlist indexexpandeval u collect reval1(p,v);
44
45
46symbolic procedure indexexpandeval u;
47   % u:list of prefix -> indexexpandeval:prefix list
48   if length u neq 1 then
49      rerror(eds,000,"Wrong number of arguments to index_expand")
50   else if rlistp(u := reval car u) then
51      makelist purge foreach x in cdr u join cdr indexexpandeval {x}
52   else makelist indexexpand u;
53
54
55symbolic procedure indexexpand u;
56   % u:prefix -> indexexpand:list of prefix
57   % u has always been reval'd, so there is no need to expand
58   % summations.
59   if eqexpr u then indexexpandeqn u else
60   if boolexpr u then indexexpandbool u else
61   begin scalar i,v,alglist!*;
62   u := simp!* u;
63   % Expand free indices (put them into some order for safety)
64   i := idsort purge if numr u and not domainp numr u then
65                              flatindxl allind !*t2f lt numr u;
66   v := foreach j in mkaindxc(i,nil) join
67                 if numr(j := subfreeindices(numr u,pair(i,j))) then
68              {absf numr j ./ denr j};  % nprimitive too?
69   return for each q in purge v collect mk!*sq multsq(q,1 ./ denr u);
70   end;
71
72
73symbolic procedure indexexpandeqn u;
74   % u:rule|equation -> indexexpandeqn:list of rule|equation
75   begin scalar i,v,lh,rh;
76               scalar alglist!*;
77   % Expand free indices on lhs (put them into some order for safety)
78   lh := reval cadr u where subfg!* = nil; % avoid let rules
79   i := idsort purge flatindxl allindk lh;
80   rh := aeval caddr u;
81   v := foreach j in mkaindxc(i,nil) join
82      if j := subfreeindeqn({car u,lh,rh},pair(i,j)) then {j};
83   % Remove duplicates
84   i := {};
85   v := foreach r in v join
86      if not(cadr r member i) then
87               << i := cadr r . i; {r} >>;
88   return v;
89   end;
90
91
92symbolic procedure subfreeindeqn(u,l);
93   % u:rule|equation, l:alist -> subfreeindeqn:rule|equation|nil
94   % Make index substitution l in u. Only index symmetry simplifications
95   % are allowed, so the lhs can either vanish (nil returned) or acquire
96   % an overall sign (sign transferred to rhs).
97   begin scalar lh,rh;
98   lh := subfreeindk(cadr u,l);
99   if null atomf lh then lh := revop1 lh; % gets done in rule!-list
100                                          % anyway
101   lh := reval lh where subfg!* = nil; % avoid let rules;
102   if lh = 0 then return nil;
103   rh := simp!* caddr u;
104   rh := quotsq(subfreeindices(numr rh,l),subfreeindices(denr rh,l));
105   if eqcar(lh,'minus) then
106   << lh := cadr lh;
107      rh := negsq rh >>;
108   return {car u,lh,mk!*sq rh};
109   end;
110
111
112symbolic procedure boolexpr u;
113   % u:any -> boolexpr:bool
114   not atom u and flagp(car u,'boolean);
115
116
117symbolic procedure indexexpandbool u;
118   % u:prefix -> indexexpandbool:list of prefix
119   begin scalar i,v,alglist!*;
120   % Expand free indices on lhs (put them into some order for safety)
121   i := idsort purge flatindxl allindk u;
122   v := foreach j in mkaindxc(i,nil) collect
123           car u .
124              foreach a in cdr u collect reval subfreeindk(a,pair(i,j));
125   return purge v;
126   end;
127
128
129symbolic procedure subfreeindices(u,l);
130   % u:sf, l:a-list -> subfreeindices:sq
131   % Discriminates indices from variables, modified from EXCALC's
132   % subfindices to go inside operators other than EXCALC's.
133   begin scalar alglist!*;
134   return
135      if domainp u then u ./ 1
136      else
137         addsq(
138            multsq(if atom mvar u then
139                             !*p2q lpow u
140                          else if sfp mvar u then
141                      exptsq(subfreeindices(mvar u,l),ldeg u)
142                   else
143                      exptsq(simp subfreeindk(mvar u,l),ldeg u),
144                   subfreeindices(lc u,l)),
145            subfreeindices(red u,l))
146   end;
147
148
149symbolic procedure subfreeindk(u,l);
150   % u:kernel, l:a-list -> subfreeindk:kernel
151   % Extends subindk to indexed variables
152   if atom u then u
153   else if flagp(car u,'indexvar) then
154      car u . subla(l,cdr u)
155   else subindk(l,u);
156
157
158put('linear_divisors,'rtypefn,'quotelist);
159put('linear_divisors,'listfn,'lineardivisors);
160
161symbolic procedure lineardivisors(u,v);
162   % u:{prefix}, v:bool -> lineardivisors:prefix list
163   makelist foreach f in lineardivisorspf xpartitop reval car u collect
164      !*pf2a1(f,v);
165
166
167symbolic procedure lineardivisorspf f;
168   % f:pf -> lineardivisorspf:list of pf
169   begin scalar x,g,v;
170   foreach p in xpows f do x := union(wedgefax p,x);
171   foreach k in x do
172   << v := intern gensym() . v;
173      g := addpf(k .* !*k2q car v .+ nil,g)>>;
174   x := edssolve(xcoeffs wedgepf(g,f),v);
175   if length x neq 1 or caar x neq t then
176      errdhh "Bad solve result in lineardivisorspf";
177   x := cadar x;
178   v := updkordl v;
179   g := numr subf(numr !*pf2sq g,x);
180   x := {};
181   while g do
182   << x := xpartitsq(lc g ./ 1) . x;
183      g := red g >>;
184   setkorder v;
185   return reverse xautoreduce x;
186   end;
187
188
189symbolic procedure xdecomposepf f;
190   % f:pf -> xdecomposepf:list of pf
191   begin scalar x;
192   x := lineardivisorspf f;
193   if length x = degreepf f then return reverse x;
194   end;
195
196
197put('exfactors,'rtypefn,'quotelist);
198put('exfactors,'listfn,'exfactors);
199
200symbolic procedure exfactors(u,v);
201   % u:{prefix}, v:bool -> exfactors:prefix list
202   makelist foreach f in xfactorspf xpartitop reval car u collect
203      !*pf2a1(f,v);
204
205
206symbolic procedure xfactorspf f;
207   % f:pf -> xfactorspf:list of pf
208   begin scalar x;
209   x := lineardivisorspf f;
210   f := xreduce(f,foreach g in x collect
211                                  addpf(1 .* (-1 ./ 1) .+ nil,g));
212   return
213      if f = !*k2pf 1 then reverse x
214      else f . reverse x;
215   end;
216
217
218symbolic operator exact;
219symbolic procedure exact u;
220   % u:prefix -> exact:bool
221   % True if u is an exact pform kernel
222   eqcar(u,'d);
223
224flag('(exact),'boolean);
225
226
227put('derived_system,'rtypefn,'getrtypecar);
228put('derived_system,'edsfn,'deriveeds);
229put('derived_system,'listfn,'derivelist);
230
231symbolic procedure derivelist(u,v);
232   % u:{xeds|rlist}, v:bool -> derivelist:rlist
233   !*sys2a1(derive !*a2sys reval car u,v) where xvars!* = nil;
234
235
236symbolic procedure deriveeds s;
237   % s:eds -> deriveeds:eds
238   begin
239   s := copyeds s;
240   if pfaffian s then
241      eds_sys s := derive pfaffpart eds_sys s
242   else rerror(eds,000,"non-Pfaffian system in derived_system");
243   return s;
244   end;
245
246
247symbolic procedure derive s;
248   % s:sys -> derive:sys
249   begin scalar c,f;
250   if null s then s;
251   s := xautoreduce s;
252   c := for each f in s collect
253               if degreepf f = 1 then intern gensym()
254         else rerror(eds,000,"non-Pfaffian system in derived_system");
255   f := zippf(s,foreach k in c collect !*k2q k);
256   s := edssolve(xcoeffs xreduce(exdfpf f,s),c);
257   if length s neq 1 or null caar s then
258      errdhh{"Bad solve result in derive:",s};
259   s := cadr car s;
260   f := pullbackpf(f,s);
261   c := setdiff(c,foreach m in s collect car m);
262   f := xrepartit f where xvars!* = c;
263   return for each x in reverse c collect xcoeff(f,{x});
264   end;
265
266
267symbolic procedure allcoords f;
268   % f:prefix -> allcoords:list of kernel
269   % Pick up 0-form kernels in f
270   makelist purge
271      foreach k in (xvarspf xpartitop f where xvars!* = t) join
272               if xdegree k = 0 and
273            not assoc(k,depl!*) and
274             not eqcar(k,'partdf) then {k}
275         else if (xdegree k = 1) and exact k then {cadr k};
276
277endmodule;
278
279end;
280