1module sstools;
2
3create!-package('(sstools),nil);
4
5%********************************************************************
6%                                                                   *
7%  The package SSTOOLS for the study of integrability of            *
8%                    super-symmetric PDEs                           *
9%                                                                   *
10%  Author: Thomas Wolf and Eberhard Schruefer with contributions    *
11%          of Winfried Neun                                         *
12%                                                                   *
13%********************************************************************
14
15global '(!*sstools!-loaded);
16
17!*sstools!-loaded := t;
18
19%-----------------------
20% The following module for non-commutative differentiation and
21% multiplication of fermionic forms is written by Eberhard Schruefer 2007.
22
23module dop$
24
25fluid '(t_changes_parity s_changes_parity asymplis!* ncmp!*)$
26
27switch t_changes_parity,s_changes_parity$
28
29put('t_changes_parity,'simpfg,
30    '((t (rmsubs) (setq t_changes_parity t) (clear_t_changes_parity))
31         (nil (rmsubs) (setq t_changes_parity nil) (clear_t_changes_parity))))$
32
33put('s_changes_parity,'simpfg,
34    '((t (rmsubs) (setq s_changes_parity t) (clear_s_changes_parity))
35         (nil (rmsubs) (setq s_changes_parity nil) (clear_s_changes_parity))))$
36
37symbolic procedure clear_t_changes_parity()$
38   begin
39     asymplis!* := clear_t_changes_power_rules asymplis!*;
40   end$
41
42symbolic procedure clear_s_changes_parity()$
43   begin
44     asymplis!* := clear_s_changes_power_rules asymplis!*;
45   end$
46
47symbolic procedure clear_t_changes_power_rules u$
48   if null u then u
49    else if smember('t,car u) then clear_t_changes_power_rules cdr u
50    else car u . clear_t_changes_power_rules cdr u$
51
52symbolic procedure clear_s_changes_power_rules u$
53   if null u then u
54    else if smember('s,car u) then clear_s_changes_power_rules cdr u
55    else car u . clear_s_changes_power_rules cdr u$
56
57ncmp!* := t$
58
59rlistat('(fermion))$
60
61symbolic procedure fermion u$ fermion1 u$
62
63symbolic procedure fermion1 u$
64   <<flag(u,'fermionic);
65     for each j in u
66       do <<flag({j},'full); put(j,'simpfn,'simpfermion)>>;
67     for each j in u do noncom1 j>>$
68
69symbolic procedure simpfermion u$
70   begin scalar x,y,z;
71     if x := opmtch u then return simp x;
72     u := car u . for each j in cdr u collect aeval j;
73     y := mksq(u,1);
74     if domainp numr y then return y;
75     z := mvar numr y;
76     if atom z or null(car z eq car u) then return y;
77     if null atsoc(z,asymplis!*)
78        then asymplis!* := (z . 2) . asymplis!*;
79     return y
80   end$
81
82rlistat('(boson))$
83
84symbolic procedure boson u$ boson1 u$
85
86symbolic procedure boson1 u$
87   for each j in u do
88     <<flag({j},'bosonic); mkop j>>$
89
90put('d,'simpfn,'simpdop)$
91
92noncom1 'd$
93
94symbolic procedure simpdop u$
95   begin scalar x;
96     u := aeval car u . cdr u;
97     if x := opmtch ('d . u) then return simp x;
98%     if x := opmtch u then return simp u;
99     x := simp cadr u;
100     if null numr x then return x;
101     if domainp denr x or null fermionicp mvar denr x
102        then return multsq(dopf(car u, numr x), 1 ./ denr x);
103     rederr "encountered non constant den in simpdop ...";
104   end$
105
106symbolic procedure dopf(n,u)$
107   if domainp u then nil ./ 1
108    else addsq(addsq(multpq(lpow u,dopf(n,if fermionicp mvar u
109                                             then negf lc u
110                                           else lc u)),
111                     multsq(dopp(n, lpow u),lc u ./ 1)),dopf(n, red u))$
112
113symbolic procedure dopp(n, u)$
114   % the following test that car u is a constant wrt d should be checked if
115   % it is complete.
116   if atom car u or null(flagp(caar u,'fermionic) or flagp(caar u,'bosonic) or
117                         (caar u eq 'd) or (caar u eq 'df))
118      then nil ./ 1
119    else if cdr u > 1
120            then multsq(cdr u ./ 1,multsq(simpexpt {car u,cdr u-1},
121                                          simpdop {n,car u}))
122    else if caar u eq 'd
123            then if n = cadr car u
124                    then simpdf {caddr car u,'x}
125                  else if n > cadr car u
126                    then negsq simpdop({cadr car u,{'d,n,caddr car u}})
127         else mkdsq({'d,n,car u},1)
128    else mkdsq({'d,n,car u},1)$
129
130symbolic procedure mkdsq(u,n)$
131   begin scalar x,z;
132    if x := opmtch u then return simp x;
133     x := mksq(u,n);
134     if null numr x then return x;
135     z := mvar numr x;
136     if fermionicp z and null atsoc(z,asymplis!*)
137        then asymplis!* := (z . 2) . asymplis!*;
138     return x
139   end$
140
141symbolic procedure dfdp(u,v,n)$
142   begin scalar x;
143     x := simpdf {caddr u,v,n};
144     if null numr x then return x;
145     if ((v eq 't) and t_changes_parity) or
146        ((v eq 's) and s_changes_parity)
147        then x := negsq x;
148     if domainp denr x or
149                null fermionicp mvar denr x
150        then return multsq(dopf(cadr u, numr x), 1 ./ denr x)
151      else rederr "fermioic den encountered in dfdp"
152   end$
153
154put('d,'dfform,'dfdp)$
155
156symbolic procedure fermionicp u$
157   null atom u and (flagp(car u,'fermionic) or
158                    (car u eq 'd and null fermionicp caddr u) or
159                    (car u eq 'df and fermionicp caddr u) or
160                    (car u eq 'df and
161                     if fermionicp cadr u
162                        then if memq('s,cddr u) then not s_changes_parity
163                              else if memq('t,cddr u) then not t_changes_parity
164                              else t
165                      else if memq('s,cddr u) then s_changes_parity
166                            else if memq('t,cddr u) then t_changes_parity
167                            else nil))$
168
169
170symbolic procedure sstools!-multfnc(u,v)$
171   % Returns canonical product of U and V, with both main vars non-
172   % commutative.
173   begin scalar x,y;
174      x := multf(lc u,!*t2f lt v);
175      if null x then nil
176       else if not domainp x and mvar x eq mvar u
177        then x := addf(if null (y := mkspm(mvar u,ldeg u+ldeg x))
178                         then nil
179                        else if y = 1 then lc x
180                        else !*t2f(y .* lc x),
181                       multf(!*p2f lpow u,red x))
182       else if noncomp mvar u
183               then if ordop(mvar u,mvar x)
184                       then x := !*t2f(lpow u .* x)
185                     else if null(y := multf(!*p2f lpow u,lc x))
186                             then x := multf(!*p2f lpow u,red x)
187                     else x := addf(!*t2f(lpow x .*
188                                    if fermionicp mvar u and
189                                       fermionicp mvar x then negf y
190                                                         else y),
191                                    multf(!*p2f lpow u,red x))
192       else x := multf(!*p2f lpow u,x) where !*!*processed=t;
193      return addf(x,addf(multf(red u,v),multf(!*t2f lt u,red v)))
194   end$
195
196symbolic procedure difff(u,v)$
197   %U is a standard form, V a kernel.
198   %Value is the standard quotient derivative of U wrt V.
199   % Allow for differentiable domains.
200   if atom u then nil ./ 1
201    else if atom car u
202     then (if diff!-fn then apply2(diff!-fn,u,v) else nil ./ 1)
203        where (diff!-fn =  get(car u,'domain!-diff!-fn))
204    else addsq(addsq(multpq(lpow u,difff(if (fermionicp mvar u and
205                                             (((v eq 's) and s_changes_parity) or
206                                              ((v eq 't) and t_changes_parity) or
207                                              fermionicp v
208                                             )
209                                            )
210                                            then negf lc u
211                                            else      lc u,v)),
212                        multsq(if cdr lpow u = 1 and (((v eq 's) and s_changes_parity) or
213                                                      ((v eq 't) and t_changes_parity) or
214                                  fermionicp car lpow u) then diffdp(lpow u,v)
215                                else diffp(lpow u,v),lc u ./ 1)),
216               difff(red u,v))$
217
218symbolic procedure diffdp(u,v)$
219  begin scalar x,y,z,w;
220    u := car u;
221    if null atom u and (x := get(car u,'dfform)) then return apply3(x,u,v,1);
222    if null depends(u,v) then return nil ./ 1;
223    if u eq v then return 1 ./ 1;
224    if eqcar(u,'df) and v eq 's
225       and t_changes_parity and s_changes_parity
226       and uneven_t_p cddr u
227       then z := -1 ./ 1 else z := 1 ./ 1;
228    if car u eq 'df
229      then if (x := find_sub_df(w := cadr u . merge!-ind!-vars(u,v),
230                    get('df,'kvalue)))
231           then <<w := simp car x;
232                  for each el in cdr x do
233                  for i := 1:cdr el do
234                  w := diffsq(w,car el);
235                  return w>>
236           else u := 'df . w
237      else u := {'df,u,v};
238      w := 1 ./ 1;
239      if x := opmtch u then return multsq(z,simp x);
240      x := mksq(u,1);
241      if null numr x then return x;
242      y := mvar numr x;
243    if fermionicp y and null atsoc(y,asymplis!*)
244    then asymplis!* := (y . 2) . asymplis!*;
245 return multsq(z,x)
246end$
247
248symbolic procedure uneven_t_p u$
249   if car u eq 't
250      then if cdr u then null(fixp cadr u and evenp cadr u)
251            else t
252    else cdr u and uneven_t_p cdr u$
253
254bothtimes put('is_fermionic,'boolfn,'evalfermionicp)$
255
256symbolic procedure evalfermionicp u$
257fermionicfp numr simp u$
258% begin scalar x;
259%  x := numr simp u;
260%  return fermionicfp x
261% end;
262
263symbolic procedure fermionicfp u$
264if domainp u then nil
265%else (fermionicp mvar u and fermionicfp lc u) or fermionicfp red u;
266             else
267if fermionicp mvar u then not fermionicfp lc u
268                     else     fermionicfp lc u$
269
270% The following function makes reordering of sstools expressions possible
271% (as needed by coeffn for example). Unfortunately some redefinitions
272% of standard REDUCE functions are necessary.
273
274symbolic procedure reordablep(u,v)$
275   atom u or atom v or reordablekp(u,v)$
276
277symbolic procedure reordablekp(u,v)$
278   sstools_kernelp u or sstools_kernelp v$
279
280symbolic procedure sstools_kernelp u$
281   car u eq 'd or flagp(car u,'fermionic) or flagp(car u,'bosonic)$
282
283symbolic procedure reordop(u,v)$
284   if reordablep(u,v) then ordop(u,v)
285      else (!*ncmp and noncomp1 u and noncomp1 v) or ordop(u,v)$
286
287symbolic procedure rmultpf(u,v); !*q2f simp!* prepf (u .* v .+ nil)$
288
289
290endmodule$
291
292%----------------------- further declarations and initializations
293
294%algebraic$
295
296fermion1 '(f th)$
297boson1 '(b)$
298
299algebraic$
300depend  {f,b},t,x$
301
302!#if (memq 'psl lispsystem!*)
303set_bndstk_size 50000$
304!#endif
305
306lisp <<simplimit!* := 100000$
307       !*noarg := nil  % not done as 'off noarg$' to prevent compiler message
308     >>$
309on dfprint$
310
311lisp(!*allfac := nil)$ % instead of OFF ALLFAC in order to avoid difference
312                       % between compiled and uncompiled version
313
314symbolic fluid '(x flist blist fblist sysfbl N_ fname_ fname_list nfct_
315                 use_new_crackout print_ print_more collect_sol sol_list
316                 proc_list_ max_factor flin_ subst_1 pdelimit_1 old_history
317                 homogen_ max_gc_short max_gc_red_len record_hist
318                 !*time lin_test_const lines_written ineq_ html_out
319                 !*t_changes_parity !*s_changes_parity)$
320
321lrule1:={D(~n,df(~h,x   )) => lin_test_const**3,
322         D(~n,df(~h,x,~m)) => lin_test_const**(2*m+1)}$
323lrule2:={D(~n,~h)          => lin_test_const,
324         df(~h,x)          => lin_test_const**2,
325         df(~h,x,~n)       => lin_test_const**(2*n)}$
326
327subli:={}$       % the list of substitutions coming with the call of ssym
328symansatz:={}$   % the ansatz for equations and symmetry
329
330%----------------------- procedures
331
332if lisp(null(getd 'redfront_color) ) then
333symbolic procedure redfront_color(a)$ a$
334
335%-------
336
337symbolic operator is_const$
338symbolic procedure is_const(a)$
339freeof(a,'f) and freeof(a,'b)$
340
341%-------
342
343symbolic operator is_fermion$
344symbolic procedure is_fermion(a)$
345if pairp a then
346if car a = 'D  then is_boson caddr a
347               else
348if car a = 'DF then if is_fermion cadr a then
349   if caddr a ='s then not !*s_changes_parity else
350   if caddr a ='t then not !*t_changes_parity else t
351                                         else
352   if caddr a ='s then !*s_changes_parity else
353   if caddr a ='t then !*t_changes_parity else nil
354               else
355if car a = 'EXPT then
356   if is_boson cadr a then nil else
357   if (caddr a = 1) then t else
358   if fixp caddr a then nil else
359   write"###### ERROR: UNDEFINED PARITY!! ######"
360                 else
361if car a = 'TIMES then begin scalar n,h;
362                        n:=0;
363                        for each h in cdr a do
364                        if is_fermion h then n:=add1 n;
365                        return not evenp n
366                       end
367                  else
368if (car a = 'MINUS) or
369   (car a = 'QUOTIENT) then is_fermion cadr a else
370   % we assume, the denominator can not be fermionic
371if (car a = 'PLUS) then begin scalar h;
372                         h:=cdr a;
373                         while h and zerop car h do h:=cdr h;
374                         return if h then is_fermion car h
375                                     else nil
376                        end
377                   else
378if not freeof(a,'f) then t else nil
379           else if not freeof(a,'f) then t else nil$
380
381symbolic operator is_boson$
382symbolic procedure is_boson(a)$
383if pairp a then
384if car a = 'D  then is_fermion caddr a
385               else
386if car a = 'DF then if is_boson cadr a then
387   if caddr a ='s then not !*s_changes_parity else
388   if caddr a ='t then not !*t_changes_parity else t
389                                       else
390   if caddr a ='s then !*s_changes_parity else
391   if caddr a ='t then !*t_changes_parity else nil
392               else
393if car a = 'EXPT then
394   if is_boson cadr a then t else
395   if (caddr a = 1) then nil else
396   if fixp caddr a then t else
397   write"###### ERROR: UNDEFINED PARITY!! ######"
398                 else
399if (car a = 'MINUS) or
400   (car a = 'QUOTIENT) then is_boson(cadr a) else
401if car a = 'TIMES then begin scalar n,h;
402                        n:=0;
403                        for each h in cdr a do
404                        if is_fermion h then n:=add1 n;
405                        return evenp n
406                       end
407                  else
408if (car a = 'PLUS) then begin scalar h;
409                         h:=cdr a;
410                         while h and zerop car h do h:=cdr h;
411                         return if h then is_boson car h
412                                     else t
413                        end
414                   else
415if not freeof(a,'b) then t else nil
416           else if not freeof(a,'b) then t else nil$
417
418%-------
419
420symbolic operator ssini$
421symbolic procedure ssini(N,nf,nb)$
422% initializes N_, fblist and selects printing routines for f(i) and b(j)
423begin scalar j;
424 N_:=N;   % assignung the global number of super symm generators N_
425          % For customized printing:"$
426 if nf=1 then put('f,'prifn,'myfpri)  % to suppress (1) when printing f(1)
427         else put('f,'prifn,nil)$     % to print index of all f(i)
428 if nb=1 then put('b,'prifn,'myfpri)  % to suppress (1) when printing b(1)
429         else put('b,'prifn,nil)$     % to print index of all b(i)
430 fblist:=append(for j:=1:nf collect {'f,j},
431                for j:=1:nb collect {'b,j} )$
432end$
433
434%-------
435
436symbolic procedure presimplify(es,fl)$
437% es,fl are algebraic lists of equations and functions
438begin scalar w,ineq_bak,nopowersbak,eqns,m,k,cpu,p,len_es$
439 % Simplifying equations and dropping multiple ones
440 len_es:=sub1 length es$
441 write len_es," equations result"$terpri()$
442 cpu:=time();
443 % write((k:=time())-cpu)/1000,
444 %      " s: Now simplifying equations and dropping multiple ones"$terpri()$
445 w:=nil;   % list of equations consisting of only a coefficient
446 ineq_bak:=ineq_$
447 ineq_:=nil$ % for simplifyterm
448 nopowersbak:=!*nopowers;
449 algebraic(off nopowers);
450 eqns:=nil$
451 fl:=cdr reval fl;
452 m:=nil; % list of new equations consisting of only a coefficient
453 repeat <<
454  if m then <<
455   es:=algebraic(sub(lisp cons('LIST,for each k in m collect {'EQUAL,k,0}),lisp es))$
456   % w:=append(m,w)$
457   w:=nconc(m,w)$
458   m:=nil;
459  >>$
460  es:=cdr es$         % to get rid of the leading 'LIST
461  while es do
462  if zerop car es then es:=cdr es else <<
463   k:=algebraic(factorize lisp car es); es:=cdr es$
464   k:=cdr k;
465   p:=nil;
466   while k do <<
467    if not freeoflist(cadr car k,fl) then p:=cons(cadr car k,p);
468    k:=cdr k;
469   >>$
470   if null p then eqns:=cons(1,eqns) else <<
471    if cdr p then p:=cons('TIMES,p)
472             else p:=simplifyterm(reval car p,fl);
473    if atom p then m   :=union({p},   m)
474              else eqns:=union({p},eqns)
475   >>
476  >>$
477  es:=cons('LIST,eqns); eqns:=nil
478 >> until null m;
479
480 write length w," coefficients found to be zero: "$listprint w$ terpri()$
481 eqns:=nconc(es,w)$
482 write reval (len_es+1-length eqns)," redundant equations have been droped."$
483 terpri()$
484
485 ineq_:=ineq_bak$
486 if nopowersbak then algebraic(on nopowers)$
487 write(time()-cpu)/1000," s for pre-simplification"$terpri()$
488
489 return eqns
490end$
491
492%-------
493
494symbolic procedure input_consistency_test(afwlist,abwlist)$
495begin
496 if (length afwlist < 2) and
497    (length abwlist < 2) then
498    rederr "flist AND blist are both not assigned."$
499
500 afwlist:=cdr afwlist;
501 while afwlist and fixp car afwlist and car afwlist>=0 do afwlist:=cdr afwlist;
502 if afwlist then if not fixp car afwlist then
503 rederr "The fermionic weight list contains not only numbers!" else % <0
504 if not fixp reval algebraic max_deg then
505 rederr "If a fermionic weight is < 0 then max_deg must be set!"$
506
507 abwlist:=cdr abwlist;
508 while abwlist and fixp car abwlist and car abwlist>0 do abwlist:=cdr abwlist;
509 if abwlist then if not fixp car abwlist then
510 rederr "The bosonic weight list contains not only numbers!" else % <=0
511 if not fixp reval algebraic max_deg then
512 rederr "If a bosonic weight is < 1 then max_deg must be set!"
513end$
514
515%-------
516
517symbolic procedure crossprodu(seta,setb,wgt,endprod)$
518begin scalar g,k,setbcp;
519 if setb then
520 for each g in seta do <<
521  setbcp:=setb;
522  k:=car g + caar setbcp;
523  while k<=wgt do <<
524   if k=wgt then endprod:=cons(cons(k,append(cdr g,cdar setbcp)),endprod);
525   setbcp:=cdr setbcp;
526   if setbcp then k:=car g + caar setbcp
527             else k:=10000000
528  >>
529 >>$
530 return endprod
531end$
532
533%-------
534
535symbolic procedure rhs_term_list(N,allf,allb,maxwgt,linsub,forbid,verbose)$
536% generation of a list of all fermionic terms and a list of all
537% bosonic terms, all of weight maxwgt
538begin scalar g,h,j,k,minwgt,maxdeg,maxtermwgt,allv,newv,maxpow,maxwgtm2,
539             fprod,bprod,allprod,newprod,linonly,p,q,r,s$
540
541 % As weights can be negative, at first we compute the most negative weight
542 % of a partial product. For that we do not need to consider derivatives as
543 % these have all higher weights than the undifferentiated variables.
544
545 minwgt:=0;  maxdeg:=algebraic max_deg;  allv:=append(allf,allb);
546 for each g in allv do
547 if car g<0 then minwgt:=minwgt+maxdeg*car g;
548 maxtermwgt:=maxwgt-minwgt; % >=maxwgt as minwgt<=0
549
550 % Add to allv for every field its x-derivatives as high as possible
551 % and sort the entries by their weight
552
553 maxwgtm2:=maxtermwgt-2;
554 % newv:=nil; through initialization
555 for each g in allv do <<
556  while (car g <= maxwgtm2) and not member(reval {'DF,cadr g,'x},forbid) do <<
557   h:=2 + car g;
558   g:={h,{'DF,cadr g,'x}};    % not:  reval {'DF,cadr g,'x}
559   newv:=cons(g,newv)
560  >>
561 >>$
562 % allv:=idx_sort append(allv,newv);
563 allv:=append(allv,newv);
564
565 % Add to allv for each of its elements all possible combinations
566 % of d(i,..) derivatives
567 for j:=1:N do <<
568  newv:=nil;
569  for each g in allv do
570  if car g < maxtermwgt then
571  newv:=cons({add1 car g,{'D,j,cadr g}},newv);
572  % allv:=idx_sort append(allv,newv);
573  allv:=append(allv,newv);
574 >>$
575
576 % generation of all possible products of elements of allv
577 allprod:={{0,nil}};
578 if verbose then <<
579  write"maxwgt=",maxwgt$terpri()$
580  q:=length allv;
581  p:=0;
582  s:=0;
583 >>$
584 for each g in allv do <<
585  if verbose then <<
586   p:=add1 p;
587   write p,"(",q,"): ",g$terpri()
588  >>$
589
590  if null linsub then linonly:=nil else
591  if smemberl(linsub,cadr g) then linonly:=t
592                             else linonly:=nil;
593  maxpow:=if linonly or is_fermion cadr g then 1 else
594          if car g < 1 then maxdeg else 10000000;
595  newprod:=nil;
596  if verbose then r:=0;
597  for each h in allprod do
598  if null linonly or null cadr h then <<
599   % i.e. either the new factor cadr g is not linonly and
600   % can occur arbitrarily often in arbitrary products, or
601   % the factors do not contain a linonly factor yet
602   k:=car h + car g;
603   j:=1$
604   while (j <= maxpow) and (k <= maxtermwgt) do <<
605    h:=cons(k,cons(linonly or cadr h,cons(cadr g,cddr h)));
606    newprod:=cons(h,newprod);
607    k:=car h + car g; j:=add1 j
608   >>;
609   if verbose then r:=r+j-1
610  >>$
611  if verbose then  s:=s+r;
612  allprod:=nconc(allprod,newprod)
613 >>$
614
615 % sorting of the products into fermionic and bosonic products
616 allprod:=cdr allprod;  % drop of {{0,nil}}
617 while allprod do <<
618  g:=car allprod; allprod:=cdr allprod;
619  if (car g = maxwgt) and (null linsub or cadr g) then
620  if is_fermion cons('TIMES,cdr g) then
621  fprod:=cons(if null cdddr g then caddr g
622                              else cons('TIMES,cddr g),fprod)
623                                   else
624  bprod:=cons(if null cdddr g then caddr g
625                              else cons('TIMES,cddr g),bprod)
626 >>$
627
628 if verbose then <<write length fprod," fermionic and ",
629                         length bprod," bosonic monomials"$terpri()>>$
630 return {fprod,bprod}
631
632end$
633
634%-------
635
636symbolic procedure gen_eqn(awlist,rhslist,eqn_list,fl,fermionic,tname)$
637begin scalar h,w,rlcp,rhs,nf,print_bak,svar$
638 svar:=reval tname$
639 while awlist do <<
640  % an equation df(h,t)=... is generated for each car of element h of awlist
641  h:=car awlist; awlist:=cdr awlist;
642
643  % w is the weight of h
644  w:=car h; h:=cadr h;
645
646  % find the proper rhs for the weight w of h
647  rlcp:=rhslist;
648  while caar rlcp neq w do rlcp:=cdr rlcp;
649  rhs:=if fermionic then cadar  rlcp
650                    else caddar rlcp;
651  % add a constant factor to each term in rhs
652  print_bak:=print_; print_:=nil; % to avoid printing the new constants
653  rhs:=for each r in rhs collect <<nf   :=newfct(fname_,nil,nfct_)$
654                                   fl   :=cons(nf,fl);
655                                   nfct_:=add1 nfct_$
656                                   {'TIMES,nf,r}>>$
657  print_:=print_bak;
658  rhs:=if null rhs then 0
659                   else if cdr rhs then cons('PLUS,rhs)
660                                   else car rhs;
661  eqn_list:=cons({'EQUAL,{'DF,h,svar},rhs},eqn_list)
662 >>$
663
664 return {eqn_list,fl}
665end$
666
667%-------
668
669symbolic procedure sspol(N,sysfbl,afwlist,abwlist,difforder,tname,
670                         linsub,forbid,flip_par,verbose)$
671% returns {eqn_list,fl} where
672% eqn_list is the algebraic list of equations/symmetries and
673% fl is the algebraic list of undetermined constants
674% tname is the lhs differentiation variable
675begin scalar g,h,i,fwlist,bwlist,allf,allb,w_list,rhs_list,eqn_list,fl,awlist$
676
677 fwlist:=cdr reval afwlist$
678 bwlist:=cdr reval abwlist$
679
680 % Add to each field its weight as a prefix to create allf and allb
681 for i:=1:(length flist) do allf:=cons({nth(fwlist,i),nth(flist,i)},allf);
682 for i:=1:(length blist) do allb:=cons({nth(bwlist,i),nth(blist,i)},allb);
683
684 % Generation of the final form of the equations
685
686 % only those weights for which there is an eqn to generate
687 % making a list of all different weights
688 w_list:=union(fwlist,union(bwlist,nil))$
689
690 % generating the right-hand-side term-list for the different weights
691 % Because expressions are generated according to their weight and are
692 % partitioned only afterwards into fermionic and bosonic, we would not save
693 % computation when generating only fermionic expressions of specific weight
694 % and afterwards the bosonic expressions with specific weight
695 rhs_list:=for each w in w_list collect
696           cons(w,rhs_term_list(N,allf,allb,w+difforder,linsub,forbid,verbose))$
697                                   % new weight(t) based on weight(x)=2
698
699 % generating all bosonic equations
700 for each h in sysfbl do
701 if car h = 'b then << g:=allb; while g and cadar g neq h do g:=cdr g$
702  awlist:=cons(car g,awlist)
703 >>$
704 nfct_:=1$
705 % whether difforder is odd or even has no influence on parity (ferm/bos)
706 h:=gen_eqn(awlist,rhs_list,eqn_list,fl,flip_par,tname);
707
708 % generating all fermionic equations
709 awlist:=nil$
710 for each h in sysfbl do
711 if car h = 'f then << g:=allf; while g and cadar g neq h do g:=cdr g$
712  awlist:=cons(car g,awlist)
713 >>$
714 % whether difforder is odd or even has no influence on parity (ferm/bos)
715 h:=gen_eqn(awlist,rhs_list,car h,cadr h,not flip_par,tname);
716
717 return {cons('LIST,car h),cadr h}
718end$
719
720%-------
721
722symbolic procedure print_list_of_lists(li)$
723begin scalar k$
724 for each k in li do <<
725  write"{"$
726  while k do <<write car k$k:=cdr k$if k then write",">>$
727  write"}"$terpri()$
728 >>$
729end$
730
731%-------
732
733symbolic operator listine$
734symbolic procedure listine(N,sys,sym,fl,ineql,non_lin_test,save_lists)$
735% generation of conditions for the unknown coefficients
736begin scalar eqn,ineql,cnd,h,k,fb,flcp,nfl,evolist,rs$
737 fl:=cdr fl$
738 ineql:=for each h in cdr ineql collect cdr h$
739
740 evolist:=for each h in append(cdr sys,cdr sym) collect caddr h$
741 k:=gensym()$
742
743 % The online data base was computed with:
744 % % none of the equations and symmetries should vanish identically
745 % for each eqn in evolist do <<
746 %  if null (h:=smemberl(fl,eqn)) then h:={0};
747 %  ineql:=cons(h,ineql)
748 % >>$
749 % %write"ineql1=",ineql$terpri()$
750
751 % none of the equations should vanish identically
752 for each eqn in cdr sys do <<
753  if null (h:=smemberl(fl,caddr eqn)) then h:={0};
754  ineql:=cons(h,ineql)
755 >>$
756
757 % at least one of the symmetries must not vanish identically
758 h:=for each h in cdr sym collect caddr h$
759 if null (h:=smemberl(fl,h)) then h:={0};
760 ineql:=cons(h,ineql)$
761
762 % None of the evolution equations should involve only one field,
763 % the system should not be triangular.
764 if cdr fblist then  % if more than one field
765 for each eqn in cdr sys do <<
766  rs:=caddr eqn; % rhs of eqn
767  fb:=delete(cadr cadr eqn,fblist); % one of the functions fb must occur in rs
768  if pairp rs and (car rs = '!*SQ)
769  then rs:={'!*SQ, subsq(cadr rs, for each h in fb collect (h . 0)), t}
770  else for each h in fb do rs:=subst(0,h,rs);
771  h:=setdiff(smemberl(fl,caddr eqn),smemberl(fl,reval rs));
772  if null h then ineql:=cons({0},ineql)
773            else ineql:=cons(h,ineql)
774 >>$
775
776 % At least one evolutionary equation should be non-linear.
777 if non_lin_test then <<
778  nfl:=nil;   % we drop all constants from fl which are coefficient to
779              % a linear term
780  for each eqn in evolist do <<
781   flcp:=smemberl(fl,eqn);
782   for each h in fblist do eqn:=subst({'TIMES,k,h},h,eqn);
783   h:=smemberl(flcp,coeffn(reval eqn,k,1));
784   flcp:=setdiff(flcp,h);
785   nfl:=append(nfl,flcp)
786  >>$
787  if null nfl then nfl:={0};    % contradiction issued
788  ineql:=cons(nfl,ineql);
789 >>$
790
791 % at least one equation should contain a fermion field or a D_i derivative
792 if null flist then << % otherwise no chance to get a purely classical equation
793  flcp:=fl;   % we drop all constants from fl which are coefficient to
794              % a term without D(i,..) and without fermion
795  algebraic <<
796   drop_fd_rules:={d(~i,~j) => 0, f(~i) => 0};
797   let drop_fd_rules;
798  >>$
799  for each eqn in evolist do <<
800   h:=smemberl(fl,reval eqn);
801   flcp:=setdiff(flcp,h)
802  >>$
803  algebraic ( clearrules drop_fd_rules )$
804  if null flcp then flcp:={0};    % contradiction issued
805  ineql:=cons(flcp,ineql);
806  %write"ineql5=",ineql$terpri()$
807 >>$
808
809 %  N>1 makes only sense if all D_1,..,D_N occur which is required next
810 if N>1 then for h:=1:N do <<
811  flcp:=fl;   % we drop all constants from fl which are coefficient to
812              % a term without D(h,..)
813  algebraic <<k:=lisp h;drop_fd_rules:={d(k,~j) => 0}; let drop_fd_rules>>;
814  for each eqn in evolist do <<
815   k:=smemberl(fl,reval eqn);
816   flcp:=setdiff(flcp,k)
817  >>$
818  algebraic ( clearrules drop_fd_rules )$
819  if null flcp then flcp:={0};    % contradiction issued
820  ineql:=cons(flcp,ineql);
821 >>$
822
823 % dropping redundand conditions
824 repeat <<
825  rs:=ineql;
826  while rs do <<
827    k:=delete(car rs,ineql);
828    while k and not_included(car rs,car k) do k:=cdr k;
829    if k then << % i.e. cnd is fully included in car k
830     ineql:=delete(car k,ineql);
831     rs:=ineql
832    >>   else rs:=cdr rs
833  >>
834 >> until null rs$
835
836 if ineql then <<
837  terpri()$
838  write"From each of the following lists at least one element must be non-zero:"$
839  terpri()$terpri()$
840  print_list_of_lists(ineql)$
841  terpri()$
842 >>$
843 if save_lists then <<
844  out  "inelist"$
845  print_list_of_lists(ineql)$
846  shut "inelist"$
847 >>$
848
849 ineql:=for each cnd in ineql collect
850 if null cdr cnd then car cnd
851                else cons('PLUS,for each h in cnd collect {'TIMES,h,gensym()})$
852
853 return cons('LIST,ineql)
854
855end$
856
857%-------
858
859algebraic procedure power_parti(eqns)$
860begin scalar w,afblist,equa,ff,ex,a,p,parti,maxpow,h$
861 afblist:=lisp(cons('LIST,fblist))$
862 w:=lisp lin_test_const$
863 maxpow:=0$
864 h:=for each equa in eqns collect <<
865  ff:=part(equa,1);ex:=part(equa,2);
866  for each a in afblist do ex:=sub(a=a*w,ex)$
867  parti:=rest coeff(ex,w);
868  if hipow!*>maxpow then maxpow:=hipow!*;
869  ff=parti
870 >>;
871
872 % Filling up each element of h with zero's so that all have same length
873 if length h > 1 then
874 h:=for each equa in h collect part(equa,1)=<<
875  le:=length part(equa,2);
876  if le=maxpow then part(equa,2)
877               else <<ex:={};for g:=1:(maxpow-le) do ex:=cons(0,ex)$
878                      append(part(equa,2),ex)>>
879 >>$
880 return {maxpow,h}
881end$
882
883%-------
884
885symbolic procedure add_df_rules_to_D_rule(subl)$
886% subl is an algebraic list of rules
887% adds to each substitution rule D(i,A)=>.. all differential consequences
888if null cdr subl then subl else
889begin scalar g,h$
890 g:=cdr subl;
891 for each h in g do
892 if pairp h and car h='REPLACEBY and pairp cadr h and caadr h='D and
893    pairp caddr cadr h and ((caaddr cadr h='f) or (caaddr cadr h='b)) then <<
894  % In algebraic mode rules are not properly printed, i.e. without brackets
895  % around sums to be differentiated, so they may look wrong when being right
896  subl:=cons('LIST,union(
897      {{'REPLACEBY,{'DF,caddr cadr h,'x},reval {'D,cadr cadr h,caddr h}},
898       {'REPLACEBY,{'DF,caddr cadr h,'x,{'~,'n}},
899                   {'D,cadr cadr h,{'DF,caddr h,'x,{'DIFFERENCE,'n,1}}}},
900       {'REPLACEBY,{'D,cadr cadr h,{'DF,caddr cadr h,{'~,'x}}},
901                   reval {'DF,caddr h,'x}},
902       {'REPLACEBY,{'D,cadr cadr h,{'DF,caddr cadr h,{'~,'x},{'~,'n}}},
903                   {'DF,caddr h,'x,'n}},
904       {'REPLACEBY,{'D,cadr cadr h,{'DF,caddr cadr h,{'~,'t},{'~,'x},{'~,'n}}},
905                   {'DF,caddr h,'t,'x,'n}}},
906      cdr subl))
907 >>$
908 return subl
909end$
910
911%-------
912
913algebraic procedure sieve(inpu,wgl)$
914% inpu={sym,cl}
915% sym is an algebraic list of equations df(f(i),s)=.., df(b(j),s=..,
916% cl is an algebraic list of coefficients in sym
917% wgl is a list  {sw,afwlist,abwlist}
918% It returns a list (new symmetry, list of coefficients}
919begin scalar sym,cl,sw,fli,bli,su,n,h,newsy,to_drop,sy$
920 cl:=second inpu;
921 sw:=first wgl;
922 fli:=second wgl;
923 bli:=third wgl;
924
925 % creating a substitution rule
926 su:={};
927 n:=0$for each h in fli do <<n:=n+1;su:=cons(f(n)=f(n)*lin_test_const**h,su)>>$
928 n:=0$for each h in bli do <<n:=n+1;su:=cons(b(n)=b(n)*lin_test_const**h,su)>>$
929 sym:=sub(su,first inpu);
930 let lrule1; sym:=sym; clearrules lrule1;
931 let lrule2; sym:=sym; clearrules lrule2;
932
933 newsy:={};
934 to_drop:={}$
935 for each sy in sym do <<
936  coeff(lhs sy,lin_test_const);
937  n:=hipow!*+sw;
938  if n<0 then <<h:=lin_test_const**(-n)*rhs sy;
939                newsy:=cons(coeffn(h,lin_test_const,n),newsy)>>
940         else newsy:=cons(coeffn(rhs sy,lin_test_const,n),newsy);
941 >>;
942 su:={};
943 for each h in cl do if freeof(newsy,h) then <<
944  su:=cons(h=0,su);
945  to_drop:=cons(h,to_drop)
946 >>$
947
948 return {sub(su,first inpu),lisp cons('LIST,setdiff(cdr cl,cdr to_drop))}
949end$
950
951%-------
952
953symbolic procedure non_t_lhs_dvs(pdes)$
954% listing lhs's of equations and substitutions which are no t-derivatives
955begin scalar h,forbid$
956 for each h in pdes do
957 if pairp h and
958    ((car h = 'EQUAL) or (car h = 'REPLACE)) and
959    (pairp cadr h) and
960    ((caadr h = 'DF) or (caadr h = 'D) or (caadr h = 'F) or (caadr h = 'B)) and
961    freeof(cadr h,t) then <<
962  forbid:=cons(cadr h,forbid);
963  if caadr h = 'D  then forbid:=cons({'DF,caddr cadr h,'x},forbid) else
964  if caadr h = 'df then for g:=1:N_ do forbid:=cons({'D,g,cadr h},forbid)
965 >>$
966 return forbid
967end$
968
969%-------
970
971algebraic procedure ssym(N,tw,sw,afwlist,abwlist,eqnlist,fl,inelist,mode)$
972% uses the global algebraic variable x
973
974% N       .. the number of superfields theta_i
975% tw      .. 2 times the differential order of the equation = weight(t)
976% sw      .. 2 times the differential order of the symmetry = weight(s)
977% afwlist .. (algebraic) list of weights of the fermion fields
978%            f(1), f(2), ...
979% abwlist .. (algebraic) list of weights of the boson fields
980%            b(1), b(2), ...
981% The number of elements in the last two lists determines the
982% number of fermion and boson fields.
983% eqnlist .. in the nonlinear case a list of extra conditions on the
984%            undetermined coefficients where conditions a3=.. are executed
985%            instantly, these and any, expressions are added to equations
986%            when calling crack
987%            in the linear case the system in form of replacements and the
988%            linearized system in form of equations
989% fl      .. extra unknowns in eqnlist to be determined
990% inelist .. a list, each element of it is a list with at least one of
991%            its elements being non-zero
992% mode    .. list of flags:
993%            init : only initialization of global data
994%            plain_com : direct computation of the commutator (for tests)
995%            power_split_com : alternative power splitting of commutator,
996%                              (not if substitutions '='> are present)
997%            % const_coeff_in_eqn : have only constant coefficients in eqn.
998%            % const_coeff_in_sym : have only constant coefficients in sym.
999%            zerocoeff : all coefficients = 0 which are not listed in ine
1000%            % verbose : extra output for tracing and debugging
1001%            interactive : whether the solution process is interactive
1002%            % no_t_in_eqn : no explicit t-derivatives occur in equations
1003%            % no_t_in_sym : no explicit t-derivatives occur in symmetries
1004%            filter: extra homogeneity weights as given in the list"$
1005%                    hom_wei have to be satisfied by the symmetry"$
1006%            lin: find symmetry for all linear fields,
1007%                 if null spar (bosonic s) then
1008%                    bosonic   lin fields are b(nb+1)..b(2*nb)
1009%                    fermionic lin fields are f(nf+1)..b(2*nf)
1010%                                          else
1011%                    bosonic   lin fields are b(nb+1)..b(nf+nb)
1012%                    fermionic lin fields are f(nf+1)..b(nf+nb)
1013%            tpar: t/nil, whether time variable t changes parity"$
1014%            spar: t/nil, whether symmetry variable s changes parity"$
1015%            log : t/nil, whether files drvlist, evolist, unolist, inelist
1016%                         are generated to hold data, for example, for
1017%                         automatic web page generation
1018%
1019begin scalar g,h,k,cpu,gti,fbno,psys,psym,msysp,msymp,totpow,syspow,sympow,
1020             newcd,afblist,sublist,zerocoeff,fl_e,non_lin_test,do_ine_test,
1021             init,subl,subl2,sys,sym,ls,rs,linsub,filter,forbid,rhssyl,sycon,
1022             nw,w,np,p,nq,q,psycon,cn,lp,verbose,plain_com,power_split_com,
1023             msgbak,interactive,nfr,nfe,nbr,nbe,nf,nb,lhslist,save_lists$
1024 backup_reduce_flags()$
1025 lisp <<
1026  record_hist:=nil;
1027  if !*time then <<cpu:=time()$ gti:=gctime()>>$
1028
1029  eqnlist:=cdr eqnlist;
1030  fl     :=cdr fl;       if fl then homogen_:=nil else homogen_:=t;
1031  mode   :=cdr mode;
1032
1033  if member('init,mode) then init:=t else init:=nil;
1034  % if member('const_coeff_in_eqn,mode) then cc_eqn:=t else cc_eqn:=nil;
1035  % if member('const_coeff_in_sym,mode) then cc_sym:=t else cc_sym:=nil;
1036  if member('zerocoeff,mode) then zerocoeff:=t else zerocoeff:=nil;
1037  if member('plain_com,mode) then plain_com:=t else plain_com:=nil;
1038  if member('power_split_com,mode) then power_split_com:=t
1039                                   else power_split_com:=nil;
1040  % if member('term_split_com,mode) then term_split_com:=t
1041  %                                 else term_split_com:=nil;
1042  if member('verbose,mode) then verbose:=t else verbose:=nil;
1043  if member('interactive,mode) then interactive:=t else interactive:=nil;
1044  % if member('no_t_in_eqn,mode) then no_t_in_eqn:=t else no_t_in_eqn:=nil;
1045  % if member('no_t_in_sym,mode) then no_t_in_sym:=t else no_t_in_sym:=nil;
1046  if member('filter,mode) then filter:=t else filter:=nil;
1047  if member('tpar,mode) then algebraic(on  t_changes_parity)
1048                        else algebraic(off t_changes_parity);
1049  if member('spar,mode) then algebraic(on  s_changes_parity)
1050                        else algebraic(off s_changes_parity);
1051  if member('log,mode) then save_lists:=t$
1052
1053  if member('lin,mode) then << % @#@#
1054   nf:=length afwlist - 1;
1055   nb:=length abwlist - 1;
1056   if member('spar,mode) then
1057   if nf neq nb then <<write"In the case of spar #(boson fields)=#(fermion fields)",
1058                             " which is not the case."; forbid:=t>>
1059                else % to do more tests one would have to know/find the numbers
1060                     % of b/f fields before linearization
1061                         else << % null member('spar,mode)
1062    if null evenp nf then <<write"If no spar then #(fermion fields) needs to be even"$
1063                      forbid:=t>>;
1064    if null evenp nb then <<write"If no spar then #(boson fields) needs to be even"$
1065                      forbid:=t>>;
1066    for each g in eqnlist do if pairp g then lhslist:=cons(cadr g,lhslist);
1067    nf:=nf/2; nb:=nb/2;
1068    for each g in lhslist do
1069
1070if car g = 'df then
1071if caadr g = 'f then if ((cadadr g <= nf) and not member({'df,{'f,cadadr g + nf},caddr g},lhslist)) or
1072                        ((cadadr g  > nf) and not member({'df,{'f,cadadr g - nf},caddr g},lhslist)) then <<
1073                       write"The counterpart of ",g," is missing on a left hand side."$ forbid:=t>> else
1074                else
1075if caadr g = 'b then if ((cadadr g <= nb) and not member({'df,{'b,cadadr g + nb},caddr g},lhslist)) or
1076                        ((cadadr g  > nb) and not member({'df,{'b,cadadr g - nb},caddr g},lhslist)) then <<
1077                       write"The counterpart of ",g," is missing on a left hand side."$ forbid:=t>> else
1078                else
1079               else
1080if car g = 'd  then
1081if caaddr g = 'f then if ((car cdaddr g <= nf) and not member({'d,cadr g,{'f,car cdaddr g + nf}},lhslist)) or
1082                         ((car cdaddr g  > nf) and not member({'d,cadr g,{'f,car cdaddr g - nf}},lhslist)) then <<
1083                        write"The counterpart of ",g," is missing on a left hand side."$ forbid:=t>> else
1084                 else
1085if caaddr g = 'b then if ((car cdaddr g <= nb) and not member({'d,cadr g,{'b,car cdaddr g + nb}},lhslist)) or
1086                         ((car cdaddr g  > nb) and not member({'d,cadr g,{'b,car cdaddr g - nb}},lhslist)) then <<
1087                        write"The counterpart of ",g," is missing on a left hand side."$ forbid:=t>> else
1088                 else
1089               else
1090   >>; % of null member('spar,mode)
1091   if forbid then <<forbid:=nil;rederr "">>$
1092  >>$
1093
1094  % In the linearized system not half of the equations need to
1095  % be = and half => because there may be nonlocal potential variables
1096  % introduced through a => relation
1097  if nil then
1098  if member('lin,mode) then << % @#@#
1099   % determine nf and nb
1100   nfr:=0;  nbr:=0;  % number of f(i) and b(i) in replaceby
1101   nfe:=0;  nbe:=0;  % number of f(i) and b(i) in equal
1102
1103   for each g in eqnlist do
1104   if pairp g then
1105   if (car g='equal) then
1106   if not freeof(cadr g,'f) then nfe:=nfe+1 else
1107   if not freeof(cadr g,'b) then nbe:=nbe+1 else
1108                     else
1109   if (car g='replaceby) then
1110   if not freeof(cadr g,'f) then nfr:=nfr+1 else
1111   if not freeof(cadr g,'b) then nbr:=nbr+1 else
1112                         else $
1113
1114%check nfr+nfe = length afwlist - 1
1115%check nbr+nbe = length abwlist - 1
1116
1117write"nfr= ",nfr," nfe=",nfe,
1118     "nbr= ",nbr," nbe=",nbe$
1119
1120   if (    member('spar,mode) and (nfr=nbe) and (nbr=nfe)) or
1121      (not member('spar,mode) and (nfr=nfe) and (nbr=nbe)) then
1122                                                           else <<
1123    if member('spar,mode) then <<
1124     if nfr neq nbe then <<write "For spar
1125      the number of '=>' relations with a fermion on the lhs should be equal
1126      the number of '='  relations with a boson   on the lhs."$ terpri()
1127     >>$
1128     if nbr neq nfe then <<write "For spar
1129      the number of '=>' relations with a boson   on the lhs should be equal
1130      the number of '='  relations with a fermion on the lhs."$ terpri()
1131     >>
1132    >>                    else
1133
1134    if not member('spar,mode) then <<
1135     if nfr neq nfe then <<write "For missing spar
1136      the number of '=>' relations with a fermion on the lhs should be equal
1137      the number of '='  relations with a fermion on the lhs."$ terpri()
1138     >>$
1139     if nbr neq nbe then <<write "For missing spar
1140      the number of '=>' relations with a boson   on the lhs should be equal
1141      the number of '='  relations with a boson   on the lhs."$ terpri()
1142     >>
1143    >>$
1144    rederr ""
1145   >>
1146  >>;
1147
1148  % Initialization of fermionic and bosonic fields + their printing
1149  input_consistency_test(afwlist,abwlist)$
1150  N_:=N;  % to avoid printing the index of D in crack_out for N=1
1151  flist := for g:=1:(length afwlist-1) collect {'f,g};
1152  if flist and length flist=1 then put('f,'prifn,'myfpri)
1153                              else put('f,'prifn,nil);
1154  blist := for g:=1:(length abwlist-1) collect {'b,g};
1155  if blist and length blist=1 then put('b,'prifn,'myfpri)
1156                              else put('b,'prifn,nil);
1157  fblist:=append(flist,blist)
1158 >>$
1159
1160 afblist:=lisp(cons('LIST,fblist))$
1161 for each g in afblist do depend g,x,t$
1162
1163 lisp <<
1164
1165  % listing lhs's of equations and substitutions which are no t-derivatives
1166  forbid:=non_t_lhs_dvs(eqnlist)$
1167
1168  terpri()$
1169  write"#######  This is the case N",N,"f",length flist,"b",
1170       length blist,"t",tw,"s",sw,"w"$
1171  for each g in append(cdr afwlist,cdr abwlist) do write g;
1172  write".  #######"$terpri()$
1173  % specifying the names of constants in the equations
1174  if fname_list then <<fname_:=car fname_list;
1175                       fname_list:=cdr fname_list>>
1176                else   fname_:='p$
1177
1178  % generating the equations
1179  % Has an ansatz for the system been made?
1180  if eqnlist                                    and
1181     (pairp car eqnlist                       ) and
1182     (((caar eqnlist='EQUAL       ) and
1183       (pairp cadar eqnlist       ) and
1184       (caadar eqnlist = 'DF      )      ) or
1185      ((caar eqnlist='REPLACEBY   ) and
1186       (pairp cadar eqnlist       ) and
1187       ((caadar eqnlist = 'DF) or
1188        (caadar eqnlist = 'D )    )      )    ) then <<
1189   g:=nil;
1190   for each h in eqnlist do if pairp h then
1191   if pairp h and
1192      car h='EQUAL and
1193      pairp cadr h and
1194      caadr h='DF and
1195      reval caddr cadr h ='t then <<
1196    sys:=cons(h,sys);
1197    subl2:=cons({'REPLACEBY,cadr h,caddr h},subl2) % substitutions based on sys
1198   >>                        else
1199   if car h='REPLACEBY and
1200      pairp cadr h and
1201      ((caadr h='DF) or
1202       (caadr h='D)     ) then subl:=cons(h,subl)
1203                          else g:=cons(h,g)$
1204   eqnlist:=g$
1205
1206   sysfbl:=reverse for each h in sys collect cadadr h;
1207
1208   % simplifying the rhs of substitutions
1209   if verbose then <<write"Simplifying the substitutions."$terpri()>>$
1210   subl2:=cons('LIST,append(subl,subl2));
1211   algebraic(let subl2);
1212   subl:=for each h in subl collect {'REPLACEBY,cadr h,algebraic(lisp(caddr h))}}$
1213   subl:=cons('LIST,reverse subl);
1214   algebraic(clearrules subl2);
1215   subl2:=nil;
1216
1217   sys:=cons('LIST,reverse sys);
1218   fl_e:=fl$
1219   fl:=nil
1220  >>                                     else <<
1221   do_ine_test:=t$
1222   if null eqnlist then non_lin_test:=t$
1223   sysfbl:=append(flist,blist);
1224   subl2:=cons('LIST,append(subl,subl2));
1225   if verbose then <<write"Formulating rhs's of the system"$terpri()>>$
1226   h:=sspol(N,sysfbl,afwlist,abwlist,tw,algebraic t,linsub,
1227            forbid,!*t_changes_parity,verbose)$
1228   sys:=car h;
1229   fl_e:=cadr h;
1230   subl:={'LIST}
1231  >>$
1232
1233  algebraic<<nodepend f,s; nodepend b,s>>$
1234  for each h in fblist do algebraic(nodepend h,s);
1235  for each h in sysfbl do algebraic(  depend h,s);
1236
1237  % add to each substitution rule D(i,A)=>.. all differential consequences
1238  subl:=add_df_rules_to_D_rule subl$
1239
1240  % specifying the names of constants in the symmetries
1241  if fname_list then <<fname_:=car fname_list;
1242                       fname_list:=cdr fname_list>>
1243                else   fname_:='q$
1244  % generating the symmetry
1245  if verbose then <<write"Formulating rhs's of the symmetry"$terpri()>>$
1246  h:=sspol(N,sysfbl,afwlist,abwlist,sw,algebraic s,linsub,
1247           forbid,!*s_changes_parity,verbose)$
1248
1249  % applying a homogeneity filter
1250  if filter and pairp algebraic hom_wei
1251            and car (algebraic hom_wei) = 'LIST then <<
1252   if verbose then <<write"Applying an extra homogeneity filter"$terpri()>>$
1253   h:={'LIST,car h,cons('LIST,cadr h)};
1254   algebraic(for each g in hom_wei do h:=sieve(h,g))$
1255   h:={cadr h,cdaddr h} % converting to lisp list of two lisp lists
1256  >>$
1257  sym:=car h$
1258  flin_:=cadr h$  h:=nil;
1259
1260  % Extracting substitutions from eqnlist to be done instantly
1261  g:=eqnlist; eqnlist:=nil;
1262  for each h in g do
1263  if pairp h and car h = 'EQUAL then <<
1264   sublist:=cons(h,sublist)$
1265   eqnlist:=cons({'DIFFERENCE,cadr h,caddr h},eqnlist)
1266  >>                            else eqnlist:=cons(h,eqnlist)$
1267
1268  if zerocoeff then
1269  for each h in append(fl_e,flin_) do
1270  if freeof(inelist,h) then sublist:=cons({'EQUAL,h,0},sublist)
1271                       else fl:=cons(h,fl)
1272               else fl:=append(fl_e,flin_);
1273  fl_e:=nil$
1274  fl:=cons('LIST,fl)$
1275  sublist:=cons('LIST,sublist)$
1276  eqnlist:=cons('LIST,eqnlist)$
1277 >>$
1278
1279 if sublist then <<
1280  sys:=sub(sublist,sys)$
1281  sym:=sub(sublist,sym)
1282 >>$
1283 if lisp(save_lists) then <<
1284  off nat$
1285  out  "drvlist"$
1286  for each g in sys do write lhs g," := ",rhs g$
1287  for each g in sym do write lhs g," := ",rhs g$
1288  write"end$"$
1289  shut "drvlist"$
1290 >>$
1291 on nat$
1292 on dfprint$
1293 off noarg$
1294 if lisp(save_lists) then <<
1295  out  "evolist"$
1296  for each g in sys do write lhs g," := ",rhs g$
1297  lisp <<terpri()$
1298         write"</pre> with symmetries"$ terpri()$
1299         write"<pre>">>$
1300  for each g in sym do write lhs g," := ",rhs g$
1301  shut "evolist"$
1302  out  "unolist"$
1303  lisp <<listprint(reverse cdr fl)$terpri()>>$
1304  shut "unolist"$
1305 >>$
1306
1307 if interactive then <<off batch_mode$ lisp(print_:= 6 )>>
1308                else <<on  batch_mode$ lisp(print_:=nil)>>$
1309 lisp <<
1310  use_new_crackout:=t;
1311  %confirm_subst:=t;
1312  collect_sol:=nil;  % no collection of solutions, i.e. use of crackout()
1313  subst_1:=11;
1314  pdelimit_1:=400;
1315  max_gc_short:=4;
1316  max_gc_red_len:=4;
1317
1318  % instead of assigning proc_list here in the form
1319  %  proc_list_:='(
1320  %   to_do
1321  %   separation
1322  %   subst_level_0
1323  %   subst_level_04
1324  %   alg_length_reduction
1325  %   find_and_use_sub_systems22
1326  %   diff_length_reduction
1327  %   factorize_to_substitute
1328  %   subst_level_35
1329  %   decoupling
1330  %   factorize_any
1331  %   subst_level_4
1332  %   alg_solve_single
1333  %   stop_batch
1334  %  )
1335  % it is done in the assignment of old_history before the call of
1336  % CRACK, in order to drop the separation step after its initial
1337  % execution
1338
1339 >>$
1340
1341 if length sys=1 then write"The equation:" else write"The system:"$
1342 for each g in sys do write g$
1343 if subl neq {} then <<
1344  if length subl=1 then write"Extra condition:" else write"Extra conditions:"$
1345  off nat;
1346  for each g in subl do write g;
1347  on nat
1348 >>$
1349 if verbose then <<
1350  write"The symmetry:"$
1351  for each g in sym do write g$
1352 >>$
1353
1354 if do_ine_test then <<
1355  inelist:=listine(N,sys,sym,fl,inelist,non_lin_test,save_lists);
1356  g:=inelist$
1357  while g neq {} and first g neq 0 do g:=rest g$
1358  if g neq {} then return <<
1359   out "invalid"$
1360   write"SYSTEM & SYMMETRY DO NOT SATISFY MINIMAL REQUIREMENTS!"$
1361   shut "invalid"$
1362   write"SYSTEM & SYMMETRY DO NOT SATISFY MINIMAL REQUIREMENTS!"$
1363  >>
1364 >>$
1365
1366 % Formulating all conditions in stages
1367
1368 if subl neq {} then let subl$% this activates the extra conditions
1369 fbno:=lisp length sysfbl$
1370 if subl neq {} then power_split_com:=nil$ % because automatic substitutions from
1371                      % subs most likely change powers by introducing new factors
1372 if power_split_com then <<
1373  lisp (lin_test_const:=gensym())$
1374  psys:=power_parti(sys)$ msysp:=first psys; psys:=second psys;
1375  psym:=power_parti(sym)$ msymp:=first psym; psym:=second psym;
1376  for totpow:=2:(msysp+msymp) do << % generating all terms with total
1377                                    % power totpow of all elements of fblist
1378   lisp<<write"Generating all terms of total degree ",totpow$terpri()>>$
1379   % initialize all partial commutators to zero
1380   newcd:={}$ for g:=1:fbno do newcd:=cons(0,newcd)$
1381   % generate all pairings giving totpow
1382   for syspow:=1:msysp do
1383   for sympow:=1:msymp do
1384   if (syspow+sympow)=totpow then <<
1385    lisp <<
1386     write"  Pairing now degree ",syspow," terms from the system"$terpri()$
1387     write"         with degree ",sympow," terms from the symmetry"$terpri()>>$
1388
1389    subl2:=append(
1390     for each g in psys collect <<ls:=lhs g;rs:=part(rhs g,syspow);ls => rs>>,
1391     for each g in psym collect <<ls:=lhs g;rs:=part(rhs g,sympow);ls => rs>>)$
1392    let subl2;
1393
1394    % generation of partial commutators and adding them to newcd
1395    for g:=1:fbno do % for the g'th evol. equation + related symmetry equation
1396    newcd:=(part(newcd,g):=part(newcd,g)    + df(part(rhs part(psys,g),syspow),s) +
1397            (if (lisp t_changes_parity) and
1398                (lisp s_changes_parity) then   df(part(rhs part(psym,g),sympow),t)
1399                                        else (-df(part(rhs part(psym,g),sympow),t))));
1400    % clearing the t- and s-derivatives of all field variables
1401    clearrules subl2
1402
1403   >>$
1404   eqnlist:=append(newcd,eqnlist)
1405  >>$
1406  psys:=psym:=0
1407 >>$
1408
1409 % Assign t- and s-derivatives for all members of fblist for use in crack_out.
1410 subl2:=append(for each g in sys collect <<ls:=lhs g;rs:=rhs g; ls => rs>>,
1411               for each g in sym collect <<ls:=lhs g;rs:=rhs g; ls => rs>> )$
1412 let subl2; % This is the ansatz for the system and symmetry to be used
1413            % for the commutator.
1414 if verbose then <<
1415  write"sys = ",sys$
1416  write"sym = ",sym$
1417  write"df(rhs part(sys,1),s) = ",df(rhs part(sys,1),s)$
1418  write"df(rhs part(sym,1),t) = ",df(rhs part(sym,1),t)$
1419 >>$
1420
1421 if not power_split_com and plain_com then
1422 for g:=1:fbno do
1423 eqnlist:=cons(df(rhs part(sys,g),s) +
1424               (if (lisp t_changes_parity) and
1425                   (lisp s_changes_parity) then   df(rhs part(sym,g),t)
1426                                           else (-df(rhs part(sym,g),t))),
1427               eqnlist)
1428                                      else
1429 if not power_split_com then lisp <<
1430
1431  sym:=reval sym;
1432  eqnlist:=cdr eqnlist;   % to get rid of 'LIST
1433
1434  % for each symmetry h collect the list of terms of the rhs caddr h
1435  rhssyl:= for each h in cdr sym collect
1436           if pairp caddr h and caaddr h='PLUS then cdaddr h
1437                                               else {caddr h};
1438  nw:=0$
1439  for each w in sysfbl do <<
1440   if verbose then write"Integrability conditions are computed for ",w," :"$ terpri()$
1441   % For each function w of sysfbl generate the full condition
1442   nw:=add1 nw$ % the index of the function w in sysfbl
1443   sycon:=0$    % the condition for this function
1444   np:=0$       % the index of the symmetry of which one term is kept non-zero
1445
1446   % Initialize the s-derivative of all functions to zero
1447   msgbak:=!*msg; !*msg:=nil;
1448   for each h in sysfbl do algebraic <<clear df(lisp h,s);df(lisp h,s):=0>>$
1449   !*msg:=msgbak;
1450   cn:=0$                                 % counter of computed partial commutators
1451   lp:=for each p in rhssyl sum length p$ % counter of all partial commutators
1452
1453   psycon:=nil$   % a partial condition for this function
1454   for each p in rhssyl do << % for each rhs of a symmetry
1455
1456    % at first set previous s-deriv again to zero.
1457    msgbak:=!*msg; !*msg:=nil;
1458    if np>0 then <<h:=nth(sysfbl,np)$
1459                   algebraic <<clear df(h,s);df(h,s):=0>> >>$
1460    !*msg:=msgbak;
1461
1462    % Now we take the next rhs of which one term is not vanishing
1463    np:=add1 np;
1464    h:=nth(sysfbl,np)$
1465    nq:=0;      % the index of the term in p which is kept non-zero
1466    for each q in p do <<    % compute partial condition and add to psycon
1467     nq:=add1 nq$
1468     cn:=add1 cn$
1469     if verbose then algebraic write "   ",cn,"(",lp,")/",nw,"(",fbno,
1470                                     "). Currently non-vanishing: ",q$
1471     msgbak:=!*msg; !*msg:=nil;
1472     algebraic <<clear df(lisp h,s);
1473                 df(lisp h,s):=lisp nth(p,nq)>>$
1474     !*msg:=msgbak;
1475     psycon:=cons(algebraic <<
1476                   k:=df(w,t);g:=df(w,s);
1477                   df(k,s) + (if (lisp t_changes_parity) and
1478                                 (lisp s_changes_parity) then    df(g,t)
1479                                                         else (- df(g,t)))
1480                  >>,psycon)
1481    >>
1482   >>;
1483   algebraic <<
1484    sycon:=num lisp cons('PLUS,psycon);
1485    h:=length sycon;
1486    if verbose then
1487    if h=1 then write"The symmetry condition for ",w," has ",h," term."
1488           else write"The symmetry condition for ",w," has ",h," terms."
1489   >>$
1490   eqnlist:=cons(sycon,eqnlist)
1491  >>$
1492  eqnlist:=cons('LIST,eqnlist)$
1493
1494  % cleaning up assignments for s-derivatives
1495  msgbak:=!*msg; !*msg:=nil;
1496  for each h in sysfbl do algebraic clear df(lisp h,s);
1497  !*msg:=msgbak;
1498
1499 >>;
1500
1501 lisp<<msgbak:=!*msg; !*msg:=nil>>;
1502 clearrules subl2$
1503 clearrules subl$
1504 lisp (!*msg:=msgbak)$
1505
1506 subli:=subl$
1507 symansatz:=subl2$
1508
1509 sys:=sym:=0$
1510
1511 if not init then <<
1512  lisp(if !*time then <<terpri()$
1513   write "CPU time so far: ",   time() - cpu," ms  ",
1514         " GC time so far: ", gctime() - gti," ms"$
1515   terpri()
1516  >>);
1517
1518  lisp <<
1519
1520   % Listing the splitting `variables'
1521   h:='(t s x)$
1522   eqnlist:=split_simplify({eqnlist,{'LIST},fl,cons('LIST,h),t})$
1523   if print_ then <<write"Now crack is called."$terpri()>>$
1524   old_history:='(
1525
1526    cp (!*comma!* 1 3 45 11 8 20 27 30 47 21 38)
1527           % to specify an automatic solving strategy with priorities:
1528           %
1529           %  1 : Hot list of urgent steps
1530           %  3 : Substitution of <=2 terms, only fcts. of less vars., no cases
1531           % 45 : Substitution of <=8 terms in <=1000 terms,
1532           %      alg. expressions, no cases
1533           % 11 : Algebraic length reduction of equations
1534           % 53 : Find sub-systems with 2 flin_ functions % taken out,
1535           %                                % may generate too many equations
1536           % 27 : Length reducing decoupling steps
1537           %  8 : Factorization to subcases leading to substitutions
1538           % 20 : Substitution of <=1000 terms, no cases
1539           % 30 : Do one decoupling step
1540           % 47 : Any factorization
1541           % 21 : Substitution of <=1000 terms
1542           % 34 : Solving an algebraic equation. % taken out as too risky
1543           % 38 : Stop batch mode
1544    s      % to display a statistics for the original overdetermined system
1545    % w nil n "w"   % to store the system in ASCII for later web-page generation
1546    %               This has been disabled for the online demo.
1547    dp     % to disable parallellism
1548    % sb "0" % to backup the original system
1549    %               This has been disabled for the online demo.
1550   );
1551  >>;
1552  % on time$
1553  %  eqnlist:=lisp presimplify(eqnlist,fl)$   ############
1554  crack(eqnlist,inelist,fl,{})$ % ============= CRACK ==============
1555  lisp <<terpri()$
1556         write if null sol_list then "No" else length sol_list,
1557               if length sol_list=1 then " solution was" else " solutions were",
1558               " found."$ terpri();
1559       >>$
1560 >>
1561end$ % of ssym
1562
1563%-------
1564
1565algebraic procedure samelists(a,b)$
1566if a={} then if b={} then true
1567                     else nil
1568        else if b={} then nil
1569                     else
1570if first a neq first b then nil
1571                       else samelists(rest a,rest b)$
1572
1573%-------
1574
1575algebraic procedure crack_out(eqns,assigns,freef,ineq,sol_count)$
1576% to be used by crack called in ssym()
1577% lines marked with %#### have to be commented when ssym is called from script
1578% uses global variablse fblist,t,s,html_out
1579% assumes that t- and s-derivatives for all members of fblist are assigned
1580if lisp use_new_crackout then
1581begin scalar g,h,tm,ratbak,allbak,eqlist,eqcp1,eqcp2,afblist,ff,symlist,frsym,
1582             all,msgbak,sol_count2;
1583             %,natbak
1584 allbak:=lisp !*allfac$    on allfac$
1585 ratbak:=lisp !*rational$  on rational$
1586 afblist:=lisp(cons('LIST,sysfbl))$
1587 let symansatz$
1588% let ssrules;
1589 html_out:=nil$ % could be made a flag in ssym()
1590
1591 % At first printing the particular solutions correspponding to each
1592 % individual free parameter
1593
1594 eqlist:=for each h in afblist collect sub(assigns,df(h,t))$
1595 eqlist:=append(eqlist,for each h in afblist collect sub(assigns,df(h,s)));
1596 symlist:=for each h in afblist collect df(h,s);
1597 frsym:={}$
1598 for each ff in freef do if not freeof(symlist,ff) then frsym:=cons(ff,frsym);
1599 all:={}$
1600 sol_count2:=0;
1601
1602 write "All symmetries of solution ",sol_count," printed separately:"$
1603 for each ff in frsym do
1604 if not freeof(eqlist,ff) then <<   % which really appear             %####
1605
1606  eqcp1:=eqlist;
1607  for each h in frsym do if h neq ff then eqcp1:=sub(h=0,eqcp1);
1608
1609  % For the remaining constant ff we take a value that normalizes the numerical
1610  % coefficient of the first term of the first non-vanishing rhs of this
1611  % symmetry to one.
1612  g:=eqcp1;
1613  while (g neq {}) and freeof(first g,ff) do g:=rest g;
1614  if g neq {} then <<
1615   if my_freeof(den first g,x) then h:=den first g else h:=1$
1616   tm:=first g$
1617   h:=h/(numcoeff tm)$
1618   eqcp1:=sub(ff=h,eqcp1)
1619  >>$
1620
1621  % Has this solution already been found from other constants?
1622  g:=all$
1623  while (g neq {}) and not samelists(eqcp1,first g) do g:=rest g;
1624  if g={} then << % a new solution
1625   g:={};
1626   all:=cons(eqcp1,all);
1627   sol_count2:=sol_count2+1;
1628   eqcp2:=eqcp1$
1629
1630   lisp <<                                                             %####
1631    terpri()$
1632    if html_out then write"</pre>"$
1633    write "Solution ",sol_count,", symmetry ",sol_count2,". :"$ terpri()$terpri()$
1634    if length afblist=2 then write"The equation: "
1635                        else write"The system: "$  terpri()$
1636    if html_out then <<write"<pre>"$terpri()>>
1637
1638   >>$
1639   for each h in afblist do <<
1640    g:=df(h,t);
1641    msgbak:=!*msg; !*msg:=nil;
1642    clear df(h,t);
1643    !*msg:=msgbak;
1644    write df(h,t),"=",first eqcp1;
1645    df(h,t):=g;
1646    eqcp1:=rest eqcp1
1647   >>$
1648
1649   % if (lisp print_more) then
1650   <<
1651    lisp <<
1652     if html_out then write"</pre>"$
1653     terpri()$
1654     write"The symmetry: "$ terpri()$
1655     if html_out then write"<pre>"$
1656    >>$
1657    for each h in afblist do <<
1658     g:=df(h,s);
1659     msgbak:=!*msg; !*msg:=nil;
1660     clear df(h,s);
1661     !*msg:=msgbak;
1662     write df(h,s),"=",first eqcp1;
1663     df(h,s):=g;
1664     eqcp1:=rest eqcp1
1665    >>
1666   >>$
1667   lisp <<
1668    if html_out then write"</pre>"$
1669    terpri()$
1670    write"And now in machine readable form: "$ terpri()$
1671    if html_out then write"<p>"$
1672    if length afblist=2 then write"The equation: "
1673                        else write"The system: "$  terpri()$
1674    if html_out then write"<pre>"
1675   >>$
1676   off nat$
1677   for each h in afblist do <<
1678    write"df(",h,",t)=",first eqcp2;
1679    eqcp2:=rest eqcp2
1680   >>$
1681   on nat$
1682   % if (lisp print_more) then
1683   <<
1684    lisp <<
1685     if html_out then write"</pre>"$
1686     terpri()$
1687     write"The symmetry: "$terpri()$
1688     if html_out then write"<pre>"
1689    >>$
1690    off nat$
1691    for each h in afblist do <<
1692     write"df(",h,",s)=",first eqcp2;
1693     eqcp2:=rest eqcp2
1694    >>$
1695    on nat$
1696   >>
1697  >>   % of `g neq {}'
1698 >>$  % of `for each ff'
1699
1700 % If there are 2 or more symmetries then printing again the solution
1701 % with a linear combination of its symmetries:
1702
1703 if length(sfrsym) > 1 then <<
1704
1705  lisp <<                                                             %####
1706   terpri()$
1707   if html_out then write"</pre>"$
1708   terpri()$
1709   write "Again solution ",sol_count,", now with ALL its parameters (symmetries):"$ terpri()$
1710   if length afblist=2 then write"The equation: "
1711                       else write"The system: "$  terpri()$
1712   if html_out then <<write"<pre>"$terpri()>>
1713  >>$
1714  for each h in afblist do <<
1715   g:=df(h,t);
1716   msgbak:=!*msg; !*msg:=nil;
1717   clear df(h,t);
1718   !*msg:=msgbak;
1719   write df(h,t),"=",sub(assigns,g)$
1720   df(h,t):=g
1721  >>$
1722
1723  % if (lisp print_more) then
1724  <<
1725   lisp <<
1726    if html_out then write"</pre>"$
1727    terpri()$
1728    write"The symmetry: "$ terpri()$
1729    if html_out then write"<pre>"$
1730   >>$
1731   for each h in afblist do <<
1732    g:=df(h,s);
1733    msgbak:=!*msg; !*msg:=nil;
1734    clear df(h,s);
1735    !*msg:=msgbak;
1736    write df(h,s),"=",sub(assigns,g)$
1737    df(h,s):=g;
1738   >>
1739  >>$
1740
1741  % And now in machine readable form:
1742
1743  lisp <<
1744   if html_out then write"</pre>"$
1745   terpri()$
1746   write"And now in machine readable form: "$ terpri()$
1747   if html_out then write"<p>"$
1748   if length afblist=2 then write"The equation: "
1749                       else write"The system: "$  terpri()$
1750   if html_out then write"<pre>"
1751  >>$
1752  off nat$
1753  for each h in afblist do
1754  write"df(",h,",t)=",sub(assigns,df(h,t));
1755  on nat$
1756
1757  % if (lisp print_more) then
1758  <<
1759   lisp <<
1760    if html_out then write"</pre>"$
1761    terpri()$
1762    write"The symmetry: "$terpri()$
1763    if html_out then write"<pre>"
1764   >>$
1765   off nat$
1766   for each h in afblist do
1767   write"df(",h,",s)=",sub(assigns,df(h,s));
1768   on nat$
1769  >>$
1770 >>$  % if more than 1 symmetry
1771
1772 lisp <<write"-----------------------------------"$ terpri()>>$  %####
1773 clearrules symansatz$
1774 % ssrules should stay on
1775 if allbak=nil then off allfac$
1776 if ratbak=nil then off rational$
1777
1778end$ % of crack_out
1779
1780%-------
1781
1782symbolic procedure bigdpri l$
1783begin scalar dd,f;
1784 if not !*nat or !*fort then return 'failed;
1785 f:=cadr l; dd:=caddr l;
1786 prin2!* '!D;
1787 ycoord!* := ycoord!*-1;
1788 if ycoord!* < ymin!* then ymin!*:=ycoord!*;
1789 if N_ neq 1 then prin2!* f;
1790 ycoord!* := ycoord!*+1;
1791 if ycoord!* < ymin!* then ymin!*:=ycoord!*;
1792 if dd then
1793 << % prin2!* "(";    % to print brackets
1794    maprin dd;
1795    % prin2!* ")";
1796 >>;
1797 return l;
1798end$
1799
1800%-------
1801
1802put('!d,'prifn,'bigdpri)$
1803
1804put('bigdpri,'expt,'inbrackets)$
1805
1806% put('f,'prifn,'myfpri);  % done in ssym() if only one fermion field
1807% put('b,'prifn,'myfpri);  % done in ssym() if only one boson   field
1808
1809symbolic procedure myfpri(x)$
1810%  if not !*nat or !*fort or null !*no_index then  'failed
1811if not !*nat or !*fort then 'failed
1812                       else << prin2!* car x;x>>$
1813
1814%-------
1815%>>>>>>> The 2001 patch, now in the Reduce development version and in Reduce 3.8
1816
1817%    symbolic procedure noncomp u$
1818%       !*ncmp and noncomp1 u$
1819
1820%-------
1821
1822%   symbolic procedure noncomp1 u$
1823%      if null pairp u then nil
1824%       else if pairp car u then noncomfp u
1825%       else flagp(car u,'noncom) or noncomlistp cdr u$
1826
1827%-------
1828
1829%    symbolic procedure noncomlistp u$
1830%       pairp u and (noncomp1 car u or noncomlistp cdr u)$
1831
1832%-------
1833
1834%    symbolic procedure mchcomb(u,v,op)$
1835%       begin integer n;
1836%	  n := length u - length v +1;
1837%	  if n<1 then return nil
1838%	   else if n=1 then return mchsarg(u,v,op)
1839%	   else if not smemqlp(frlis!*,v) then return nil;
1840%	  return for each x in comb(u,n) join
1841%	       if null ncmp!* then mchsarg((op . x) . setdiff(u,x),v,op)
1842%		else (if null y then nil
1843%		       else mchsarg((op . x) . car y,
1844%				    if cdr y then reverse v else v,op))
1845%		where y = mchcomb2(x,u,nil,nil,nil)
1846%       end$
1847
1848%-------
1849
1850%    symbolic procedure mchcomb2(u,v,w,bool1,bool2)$
1851%       if null u then nconc(reversip w,v) . bool2
1852%	else if car u = car v
1853%	       then if noncomp car u then mchcomb2(cdr u,cdr v,w,t,bool2)
1854%		     else mchcomb2(cdr u,cdr v,w,bool1,bool2)
1855%	else if noncomp car u
1856%	 then if bool1 then nil
1857%	       else mchcomb2(u,cdr v,car v . w,t,if bool2 then bool2 else t)
1858%	else mchcomb2(u,cdr v,car v . w,bool1,bool2)$
1859
1860%-------
1861
1862symbolic procedure search_li3(l,carli)$
1863% Find all sublists of l which start with an element of carli (no nesting)
1864begin scalar b,res$
1865
1866 if pairp l then <<
1867  b:=carli;
1868  while b and null res do if car l = car b then res:=list l
1869                                           else b:=cdr b
1870 >>         else if l and member(l,carli) then res:=list l;
1871 if null res then
1872 while pairp l do <<
1873  if b:=search_li3(car l,carli) then res:=union(b,res);
1874  l:=cdr l
1875 >>$
1876 return res
1877end$
1878
1879%-------
1880
1881%symbolic procedure itercoeff(sf,splitvar)$                 %#################
1882%if not pairp sf then {!*f2a sf} else                       %#################
1883%if not memq(mvar sf,splitvar) then {!*f2a sf} else         %#################
1884%nconc(itercoeff(lc sf,splitvar),itercoeff(red sf,splitvar))$   %#################
1885
1886%-------
1887
1888symbolic procedure add_terms(h)$
1889% h is a lisp list of terms that each has to get a constant coefficient
1890% and the terms have to be added and returned
1891begin scalar fl,nf$
1892 h:=for each r in h collect <<nf   :=newfct(fname_,nil,nfct_)$
1893                              fl   :=cons(nf,fl);
1894                              nfct_:=add1 nfct_$
1895                              {'TIMES,nf,r}>>$
1896 return {if null h then 0 else
1897         if null cdr h then reval car h else reval cons('PLUS,h),
1898         fl}
1899end$
1900
1901%-------
1902
1903algebraic procedure ssconlaw(N,tw,cw,afwlist,abwlist,pdes,fermi)$
1904 % N       .. the number of superfields theta_i
1905 % tw      .. 2 times the differential order of the equation = weight(t)
1906 % cw      .. 2 times the differential order of the conservation law
1907 % afwlist .. (algebraic) list of weights of the fermion fields
1908 %            f(1), f(2), ...
1909 % abwlist .. (algebraic) list of weights of the boson fields
1910 %            b(1), b(2), ...
1911 % The number of elements in the last two lists determines the
1912 % number of fermion and boson fields.
1913 % pdes    .. an algebraic list of the pde(s) for which a conservation
1914 %            law is to be found, for example
1915 %            {df(f(1),t)=df(f(1),x)*b(1)*p9,
1916 %             df(b(1),t)=b(1)**3*p3 + d(1,f(1))**2*p2 +
1917 %                        df(b(1),x)*b(1)*p9 + df(f(1),x)*f(1)*p4}$
1918 % fermi   .. if =t then conserved flow is fermionic, if =nil then bosonic
1919begin
1920 scalar cpu,gti,sw,fwlist,bwlist,allf,allb,g,h,k,l,m,Pt,Px,Pdl,fl,nf,fbno,eqn,
1921        Ptcp,Pxcp,Pdlcp,sol,spezsol,Ql,terms,te,minu,factors,prefac,fa,dli,
1922        dfli,inte,delta,dlirevcp,tr_cl,pdes2,pdes3,fno,bno,syli,deno,allcl,
1923        nontriv,forbid,Pt1,Px1,Pdl1,verbose$
1924
1925% tr_cl:=t$
1926
1927 backup_reduce_flags()$
1928
1929 lisp <<
1930  record_hist:=nil;homogen_:=t;
1931  if !*time then <<cpu:=time()$ gti:=gctime()>>$
1932  input_consistency_test(afwlist,abwlist)$
1933  N_:=N;  % to avoid printing the index of D in crack_out for N=1
1934  fno:=length afwlist-1;
1935  bno:=length abwlist-1;
1936  flist := for g:=1:fno collect {'f,g};
1937  if flist and length flist=1 and tr_cl=nil then put('f,'prifn,'myfpri)
1938                                            else put('f,'prifn,nil);
1939  blist := for g:=1:bno collect {'b,g};
1940  if blist and length blist=1 and tr_cl=nil then put('b,'prifn,'myfpri)
1941                                            else put('b,'prifn,nil);
1942  fblist:=append(flist,blist)
1943 >>$
1944 afblist:=lisp(cons('LIST,fblist))$
1945 for each g in afblist do depend g,x,t$
1946
1947 lisp <<
1948
1949  % creating more substitutions based on pdes
1950  for each g in cdr pdes do <<
1951   if is_fermion cadr g then <<fno:=add1 fno; h:={'f,fno}>>
1952                        else <<bno:=add1 bno; h:={'b,bno}>>;
1953   algebraic(depend lisp(h),x,t);
1954   pdes2:=cons({'REPLACEBY,cadr g,{'PLUS,caddr g,h}},pdes2);
1955   pdes3:=cons({h,0,{'DIFFERENCE,cadr g,caddr g}},pdes3);
1956   syli:=cons(h,syli)
1957  >>$
1958  % pde2 is a substitution list introducing a symbol for each equation
1959  % pde3 is an assoc list with an entry for each equation, eg.
1960  % { {symbol_for_eqn, Q(1), df(f(1),t)-...},
1961  %   {symbol_for_eqn, Q(2), df(b(1),t)-...}, ... }
1962
1963  % listing lhs's of equations and substitutions which are no t-derivatives
1964  forbid:=non_t_lhs_dvs(cdr pdes)$
1965
1966  pdes2:=cons('LIST,pdes2);
1967  if tr_cl then <<
1968   algebraic(write"pdes2=",lisp pdes2)$
1969   write"pdes3="$prettyprint pdes3$ terpri()$
1970   write"syli=",syli$terpri()
1971  >>$
1972
1973  % specifying the names of constants in the equations
1974  if fname_list then <<fname_:=car fname_list;
1975                       fname_list:=cdr fname_list>>
1976                else   fname_:='r$
1977
1978  fwlist:=cdr reval afwlist$
1979  bwlist:=cdr reval abwlist$
1980
1981  % Add to each field its weight as a prefix to create allf and allb
1982  for i:=1:(length flist) do allf:=cons({nth(fwlist,i),nth(flist,i)},allf);
1983  for i:=1:(length blist) do allb:=cons({nth(bwlist,i),nth(blist,i)},allb);
1984
1985  nfct_:=1$ print_:=nil$
1986
1987  h:=rhs_term_list(N,allf,allb,cw-tw,nil,forbid,verbose);
1988  if fermi then h:= car h else h:=cadr h;
1989  h:=add_terms(h)$ Pt:=car h$ fl:=append(cadr h,fl)$
1990  eqn:={{'DF,Pt,'t}}$
1991
1992  if N=0 then <<
1993   h:=rhs_term_list(N,allf,allb,cw-2,nil,forbid,verbose);
1994   if fermi then h:=car h else h:=cadr h;
1995   h:=add_terms(h)$ Px:=car h$ fl:=append(cadr h,fl)$
1996   eqn:=cons({'DF,Px,'x},eqn)$
1997   Pdl:=nil$
1998  >>     else <<
1999   Px:=0;    % which is no loss of generality
2000
2001   h:=rhs_term_list(N,allf,allb,cw-1,nil,forbid,verbose);
2002   if !*t_changes_parity                  then
2003   if fermi then h:= car h else h:=cadr h else
2004   if fermi then h:=cadr h else h:= car h;
2005   Pdl:=nil$
2006   for g:=N step -1 until 1 do << % step -1 is crucial for backintegration
2007    k:=add_terms(h)$
2008    Pdl:=cons({'LIST,g,car k},Pdl)$
2009    fl:=append(cadr k,fl)$
2010    eqn:=cons({'D,g,car k},eqn)$
2011   >>
2012  >>$
2013
2014  eqn:=cons('PLUS,eqn)$
2015  Pdl:=cons('LIST,Pdl)$
2016  fl:=cons('LIST,fl)
2017
2018 >>$
2019
2020 if tr_cl then <<
2021  write"Ansatz for Pt: ",Pt$
2022  write"Ansatz for Px: ",Px$
2023  write"Ansatz for Pdl: ",Pdl$
2024  write"list of unknown coefficients: ",fl$
2025 >>$
2026
2027 if fl={} then return <<
2028  write "No valid ansatz in this case."$
2029  nil
2030 >>$
2031
2032 let pdes; % initializing all substitutions based on equations to generate eqn.
2033 eqn:=eqn; % only now after `let pdes'
2034 clearrules pdes;
2035
2036 lisp << % Formulating the conservation law conditions
2037
2038  % Listing the splitting `variables'
2039  h:=search_li3(reval eqn,{'DF,'F,'B,'D,'X,'T});
2040  if !*time then <<terpri()$
2041   write "CPU time so far: ",   time() - cpu," ms  ",
2042         " GC time so far: ", gctime() - gti," ms"$
2043   terpri()$
2044  >>$
2045
2046  % Splitting
2047  k:=setkorder h$
2048  eqn:=cons('LIST,
2049            for each g in itercoeff(if not pairp eqn then eqn else
2050                                    if car eqn = '!*sq then reorder numr cadr eqn
2051                                                       else numr simp eqn,h)
2052            collect {'!*SQ,(g . 1),t})$
2053  setkorder k$
2054
2055  !!arbint:=0$
2056  terpri()$
2057  write length cdr eqn," conditions for ",length cdr fl," unknowns"$
2058  terpri()
2059
2060 >>$
2061
2062 % Solving the conservation law conditions
2063 sol:=first solve(eqn,fl)$
2064
2065 lisp <<
2066  % terpri()$
2067  write if 0=!!arbint then "No" else !!arbint,
2068        if 0 neq !!arbint then " (possibly trivial) " else " ",
2069        if fermi then "fermionic" else "bosonic"," conservation ",
2070        if 1=!!arbint then "law" else "laws"," of weight ",cw,
2071        if 0=!!arbint then "." else ":"$
2072  terpri()$
2073
2074 >>$
2075
2076 nontriv:=0$
2077 for g:=1:lisp !!arbint do <<
2078  spezsol:=sub(arbcomplex(g)=1,sol);
2079  for h:=1:lisp !!arbint do
2080  if h neq g then spezsol:=sub(arbcomplex(h)=0,spezsol);
2081
2082  Ptcp :=sub(spezsol,Pt)$
2083  Pxcp :=sub(spezsol,Px)$
2084  Pdlcp:=sub(spezsol,Pdl)$
2085
2086  % Counting how many P's are nonzero
2087  if Ptcp=0 then k:=0 else k:=1;
2088  if Pxcp neq 0 then k:=k+1;
2089  h:=Pdlcp$
2090  while h neq {} do << if second first h neq 0 then k:=k+1;
2091                       h:=rest h >>$
2092
2093  % Drop conservation laws with only one non-vanishing component
2094  if k=0 then write "We drop a conservation law with Pt, Px, Pd[i] ",
2095                     "all being zero"
2096         else
2097  if k=1 then write "We drop a conservation law with only ",
2098                     if Ptcp neq 0 then "Pt" else "one Pd[i]"," nonzero."
2099         else <<
2100
2101   Pt1:=Ptcp$
2102   Px1:=Pxcp$
2103   Pdl1:=Pdlcp$
2104
2105   if tr_cl then <<
2106    write ">>>>> ",g,". Conservation law: "$
2107    write"Pt = ",Ptcp$
2108    write"Px = ",Pxcp$
2109    h:=Pdlcp$
2110    while h neq {} do << write"Pd[",first first h,"] = ",second first h;
2111                         h:=rest h >>
2112   >>$
2113
2114   % Determination of the characteristic functions
2115   eqn:=df(Ptcp,t)+df(Pxcp,x)+
2116        for l:=1:N sum <<m:=first part(Pdlcp,l);D(m,second part(Pdlcp,l))>>$
2117   if tr_cl then write"rhs = ",eqn$
2118
2119   if eqn=0 then % write "This conservation law is trivial!"
2120            else lisp <<
2121
2122    Ql:=pdes3$
2123
2124    pdes2:=add_df_rules_to_D_rule pdes2$
2125    algebraic(let pdes2);
2126    eqn:=reval eqn;
2127    algebraic(clearrules pdes2);
2128
2129    if tr_cl then <<write"rhs in terms of equations = "$
2130     % by now eqn should not involve any lhs or derivative
2131     % of any lhs of any equation in pdes.
2132     prettyprint eqn;
2133     terpri()
2134    >>$
2135
2136    % taking care of quotients
2137    deno:=nil;
2138    if pairp eqn and car eqn = 'QUOTIENT then
2139    if is_const caddr eqn then <<deno:=caddr eqn; eqn:=cadr eqn>>
2140                          else rederr" Something is wrong!! (0) *****"$
2141    % converting to a list of terms
2142    if not pairp eqn or car eqn neq 'PLUS then terms:={eqn} else terms:=cdr eqn;
2143    for each te in terms do <<
2144
2145     % taking care of the minus sign
2146     minu:=nil; if pairp te and car te = 'MINUS then <<minu:=t; te:=cadr te>>$
2147     % checking each factor of which only one should have a t-derivative
2148     if not pairp te or car te neq 'TIMES then factors:={te}
2149                                          else factors:=cdr te;
2150     prefac:=nil;
2151     while factors do << % search until the first factor including an
2152                         % element of syli is found
2153      fa:=car factors;  factors:=cdr factors;
2154      if freeoflist(fa,syli) then prefac:=cons(fa,prefac)
2155                             else <<
2156       if tr_cl then <<
2157        write" Next term to be partially integrated:";terpri()$
2158        prettyprint prefac$ terpri()$ write" times "$terpri()$
2159        prettyprint fa$     terpri()$ write" times "$terpri()$
2160        prettyprint factors$terpri()$
2161        write"minu=",minu$  terpri()
2162       >>$
2163
2164       % In the case of an exponent, a single factor is good enough for us
2165       if pairp fa and car fa='EXPT then <<
2166        factors:=cons(if caddr fa=2 then cadr fa
2167                                    else {'EXPT,cadr fa,sub1 caddr fa},factors);
2168        fa:=cadr fa
2169       >>$
2170
2171       if null factors then factors:=1 else
2172       if null cdr factors then factors:=car factors
2173                           else factors:=cons('TIMES,factors);
2174
2175       if null prefac then prefac:=1 else
2176       if null cdr prefac then prefac:=car prefac
2177                          else prefac:=cons('TIMES,reverse prefac);
2178
2179       % Merging prefac and factors, i.e. moving fa to be last factor
2180       if is_fermion fa and is_fermion factors then minu:=not minu;
2181       prefac:={'TIMES,prefac,factors};
2182       factors:=nil; % i.e. stop of while loop
2183       if minu then prefac:={'MINUS,prefac}; minu:=nil;
2184
2185       dli:=nil;
2186       while pairp fa and car fa='D do <<dli:=cons(cadr fa,dli);fa:=caddr fa>>$
2187       dli:=reverse dli;
2188
2189       if pairp fa and car fa = 'DF then <<
2190        delta:=cadr fa;
2191        dfli:=cddr fa;
2192        h:=nil;
2193        while dfli do <<
2194         if fixp car dfli then for k:=2:(car dfli) do h:=cons(car h,h)
2195                          else h:=cons(car dfli,h);
2196         dfli:=cdr dfli
2197        >>$
2198        dfli:=reverse h
2199       >>          else <<delta:=fa;dfli:=nil>>;
2200       k:=assoc(delta,Ql);  % k has form: {symb for eqn, Q_f(i), df(f(i),t)-rhs}
2201       Ql:=delete(k,Ql);
2202       delta:=caddr k;
2203       if deno then delta:={'QUOTIENT,delta,deno};
2204
2205       if tr_cl then <<write"delta=",delta,"  dfli=",dfli$ terpri()>>$
2206
2207       % At first partial integration of the D-derivatives
2208       while dli do <<
2209        h:=car dli; dli:=cdr dli; % h in index of supersym. generator
2210        if is_fermion prefac then prefac:={'MINUS,prefac};
2211        inte:=delta;
2212        if dfli then inte:=cons('DF,cons(inte,dfli));
2213
2214        dlirevcp:=reverse dli;
2215        while dlirevcp do <<
2216         inte:={'D,car dlirevcp,inte};
2217         dlirevcp:=cdr dlirevcp
2218        >>$
2219        algebraic(if tr_cl then write"Pdlcp before = ",Pdlcp)$
2220        algebraic(Pdlcp:=(part(Pdlcp,h):={first part(Pdlcp,h),
2221                                          second part(Pdlcp,h) -
2222                                          lisp({'TIMES,prefac,inte})}))$
2223        algebraic(if tr_cl then write"Pdlcp after = ",Pdlcp)$
2224        prefac:=reval {'MINUS,{'D,h,prefac}}
2225       >>$
2226       % Then partial integration of the x,t-derivatives
2227       while dfli and not zerop prefac do <<
2228        if tr_cl then <<write"Partial ",car dfli,"-integration:"$ terpri()>>$
2229        inte:=delta;
2230        if cdr dfli then inte:=cons('DF,cons(inte,cdr dfli));
2231        if car dfli='x then <<
2232         algebraic(if tr_cl then write"Pxcp before = ",Pxcp)$
2233         Pxcp:={'DIFFERENCE,Pxcp,{'TIMES,prefac,inte}};
2234         algebraic(if tr_cl then write"Pxcp after = ",Pxcp)$
2235        >>             else <<
2236         algebraic(if tr_cl then write"Ptcp before = ",Ptcp)$
2237         Ptcp:={'DIFFERENCE,Ptcp,{'TIMES,prefac,inte}};
2238         algebraic(if tr_cl then write"Ptcp after = ",Ptcp)$
2239        >>$
2240
2241        prefac:=reval {'MINUS,{'DF,prefac,car dfli}};
2242        dfli:=cdr dfli
2243       >>$
2244
2245       if deno then prefac:={'QUOTIENT,prefac,deno};
2246       Ql:=cons({car k,{'PLUS,prefac,cadr k},caddr k},Ql);
2247      >>  % of not freeoflist(fa,syli)
2248     >>  % while factors do
2249    >>; % for each te in terms do
2250
2251    % The integrating factors as well as the conserved densities may now
2252    % involve the symbols, like f(4), b(3) introduced for the equations.
2253    % These symbols have to be replaced by the equations now.
2254    nf:=cons('LIST,for each g in pdes3 collect {'EQUAL,car g,caddr g}) ;
2255    Ql:=for each h in Ql collect {car h,algebraic(sub(nf,lisp cadr h)),caddr h};
2256
2257    % Is the conservation law trivial?
2258    h:=Ql$
2259    while h and zerop cadar h do h:=cdr h;
2260    if tr_cl then <<
2261     if null h then write"CL is TRIVIAL!"
2262               else write"CL is GENUINE!"$
2263     terpri()$
2264    >>$
2265
2266    if h then algebraic << % CL is not trivial
2267
2268     Ptcp :=sub(nf,Ptcp)$
2269     Pxcp :=sub(nf,Pxcp)$
2270     if N>0 then <<
2271      % We use D_x(Pxcp) = D(m,D(m,Pxcp)), set Pxcp=0 and add D(m,Pxcp) to
2272      % Pd(m) where m is one index in 1..N (probably 1).
2273      m:=first first Pdlcp;
2274      if tr_cl then write"We add Pxcp=",Pxcp," to Pd[",m,"]."$
2275      Pdlcp:=cons({m, D(m,Pxcp)+second first Pdlcp},
2276                  rest Pdlcp);
2277      Pxcp:=0$
2278      Pdlcp:=sub(nf,Pdlcp)$
2279     >>$
2280
2281     % Test of correctness
2282     lisp(h:=if length Ql=1 then {'TIMES,cadar Ql,caddar Ql}
2283                            else cons('PLUS,for each h in Ql collect
2284                                            {'TIMES,cadr h,caddr h}));
2285
2286     k:=df(Ptcp,t)+df(Pxcp,x) - (lisp h) +
2287        for g:=1:N sum <<m:=first part(Pdlcp,g);D(m,second part(Pdlcp,g))>>$
2288
2289     % Making Ql an algebraic list
2290     lisp (Ql:=cons('LIST,for each h in Ql collect {'LIST,caddr h,cadr h}));
2291
2292     % Is conservation law new? --> compare with CL in allcl.
2293     % structure of allcl: {cl1,cl2,...} where each cli is a conservation law
2294     % and has the form {Ql,Pdlcp,Pxcp,Ptcp,no_of_terms_in_all_P}
2295     %    h:=allcl$
2296     %    while h neq {} do <<
2297     %
2298     %    >>$
2299     %. . . . . . .
2300
2301     if t then << % Print new conservation law
2302      lisp <<
2303       nontriv:=add1 nontriv;
2304       terpri()$
2305       write ">>>>> ",nontriv,". Non-trivial conservation law: "$
2306       terpri()$terpri()$
2307       write"At first in the form of a conserved current vanishing mod eqn.s:"$
2308       terpri()$terpri()$
2309      >>$
2310      write"Pt = ",Pt1$
2311      write"Px = ",Px1$
2312      % Printing of Pdl1
2313      h:=Pdl1$
2314      while h neq {} do <<write"Pd[",first first h,"] = ",second first h;
2315                          h:=rest h >>$
2316
2317      write"----- and in machine readable form:"$
2318      off nat$
2319      write"Pt = ",Pt1$
2320      write"Px = ",Px1$
2321      h:=Pdl1$
2322      while h neq {} do << write"Pd[",first first h,"] = ",second first h;
2323                           h:=rest h >>$
2324      on nat;
2325      write"Now as conserved current with characteristic functions:"$
2326
2327      write"Pt = ",Ptcp$
2328      write"Px = ",Pxcp$
2329      % Printing of Pdlcp
2330      h:=Pdlcp$
2331      while h neq {} do <<write"Pd[",first first h,"] = ",second first h;
2332                          h:=rest h >>$
2333      for each h in Ql do write"Q[",first h,"] = ",second h;
2334      write"----- and in machine readable form:"$
2335      off nat$
2336      write"Pt = ",Ptcp$
2337      write"Px = ",Pxcp$
2338      h:=Pdlcp$
2339      while h neq {} do << write"Pd[",first first h,"] = ",second first h;
2340                           h:=rest h >>$
2341      for each h in Ql do write"Q[",first h,"] = ",second h;
2342      on nat;
2343     >>;
2344
2345     % We print contradiction message after the values to be able to check
2346     if k neq 0 then <<
2347      write"***** ERROR: "$
2348      write" 0 <> ",k$
2349      rederr" A test shows a contradiction! *****"$
2350%      write" A test shows a contradiction! *****"$    %@@@@@@@@@@@@@@@@
2351     >>
2352
2353    >>   %  there is a non-zero Q after substitution of equations
2354   >>   %  divergence does not vanish identically
2355  >>   %  if k>1 then <<  (more than one non-vanishing component)
2356 >>$  %  for g:=1:lisp !!arbint do <<
2357 %#1#
2358end$
2359
2360%-------
2361
2362algebraic procedure ssconl(N,tw,mincw,maxcw,afwlist,abwlist,pdes)$
2363<<lisp <<
2364  terpri()$
2365  write"##### ssconl: This is the case N",N,"f",length cdr afwlist,"b",
2366       length cdr abwlist,"t",tw,"w"$
2367  for each g in append(cdr afwlist,cdr abwlist) do write g;
2368  write". #####"$terpri()$
2369 >>$
2370 write"The ",if length pdes=1 then "equation" else "system",
2371      " to be investigated:"$
2372 for each h in pdes do write h$
2373 pdes:= for each h in pdes collect
2374        lisp {'REPLACEBY,cadr algebraic h,caddr algebraic h}$
2375 lisp <<
2376  terpri()$
2377  write"Each CL has the form:  "$terpri()$
2378  if N=0 then write"    Dt(Pt) + Dx(Px) = Q1*PDE1 + .."
2379         else write"    Dt(Pt) + sum_i Di(Pd(i)) = Q1*PDE1 + .."$
2380  terpri()$
2381  write"where, e.g. PDE1 is (left hand side) - (right hand side) of PDE 1."$
2382  terpri()
2383 >>$
2384 for h:=mincw:maxcw do <<
2385  write"NEXT: FERMIONIC CONSERVATION LAWS OF WEIGHT ",h$
2386  ssconlaw(N,tw,h,afwlist,abwlist,pdes,t)$
2387  write"NEXT: BOSONIC CONSERVATION LAWS OF WEIGHT ",h$
2388  ssconlaw(N,tw,h,afwlist,abwlist,pdes,nil)
2389 >>
2390>>$
2391
2392%-------
2393
2394symbolic procedure lofl(n,mwt)$
2395% produces lists with n elements, each having a value 0...mwt
2396if n=1 then for i:=0:mwt collect {i,{i}} else
2397for each l in lofl(n-1,mwt) join
2398for j:=0:(mwt-car l) collect {j+car l,cons(j,cadr l)}$
2399
2400%-------
2401
2402symbolic procedure wgtof(h,ali,wl)$
2403% subroutine for FindSSWeights below, determines weight of a kernel
2404begin scalar w,wt,m$
2405 return
2406 if w:=assoc(h,ali) then {cdr w,ali,wl} else
2407 if pairp h and car h='D then <<
2408  w:=wgtof(caddr h,ali,wl);
2409  cons({'PLUS,car w,1},cdr w)
2410 >>                      else
2411 if pairp h and car h='DF then <<
2412  w:=wgtof(cadr h,ali,wl);
2413  wt:=car w;
2414  ali:=cadr w;
2415  wl:=caddr w;
2416  h:=cddr h; % h is now a list of differentiations
2417  while h do <<
2418   w:=wgtof(car h,ali,wl);
2419   h:=cdr h;
2420   ali:=cadr w;
2421   wl:=caddr w;
2422   if null h or not fixp car h then m:=1
2423                               else <<m:=car h; h:=cdr h>>;
2424   wt:={'DIFFERENCE,wt,{'TIMES,car w,m}}
2425  >>$
2426  {wt,ali,wl}
2427 >>                       else <<
2428  wt:=mkid('w_,if pairp h then gensym()
2429                          else h       );
2430  ali:=cons((h . wt), ali);
2431  wl:=cons(wt,wl);
2432  {wt,ali,wl}
2433 >>
2434end$
2435
2436%-------
2437
2438symbolic operator FindSSWeights$
2439symbolic procedure FindSSWeights(N,nf,nb,exli,zerowei,verbose)$
2440% This procedure computes all homogeneities of a list of
2441% expressios exli which can but need not be equations.
2442% zerowei is a list of identifiers which should get zero weight.
2443begin
2444 scalar j,k,s,bsn,w,jmax,h,ex,fwli,bwli,sc,fh,bh,eh,posi,zrop,posh,negh,
2445        bs,    % list of all 2^N sequences of 0's and 1's
2446        ali,   % association list of (kernel . weight_name)
2447        allali,% list of association lists for all homogeneities
2448        wl,    % list of all weight names (needed for solveeval)
2449        sf,    % numerator of an equation (minus first terms..)
2450        tf,    % (first) term of sf
2451        wtli,  % list of conditions resulting from one equation
2452        wt,    % weight of one term in one equation
2453        eli,   % list of all conditions on weights
2454        ewli,  % list of weights of expressions/equations in exli
2455        hlist  % list of all homogeneities
2456$
2457
2458 % Consistency of input:
2459 if not pairp exli or car exli neq 'LIST then return <<
2460  write"The input expression is not a list  {  } ."$terpri()$
2461  write"Try again."$terpri()
2462 >>$
2463
2464 % Generation of bs
2465 bs:={nil};
2466 for j:=1:N do <<
2467  for each s in bs do <<
2468   bsn:=cons(cons(0,s),bsn);
2469   bsn:=cons(cons(1,s),bsn)
2470  >>;
2471  bs:=bsn; bsn:=nil
2472 >>$
2473 bs:=for each s in bs collect cons(for each j in s sum j,s)$
2474
2475 % initialization of ali
2476 ali:={('x . -2)}$
2477
2478 % identifiers with zero weight
2479 for each h in cdr reval zerowei do
2480 ali:=cons((h . 0), ali);
2481
2482 % the weights of further differentiation variables will be
2483 % added dynamically as they occur
2484
2485 % the weights of all th-variables
2486 for j:=1:N do ali:=cons(({'th,j} . -1), ali);
2487
2488 % the weights of all fermionic fields
2489 for j:=1:nf do <<
2490
2491  % adding f(j):
2492  w:=mkid('w_,mkid('f,j));
2493  wl:=cons(w,wl);
2494  ali:=cons(({'f,j} . w), ali);
2495
2496  % now all coefficients f(j,..), b(j,..) of all terms in the
2497  % th(k) power expansion of f(j):
2498  for each s in bs do
2499  ali:=cons(((if evenp car s then append({'f,j},cdr s)
2500                             else append({'b,j},cdr s)) .
2501             (if zerop car s then w
2502                             else {'PLUS,w,car s})        ), ali)
2503 >>$
2504
2505 % the weights of all bosonic fields
2506 for j:=1:nb do <<
2507
2508  % adding b(j):
2509  w:=mkid('w_,mkid('b,j));
2510  wl:=cons(w,wl);
2511  ali:=cons(({'b,j} . w), ali);
2512
2513  % now all coefficients f(j,..), b(j,..) of all terms in the
2514  % th(k) power expansion of b(j):
2515  for each s in bs do
2516  ali:=cons(((if evenp car s then append({'b,j},cdr s)
2517                             else append({'f,j},cdr s)) .
2518             (if zerop car s then w
2519                             else {'PLUS,w,car s})        ), ali)
2520 >>$
2521
2522 for each ex in cdr exli do <<
2523
2524  sf:=numr (if pairp ex and (car ex = 'EQUAL) then
2525  subtrsq(simp cadr ex,simp caddr ex)         else simp ex);
2526  % sf is the numerator of the expression
2527
2528  wtli:=nil;
2529  while sf do <<
2530   tf:=first_term_SF sf; sf:=subtrf(sf,tf); % tf is the first term
2531   wt:=nil;
2532   while tf and not domainp tf do <<
2533
2534    w:=wgtof(mvar tf,ali,wl);
2535    wt:=cons({'TIMES,ldeg tf,car w},wt);
2536    ali:=cadr w;
2537    wl:=caddr w;
2538    tf:=lc tf
2539   >>;
2540   wt:=if null wt then 0              else
2541       if cdr  wt then cons('PLUS,wt) else
2542		       car wt$
2543   wtli:=cons(reval wt,wtli)
2544  >>;
2545
2546  if wtli and cdr wtli then % ie. if there are at least 2 terms
2547  % then formulate a weight condition for each term to have
2548  % the same weight as the first term, eli is the resulting
2549  % list of conditions
2550  for each w in cdr wtli do
2551  eli:=cons(reval {'DIFFERENCE,car wtli,w}, eli)$
2552  ewli:=cons(car wtli,ewli)
2553 >>;
2554
2555 % solving the system of conditions
2556 !!arbint:=0;
2557 if null eli then
2558 if null wl then s:=nil
2559            else <<
2560  s:=nil;
2561  for each w in wl do <<
2562   !!arbint:=add1 !!arbint;
2563   s:=cons({'EQUAL,w,{'ARBCOMPLEX,!!arbint}},s)
2564  >>$
2565  s:=cons('LIST,s)
2566 >>          else <<
2567  s:=solveeval {cons('LIST,eli),cons('LIST,wl)};
2568  if s then s:=cdr s;  % to get rid of 'LIST
2569  if s then s:=car s;  % This was a linear homogeneous problem with
2570  % only one solution which we take (possibly with free parameters).
2571
2572  % If there is just a single unknown then the solution will consist
2573  % only on one {'EQUAL,w_..,..} and not a list --> we make it a list
2574  if pairp s and car s = 'EQUAL then s:={'LIST,s}
2575
2576 >>;
2577
2578 if s then <<
2579
2580  fwli:=cons('LIST,for j:=1:nf collect cdr assoc({'f,j},ali));
2581  bwli:=cons('LIST,for j:=1:nb collect cdr assoc({'b,j},ali));
2582  ewli:=cons('LIST,ewli);
2583
2584
2585  jmax:=if !!arbint=0 then 1 else !!arbint$
2586  for j:=1:jmax do <<
2587   sc:=s;
2588   for k:=1:!!arbint do
2589   if k neq j then sc:=algebraic(sub(arbcomplex(lisp k)=0,sc));
2590
2591   % storing all homogeneities
2592   k:=ali;
2593   for each h in cdr sc do
2594   k:=subst(caddr h,reval cadr h,k);
2595   allali:=cons(k,allali);
2596
2597   if j neq 0 then
2598   sc:=algebraic(sub(arbcomplex(lisp j)=1,sc));
2599
2600   % storing the homogeneity weights of f(1),f(2),.. b(1),b(2),..
2601   % in sorted lists:
2602   fh:=algebraic(sub(lisp sc,lisp fwli));
2603   bh:=algebraic(sub(lisp sc,lisp bwli));
2604   eh:=algebraic(sub(lisp sc,lisp ewli));
2605
2606   % All homogeneities with strictly positive weights will be
2607   % collected in posh and all others in negh. Then they are
2608   % appended with posh coming first.
2609   posi:=t; zrop:=t;
2610   for each k in cdr fh do <<
2611    if k<=0 then posi:=nil;
2612    if k neq 0 then zrop:=nil
2613   >>;
2614   for each k in cdr bh do <<
2615    if k<=0 then posi:=nil;
2616    if k neq 0 then zrop:=nil
2617   >>$
2618   if null zrop then <<
2619    if posi then posh:=cons({'LIST,fh,bh,eh},posh)
2620	    else negh:=cons({'LIST,fh,bh,eh},negh);
2621   >>
2622  >>
2623 >>;
2624 hlist:=cons('LIST,append(posh,negh));
2625
2626 % Printing all homogeneities
2627 if verbose then <<
2628  if hlist={'LIST} then write"This system is not homogeneous."
2629		   else <<
2630   if !!arbint=0 then
2631   write"This system has the following homogeneity:" else
2632   write"This system has the following homogeneities:"$
2633   terpri()$
2634   for each ali in allali do <<
2635    for each w in ali do
2636    algebraic(write"W[",lisp car w,"] = ",lisp reval cdr w)$
2637    write"================================="$terpri()$
2638   >>$
2639
2640   write"The program returns a list of "$
2641   h:=length hlist - 1;
2642   if h=1 then <<write"one homogeneity"$terpri()$
2643    write"which is a list of"
2644   >>     else <<write h," homogeneities,"$ terpri()$
2645    write"each homogeneity being a list of"$
2646   >>$
2647   terpri()$
2648   write"a list of all f(1),b(2),.. weights and"$terpri()$
2649   write"a list of all b(1),b(2),.. weights and"$terpri()$
2650   write"a list of the weights of equations/expressions"$terpri()$
2651   write"in the input. "$
2652   if !!arbint neq 0 then <<
2653    write"In this output free parameters"$terpri()$
2654    write"arbcomplex(i) are replaced by 1:"$
2655   >>$
2656   terpri()$
2657   mathprint hlist$
2658  >>
2659 >>$
2660
2661 return hlist
2662
2663end$
2664
2665%-------
2666
2667algebraic procedure linearize(pdes,nf,nb,tpar,spar)$
2668% pdes   .. a list of equations with a t-derivative on the left hand side
2669% nf     .. number of fermion fields f(1) .. f(nf)
2670% nb     .. number of boson   fields b(1) .. b(nb)
2671% tpar   .. whether the variable t is parity changing
2672% spar   .. whether the symmetry parameter s of the symmetry that corresponds
2673%           to this linearization is parity changing (t) or not (nil)
2674% linearize returns
2675% {list of relations that define the introduced fields like df(f(3),s)=f(6),
2676%  list of linearized equations}
2677% Newly introduced fields depend on x,t,s.
2678begin scalar spar_bak,tpar_bak,n,m,p,npdes,rel$
2679
2680 spar_bak:=lisp s_changes_parity$
2681 tpar_bak:=lisp t_changes_parity$
2682 if tpar then on  t_changes_parity
2683         else off t_changes_parity$
2684
2685 % generation of a substitution list f --> F, b --> B for spar=nil
2686 %                                or f --> B, f --> F got spar=t
2687 n:=length pdes$
2688 lisp << terpri()$
2689  write"The following linearization generates ",
2690       if n = 1 then "a condition" else "conditions",
2691       " for the right hand side",
2692       if n = 1 then " " else "s "$terpri()$
2693  write"of the symmetry: "$
2694  terpri()$terpri()
2695 >>$
2696 rel:={}$
2697 if spar then <<
2698  on  s_changes_parity$
2699  for n:=1:nf do <<m:=nb+n;depend b(m),x,t;depend f(n),s;
2700                   rel:=cons(df(f(n),s)=b(m),rel);df(f(n),s):=b(m);
2701                   lisp write "   D_s f(",n,") = b(",m,")">>$
2702  for n:=1:nb do <<m:=nf+n;depend f(m),x,t;depend b(n),s;
2703                   rel:=cons(df(b(n),s)=f(m),rel);df(b(n),s):=f(m);
2704                   lisp write "   D_s b(",n,") = f(",m,")">>
2705 >>      else <<
2706  off s_changes_parity$
2707  for n:=1:nf do <<m:=nf+n;depend b(m),x,t;depend f(n),s;
2708                   rel:=cons(df(f(n),s)=f(m),rel);df(f(n),s):=f(m);
2709                   lisp write "   D_s f(",n,") = f(",m,")">>$
2710  for n:=1:nb do <<m:=nb+n;depend f(m),x,t;depend b(n),s;
2711                   rel:=cons(df(b(n),s)=b(m),rel);df(b(n),s):=b(m);
2712                   lisp write "   D_s b(",n,") = b(",m,")">>
2713 >>$
2714 lisp terpri()$
2715
2716 write"The original system:"$
2717 for each p in pdes do write p$
2718
2719 npdes:={};
2720 for each p in pdes do <<
2721  m:=df(lhs p,s)$
2722  if part(m,0)=minus then n:=t else n:=nil$
2723  npdes:=sqcons(if tpar and spar then if n then (-m =  df(rhs p,s))
2724                                           else ( m = -df(rhs p,s))
2725                                 else if n then (-m = -df(rhs p,s))
2726                                           else ( m =  df(rhs p,s)), npdes)
2727 >>$
2728
2729 npdes:=reverse npdes;
2730 write"The linearized system: "$
2731 for each p in npdes do write p$
2732 write"and in machine readable form: "$ off nat$
2733 for each p in npdes do write p$
2734 on  nat$
2735
2736 if lisp(spar_bak neq s_changes_parity) then
2737 if spar_bak then on  s_changes_parity
2738             else off s_changes_parity$
2739
2740 if lisp(tpar_bak neq t_changes_parity) then
2741 if tpar_bak then on  t_changes_parity
2742             else off t_changes_parity$
2743
2744 % for n:=1:nf do clear df(f(n),s)$
2745 % for n:=1:nb do clear df(b(n),s)$
2746
2747 return {rel,npdes}
2748end$
2749
2750%-------
2751
2752symbolic operator  GenSSPoly$
2753symbolic procedure GenSSPoly(N,wgtlist,cname,mode)$
2754begin
2755 scalar linsub,g,h,N,flist,blist,afblist,forbid,non0coeff,fname_bak,
2756        fl,allf,allb,fwlist,bwlist,nf,nb,verbose,afwlist,abwlist,pol,
2757        fonly,bonly,wgt,newf,newb,print_bak$
2758
2759 wgtlist:=cdr wgtlist; % to drop 'LIST
2760 afwlist:=cadr   car wgtlist;
2761 abwlist:=caddr  car wgtlist;
2762 wgt    :=cadddr car wgtlist;
2763 wgtlist:=cdr wgtlist;
2764
2765 nf:=length afwlist - 1; fwlist:=cdr reval afwlist$
2766 nb:=length abwlist - 1; bwlist:=cdr reval abwlist$
2767
2768 % Initialization of flags
2769 mode:=cdr mode;       % to drop 'LIST
2770 while mode do <<
2771  if car mode = 'lin then <<
2772   if not evenp nf then
2773   rederr"The flag `lin' can not be run with odd many fermions"$
2774   linsub:=for h:=(nf/2+1):nf collect {'f,h}$
2775   if not evenp nb then
2776   rederr"The flag `lin' can not be run with odd many bosons"$
2777   linsub:=cons('LIST,append(linsub,for h:=(nb/2+1):nb collect {'b,h}))
2778  >>                 else
2779  if car mode = 'verbose then verbose:=t else
2780  if car mode = 'fonly then fonly:=t else
2781  if car mode = 'bonly then bonly:=t else
2782  if pairp car mode and cadar mode = 'forbid    then forbid   :=cddar mode else
2783  if pairp car mode and cadar mode = 'non0coeff then non0coeff:=cddar mode$
2784  mode:=cdr mode
2785 >>$
2786
2787 % Initialization of fermionic and bosonic fields + their printing
2788 input_consistency_test(afwlist,abwlist)$
2789 N_:=N;  % to avoid printing the index of D in crack_out for N=1
2790
2791 flist := for g:=1:nf collect {'f,g};
2792 % if nf=1 then put('f,'prifn,'myfpri) else
2793 % if nf>1 then put('f,'prifn,nil);
2794 % This was commented out to avoid confusion when nf=1 and using th(i)
2795 % expansions and coefficients f(1,1,0,0,1,..).
2796
2797 blist := for g:=1:nb collect {'b,g};
2798 % if nb=1 then put('b,'prifn,'myfpri) else
2799 % if nb>1 then put('b,'prifn,nil);
2800
2801 % Making f(i) and b(i) x,t-dependent
2802 algebraic <<
2803  afblist:=lisp(cons('LIST,append(flist,blist)))$
2804  for each g in afblist do depend g,x,t$  %<############## t-dependence?
2805 >>$
2806
2807 % specifying the names of constants in the symmetries
2808 fname_bak:=fname_; fname_:=reval cname;
2809
2810 % Add to each field its weight as a prefix to create allf and allb
2811 for i:=1:nf do allf:=cons({nth(fwlist,i),nth(flist,i)},allf);
2812 for i:=1:nb do allb:=cons({nth(bwlist,i),nth(blist,i)},allb);
2813
2814 % generation of a list of all terms
2815 pol:=rhs_term_list(N,allf,allb,wgt,linsub,forbid,verbose)$ % returns {fprod,bprod}
2816 pol:={'LIST,cons('LIST,car pol),cons('LIST,cadr pol)}$
2817
2818 %>>>>>> add to avoid all forbid factors also if they are f(i) or b(i) or d(i,..)
2819
2820 % applying the remaining homogeneities
2821 while wgtlist do <<
2822  afwlist:=cadr   car wgtlist;
2823  abwlist:=caddr  car wgtlist;
2824  wgt    :=cadddr car wgtlist;
2825  wgtlist:=cdr wgtlist;
2826
2827  algebraic <<
2828   % creating a substitution rule
2829   su:={};
2830   g:=0$for each h in afwlist do <<g:=g+1;su:=cons(f(g)=f(g)*lin_test_const**h,su)>>$
2831   g:=0$for each h in abwlist do <<g:=g+1;su:=cons(b(g)=b(g)*lin_test_const**h,su)>>$
2832   pol:=sub(su,pol);
2833   let lrule1; pol:=pol; clearrules lrule1;
2834   let lrule2; pol:=pol; clearrules lrule2;
2835
2836   newf:={}$
2837   for each h in first pol do <<
2838    h:=h/lin_test_const**wgt;
2839    if freeof(h,lin_test_const) then newf:=cons(h,newf)
2840   >>$
2841
2842   newb:={}$
2843   for each h in second pol do <<
2844    h:=h/lin_test_const**wgt;
2845    if freeof(h,lin_test_const) then newb:=cons(h,newb)
2846   >>$
2847
2848   pol:={newf,newb}
2849  >>
2850 >>$
2851
2852 print_bak:=print_; print_:=nil;
2853 newf:=if bonly then 0
2854                else <<
2855  h:=for each g in cdadr pol collect <<h    :=newfct(fname_,nil,nfct_)$
2856	 			       fl   :=cons(h,fl);
2857				       nfct_:=add1 nfct_$
2858				       {'TIMES,h,g}>>$
2859  if h then if cdr h then cons('PLUS,h)
2860                     else car h
2861       else 0
2862 >>$
2863 newb:=if fonly then 0
2864                else <<
2865  h:=for each g in cdaddr pol collect <<h    :=newfct(fname_,nil,nfct_)$
2866	 		   	        fl   :=cons(h,fl);
2867				        nfct_:=add1 nfct_$
2868				        {'TIMES,h,g}>>$
2869  if h then if cdr h then cons('PLUS,h)
2870                     else car h
2871       else 0
2872 >>$
2873 print_:=print_bak;
2874 fname_:=fname_bak;
2875 flin_:=fl$
2876
2877 return {'LIST,newf,newb,cons('LIST,fl)}
2878end$
2879
2880%-------
2881
2882algebraic procedure thgen(k)$
2883% generates a list with 2^N elements, each element generating
2884% a term in the th-power series and has the form:
2885% {product_of_th(i),
2886%  number_of_factors_th(i),
2887%  {0's and 1's, i.e. a 1 for each th() otherwise zeros} }
2888if k=0 then {{1,0,{}}}
2889       else begin
2890 scalar a,b$
2891 a:= thgen(k-1)$
2892 return for each b in a
2893        join {{first b,       second b,   append (third b,{0})},
2894              {th(k)*first(b),second(b)+1,append (third b,{1})} }
2895end$
2896
2897Drule:={D(~i,~a)=>df(a,th(i))+th(i)*df(a,x)}$ % can not be defined within tocoo
2898
2899algebraic procedure tocoo(N,nf,nb,ex)$
2900% nf: number of fermion fields
2901% nb: number of boson fields
2902% ex: polynomial expression in field form to be transformed into coordinate form
2903begin
2904 scalar thli,subfnb,i,h$
2905 lisp(N_:=N);
2906 put('f,'prifn,nil);
2907 put('b,'prifn,nil);
2908 thli:=thgen(N)$
2909
2910 % Generate the coordinate form of f(i)
2911 subfnb:={};
2912 for i:=1:nf do
2913 subfnb:=cons(f(i)=for each h in thli sum
2914   if evenp second h then lisp(cons('f,cdr algebraic cons(i,third h)))*first h
2915	             else lisp(cons('b,cdr algebraic cons(i,third h)))*first h,
2916              subfnb);
2917
2918 % Generate the coordinate form of b(i)
2919 for i:=1:nb do
2920 subfnb:=cons(b(i)=for each h in thli sum
2921   if evenp second h then lisp(cons('b,cdr algebraic cons(i,third h)))*first h
2922                     else lisp(cons('f,cdr algebraic cons(i,third h)))*first h,
2923              subfnb);
2924
2925 % Do the substitution
2926 ex:=sub(subfnb,ex)$
2927 % Replace D(i,..) which has to be done AFTER the above substitution
2928 let Drule;
2929 ex:=ex;
2930 clearrules Drule;
2931
2932 return ex
2933
2934end$
2935
2936%-------
2937
2938algebraic procedure tofield(N,nf,nb,ex,zerowei)$
2939% tw the total weight of the input expression ex which is
2940% expected to be either purely fermionic or purely bosonic.
2941
2942% In a call to GenSSPoly() an ansatz (ssp) in field form with the proper
2943% homogeneity weights is generated, then converted to coordinate form (ans)
2944% in a call to tocoo() and compared to the input in coordinate form (ex).
2945% The difference is splitted through a call to split_simp() and the
2946% resulting system of conditions for the undetermined coefficients of
2947% the ansatz has to vanish.
2948
2949if lisp (pairp algebraic ex and car ex = 'LIST)            then
2950write"The input should be a polynomial, not an equations!" else
2951begin scalar ssp,ans,fieldans,rlist,indepv,h,eqnli,solu,wflist,wblist,tw,
2952 gs,k,ft,sb,notnum,notpos$
2953
2954 w:=FindSSWeights(N,nf,nb,{ex},zerowei,t)$ % e.g.  w={{{1,3},{},{5}}}
2955 if w={} then return lisp <<
2956  write"The input expression does not have a homogeneity"$ terpri()$
2957  write"which could be used to compute a field form."$     terpri()
2958 >>$
2959
2960 w:=first w$ % We just use the fist homogeneity because if there are
2961             % several and there is one with strictly positive weights
2962             % then this comes first.
2963
2964 wflist:=first w;
2965 wblist:=second w;
2966 tw:=first third w;
2967
2968 for each h in append(wflist,wblist) do
2969 if not fixp h then notnum:=t else
2970 if h <= 0 then notpos:=t$
2971
2972 if notnum then lisp <<
2973  write"#####   Strange error: returned weights returned"$terpri()$
2974  write"        from FindSSWeights() are not integers!   #####"$terpri()
2975 >>        else
2976 if notpos then lisp <<
2977  gs:=gensym()$
2978  k:=setkorder list gs$
2979  ft:=append(search_li2(algebraic ex,'f),
2980             search_li2(algebraic ex,'b) )$   % all f(..), b(..) functions
2981  for each h in ft do sb:=cons((h . {'TIMES,h,gs}), sb)$
2982  h:=numr simp {'!*sq,subsq(simp algebraic ex,sb),nil}$
2983  algebraic(max_deg):=if mvar h neq gs then 0 else
2984                      ldeg numr simp {'!*sq,subsq(simp algebraic ex,sb),nil}$
2985  setkorder k
2986 >>$
2987
2988 % generating an ansatz in field form (ssp)
2989 % and convert it back to coordinate form (ans)
2990 % Next line works so far only if not compilied.
2991 if is_fermionic(lisp{'!*SQ,((first_term_SF numr simp algebraic ex) . 1),t}) then <<
2992 % if (lisp{'!*SQ,((first_term_SF numr simp algebraic ex) . 1),t})**2=0 then <<
2993  ssp:=gensspoly(N,{{wflist,wblist,tw}},r,{fonly});
2994  fieldans:=first(ssp);
2995 >>                                                                          else <<
2996  ssp:=gensspoly(N,{{wflist,wblist,tw}},r,{bonly});
2997  fieldans:=second(ssp);
2998 >>;
2999
3000 ans:=tocoo(N,nf,nb,fieldans)$
3001
3002 %spliting the equation
3003 rlist:=third(ssp); % a list of undetermined coefficients r1,r2, ...
3004 indepv:=append({x,t},for h:=1:N collect th(h))$ % all independent variables
3005 lisp(print_:=nil)$ % to suppress printing about the splitting
3006 eqnli:=split_simp({ans-ex},{},rlist,indepv,t)$
3007
3008 %applying crack to find solutions
3009 solu:=crack(eqnli,{},rlist,{})$
3010
3011 return if solu={} then lisp <<
3012  write"The expression to be transformed has a homogeneity"$terpri()$
3013  write"but a field form does not exist."$terpri()
3014 >>                else <<solu:=second first (solu);
3015                          sub(solu,fieldans)>>
3016
3017end$
3018
3019%-------
3020
3021symbolic procedure restore_interactive_prompt$
3022  begin scalar !*redefmsg,!*usermode$
3023    copyd('update_prompt,'restore_update_prompt)
3024  end$
3025
3026%-------
3027
3028symbolic procedure change_prompt$
3029  begin scalar !*usermode$
3030    if null promptstring!* then promptstring!* := "";
3031    setpchar promptstring!*;
3032    promptexp!* := promptstring!*
3033  end$
3034
3035%-------
3036
3037symbolic procedure change_prompt_to u$
3038   begin scalar oldprompt,!*redefmsg,!*usermode$
3039     oldprompt := promptstring!*;
3040     promptstring!* := u;
3041     copyd('restore_update_prompt,'update_prompt);
3042     copyd('update_prompt,'change_prompt);
3043     update_prompt();
3044     restore_interactive_prompt()$
3045     return oldprompt
3046   end$
3047
3048%-------
3049
3050symbolic procedure cnt_l_$
3051if lines_written geq 10 then <<
3052  change_prompt_to ""$
3053  write"                                   Press Enter to continue "$
3054  restore_interactive_prompt()$
3055  rds nil; wrs nil$         % Switch I/O to terminal
3056  readch()$
3057  if ifl!* then rds cadr ifl!*$  %  Resets I/O streams
3058  if ofl!* then wrs cdr ofl!*$
3059  lines_written:=0
3060>>$
3061
3062%-------
3063
3064symbolic procedure wl1(l)$
3065<<write l$ terpri()$ cnt_l_()$ lines_written:=lines_written+1>>$
3066symbolic procedure wl2(l)$
3067<<write l$ terpri()$ terpri()$ cnt_l_()$ lines_written:=lines_written+2>>$
3068
3069%-------
3070
3071symbolic procedure sshelp1$
3072<<wl2"Purpose: "$
3073
3074 wl1"to determine evolutionary supersymmetric PDEs with higher"$
3075 wl1"symmetries, to compute symmetries and conservation laws for"$
3076 wl1"given supersymmetric PDEs, or to do any other algebra or"$
3077 wl2"differentiations of polynomial supersymmetric expressions."$
3078>>$
3079
3080%-------
3081
3082symbolic procedure sshelp2$
3083<<wl2"Notation: "$
3084
3085 wl1"Fields f(1),f(2),.. are treated as fermionic, and b(1),.. are"$
3086 wl1"treated as bosonic. The independent variables are t, x and in"$
3087 wl1"symmetry computations the symmetry `time' variable is s. "$
3088 wl2"Further fermionic variables can be defined, for example, theta."$
3089
3090 wl1"The super derivatives D_i which satisfy (D_i)**2 = d/dx "$
3091 wl2"are implemented as d(i, .. ) .  "$
3092
3093 wl1"When refering to homogeneity weights, all weights are scaled"$
3094 wl1"such that the weight of d(i, .. ) is 1 and thus the weight"$
3095 wl2"of df( .. , x) is two, i.e. twice the usual value. "$
3096>>$
3097
3098%-------
3099
3100symbolic procedure sshelp3$
3101<<wl2"Initializations: "$
3102
3103 wl1"Before starting super-computations the number N of superfields theta_i"$
3104 wl1"needs to be known. Also the number nb of boson fields and the nf of"$
3105 wl1"fermion fields is needed to have compact printing on the screen (to"$
3106 wl1"avoid printing the index 1 if nf=1 or nb=1). The input of N,nf,nb is"$
3107 wl1"done either through a call of the procedure ssym (to compute"$
3108 wl1"symmetries) or of the procedure ssconl (to compute conservation laws)"$
3109 wl1"or with a call of procedure ssini. This help item describes ssini."$
3110 wl2"It's format is:"$
3111
3112 wl2"  ssini(N,nf,nb)$"$
3113
3114 wl2"where"$
3115
3116 wl1" N      .. the number of superfields theta_i"$
3117 wl1" nf     .. number of fermion fields f(1) .. f(nf)"$
3118 wl2" nb     .. number of boson   fields b(1) .. b(nb)"$
3119
3120 wl2"-----------------"$
3121
3122 wl2"Example: For N=1 and 2 fermion and 2 boson fields we do"$
3123
3124 wl2"ssini(1,2,2)$ "$
3125
3126 wl2"and verify easily that fermions anticommute: "$
3127
3128 wl1"f(1)**2; "$
3129 wl2"f(1)*f(2)+f(2)*f(1);"$
3130
3131 wl2"that bosons commute:"$
3132
3133 wl1"b(1)**3;"$
3134 wl2"b(1)*b(2)-b(2)*b(1);"$
3135
3136 wl2"that x,t-differentiations by default do not change parity"$
3137
3138 wl1"df(b(1),t)**2;"$
3139 wl2"df(b(1),x)**2;"$
3140
3141 wl2"unless one explicitly defines d/dt to be parity changing through"$
3142
3143 wl2"on t_changes_parity$"$
3144
3145 wl2"demonstrated by"$
3146
3147 wl2"df(b(1),t)**2;"$
3148
3149 wl1"which is also possible to define for the symmetry variable s "$
3150 wl1"(through       on s_changes_parity$             , see below)"$
3151 wl2"but not for x. For the remainder of this demo we change back"$
3152
3153 wl1"off t_changes_parity$"$
3154 wl2"df(b(1),t)**2;"$
3155
3156 wl1"Differentiations d/dt, d/dx, D_i follow the supersymmetric"$
3157 wl2"product rule:"$
3158
3159 wl1"df(b(1)*f(2),x);"$
3160 wl1"d(1,b(1)*f(2));"$
3161 wl1"df(f(1)*f(2),x);"$
3162 wl2"d(1,f(1)*f(2));"$
3163
3164 wl2"furthermore D_i changes parity:"$
3165
3166 wl1"d(1,f(1))**2;"$
3167 wl2"d(1,b(1))**2;"$
3168
3169 wl2"and satisfies  D^2 = d/dx:"$
3170
3171 wl2"d(1,d(1,b(1)));"$
3172
3173 wl1"To compactify output, in the special case N=1 the index of D_i"$
3174 wl1"is not printed as seen above. Similarly, indices of f(1), b(1)"$
3175 wl2"are not printed if we have only one fermion and/or one boson field:"$
3176
3177 wl1"ssini(2,1,2)$"$
3178 wl2"f(1)*b(1)*d(1,b(2));"$
3179
3180 wl1"Any initializations made by ssini are overwritten by a call to the"$
3181 wl2"procedures ssym or ssconl."$
3182
3183>>$
3184
3185%-------
3186
3187symbolic procedure sshelp4$
3188<<wl2"Coefficients: "$
3189
3190 wl1"coeffn is a standard REDUCE command although its implementation"$
3191 wl1"(at least in REDUCE distributions up to version 3.8) does not work"$
3192 wl1"with non-commuting variables. By loading SsTools some of the REDUCE"$
3193 wl1"routines are re-defined so that coeffn can be used for supersymmetric"$
3194 wl2"expressions. The call is unchanged from the standard REDUCE command:"$
3195
3196 wl2"coeffn(ex,h,n)$"$
3197
3198 wl1"which computes the coefficient of h**n in the expression ex which"$
3199 wl1"is supposed to be polynomial in h. The returned result is the"$
3200 wl1"coefficient AFTER h**n has been moved to be the first factor in"$
3201 wl2"each term of ex. This matters if h is fermionic and thus n=1."$
3202
3203 wl2"Example:"$
3204
3205 wl1"SSIni(2,0,2)$"$
3206 wl1"ex:=D(2,b(1))*D(1,df(b(2),x,2));"$
3207 wl2"coeffn(ex,D(1,df(b(2),x,2)),1);"$
3208
3209 wl1"Because D_2 b(1) and D_1 D_x D_x b(2) are fermionic a minus sign is"$
3210 wl1"introduced when making D_1 D_x D_x b(2) the first factor in the"$
3211 wl2"product. This explains the minus in the result of the coeffn call."$
3212
3213 wl1"NOTE: When using SsTools the REDUCE command noncom has no effect,"$
3214 wl1"i.e. it is not possible to define additional noncommuting quantities"$
3215 wl2"besides fermionic ones which are introduced by the command fermionic."$
3216
3217>>$
3218
3219%-------
3220
3221symbolic procedure sshelp5$
3222<<% special use comments:
3223 % in drvlist$ % loads the ansatz of the previous ssym run
3224 % use_new_crackout:=t$ % needed to use the ssym specific crack_out()
3225
3226 wl2"Computing symmetries:"$
3227
3228 wl1"To compute PDEs together with higher symmetries, or to"$
3229 wl2"compute higher symmetries for a given PDE(-system) the call is:"$
3230
3231 wl2"  ssym(N,tw,sw,afwlist,abwlist,eqnlist,fl,inelist,mode)$"$
3232
3233 wl2"where"$
3234
3235 wl1" N       .. the number of superfields theta_i"$
3236 wl1" tw      .. 2 times the differential order of the equation = weight(t)"$
3237 wl1" sw      .. 2 times the differential order of the symmetry = weight(s)"$
3238 wl1" afwlist .. (algebraic mode) list of weights of the fermion fields"$
3239 wl1"            f(1), f(2), ... . The number of elements of this list"$
3240 wl1"            determines the number of fermion fields. "$
3241 wl1" abwlist .. (algebraic mode) list of weights of the boson fields"$
3242 wl1"            b(1), b(2), ... . The number of elements of this list"$
3243 wl1"            determines the number of boson fields. "$
3244 wl1" eqnlist .. - in the nonlinear case a list of extra conditions on the"$
3245 wl1"              undetermined coefficients where conditions a3=.. are executed"$
3246 wl1"              instantly. These and any expressions are added to equations"$
3247 wl1"              when calling crack"$
3248 wl1"            - in the linear case the system in form of replacements => and"$
3249 wl1"              the linearized system in form of equations = ."$
3250 wl1" fl      .. extra unknowns in eqnlist to be determined"$
3251 wl1" inelist .. a list, each element of it is a list with at least one of"$
3252 wl1"            its elements being non-zero"$
3253 wl1" mode    .. list of flags: "$
3254 wl1"            init: only initialization of global data"$
3255 wl1"            plain_com : direct computation of the commutator (for tests)"$
3256 wl1"            power_split_com : alternatice power splitting of commutator,"$
3257 wl1"                             (not if substitutions '=>' are present,"$
3258 wl1"                              see below)"$
3259 wl1"            zerocoeff: all coefficients = 0 which do not appear in inelist"$
3260 wl1"            filter: extra homogeneity weights as given in the list"$
3261 wl1"                    hom_wei have to be satisfied by the symmetry"$
3262 wl1"            lin: find symmetry that is linear homogeneous in all fields "$
3263 wl1"                 with (index)>(maxindex/2), i.e. if lin then the"$
3264 wl1"                 number of fermions and boson fields must both be even"$
3265 wl1"            tpar: t/nil, whether time variable t changes parity"$
3266 wl2"            spar: t/nil, whether symmetry variable s changes parity"$
3267
3268 wl2"-----------------"$
3269
3270 wl1"The question whether ssym is to be used to compute a PDE(-system) with"$
3271 wl1"symmetries, or to compute symmetries for a `given' PDE(-system) is"$
3272 wl1"decided based on the form of the input eqnlist. "$
3273 wl1"If symmetries of a specific weight for a given PDE(-system) are to"$
3274 wl1"be computed then eqnlist has to have the form"$
3275 wl1"{df(f(1),t)=..., df(b(1),t)=...} with as many t-derivatives"$
3276 wl1"of fermion fields as there are elements in afwlist and as many"$
3277 wl2"t-derivatives of boson fields as there are elements in abwlist."$
3278
3279 wl1"If the right hand sides of the PDE(s) in eqnlist involve any "$
3280 wl1"constants which are to be computed then these constants have to "$
3281 wl1"be listed in the input list fl, like {p1,p2,...}. If these constant "$
3282 wl1"coefficients are not given in the list fl then they are treated as "$
3283 wl1"independent parameters and symmetries are determined which are "$
3284 wl1"symmetries for generic values of these constant coefficients."$
3285 wl1"Also, when computing symmetries for a given PDE(-system) then "$
3286 wl1"dependencies of all the fields on t,x has to be declared"$
3287 wl1"beforehand in order to be able to input derivatives, like df(f(1),t) "$
3288 wl2"which otherwise would be automatically evaluated to zero."$
3289
3290 wl1"If the input list eqnlist does not start with {df(..)=..} then"$
3291 wl1"an ansatz for the PDE(-system) is generated and this system and"$
3292 wl1"its symmetries are computed. In this case eqnlist can contain"$
3293 wl1"substitutions, like p2=0, or expressions which are to be set to"$
3294 wl1"zero, like p3*r4+p2*r3. Typically one would run the program first "$
3295 wl1"to see which ansatz for the PDE(-system) and the symmetry is "$
3296 wl1"generated and then start ssym again with the intended extra"$
3297 wl2"conditions."$
3298
3299 wl1"If a boson weight is non-positive or a fermion weight is negative"$
3300 wl1"then the global (algebraic mode) variable max_deg must have a positive"$
3301 wl2"integer value which is the highest degree of such a field in any ansatz."$
3302
3303 wl1"When determining symmetries (either for a given system of equations,"$
3304 wl1"or when determining equations and symmetries in one run) one has two"$
3305 wl1"cases: whether the symmetry variable s flips parity, or not, in other"$
3306 wl1"words, whether s itself is fermionic or bosonic. This has nothing to do"$
3307 wl1"with the question whether s has an odd or even weight."$
3308 wl1"Similarly, when determining equations and symmetries in one computation"$
3309 wl1"one can consider the two cases whether the `time' derivative changes"$
3310 wl1"parity or not. A parity changing t is specified through the flag `tpar'"$
3311 wl2"and a parity changing s through the flag `spar'."$
3312
3313 wl2"-----------------"$
3314
3315 wl2"Example for determining PDEs + symmetries:"$
3316
3317 wl2"ssym(1,4,5,{2},{2},{},{},{},{})$ "$
3318
3319 wl2"Example for determining the symmetries of a given PDE-system:"$
3320
3321 wl1"ssym(1,4,5,{2},{2},"$
3322 wl1"     {df(f(1),t)=df(f(1),x)*b(1)*p9,"$
3323 wl1"      df(b(1),t)=df(b(1),x)*b(1)*p9 + df(f(1),x)*f(1)*p4},"$
3324 wl2"     {},{},{})$ "$
3325
3326 wl2"-----------------"$
3327
3328 wl1"Example for determining symmetries of given PDEs in the"$
3329 wl1"presence of extra substitution rules which are given"$
3330 wl1"using `=>' instead of `=' for the equations for which"$
3331 wl1"matching symmetries are to be found:"$
3332
3333 wl1"lisp put('f,'prifn,nil)$  % from now on more than one fermion "$
3334 wl1"lisp put('b,'prifn,nil)$  % from now on more than one boson "$
3335 wl1"ssym(1,1,4,{1,1},{1,3,1,3}, "$
3336 wl1"     {df(f(1),t)=>-2*f(1)*b(1)*p1,  "$
3337 wl1"      df(b(1),t)=>b(1)**2*p1+d(1,f(1))*p2,  "$
3338 wl1"      d(1,b(2)) =>-df(b(1),x)*f(1),  "$
3339 wl1"      df(b(2),t)=>-2*d(1,b(1))*f(1)*b(1)*p1+d(1,df(b(1),t))*f(1)"$
3340 wl1"                  -d(1,f(1))**2*p2/2-d(1,f(1))*b(1)**2*p1"$
3341 wl1"                  +d(1,f(1))*df(b(1),t),"$
3342 wl1"      df(f(2),t)=-2*f(2)*b(1)*p1-2*f(1)*b(3)*p1,  "$
3343 wl1"      df(b(3),t)=2*b(1)*b(3)*p1+d(1,f(2))*p2,  "$
3344 wl1"      d(1,b(4))=>-df(b(3),x)*f(1)-df(b(1),x)*f(2),  "$
3345 wl1"      df(b(4),t)=>-2*d(1,b(3))*f(1)*b(1)*p1-2*d(1,b(1))*f(2)*b(1)*p1"$
3346 wl1"                  -2*d(1,b(1))*f(1)*b(3)*p1+d(1,df(b(3),t))*f(1) "$
3347 wl1"                  +d(1,df(b(1),t))*f(2)-d(1,f(1))*d(1,f(2))*p2 "$
3348 wl1"                  -d(1,f(2))*b(1)**2*p1-2*d(1,f(1))*b(3)*b(1)*p1"$
3349 wl1"                  +d(1,f(2))*df(b(1),t)+d(1,f(1))*df(b(3),t)},"$
3350 wl2"     {},{},{});  "$
3351
3352 wl2"-----------------"$
3353
3354 wl1"The same example again but now with the flag `lin' to generate a"$
3355 wl1"symmetry that is linear homogeneous in all fields with an index"$
3356 wl1"greater than half the maximum index for that type of field. When "$
3357 wl1"setting the flag 'lin' the number of fermions and boson fields must"$
3358 wl1"both be even. In the following example that means, the symmetry will"$
3359 wl2"be linear and homogeneous in f(2),b(3),b(4)."$
3360
3361 wl1"ssym(1,1,4,{1,1},{1,3,1,3},  "$
3362 wl1"     {df(f(1),t)=>-2*f(1)*b(1)*p1,  "$
3363 wl1"      df(b(1),t)=>b(1)**2*p1+d(1,f(1))*p2,  "$
3364 wl1"      d(1,b(2)) =>-df(b(1),x)*f(1),  "$
3365 wl1"      df(b(2),t)=>-2*d(1,b(1))*f(1)*b(1)*p1+d(1,df(b(1),t))*f(1)"$
3366 wl1"                  -d(1,f(1))**2*p2/2-d(1,f(1))*b(1)**2*p1"$
3367 wl1"                  +d(1,f(1))*df(b(1),t),"$
3368 wl1"      df(f(2),t)=-2*f(2)*b(1)*p1-2*f(1)*b(3)*p1,  "$
3369 wl1"      df(b(3),t)=2*b(1)*b(3)*p1+d(1,f(2))*p2,  "$
3370 wl1"      d(1,b(4)) =>-df(b(3),x)*f(1)-df(b(1),x)*f(2),  "$
3371 wl1"      df(b(4),t)=>-2*d(1,b(3))*f(1)*b(1)*p1-2*d(1,b(1))*f(2)*b(1)*p1"$
3372 wl1"                  -2*d(1,b(1))*f(1)*b(3)*p1+d(1,df(b(3),t))*f(1) "$
3373 wl1"                  +d(1,df(b(1),t))*f(2)-d(1,f(1))*d(1,f(2))*p2"$
3374 wl1"                  -d(1,f(2))*b(1)**2*p1-2*d(1,f(1))*b(3)*b(1)*p1"$
3375 wl1"                  +d(1,f(2))*df(b(1),t)+d(1,f(1))*df(b(3),t)},"$
3376 wl2"     {},{},{lin});  "$
3377
3378 wl2"-----------------"$
3379
3380 wl1"If one wants to find for a given system of equations df(..,t)=.. a "$
3381 wl1"symmetry which satisfies in addition to a symmetry weight sw, a list "$
3382 wl1"of fermion weights afwlist and list of boson weights abwlist other "$
3383 wl1"homogeneities, then in the last parameter of the ssym call one adds the "$
3384 wl1"flag `filter' and has a global variable hom_wei which is a list of "$
3385 wl1"lists {sw,afwlist,abwlist} each being an additional set of homogeneity"$
3386 wl2"weights."$
3387
3388 wl2"Example: "$
3389
3390 wl2"hom_wei:={{10,{3,3},{2,2}}}$"$
3391
3392 wl1"ssym(1,1,6,{1,1},{1,1},"$
3393 wl1"     {df(f(1),t)=>-2*f(1)*b(1)*p1,"$
3394 wl1"      df(b(1),t)=>b(1)*b(1)*p1+d(1,f(1))*p2,"$
3395 wl1"      df(f(2),t)= - 2*f(2)*b(1)*p1 - 2*f(1)*b(2)*p1,"$
3396 wl1"      df(b(2),t)=2*b(2)*b(1)*p1 + d(1,f(2))*p2},"$
3397 wl2"     {},{},{lin,filter});"$
3398
3399 wl2"-----------------"$
3400
3401 wl1"After a run of ssym() the substitution rules, like df(f(1),t)=> ..  "$
3402 wl1"are not active. If one wanted to check symmetries interactively  "$
3403 wl1"one would have to activate these rules which are stored under the"$
3404 wl1"name `subli' by:                                    let subli$"$
3405 wl1"and if one would want to de-activate them then do:  clearrules subli$ "$
3406 wl2"For example, to check symmetries after the above run one could do  "$
3407
3408 wl2" let subli$  "$
3409
3410 wl1" df(b(3),t):=2*b(3)*b(1)*p1 + d(1,f(2))*p2$  "$
3411 wl2" df(f(2),t):=- 2*f(2)*b(1)*p1 - 2*f(1)*b(3)*p1$  "$
3412
3413 wl1" b3t:=df(b(3),t)$  "$
3414 wl2" f2t:=df(f(2),t)$  "$
3415
3416 wl1" df(b(3),s):=(1/2*d(1,b(3))*f(1)*b(1)**2*p1*p2 + 4/11*d(1,b(1))*f(2)  "$
3417 wl1"  *b(1)**2*p1*p2 + d(1,b(1))*f(1)*b(3)*b(1)*p1*p2 + 7/22*d(1,f(1))**2 "$
3418 wl1"  *b(3)*p2**2 + 7/11*d(1,f(1))*b(3)*b(1)**2*p1*p2 + 1/4*d(1,f(1))*  "$
3419 wl1"  d(1,b(3))*f(1)*p2**2 + 5/44*d(1,f(1))*d(1,b(1))*f(2)*p2**2 + 1/44*  "$
3420 wl1"  df(b(1),x)*f(2)*f(1)*p2**2 + 1/22*df(f(1),x)*f(2)*b(1)*p2**2  "$
3421 wl1"  + 5/22*df(f(1),x)*f(1)*b(3)*p2**2)/(p1*p2)$  "$
3422 wl1" df(f(2),s):=(1/4*d(1,f(2))*d(1,f(1))*f(1)*p2**2 + 1/2*d(1,f(2))  "$
3423 wl1"  *f(1)*b(1)**2*p1*p2 + 3/44*d(1,f(1))**2*f(2)*p2**2 + 3/22*d(1,f(1)) "$
3424 wl2"  *f(2)*b(1)**2*p1*p2 + 1/44*df(f(1),x)*f(2)*f(1)*p2**2)/(p1*p2)$  "$
3425
3426 wl1" b3s:=df(b(3),s)$  "$
3427 wl2" f2s:=df(f(2),s)$  "$
3428
3429 wl1" write ""ZERO = "",df(b3t,s)-df(b3s,t)$  "$
3430 wl2" write ""ZERO = "",df(f2t,s)-df(f2s,t)$  "$
3431
3432 wl2" clearrules subli$ "$
3433>>$
3434
3435%-------
3436
3437symbolic procedure sshelp6$
3438<<wl2"To compute conservation laws the call is:"$
3439
3440 wl2"  ssconl(N,tw,mincw,maxcw,afwlist,abwlist,pdes)$"$
3441
3442 wl2"where"$
3443
3444 wl1" N       .. the number of superfields theta_i"$
3445 wl1" tw      .. weight(d_t) = 2 times the differential order of the equation"$
3446 wl1" mincw   .. min weight of the conservation law"$
3447 wl1" maxcw   .. max weight of the conservation law"$
3448 wl1" afwlist .. (algebraic) list of weights of the fermion fields"$
3449 wl1"            f(1), f(2), ... "$
3450 wl1" abwlist .. (algebraic) list of weights of the boson fields"$
3451 wl1"            b(1), b(2), ..."$
3452 wl1" pdes    .. an algebraic list of the pde(s) for which a conservation"$
3453 wl1"            law is to be found, for example"$
3454 wl1"            {df(f(1),t)=df(f(1),x)*b(1)*p9,"$
3455 wl1"             df(b(1),t)=b(1)**3*p3 + d(1,f(1))**2*p2 + "$
3456 wl2"                        df(b(1),x)*b(1)*p9 + df(f(1),x)*f(1)*p4}$"$
3457
3458 wl1"The number of elements in afwlist and abwlist "$
3459 wl1"determines the number of fermion and boson fields. All weights are "$
3460 wl1"based on weight(x)=2."$
3461 wl1"Important: All fields must be made dependent on x and t before."$
3462 wl2"Fermion fields are called f(1),f(2),.. and boson fields b(1),.. ."$
3463
3464 wl1"If a boson weight is non-positive or a fermion weight is negative"$
3465 wl1"then the global (algebraic mode) variable max_deg must have a positive"$
3466 wl1"integer value which is the highest degree of such a field or any of"$
3467 wl2"its derivatives in any ansatz."$
3468
3469 wl2"-----------------"$
3470
3471 wl2"Example for determining conservation laws for a given PDE-system:"$
3472
3473 wl1"ssconl(1,4,10,15,{2},{2},"$
3474 wl1"       {df(f(1),t)=df(f(1),x)*b(1)*p9,"$
3475 wl1"        df(b(1),t)=b(1)**3*p3 + d(1,f(1))**2*p2 + df(b(1),x)*b(1)*p9 "$
3476 wl1"                   + df(f(1),x)*f(1)*p4}"$
3477 wl2"      );"$
3478>>$
3479
3480%-------
3481
3482symbolic procedure sshelp7$
3483<<wl1"To compute all possible homogeneities of a list of expressions"$
3484 wl1"or equations, i.e. to determine for each homogeneity the weight"$
3485 wl2"of d_t, of fields f(i),b(j) and other constants, the call is:"$
3486
3487 wl2"  FindSSWeights(N,nf,nb,exli,zerowei,verbose)$"$
3488
3489 wl2"where"$
3490
3491 wl1" N       .. the number of superfields theta_i"$
3492 wl1" nf      .. number of fermion fields f(1) .. f(nf)"$
3493 wl1" nb      .. number of boson   fields b(1) .. b(nb)"$
3494 wl1" exli    .. list of equations or expressions"$
3495 wl1" zerowei .. list of constants or other kernels that"$
3496 wl1"            should have zero weight"$
3497 wl2" verbose .. =t/nil whether detailed comments shall be made"$
3498
3499 wl1"Weights are scaled such that weight of d_x is 2, i.e. the weight of"$
3500 wl2"any D is 1. Input expressions can be in field form or coordinate form."$
3501
3502 wl1"The program returns a list of homogeneities,"$
3503 wl1"each homogeneity being a list of"$
3504 wl1"a list of the weights of f(1),b(2),.. and"$
3505 wl1"a list of the weights of b(1),b(2),.. and"$
3506 wl2"a list of the weights of equations/expressions in the input. "$
3507
3508
3509 wl2"-----------------"$
3510
3511 wl1"Example for determining all possible weights for a given"$
3512 wl2"(artificial) system of equations: "$
3513
3514 wl1"- th(i) denotes superfields theta_i,"$
3515 wl1"- b's and f's with 3 indices are the coefficients in the"$
3516 wl1"  Taylor expansions"$
3517 wl1"  f(i)=f(i,0,0)+b(i,1,0)*th(1)+b(i,0,1)*th(2)+f(i,1,1)*th(1)*th(2)"$
3518 wl1"- x,t,th(i) may occur not only as differentiation variables"$
3519 wl1"  but also explicitly polynomially,"$
3520 wl2"- p3,p5 are supposed to have zero weight."$
3521
3522 wl1"FindSSWeights(2,2,0,"$
3523 wl1"              {x*df(b(1,1,0),th(1),t)= d(1,d(2,f(2,0,0)))*p5 - "$
3524 wl1"                                       df(f(1,1,1)*th(1)*th(2),x,2)*p3,"$
3525 wl1"               df(f(2),t)=(df(f(2),x,2)*p3*p5 - df(f(1),x,3)*p3**2)/p5},"$
3526 wl1"              {p3,p5},"$
3527 wl2"              t)$"$
3528>>$
3529
3530%-------
3531
3532symbolic procedure sshelp8$
3533<<wl1"To compute a linearization of a system of evolution equations,"$
3534 wl2"the call is:"$
3535
3536 wl2"  linearize(pdes,nf,nb,tpar,spar)$"$
3537
3538 wl2"where"$
3539
3540 wl1" pdes   .. list of equations with a t-derivative on left hand side"$
3541 wl1" nf     .. number of fermion fields f(1) .. f(nf)"$
3542 wl1" nb     .. number of boson   fields b(1) .. b(nb)"$
3543 wl1" tpar   .. t/nil, whether time variable t changes parity"$
3544 wl2" spar   .. t/nil, whether symmetry variable s changes parity"$
3545
3546 wl1" linearize returns {list of relations like df(f(3),s) = f(6),"$
3547 wl1"                    list of linearized equations}"$
3548 wl2"Newly introduced fields depend on x,t,s."$
3549
3550 wl1"A linearization is associated with a symmetry D_s. The linearized"$
3551 wl1"equation depends on whether s is parity changing or not and furthermore"$
3552 wl1"if s is parity changing then the sign of the lhs of the equation (or,"$
3553 wl2"equally of the rhs) depends on whether t is parity changing or not."$
3554
3555 wl2"-----------------"$
3556
3557 wl2"An (artificial) example for a linearization:"$
3558
3559 wl1"linearize({df(f(1),t)=df(f(2),x)*b(1)**3*p1 "$
3560 wl1"                      + d(1,f(2))**5*df(f(1),x,2)*p2,"$
3561 wl1"           df(f(2),t)=2*d(1,df(b(1),x))*df(f(2),x,2)*df(f(1),x)*p3 "$
3562 wl1"                      + df(f(1),x,3)*p3**2,"$
3563 wl2"           d(1,d(2,b(1)))=b(1)*d(1,d(2,f(1)))*f(1) },2,1,nil,nil)$"$
3564>>$
3565
3566%-------
3567
3568symbolic procedure sshelp9$
3569<<wl1"To make an ansatz for a polynom in a number of fermionic and/or"$
3570 wl1"bosonic variables, and their x-derivatives and D(i,..) derivatives"$
3571 wl2"The call is:"$
3572
3573 wl2"  GenSSPoly(N,wgtlist,cname,mode)$"$
3574
3575 wl2"where"$
3576
3577 wl1" N       .. the number of superfields theta_i"$
3578 wl1" wgtlist .. an algebraic mode list of at least one alg. mode list"$
3579 wl1"            {afwlist,abwlist,wgt} where"$
3580 wl1"            afwlist .. (algebraic mode) list of weights of the"$
3581 wl1"                       fermion fields f(1), f(2), ..."$
3582 wl1"            abwlist .. (algebraic mode) list of weights of the"$
3583 wl1"                       boson   fields b(1), b(2), ..."$
3584 wl1"            wgt     .. the weight of each term in the generated"$
3585 wl1"                       polynomial according to afwlist and abwlist"$
3586 wl1"            The number of elements in afwlist, abwlist determines"$
3587 wl1"            the number of fermion and boson fields."$
3588 wl1"            If the generated polynomial should have more than one"$
3589 wl1"            homogeneity symmetry then each extra one is specified"$
3590 wl1"            by another element {afwlist,abwlist,wgt} in wgtlist."$
3591 wl1" cname   .. name of the coefficients which gets added an index"$
3592 wl1" mode    .. list of flags: "$
3593 wl1"            lin: the generated polynomial has to be linear homogeneous"$
3594 wl1"                 in all fields with (index)>(maxindex/2), i.e."$
3595 wl1"                 if lin then the number of fermions and boson fields"$
3596 wl1"                 must both be even"$
3597 wl1"            fonly: only fermionic terms are generated"$
3598 wl1"            bonly: only bosonic   terms are generated"$
3599 wl1"            {forbid,d(1,f(1)),..}: a list of fermionic and bosonic"$
3600 wl1"                 variables and their derivatives that may not occur"$
3601% wl1"            {non0coeff,p7,p18,..}: an exclusive list of coefficients"$
3602% wl2"                 that may occur"$
3603
3604 wl2"The result is a list"$
3605
3606 wl2"{fermonic_polynomial, bosonic_polynomial, {coefficients}}."$
3607
3608 wl2"Comments:"$
3609
3610 wl1"- The weight of d/dx is 2, the weight of d(i,..) is 1, no other"$
3611 wl2"  derivatives occur."$
3612
3613 wl1"- If in the first homogeneity, i.e. in the first element of wgtlist"$
3614 wl1"  a boson weight is non-positive or a fermion weight is negative"$
3615 wl1"  then the global (algebraic mode) variable max_deg must have a"$
3616 wl1"  positive integer value which is the highest degree of such a"$
3617 wl1"  field in any ansatz. Therefore if there is a set of strictly"$
3618 wl1"  positive homogeneity weights, then they should come first, needing"$
3619 wl2"  no extra restriction through max_deg which might restrict generality."$
3620>>$
3621
3622%-------
3623
3624symbolic procedure sshelp10$
3625<<wl1"The N symmetry generators theta_1,..,theta_N are called th(1),..,th(N)"$
3626 wl1"in SSTools. As they are of odd parity only 2^N many products of their"$
3627 wl1"powers exist. The name convention used in SSTools for coefficients"$
3628 wl1"in the Taylor expansion wrt. theta_i becomes apparent in the following"$
3629 wl2"two expansions (here for N=2):"$
3630
3631 wl1"f(i)=f(i,0,0)+b(i,1,0)*th(1)+b(i,0,1)*th(2)+f(i,1,1)*th(1)*th(2)"$
3632 wl2"b(i)=b(i,0,0)+f(i,1,0)*th(1)+f(i,0,1)*th(2)+b(i,1,1)*th(1)*th(2)."$
3633
3634 wl2"The theta expansion of derivatives D(i,P) is"$
3635
3636 wl2"D(i,P) = df(P,th(i))+th(i)*df(P,x)} ."$
3637
3638 wl1"These expansions can be used to convert a polynomial expression in"$
3639 wl1"terms of f(i), b(j), D(k,..) (what we call `field form') into what"$
3640 wl2"we call `coordinate form'. The call is:"$
3641
3642 wl2"  tocoo(N,nf,nb,ex)$"$
3643
3644 wl2"where"$
3645
3646 wl1" N      .. the number of superfields th(1) .. th(N)"$
3647 wl1" nf     .. number of fermion fields f(1) .. f(nf)"$
3648 wl1" nb     .. number of boson   fields b(1) .. b(nb)"$
3649 wl2" ex     .. the expression to be converted"$
3650
3651 wl2"-----------------"$
3652
3653 wl2"Example: For N=2 and 2 fermion fields the call is"$
3654
3655 wl1"    hh := tocoo(2,2,0,"$
3656 wl1"                f(1)*d(1,f(2))+d(2,d(1,DF(f(2),x)))+"$
3657 wl2"                df(f(2),x,2)*p3*p5 - df(f(1),x,3)*p3**2);"$
3658>>$
3659
3660%-------
3661
3662symbolic procedure sshelp11$
3663<<wl1"The conversion back from coordinate form to field form is done"$
3664 wl1"by the function tofield. For this function to be applicable "$
3665 wl1"the expression in coordinate form has to be homogeneous with "$
3666 wl1"suitable weights and in the current form only involve x- and "$
3667 wl2"th(i)-derivatives. The call is:"$
3668
3669 wl2"  tofield(N,nf,nb,ex,zerowei)$"$
3670
3671 wl2"where"$
3672
3673 wl1" N       .. the number of superfields theta_i"$
3674 wl1" nf      .. number of fermion fields f(1) .. f(nf)"$
3675 wl1" nb      .. number of boson   fields b(1) .. b(nb)"$
3676 wl1" ex      .. the expression to be converted"$
3677 wl1" zerowei .. list of constants or other kernels that"$
3678 wl2"            should have zero weight"$
3679
3680 wl2"-----------------"$
3681
3682 wl1"Examples:"$
3683 wl2"The following call inverts the above computation:"$
3684
3685 wl2"    tofield(2,2,0,hh,{p3,p5});"$
3686
3687 wl1"The next call shows what happens when the input does"$
3688 wl2"not have a homogeneity:"$
3689
3690 wl2"    tofield(1,1,0,f(1,0)+df(f(1,0),x)+f(1,0)*th(1)*df(f(1,0),x),{});"$
3691
3692 wl1"In the next call the input has a homogeneity but the resulting"$
3693 wl1"homogeneity weight for f(1) is negative which does not prevent"$
3694 wl1"the computation but makes it more expensive but that does not"$
3695 wl2"matter in the following short example:"$
3696
3697 wl1"    tofield(1,2,0,f(2,0)*b(1,1) + f(2,0)*th(1)*df(f(1,0),x) + "$
3698 wl2"                  f(1,0) + th(1)*b(2,1)*b(1,1) + th(1)*b(1,1),{}); "$
3699
3700 wl1"In the next call the input does have a homogeneity, the resulting"$
3701 wl2"weight for f(1) is positive but still a field form does not exist:"$
3702
3703 wl2"    tofield(1,1,0,df(f(1,0),x)+f(1,0)*th(1)*df(f(1,0),x),{});"$
3704>>$
3705
3706%-------
3707
3708symbolic operator sshelp$
3709symbolic procedure sshelp$
3710begin scalar ps,s$
3711 lines_written:=0$
3712 rds nil; wrs nil$              % Switch I/O to terminal
3713 ps:=promptstring!*$  promptstring!*:=redfront_color ""$
3714 repeat <<
3715  write"To read about the following topics input the corresponding number:"$
3716  terpri()$ terpri()$
3717  write"  Purpose                        (1)"$ terpri()$
3718  write"  Notation                       (2)"$ terpri()$
3719  write"  Initializations                (3)"$ terpri()$
3720  write"  Coefficients                   (4)"$ terpri()$
3721  write"  Symmetries                     (5)"$ terpri()$
3722  write"  Conservation laws              (6)"$ terpri()$
3723  write"  Homogeneities                  (7)"$ terpri()$
3724  write"  Linearizations                 (8)"$ terpri()$
3725  write"  Generating polynomials         (9)"$ terpri()$
3726  write"  Taylor expansions             (10)"$ terpri()$
3727  write"  Inverting Taylor expansions   (11)"$ terpri()$
3728  write"  Exit help                     (12)"$ terpri()$
3729  terpri()$
3730  write"Choice: "$
3731  s:=read()$
3732  if ifl!* then rds cadr ifl!*$  %  Resets I/O streams
3733  if ofl!* then wrs cdr ofl!*$
3734  wl2"=========================================================="$
3735  if s= 1 then sshelp1() else
3736  if s= 2 then sshelp2() else
3737  if s= 3 then sshelp3() else
3738  if s= 4 then sshelp4() else
3739  if s= 5 then sshelp5() else
3740  if s= 6 then sshelp6() else
3741  if s= 7 then sshelp7() else
3742  if s= 8 then sshelp8() else
3743  if s= 9 then sshelp9() else
3744  if s=10 then sshelp10() else
3745  if s=11 then sshelp11() else
3746  if s neq 12 then write"Wrong input, try again."$
3747 >> until s=12$
3748 promptstring!*:=ps
3749end$
3750
3751%-------
3752
3753lisp <<write"For help type:  sshelp()$  "$terpri()>>$
3754
3755endmodule;
3756
3757end$
3758
3759
3760