1module hidden;
2
3% Author: James H. Davenport.
4
5% Redistribution and use in source and binary forms, with or without
6% modification, are permitted provided that the following conditions are met:
7%
8%    * Redistributions of source code must retain the relevant copyright
9%      notice, this list of conditions and the following disclaimer.
10%    * Redistributions in binary form must reproduce the above copyright
11%      notice, this list of conditions and the following disclaimer in the
12%      documentation and/or other materials provided with the distribution.
13%
14% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
15% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
16% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
17% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR
18% CONTRIBUTORS
19% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
20% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
21% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
22% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
23% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
24% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
25% POSSIBILITY OF SUCH DAMAGE.
26%
27
28
29fluid '(!*trint
30        !*trjhd
31        cancellationlist
32        equalitystack
33        hide
34        indexlist
35        indexmap
36        inequalitystack
37        lastreduction
38        nindex
39        popstack
40        rhs
41        ulist);
42
43exports hiddentermgenerators;
44imports interr,nextprime,mksp,addf,getv,!*multf,printsf,gcdf,quotf,negf;
45imports makeprim,prin2,mkvect,copyvec,printdf;
46
47
48% Package for dealing with linear subspaces defined
49% by both equalities and inequalities
50
51
52% Major datastructure -
53% An equation is stored in a vector of length (length indexlist)
54% elements 2::n are the coeffs of the various index quantities
55% element 1 is constant part
56% element 0 is number defining a pivot in the row
57%
58% all vector elements are standard forms
59
60symbolic procedure newspace;
61 begin scalar w;
62    indexmap:=nil;
63    nindex:=1;
64    w:=indexlist;
65    while w do <<
66        nindex:=nindex+1;
67        indexmap:=(car w . nindex) . indexmap;
68        w:=cdr w >>;
69    popstack:=equalitystack:=inequalitystack:=nil;
70    if !*trint then printc "New space set up"
71 end;
72
73symbolic procedure currentspace;
74% Picks up the current set of linear constraints so they can be stored
75    (equalitystack . inequalitystack) . popstack;
76
77% Restores a space saved by 'currentspace'
78symbolic procedure setlinspace a;
79 << popstack:=cdr a;
80    a:=car a;
81    equalitystack:=car a;
82    inequalitystack:=cdr a;
83    nil >>;
84
85
86symbolic procedure newequation(eqqn,notflag);
87% Declare the equation
88%    eqqn=0
89% (or if notflag is set)    eqqn ne 0
90% to hold in addition to all previously declared relationships
91% return 'inconsistent if this is not possible
92 begin scalar newrow,newpivot,q1,r;
93    if !*trint then <<
94       prin2 "New ";
95       if notflag then prin2 "in";
96       printc "equality";
97       printsf eqqn >>;
98    if inequalitystack='inconsistent then interr "double inconsistency";
99    popstack:=(equalitystack . inequalitystack) . popstack;
100    newrow:=mkvect nindex; % Space for new row in eqqnstack
101%   for i:=0:nindex do putv(newrow,i,nil); %set to zero
102    expandinto(newrow,eqqn,1);
103    makeprim(newrow,nindex); % Remove any common factors
104    mapc(reverse equalitystack,function lambda j;
105      eliminatewith(newrow,j));
106    newpivot:=0;
107    for i:=2:nindex do
108      if getv(newrow,i) neq nil then newpivot:=i;
109    if newpivot=0 then
110       if (getv(newrow,1)=nil)=notflag then go to contradict
111         else return 'redundant;
112    putv(newrow,0,newpivot);
113    if notflag then <<
114      inequalitystack:=newrow . inequalitystack;
115      if !*trint then printc "New inequality added";
116      return 'consistent! inequality >>;
117    equalitystack:=newrow . equalitystack;
118    q1:=inequalitystack;
119    inequalitystack:=nil;
120wind:
121    if q1=nil then return 'consistent! equality;
122    r:=copyvec(car q1,nindex);
123    q1:=cdr q1;
124    eliminatewith(r,newrow);
125    newpivot:=0;
126    for i:=2:nindex do
127      if getv(r,i) neq nil then newpivot:=i;
128    if newpivot neq 0 then go to ok;
129    if getv(r,1)=nil then go to contradict;
130    go to wind;
131contradict:
132    equalitystack:=inequalitystack:='inconsistent;
133    return 'inconsistent;
134ok: putv(r,0,nil);
135    makeprim(r,nindex);
136    inequalitystack:=r . inequalitystack;
137    go to wind
138 end;
139
140symbolic procedure expandinto(row,val,above);
141    if domainp val then putv(row,1,addf(getv(row,1),!*multf(val,above)))
142    else begin scalar x;
143      x:=atsoc(mvar val,indexmap); % Process leading term first
144      if x then << x:=cdr x;
145        if not (tdeg lt val = 1) then
146          interr "not linear in expandinto";
147        putv(row,x,addf(getv(row,x),!*multf(lc val,above))) >>
148    else expandinto(row,lc val,!*multf(above,((lpow val) .* 1) .+ nil));
149    return expandinto(row,red val,above)
150 end;
151
152symbolic procedure eliminatewith(newrow,oldrow);
153   begin scalar pivotx,p,q;
154    pivotx:=getv(oldrow,0);
155    p:=getv(newrow,pivotx);
156    if p=nil then return nil; %Nothing needs doing
157    q:=getv(oldrow,pivotx);
158    begin scalar gg;
159      gg:=gcdf(p,q);
160      p:=quotf(p,gg);
161      q:=negf quotf(q,gg) end;
162    for i:=1:nindex do
163      putv(newrow,i,addf(!*multf(getv(newrow,i),q),
164                         !*multf(getv(oldrow,i),p)));
165    pivotx:=getv(newrow,0);
166    putv(newrow,0,nil);
167    makeprim(newrow,nindex);
168    putv(newrow,0,pivotx);
169    return nil
170   end;
171
172symbolic procedure dropequation;
173   begin
174    if atom popstack then interr "popstack underflow";
175    equalitystack:=caar popstack;
176    inequalitystack:=cdar popstack;
177    popstack:=cdr popstack
178   end;
179
180
181symbolic procedure printlinspace();
182 begin if equalitystack='inconsistent then <<
183      printc "never happens";
184      return nil >>;
185    if equalitystack=nil then
186      printc "no equalityies active"
187    else << printc "subject to"; printeqns equalitystack >>;
188    if atom inequalitystack then return nil;
189    printc "except when";
190    printeqns inequalitystack;
191    return nil end;
192
193symbolic procedure printeqns l;
194 begin while l do <<
195      printsf vec2sf(car l,indexmap);
196      terpri();
197      l:=cdr l >>;
198    terpri();
199    return nil end;
200
201symbolic procedure vec2sf(x,indexmap);
202 begin scalar r;
203    r:=getv(x,1); % constant part
204    for i:=2:nindex do begin
205        scalar v;
206        v:=indexmap;
207        while not (i=cdar v) do v:=cdr v;
208        v:=caar v;
209        v:=(mksp(v,1) .* 1) .+ nil;
210        r:=addf(r,!*multf(v,getv(x,i))) end;
211    return r end;
212
213
214
215
216
217symbolic procedure hiddentermgenerators(rhs);
218 begin scalar p,q,reductions;
219    cancellationlist:=nil;
220    hide:=nil;
221    newspace(); %Init linear subspace package
222    eliminateconstofint();
223    reductions:=constrainedreductions rhs;
224    if !*trint then mapc(reductions,function printreduction);
225    lastreduction:=lvlast reductions;
226    p:=cdr reductions;
227    while p do <<
228        q:=reductions;
229        while not (p eq q) do <<
230            interfere(car p,car q); %May extend reductions list
231            q:=cdr q >>;
232        p:=cdr p >>;
233    return nil
234   end;
235
236
237
238
239symbolic procedure lvlast l;
240    if null l then interr "l=nil in lvlast"
241    else if null cdr l then l
242    else lvlast cdr l;
243
244%symbolic procedure newreduction a;
245%% Record a new reduction formula
246%  begin
247%    a:=list a;
248%    rplacd(lastreduction,a);
249%    lastreduction:=a
250%   end;
251
252symbolic procedure eliminateconstofint();
253% Set up an inequality that gets rid of the constant of
254% integration but no other positive terms. maybe this is
255% done by a fudge.....
256 begin scalar p,r,m,v,x;
257    x:=indexlist;
258    p:=lpow ulist;
259% (a) find largest value in this thing
260    m:=0;
261    while p do << m:=max(m,car p); p:=cdr p >>;
262    v:=0;
263    p:=lpow ulist;
264% now allocate a prime number > m to each index
265    while p do <<
266        m:=nextprime m;
267        r:=addf(r,(mksp(car x,1) .* m) .+ nil);
268        x:=cdr x;
269        v:=v-m*car p;
270        p:=cdr p >>;
271    if not(v=0) then r:=addf(r,v);
272    newequation(r,t); %Assert the inequality
273    return nil
274  end;
275
276symbolic procedure constrainedreductions rhs;
277 begin scalar rels,fg;
278top:
279    if rhs=nil then return reversewoc rels;
280    fg:=newequation(numr lc rhs,t); %non-zero leading term
281    if not(fg='inconsistent) then
282        rels:=(rhs . currentspace()) . rels;
283    dropequation();
284    fg:=newequation(numr lc rhs,nil); %now lt = 0
285    if fg='inconsistent then return reversewoc rels;
286    rhs:=red rhs;
287    go to top end;
288
289symbolic procedure printreduction a;
290 begin
291    printc "**** REDUCTION FORMULA ****";
292    printdf car a;
293    setlinspace cdr a;
294    printlinspace()
295 end;
296
297symbolic procedure interfere(a,b);
298 begin
299   if !*trjhd then printc "potential cancellation detected"
300 end;
301
302
303%**********************************************************************;
304% From here down I have code left over from the previous version
305
306
307%symbolic procedure shiftmatch(form,p1,p2);
308%    sfsublis(form,shiftassoc(p1,p2,indexlist));
309
310%symbolic procedure shiftassoc(p1,p2,l);
311%    if null l then nil
312%    else ((car l) . (caar p1 - caar p2)) .
313%                         shiftassoc(cdr p1,cdr p2,cdr l);
314
315
316%symbolic procedure sfsublis(fm,l);
317%    if domainp fm then fm
318%    else begin
319%      scalar w;
320%      w:=atsoc(mvar fm,l);
321%      if w neq nil then <<
322%        w:=cdr w;
323%        if w=0 then w:=nil >>;
324%      w:=(mksp(mvar fm,1) .* 1) .+ w;
325%      w:=!*multf(w,sfsublis(lc fm,l));
326%      return addf(w,sfsublis(red fm,l)) end;
327
328
329
330
331
332%symbolic procedure constantofintegration();
333% begin scalar n;
334%    if not null inequalitystack then return nil;
335%%To be recognized, the term here must be totally constrained
336%%by equalities, and must match the c of i already set in ulist
337%    n:=length equalitystack;
338%    if not (n=nindex-1) then return nil;
339%    return matchp(lpow ulist,equalitystack,inequalitystack) end;
340
341
342%symbolic procedure dfsublis(p,l);
343%    dfsublis1(p,l,mapcar(l,function cdr));
344
345%symbolic procedure dfsublis1(p,l1,l2);
346%    if null p then nil
347%    else (lambda coef,rest;
348%      if null coef then rest else
349%        (pluslists(lpow p,l2) .* (coef ./ denr lc p)) .+ rest)
350%      (sfsublis(numr lc p,l1),dfsublis1(red p,l1,l2));
351
352%symbolic procedure pluslists(l1,l2);
353%    if null l1 then nil
354%    else ((caar l1 + car l2) . nil) . pluslists(cdr l1,cdr l2);
355
356
357%symbolic procedure invertclashes();
358% begin
359%    cancellationlist:=nil;
360%    mapc(hide,function clash);
361%    mapc(cancellationlist,function printcancel);
362% end;
363
364%symbolic procedure printcancel l;
365% begin
366%    equalitystack:=caar l;
367%    inequalitystack:=cdar l;
368%    printc "A phantom may be needed when the leading term satisfies";
369%    printlinspace();
370%    printc "offsets =";
371%    mapc(cdr l,function printc);
372%    terpri()
373% end;
374
375
376%symbolic procedure clash(a);
377%% A is the structure ((a . b) . (equal . inequal))
378%% create a corresponding entry on calcellationlist
379% begin
380%    scalar rf1,rf2,lineq;
381%    lineq:=cdr a;
382%    rf1:=caar a;
383%    rf2:=cdar a;
384%    a:=plusdf(red rf1,red rf2);
385%% Now lt a is maybe important term left over
386%    lineq:=shifteqns(car lineq,unlist lpow a) .
387%             shifteqns(cdr lineq,unlist lpow a);
388%    a:=formdeltas(lpow a,lpow rf1);
389%    cancellationlist:=(lineq . list a) . cancellationlist;
390%    return nil
391% end;
392
393%symbolic procedure unlist l;
394%    if null l then nil
395%    else caar l . unlist cdr l;
396
397
398%symbolic procedure formdeltas(a,b);
399%    if null a then nil
400%    else (caar a - caar b) . formdeltas(cdr a,cdr b);
401
402%symbolic procedure shifteqns(l,delta);
403%    if null l then nil
404%    else shifteqn(car l,delta) . shifteqns(cdr l,delta);
405
406
407%symbolic procedure shifteqn(v,delta);
408% begin scalar w,i,new;
409%    new:=mkvect nindex;
410%    for i:=0:nindex do putv(new,i,getv(v,i));
411%    v:=new;
412%    i:=2;
413%    w:=nil;
414%    while delta do <<
415%        w:=addf(w,!*multf(getv(v,i),car delta));
416%        delta:=cdr delta >>;
417%    putv(v,1,addf(getv(v,1),negf w));
418%    return v
419% end;
420
421endmodule;
422
423end;
424
425