program coulomb;
{EatNotes+[251]}
{EatNotes+[267]}
{EatNotes+[281]}
{EatNotes+[295]}

const fa = '\frac{'; fm='}{'; fe='}'; ra='\sqrt{'; auf='\left( '; zu='\right) ';
   maxLongInteger = 2147483647;
   MachineEps	  = 2.2E-16;

type RationalNumber_T = record Z,N: integer end;

type IT  = array [1:2] of integer;
   PPT = ^PT;
   PT  = record degree: IT; {vector of degrees}
	    c	      : pointer; {to coefficient}
	    pre,post  : PPT;     {list pointers}
	 end;
   
   sparsePolynomial = record sp: ^PT;
			 level : integer
		      end;

{there are two levels of polynomials:
 1) level=0. Then the coefficients of the bivariate polynomials in X,Y are c: ^real.
 2) level=1. Then the coefficients of the bivariate polynomials are
    polynomials in x,y of level 0}

type TeXExpr_T = record
		  p0,p1,p2,pv,m0,m1,m2,mv : string;
		  vorz		          : integer
	         end;

type FctListT = ^AT;
     AT = record x,y,xx,yy: real;
                 fct,sx,sy,sxx,syy: string;
                 sp: sparsePolynomial;
                 post: FctListT
          end;
{ A variable of this type represents one on the special functions F,G,L,...
  for certain arguments.
  Fct:	Name of the function, e.g., 'F'.
  sp:	The represented function is sp(x,y;X,Y) * F(x,y;...) with the polynomial sp.
  sx:	Name of the first variable, e.g., "x".
        In the case of sx="", there is no non-constant term in the polynomial.
  x:	If sx="", x contains the value of the first variable.
  sy,y:	Similarly for the second variable.
  sxx,xx,syy,yy: Similarly for the parameters X,Y.
}

var phase,xintegration,yintegration,m12,g13,g34,  
    cem,cXY,RationalRepresentation,MaxDenominator : integer;
    msxx,msyy,fileExt,fileTex,filelog		  : string;
    protocol,ES,TX				  : text;
    Tolerance,TolRational,TreatAsRationalNumber	  : real;

procedure StartP(var p: sparsePolynomial; level: integer);
begin p.sp:=NIL; p.level:=level end; { StartP }

{ last element in the list }

function LastPT(var p: sparsePolynomial): ^PT;
var e: ^PT;
begin e:=p.sp; if e<>NIL then
   while e^.post<>NIL do e:=e^.post;
   LastPT:=e
end;

function LargestPT(var p: sparsePolynomial; xy: integer): ^PT;
{ chooses the term with maximal degree *^.degree[xy] }
var e,emax: ^PT; dmax: integer;
begin dmax:=-1; emax:=NIL; e:=p.sp;
   while e<>NIL do 
   begin if e^.degree[xy]>=dmax then begin emax:=e; dmax:=e^.degree[xy] end;
      e:=e^.post
   end;
   LargestPT:=emax
end;

{ ordering of the pairs of degrees ...................................}

function CompareI2(var t1,t2: IT): integer;
{ Result<0: t1<t2 % Result=0: t1=t2 % Result>0: t1>t2 }
var i: integer;
begin i:=t1[2]-t2[2]; if i=0 then i:=t1[1]-t2[1];
   CompareI2 := i
end; { CompareI2 }

{ search of term of degree t in the polynomial .......................}

function FindDegree_general(var p: sparsePolynomial; var t: IT;
                            var e: ^PT): integer;
{e=NIL: p is zero polynomial, i.e., p.sp=NIL (then Result<0).
 Otherwise,
 Result=0: e^.degree = t,
           i.e., e^.c contains the coefficient for degree t.
 Result<0: term for degree t does not exist.
           It can be inserted before "e".
 Result>0: term for degree t does not exist.
           It can be inserted as last element after "e" (e=NIL possible)
 }
var l: ^PT; i: integer; label 1,9;
begin l:=p.sp; i:=1; if l=NIL then goto 9;
1: i:=CompareI2(t,l^.degree);
   if i<=0 then goto 9;
   if l^.post=NIL then goto 9; {then i>0}
   l:=l^.post; goto 1;
9: e:=l; FindDegree_general:=i
end; { FindDegree_general }

function FindDegree(var p: sparsePolynomial; var t: IT): pointer;
{Result=NIL: coefficient corresponding to degree t is zero.
 Otherwise, the Result is the pointer to the coefficient  }
var e: ^PT;
begin if FindDegree_general(p,t,e)<>0 then FindDegree:=NIL else
   FindDegree:=e^.c
end; { FindDegree }

function real_from_pointer(c: pointer): real;
{ assumption: c<>NIL and c points to a real number }
var pr: ^real;
begin pr:=c; real_from_pointer:=pr^ end;

function sparsePolynomialPointer_from_pointer(c: pointer): ^sparsePolynomial;
{ assumption: c<>NIL and c points to a real number }
var pp: ^sparsePolynomial;
begin pp:=c; sparsePolynomialPointer_from_pointer:=pp end;

{ joining two list elements }
procedure join(var p: sparsePolynomial; e1,e2: ^PT);
begin if e2<>NIL then e2^.pre:=e1; 
      if e1=NIL then p.sp:=e2 else e1^.post:=e2
end;

procedure writeP(var s: string);
begin write(protocol,s) end;

procedure writelnP;
begin writeln(protocol) end;

procedure write_P(var s: string);
begin write(s); writeP(s) end;

procedure writeln_P;
begin writeln; writelnP end;

function RationalFromReal(var r: RationalNumber_T; x,eps: real): real;
{Change x -> r, so that |x-r|<=eps. Result = x-r}
var d,dabs,dmin,dabsmin: real; z,n,zmin,nmin: integer; label 8,9;
begin dabsmin:=2; dmin:=-1;
   for n:=1 to MaxDenominator do
   begin d:=x*n; if abs(d)>maxLongInteger then goto 8;
      z:=round(d); d:=x-z/n; dabs:=abs(d);
      if dabs<=eps then goto 9;
      if dabs<dabsmin then begin dabsmin:=dabs; dmin:=d; zmin:=z; nmin:=n end
   end;
8: if dmin<0 then begin n:=1; z:=0; d:=abs(x); goto 9 end;
   z:=zmin; n:=nmin; d:=dmin;
9: r.n:=n; r.z:=z; RationalFromReal:=d
end; { RationalFromReal }

function RationalFromReal0(var r: RationalNumber_T; x: real): real;
{Umwandlung x -> r. Wert = x-r}
begin RationalFromReal0:=RationalFromReal(r,x,TolRational)
end; { RationalFromReal0 }

procedure RealWithSign(var s: string; var r: real);
begin if r>=0 then writev(s,' +',r:25:22) else writev(s,' ',r:25:22) end;

procedure WriteRationalNumber(var r: RationalNumber_T);
var s: string;
begin if r.n=1 then writev(s,'(',r.z:1,')')
         else writev(s,'(',r.z:1,'/',r.n:1,')');
      write_P(s)
end; { WriteRationalNumber }

procedure WriteReal(x: real);
var s: string; r: RationalNumber_T; d: real; label 2;
begin if abs(x)<=Tolerance then x:=0;
   if RationalRepresentation=1 then
   begin if x<>0 then if abs(x) < TreatAsRationalNumber then goto 2;
      d:=RationalFromReal0(r,x); if d>TolRational then goto 2;
         s:=' + '; write_P(s); WriteRationalNumber(r)
   end else
   begin 2: RealWithSign(s,x); write_P(s) end
end;

procedure ExpSymbol(k,l,level: integer; x1,x2,c1,c2: string);
var s: string;
begin if level=0 then begin x1:=c1; x2:=c2 end;
   {writev(s,'*',x1,'^',k:1,'*',x2,'^',l:1); write_P(s);}
   if abs(k)+abs(l)>0 then begin s:='*'; write_P(s) end;
   if k=1 then write_P(x1) else
      if k<>0 then begin writev(s,x1,'^',k:1); write_P(s) end;
   if k*l<>0 then begin s:='*'; write_P(s) end;
   if l=1 then write_P(x2) else
      if l<>0 then begin writev(s,x2,'^',l:1); write_P(s) end
end;

procedure WritePolynomial(title: string; 
                         var p: sparsePolynomial; x1,x2,c1,c2: string);
var l: ^PT; p0: ^sparsePolynomial; pc: ^real; 
    r: RationalNumber_T; s: string; nr: integer;
begin nr:=0;
   if title<>'' then write_P(title);
   if p.level=1 then begin writeln_P; s:='{'; write_P(s) end;
   if p.sp=NIL then
   begin s:='ZERO'; write_P(s) end else
   begin l:=p.sp;
      while l<>NIL do
      begin {output of coefficient}
	 if p.level=0 then
	 begin WriteReal(real_from_pointer(l^.c)); nr:=nr+1 end
         else {if p.level=1 then}
	 begin p0:=l^.c;
            if nr=0 then s:='  [' else begin writeln_P; s:=' + [' end;
	    write_P(s);
            WritePolynomial('',p0^,c1,c2,c1,c2); writev(s,']');
            write_P(s); nr:=nr+1;
	 end;
	 ExpSymbol(l^.degree[1],l^.degree[2],p.level,x1,x2,c1,c2);
	 l:=l^.post
      end
   end;
   if p.level=1 then begin writev(s,' } * '); write_P(s) end
end; { WritePolynomialRational }

{ procedures for deleting list elements }

procedure DeletePolynomialTerm_direct(var p: sparsePolynomial; e: ^PT);
forward;

procedure DeletePolynomial(var p: sparsePolynomial);
begin while p.sp<>NIL do DeletePolynomialTerm_direct(p,p.sp);
end; { DeletePolynomial }

procedure DeletePolynomialPointer(p: ^sparsePolynomial);
begin if p<>NIL then begin DeletePolynomial(p^); dispose(p) end end;

procedure DeleteCoefficient(c: pointer; var level: integer);
var pp: ^sparsePolynomial; pr: ^real;
begin if c<>NIL then
   begin
      if level=1 then begin pp:=c; DeletePolynomialPointer(pp) end
      else begin pr:=c; dispose(pr) end
   end
end; { DeleteCoefficient }

procedure DeletePT(e: ^PT; level: integer);
begin if e<>NIL then DeleteCoefficient(e^.c,level); dispose(e)
end; { DeletePT }

procedure ExcludePolynomialTerm_direct(var p: sparsePolynomial; e: ^PT);
begin if e<>NIL then join(p,e^.pre,e^.post) end;

procedure DeletePolynomialTerm_direct(var p: sparsePolynomial; e: ^PT);
begin ExcludePolynomialTerm_direct(p,e); DeletePT(e,p.level) end;

procedure DeletePolynomialTerm(var p: sparsePolynomial; var t: IT);
var e: ^PT;
begin if FindDegree_general(p,t,e)=0 then DeletePolynomialTerm_direct(p,e)
end; { DeletePolynomialTerm }

procedure InsertPT_direct(var p: sparsePolynomial; e,term: ^PT; i: integer);
{Assumptions: 1) p does not contain a term for the coefficient term^.degree.
              2) term^.c<>NIL
              3) i:=FindDegree_general(p,term^.degree,e) yields i,e
                 as in the parameterlist. Note that i<>0 by assmption 1 }
begin 
   if i<0 then {insert before e}
      begin join(p,e^.pre,term); join(p,term,e) end
         else  {insert after e as last entry}
            begin join(p,e,term); join(p,term,NIL) end
end; { InsertPT_direct }

procedure DefineCoefficient(var p: sparsePolynomial; t1,t2: integer; c: pointer);
{c^ is existing coefficient corresponding to degree t=(t1,t2).
 If p already contains this term, the coefficient 
 pointer is overwritten by c (no copy made!).
 Otherwise, a new term is added to the list.
 If c=NIL, a possibly existing term is deleted}
var t: IT; i: integer; e,v,ep: ^PT; label 5;
begin t[1]:=t1; t[2]:=t2;
   i:=FindDegree_general(p,t,e);
   if i=0 then {term existing}
   begin if c<>NIL then begin DeleteCoefficient(e^.c,p.level); e^.c:=c end
         else DeletePolynomialTerm_direct(p,e)
   end else
      {term not existing}
      if c<>NIL then {insert new term}
      begin
	 new(v); v^.degree:=t; v^.c:=c;
	 InsertPT_direct(p,e,v,i)
      end;
end; { DefineCoefficient }

procedure DefineCoefficient0_real(var p: sparsePolynomial; t1,t2: integer; c: real);
{ Case of p.level=0. The coefficient is given as real number (not its address!) }
var pc: ^real;
begin new(pc); pc^:=c; DefineCoefficient(p,t1,t2,pc)
end; { DefineCoefficient0_real }

procedure DefineCoefficient1_real(var p: sparsePolynomial; t1,t2: integer; c: real);
{ Case of p.level=1. The coefficient is given as constant polynomial = c }
var pp: ^sparsePolynomial;
begin new(pp); StartP(pp^,0); 
    DefineCoefficient0_real(pp^,0,0,c); DefineCoefficient(p,t1,t2,pp)
end; { DefineCoefficient1_real }

procedure DefineCoefficient_real(var p: sparsePolynomial; t1,t2: integer; c: real);
begin if p.level=0 then 
    DefineCoefficient0_real(p,t1,t2,c) else DefineCoefficient1_real(p,t1,t2,c)
end;

function CopyPT(term: ^PT; level: integer): ^PT; forward;

function CopyPolynomial(var p: sparsePolynomial): ^sparsePolynomial;
{ CopyPolynomial^ is a copy of the polynomial p }
var q: ^sparsePolynomial; l,e,ep: ^PT;
begin new(q); q^.level:=p.level; q^.sp:=NIL; l:=p.sp; e:=NIL;
   while l<>NIL do
   begin ep:=CopyPT(l,p.level); join(q^,e,ep); ep^.post:=NIL;
         e:=ep; l:=l^.post
   end;
   CopyPolynomial:=q
end; { CopyPolynomial }

procedure CopyPolynomialsp(var pcopy,p: sparsePolynomial);
{ pcopy becomes a copy of the polynomial p }
var l,e,ep: ^PT;
begin pcopy.level:=p.level; pcopy.sp:=NIL; l:=p.sp; e:=NIL;
   while l<>NIL do
   begin ep:=CopyPT(l,p.level); join(pcopy,e,ep); ep^.post:=NIL;
         e:=ep; l:=l^.post
   end
end; { CopyPolynomialsp }

function CopyPT(term: ^PT; level: integer): ^PT;
{ term^ is copied into CopyPT^ }
var v: ^PT; pr: ^real;
begin new(v); v^.degree:=term^.degree;
   if level=0 then
   begin new(pr); pr^:=real_from_pointer(term^.c); v^.c:=pr end
   else v^.c:=CopyPolynomial(sparsePolynomialPointer_from_pointer(term^.c)^);
   CopyPT:=v
end; { CopyPT }

procedure WriteExPolynomial(var p: sparsePolynomial); forward;
function ReadExPolynomialP: ^sparsePolynomial; forward;

procedure WriteExPT(e: ^PT; level: integer);
{ assumption: e<>NIL and e^.c<>NIL }
begin writeln(ES,'degree'); writeln(ES,e^.degree[1],e^.degree[2]);
   if level=0 then writeln(ES,real_from_pointer(e^.c):25:22)
   else WriteExPolynomial(sparsePolynomialPointer_from_pointer(e^.c)^)   
end;

function ReadExPT(level: integer): ^PT;
var s: string; e: ^PT; r: ^real; p: ^sparsePolynomial; label 9;
begin readln(ES,s); if s='NIL' then begin e:=NIL; goto 9 end;
   new(e); readln(ES,e^.degree[1],e^.degree[2]);
   if level=0 then begin new(r); readln(ES,r^); e^.c:=r; goto 9 end;
   e^.c:=ReadExPolynomialP;
9: ReadExPT:=e
end;

procedure WriteExPolynomial(var p: sparsePolynomial);
var e: ^PT;
begin writeln(ES,p.level:1,'   level of polynomial'); e:=p.sp;
   while e<>NIL do begin WriteExPT(e,p.level); e:=e^.post end;
   writeln(ES,'NIL');
end;

function ReadExPolynomialP: ^sparsePolynomial;
var e,l: ^PT; s: string; p: ^sparsePolynomial; label 1,9;
begin l:=NIL; new(p); readln(ES,p^.level); p^.sp:=NIL;
1: e:=ReadExPT(p^.level); if e=NIL then goto 9;
   InsertPT_direct(p^.sp,l,e,1); l:=e; goto 1;
9: ReadExPolynomialP:=p
end;

procedure ReadExPolynomial(var p: sparsePolynomial);
var e,l: ^PT; s: string; label 1,9;
begin l:=NIL; readln(ES,p.level); p.sp:=NIL;
1: e:=ReadExPT(p.level); if e=NIL then goto 9;
   InsertPT_direct(p.sp,l,e,1); l:=e; goto 1;
9:
end;

procedure Polynomial0_input(var p: sparsePolynomial);
{ assumption: p.level=0 }
var k,l: integer; c: real; label 1,5,9;
begin if p.level<>0 then goto 9;
   writeln('##### input of coefficient polynomial ...');
1: writeln('      choose the degrees k,l   < terminate with k<0 >');
   write('      -> k = '); readln(k); if k<0 then
5: begin writeln('      value must be non-negative'); goto 9 end;
   write('      -> l = '); readln(l); if l<0 then goto 5;
   write('      -> value of coefficient[',k:1,',',l:1,'] = '); readln(c);
   DefineCoefficient0_real(p,k,l,c); goto 1;
9: writeln('##### input of coefficient polynomial finished')
end;

procedure Polynomial1_input(var p: sparsePolynomial);
{ assumption: p.level=1 }
var k,l,a: integer; c: real; q: ^sparsePolynomial; label 1,5,9;
begin if p.level<>1 then goto 9;
   writeln('### input of polynomial ...');
1: writeln('    choose the degrees k,l   < terminate with k<0 >');
   write('    -> k = '); readln(k); if k<0 then
5: begin writeln('    value must be non-negative'); goto 9 end;
   write('    -> l = '); readln(l); if l<0 then goto 5;
   writeln('    choose "0" for real coefficient (constant polynomial) or');
   write('           "1" for non-constant polynomial. Choice = ');
   readln(a); if a=0 then
   begin write('    -> value of coefficient[',k:1,',',l:1,'] = '); readln(c);
         DefineCoefficient1_real(p,k,l,c); goto 1
   end;
   new(q); StartP(q^,0); Polynomial0_input(q^);
   DefineCoefficient(p,k,l,q); goto 1;
9: writeln('### input of polynomial finished')
end;

procedure P_plus_Q(var p,q: sparsePolynomial); forward;

procedure AddTerm(var p: sparsePolynomial; term: ^PT);
{ p := p + term}
var e,term2: ^PT; i: integer; r: real;
   pr1,pr2: ^real; pp1,pp2: ^sparsePolynomial; label 9;
begin if term=NIL then goto 9; if term^.c=NIL then goto 9;
   i:=FindDegree_general(p,term^.degree,e);
   if i=0 then { e and term belong to the same degree }
   begin if p.level=0 then {sum of reals}
      begin pr1:=e^.c; pr2:=term^.c; pr1^:=pr1^+pr2^ end
      else {p.level=1, sum of polynomials}
      begin pp1:=e^.c; pp2:=term^.c; P_plus_Q(pp1^,pp2^) end
   end else { i<>0 }
   begin
      term2:=CopyPT(term,p.level); InsertPT_direct(p,e,term2,i)
   end;
9:
end; { AddTerm }

procedure P_plus_Q(var p,q: sparsePolynomial);
{ p := p + q}
var l: ^PT;
begin l:=q.sp;
   while l<>NIL do
   begin AddTerm(p,l); l:=l^.post end
end; { P_plus_Q }


procedure SparsifyPolynomial(var p: sparsePolynomial);
var e,ep: ^PT; pp: ^sparsePolynomial; label 5;
begin e:=p.sp;
   while e<>NIL do
   begin ep:=e^.post; if e^.c=NIL then goto 5;
      if p.level=0 then
      begin if abs(real_from_pointer(e^.c))<=Tolerance then
5:	 DeletePolynomialTerm_direct(p,e);
      end
      else { p.level=1 }
      begin pp:=e^.c; SparsifyPolynomial(pp^); if pp^.sp=NIL then goto 5 end;
      e:=ep
   end
end; { SparsifyPolynomial }

procedure sxP(var p: sparsePolynomial; s: real); forward;

procedure sxTerm(term: ^PT; s: real; level: integer);
{ term := term * s ; assumption: term<>NIL }
var pr: ^real;
begin if level=0 then begin pr:=term^.c; pr^:=s*pr^ end
   else {p.level=1} sxP(sparsePolynomialPointer_from_pointer(term^.c)^,s);
end;

procedure sxP(var p: sparsePolynomial; s: real);
{ p := s * p }
var l: ^PT;
begin if s=0 then
   begin DeletePolynomial(p); p.sp:=NIL end 
   else if s=1 then {no action}
   else {s<>0,<>1}
   begin l:=p.sp;
      while l<>NIL do begin sxTerm(l,s,p.level); l:=l^.post end
   end
end;

procedure P_minus_Q(var p,q: sparsePolynomial); forward;

procedure SubtractTerm(var p: sparsePolynomial; term: ^PT);
 { p := p - term}
var e,term2: ^PT; i: integer; r: real;
   pr1,pr2: ^real; pp1,pp2: ^sparsePolynomial; label 9;
begin if term=NIL then goto 9; if term^.c=NIL then goto 9;
   i:=FindDegree_general(p,term^.degree,e);
   if i=0 then { e and term belong to the same degree }
   begin if p.level=0 then {difference of reals}
      begin pr1:=e^.c; pr2:=term^.c; pr1^:=pr1^-pr2^ end
      else {p.level=1, sum of polynomials}
      begin pp1:=e^.c; pp2:=term^.c; P_minus_Q(pp1^,pp2^) end
   end else { i<>0 }
   begin term2:=CopyPT(term,p.level); sxTerm(term2,-1,p.level);
         InsertPT_direct(p,e,term2,i)
   end;
9:
end; { SubtractTerm }

procedure P_minus_Q(var p,q: sparsePolynomial);
{ p := p - q}
var l: ^PT;
begin l:=q.sp;
   while l<>NIL do
   begin SubtractTerm(p,l); l:=l^.post end
end; { P_minus_Q }

procedure ChangeTerm(term: ^PT; dk,dl: integer; s: real; level: integer);
{ Term is shifted by (dk,dl) and multiplied by s; assumption: term<>NIL }
begin sxTerm(term,s,level); with term^ do begin
   degree[1]:=degree[1]+dk; if degree[1]<0 then
      writeln('### ERROR in ChangeTerm: degree[1]=',degree[1]:1,' < 0.  dk=',dk:1);
   degree[2]:=degree[2]+dl; if degree[2]<0 then
      writeln('### ERROR in ChangeTerm: degree[2]=',degree[2]:1,' < 0.  dl=',dl:1)
end end;

procedure ChangePolynomial(var p: sparsePolynomial; dk,dl: integer; s: real);
{ All terms shifted by (dk,dl) and multiplied by s }
var l: ^PT;
begin l:=p.sp;
   while l<>NIL do begin ChangeTerm(l,dk,dl,s,p.level); l:=l^.post end
end;

procedure RedefinePolynomial(p,p_original: sparsePolynomial); forward;

procedure RedefineTerm(term,original: ^PT; level: integer);
{ term := original; assumption: term<>NIL, original<>NIL }
var pr: ^real;
begin term^.degree[1]:=original^.degree[1]; term^.degree[2]:=original^.degree[2];
   if level=0 then
   begin pr:=term^.c; pr^:=real_from_pointer(original^.c) end 
   else {p.level=1}
   RedefinePolynomial(sparsePolynomialPointer_from_pointer(term^.c)^,
                      sparsePolynomialPointer_from_pointer(original^.c)^)
end;

procedure RedefinePolynomial(p,p_original: sparsePolynomial);
{ assumption: isomorphic structures of p,p_original }
var l,ll: ^PT;
begin l:=p.sp; ll:=p_original.sp;
   while ll<>NIL do 
   begin RedefineTerm(l,ll,p.level); l:=l^.post; ll:=ll^.post end
end;

{ combined term operation }

procedure TermOp(var p: sparsePolynomial; term,term_original: ^PT; 
                 dk,dl: integer; s: real; pk,pl: integer);
{ Assumption: term = term_original at starting time.
  p := p + s*term, where
  a) the degree of term is the one of term_original shifted by (dk,dl),
  b) the polynomial coefficient is shifted by pk,pl.
  At the end, term is redefined by term_original. }
var e,term2: ^PT; i: integer; pr1,pr2: ^real; pp1,pp2: ^sparsePolynomial; label 9;
begin if term=NIL then goto 9; if s=0 then goto 9; if term^.c=NIL then goto 9;
   ChangeTerm(term,dk,dl,s,p.level);
   if (pk>0) or (pl>0) then begin pp1:=term^.c; ChangePolynomial(pp1^,pk,pl,1) end;
   i:=FindDegree_general(p,term^.degree,e);
   if i=0 then
   begin if p.level=0 then {sum of reals}
      begin pr1:=e^.c; pr2:=term^.c; pr1^:=pr1^+pr2^ end
      else {p.level=1, sum of polynomials}
      begin pp1:=e^.c; pp2:=term^.c; P_plus_Q(pp1^,pp2^) end
   end else { i<>0 }
   begin
      term2:=CopyPT(term,p.level); InsertPT_direct(p,e,term2,i)
   end;
   RedefineTerm(term,term_original,p.level);
9:
end; { TermOp }

procedure multiply(var product,p1,p2: sparsePolynomial); forward;

function multiplyTT(t1,t2: ^PT; level: integer): ^PT;
{ RESULT^ := t1^ * t2^ }
var e: ^PT; pr: ^real; pp: ^sparsePolynomial; label 9;
begin e:=NIL; if t1=NIL then goto 9; if t2=NIL then goto 9; new(e);
   e^.degree[1]:=t1^.degree[1]+t2^.degree[1];
   e^.degree[2]:=t1^.degree[2]+t2^.degree[2];
   if level=0 then {product of real numbers}
   begin new(pr);
      pr^:=real_from_pointer(t1^.c)*real_from_pointer(t2^.c);
      e^.c:=pr
   end else { level = 1 }
   begin new(pp); StartP(pp^,0);
      multiply(pp^,sparsePolynomialPointer_from_pointer(t1^.c)^,
                   sparsePolynomialPointer_from_pointer(t2^.c)^);
      e^.c:=pp
   end;
9: multiplyTT:=e
end;

procedure multiply(var product,p1,p2: sparsePolynomial);
{ product := product + p1 * p2;
  assumption: all polynomials are of the same level }
var l1,l2,e: ^PT;
begin l2:=p2.sp; while l2<>NIL do
   begin l1:=p1.sp; while l1<>NIL do
      begin e:=multiplyTT(l1,l2,p1.level);
         AddTerm(product,e); DeletePT(e,p1.level);
         l1:=l1^.post
      end;
      l2:=l2^.post
end end;

procedure multiply10(var p,q1: sparsePolynomial; d1,d2: integer; q2: ^PT);
{ p := p + q1 * q, where q = q2 * x^d1 * y^d2; q2 afterwards deleted.
  assumption: p.level=q1.level=1 and level=0 for q2 }
var q,qq: ^sparsePolynomial; e: ^PT;
begin
   new(qq); StartP(qq^,0); insertPT_direct(qq^,NIL,q2,1);
   new(q); StartP(q^,1);
   new(e); e^.degree[1]:=d1; e^.degree[2]:=d2; e^.c:=qq;
   insertPT_direct(q^,NIL,e,1);
   multiply(p,q1,q^);
   DeletePolynomialPointer(q) {deletes also qq,ee,q2}
end;

procedure replacexy(var p: sparsePolynomial; 
                    VariableNr: integer; var q: sparsePolynomial);
{ assumption: p.level=q.level, 1<=VariableNr<=2;
     The "VariableNr"-th variable (x or y) is to be replaced by the
     polynomial q, which must be constant w.r.t. the replaced variable }
var e,ec: ^PT; pe,pr: sparsePolynomial; label 1,9;
begin e:=LargestPT(q,VariableNr);
   if e<>NIL then if e^.degree[VariableNr]>0 then
   begin writeln('#### Error in REPLACE: q contains non-constant terms in the variable Nr ',
                  VariableNr:1); goto 9
   end;
   StartP(pe,p.level); StartP(pr,p.level);
1: e:=LargestPT(p,VariableNr); if e=NIL then goto 9;
   if e^.degree[VariableNr]=0 then goto 9; {variable x[VariableNr] replaced }
   DeletePolynomial(pe); DeletePolynomial(pr);
   ec:=CopyPT(e,p.level); DeletePolynomialTerm_direct(p,e);
   ec^.degree[VariableNr]:=ec^.degree[VariableNr]-1;
   InsertPT_direct(pe,NIL,ec,1); multiply(pr,pe,q); P_plus_Q(p,pr);
   goto 1;
9:
end;

procedure replace_XY(var p: sparsePolynomial; 
                     VariableNr: integer; var q: sparsePolynomial);
{ assumption: p.level=1, q.level=0, 3<=VariableNr<=4;
     The "VariableNr"-th variable (X or Y) is to be replaced by the
     polynomial q, which must be constant w.r.t. the replaced variable }
var e: ^PT;
begin VariableNr:=VariableNr-2; e:=p.sp;
   while e<>NIL do
   begin replacexy(sparsePolynomialPointer_from_pointer(e^.c)^,VariableNr,q);
         e:=e^.post
   end
end;

procedure replacexy_XY(var p: sparsePolynomial; 
                       VariableNr: integer; var q: sparsePolynomial);
{ assumption: p.level=q.level=1, 3<=VariableNr<=4;
     The "VariableNr"-th variable (X or Y) is to be replaced by the
     polynomial q, which must be constant w.r.t. the replaced variable }
var t,l: ^PT; d,v: integer; pp: ^sparsePolynomial; label 5;
begin v:=VariableNr-2;
   repeat
      d:=0; { d = largest degree w.r.t. VariableNr }
      l:=p.sp;
      while l<>NIL do
      begin pp:=l^.c; if pp=NIL then goto 5; 
         t:=LargestPT(pp^,v); if t=NIL then goto 5; if t^.degree[v]=0 then goto 5;
         ExcludePolynomialTerm_direct(pp^,t);
         t^.degree[v]:=t^.degree[v]-1; if d<t^.degree[v] then d:=t^.degree[v];
         multiply10(p,q,l^.degree[1],l^.degree[2],t); { this deletes t }
5:       l:=l^.post
      end;
      sparsifyPolynomial(p)
   until d=0
end;

procedure ReplaceByLinear(var p: sparsePolynomial; 
                          VariableNr: integer; a,c: real);
{ as "replace" but with the linear polynomial q=a*z+c, 
  where z is the other variable, i.e. y if VariableNr=1, x if VariableNr=2 }
var q: sparsePolynomial;
begin StartP(q,p.level);
   if c<>0 then DefineCoefficient_real(q,0,0,c);
   if a<>0 then if VariableNr=1 then
   DefineCoefficient_real(q,0,1,a) else DefineCoefficient_real(q,1,0,a);
   if VariableNr>2 then replace_XY(p,VariableNr,q) else replacexy(p,VariableNr,q);
   DeletePolynomial(q)
end;

procedure ReplaceByQuadratic(var p: sparsePolynomial; 
                             VariableNr: integer; a,b,c: real);
{ as "ReplaceByLinear" but with the quadratic polynomial q=a*z^2+b*z+c, 
  where z is the other variable, i.e. y if VariableNr=1, x if VariableNr=2 }
var q: sparsePolynomial; V,level: integer;
begin level:=p.level; V:=VariableNr; if V>2 then begin level:=0; V:=V-2 end;
   StartP(q,level);
   if a<>0 then
   begin if V=1 then DefineCoefficient_real(q,0,2,a)
       else DefineCoefficient_real(q,2,0,a)
   end;
   if b<>0 then
   begin if V=1 then DefineCoefficient_real(q,0,1,a)
       else DefineCoefficient_real(q,1,0,a)
   end;
   if c<>0 then DefineCoefficient_real(q,0,0,c);
   if VariableNr>2 then replace_XY(p,VariableNr,q) else replacexy(p,VariableNr,q);
   DeletePolynomial(q)
end;

procedure reinterpretationPolynomial(var p: sparsePolynomial; jxy,jj: integer);
{ assumption: 3<=jxy,jj<=4, p.level=1, p constant w.r.t. x,y.
  Each term c * X_jxy^a * X_j^b (j:=7-jxy) becomes
            c * (x-y)^a * X_jj^b }
var s: string; l1,l0: ^PT; q: ^sparsePolynomial; i: integer;
begin if jxy=jj then { interchange X <-> Y }
   begin l1:=p.sp; while l1<>NIL do
      begin q:=l1^.c; if q<>NIL then
         begin l0:=q^.sp; while l0<>NIL do
            begin i:=l0^.degree[1]; l0^.degree[1]:=l0^.degree[2]; l0^.degree[2]:=i
            end;
            l0:=l0^.post
         end;
         l1:=l1^.post
      end;
      jxy:=7-jxy
   end; { X <-> Y interchanged}
   new(q); StartP(q^,1); { q^ := x-y }
   DefineCoefficient1_real(q^,1,0,1); DefineCoefficient1_real(q^,0,1,-1);
   replacexy_XY(p,jxy,q^);
   DeletePolynomialPointer(q)
end;

function interchangePolynomial0(p: ^sparsePolynomial): ^sparsePolynomial;
{ assumption: level=0. The variables X and Y are interchanged }
var t: ^PT; d: integer; q: ^sparsePolynomial;
begin if p<>NIL then
    begin t:=p^.sp; while t<>NIL do with t^ do
      begin d:=degree[1]; degree[1]:=degree[2]; degree[2]:=d;
            t:=post
      end;
      q:=CopyPolynomial(p^); DeletePolynomialPointer(p); p:=q
    end;
    interchangePolynomial0:=p
end;

procedure interchangeXYPolynomial(var p: sparsePolynomial);
{ assumption: level=1. The variables X and Y are interchanged }
var t: ^PT; d: integer; q: ^sparsePolynomial;
begin t:=p.sp; while t<>NIL do
    begin t^.c:=interchangePolynomial0(sparsePolynomialPointer_from_pointer(t^.c));
          t:=t^.post
end end;

function constant_from_Polynomial(var p: sparsePolynomial): real;
{ p is assumed to be a constant polynomial. The constant is the return value }
var z: real; q: ^sparsePolynomial; a: ^PT; label 9;
begin z:=0; if p.sp=NIL then goto 9;
    a:=p.sp; if p.level=1 then
    begin if a^.degree[1]+a^.degree[2]>0 then
          writeln('ERROR in constant_from_Polynomial: Polynomial not constant in x,y');
          q:=a^.c; if q=NIL then goto 9; a:=q^.sp;
          if a^.degree[1]+a^.degree[2]>0 then
          writeln('ERROR in constant_from_Polynomial: Polynomial not constant in X,Y');
    end;
   if a<>NIL then z:=real_from_pointer(a^.c);
9: constant_from_Polynomial:=z
end;

function evaluatePolynomial(var p: sparsePolynomial; x,y,xx,yy: real): real;
var z: real; q,qq: ^sparsePolynomial; m: sparsePolynomial; label 9;
begin z:=0; if p.sp=NIL then goto 9;
   q:=CopyPolynomial(p); StartP(m,1);
   DefineCoefficient1_real(m,0,0,x);  replacexy(q^,1,m);  { variable "x" evaluated by x }
   DefineCoefficient1_real(m,0,0,y);  replacexy(q^,2,m);  { variable "y" evaluated by y }
   DeletePolynomial(m); StartP(m,0);
   DefineCoefficient0_real(m,0,0,xx); replace_XY(q^,3,m); { variable "X" evaluated by xx }
   DefineCoefficient0_real(m,0,0,yy); replace_XY(q^,4,m); { variable "Y" evaluated by yy }
   DeletePolynomial(m);
   z:=constant_from_Polynomial(q^); DeletePolynomialPointer(q);
9: evaluatePolynomial:=z
end;

{ quadrature rules }

procedure Fx(var G,L,FF: sparsePolynomial; ix,iy: integer; x,y: real);
{input are FF and initial values of G,L.
 ix=0: x is the value of the parameter X, otherwise X is symbolic.
 iy=0: y is the value of the parameter Y, otherwise Y is symbolic.
 The x-antiderivative of FF yields G,L-contributions which are added.
 The data FF are overwritten by the zero polynomial. }
var lF,e,ec: ^PT; k: integer; f: real; label 1,9;
begin
1: lF:=LargestPT(FF,1); if lF=NIL then goto 9;
   e:=CopyPT(lF,FF.level); ec:=CopyPT(e,FF.level);
   DeletePolynomialTerm_direct(FF,lF);
   k:=e^.degree[1];
   if k>1 then
   begin f:=(1-k)/k;
         TermOp(G, ec,e,-1,0,1/k,     0,0);
         TermOp(FF,ec,e,-1,1,1-f,     0,0);
         TermOp(FF,ec,e,-2,2,f,       0,0);
      if ix=0 then
         TermOp(FF,ec,e,-2,0,f,       2,0)  { symbolic variable X } else
      if x<>0 then
         TermOp(FF,ec,e,-2,0,f*sqr(x),0,0); { X evaluated by x }
      if iy=0 then
         TermOp(FF,ec,e,-2,0,f,       0,2)  { symbolic variable Y } else
      if y<>0 then
         TermOp(FF,ec,e,-2,0,f*sqr(y),0,0)  { Y evaluated by y }
   end else                              {cf. (5.3)}
   if k=0 then AddTerm(L,e) else         {cf. (5.1)}
   if k=1 then 
   begin TermOp(G,ec,e,-1,0,1,0,0);
         TermOp(L,ec,e,-1,1,1,0,0)
   end;                                  {cf. (5.2)}
   DeletePT(e,FF.level); DeletePT(ec,FF.level);
   goto 1;
9:
end; { Fx }

procedure Fy(var G,L,FF: sparsePolynomial; ix,iy: integer; x,y: real);
{input are FF and initial values of G,L.
 The y-antiderivative of FF yields G,L-contributions which are added.
 The data FF are overwritten by the zero polynomial }
var lF,e,ec: ^PT; ll: integer; f: real; label 1,9;
begin 
1: lF:=LargestPT(FF,2); if lF=NIL then goto 9;
   e:=CopyPT(lF,FF.level); ec:=CopyPT(e,FF.level); 
   DeletePolynomialTerm_direct(FF,lF);
   ll:=e^.degree[2]; if ll>1 then
   begin f:=(1-ll)/ll;
         TermOp(G, ec,e,0,-1,1/ll,    0,0);
         TermOp(FF,ec,e,1,-1,1-f,     0,0);
         TermOp(FF,ec,e,2,-2,f,       0,0);
      if ix=0 then
         TermOp(FF,ec,e,0,-2,f,       2,0)  { symbolic variable X } else
      if x<>0 then
         TermOp(FF,ec,e,0,-2,f*sqr(x),0,0); { X evaluated by x }
      if iy=0 then
         TermOp(FF,ec,e,0,-2,f,       0,2)  { symbolic variable Y } else
      if y<>0 then
         TermOp(FF,ec,e,0,-2,f*sqr(y),0,0)  { Y evaluated by y }
   end else                               {cf. (5.6)}
   if ll=0 then SubtractTerm(L,e) else    {cf. (5.4)}
   if ll=1 then 
   begin TermOp(G,ec,e,0,-1,1, 0,0);
         TermOp(L,ec,e,1,-1,-1,0,0)
   end;                                   {cf. (5.4)}
   DeletePT(e,FF.level); DeletePT(ec,FF.level);
   goto 1;
9:
end; { Fy }

procedure Gx_Fx(var G,Fx,GG: sparsePolynomial);
{input are GG and initial values of G,Fx.
 The x-antiderivative of GG yields G,Fx-contributions which are added.
 The data GG are overwritten by the zero polynomial }
var lF,e,ec: ^PT; fac: real;
begin while GG.sp<>NIL do
   begin lF:=GG.sp;
      e:=CopyPT(lF,GG.level); ec:=CopyPT(e,GG.level); 
      DeletePolynomialTerm_direct(GG,lF);
      fac:=1/(e^.degree[1]+1);
      TermOp(G, ec,e,1,0, fac,0,0);
      TermOp(Fx,ec,e,2,0,-fac,0,0);
      TermOp(Fx,ec,e,1,1, fac,0,0);          {cf. (5.8)}
      DeletePT(e,GG.level); DeletePT(ec,GG.level)
   end
end; { Gx_Fx }

procedure Gy_Fy(var G,Fy,GG: sparsePolynomial);
{input are GG and initial values of G,Fy.
 The y-antiderivative of GG yields G,Fy-contributions which are added.
 The data GG are overwritten by the zero polynomial }
var lF,e,ec: ^PT; fac: real;
begin while GG.sp<>NIL do
   begin lF:=GG.sp;
      e:=CopyPT(lF,GG.level); ec:=CopyPT(e,GG.level); 
      DeletePolynomialTerm_direct(GG,lF);
      fac:=1/(e^.degree[2]+1);
      TermOp(G, ec,e,0,1, fac,0,0);
      TermOp(Fy,ec,e,0,2,-fac,0,0);
      TermOp(Fy,ec,e,1,1, fac,0,0);          {cf. (5.9)}
      DeletePT(e,GG.level); DeletePT(ec,GG.level)
   end
end; { Gy_Fy }

procedure Lx_Fx(var L,Fx,LL: sparsePolynomial);
{input are LL and initial values of L,Fx.
 The x-antiderivative of LL yields L,Fx-contributions which are added.
 The data LL are overwritten by the zero polynomial }
var lF,e,ec: ^PT; fac: real;
begin while LL.sp<>NIL do
   begin lF:=LL.sp; e:=CopyPT(lF,LL.level); ec:=CopyPT(e,LL.level); 
      DeletePolynomialTerm_direct(LL,lF); fac:=1/(e^.degree[1]+1);
      TermOp(L, ec,e,1,0, fac,0,0);
      TermOp(Fx,ec,e,1,0,-fac,0,0);          {cf. (5.7)}
      DeletePT(e,LL.level); DeletePT(ec,LL.level)
   end
end; { Lx_Fx }

procedure Ly_Fy(var L,Fy,LL: sparsePolynomial);
{input are LL and initial values of L,Fy.
 The y-antiderivative of LL yields L,Fy-contributions which are added.
 The data LL are overwritten by the zero polynomial }
var lF,e,ec: ^PT; fac: real;
begin while LL.sp<>NIL do
   begin lF:=LL.sp; e:=CopyPT(lF,LL.level); ec:=CopyPT(e,LL.level); 
      DeletePolynomialTerm_direct(LL,lF); fac:=1/(e^.degree[2]+1);
      TermOp(L, ec,e,0,1,fac,0,0);
      TermOp(Fy,ec,e,0,1,fac,0,0);          {cf. (5.7)}
      DeletePT(e,LL.level); DeletePT(ec,LL.level)
   end
end; { Ly_Fy }

procedure Ax_Fx(var B,M,P,Fx,A: sparsePolynomial; ix,iz: integer; x,z: real);
{input are A and initial values of B,M,P,Fx.
 ix=0: x is the value of the parameter X, otherwise X is symbolic.
 iz=0: z is the value of the parameter Z, otherwise Z is symbolic.
 The x-antiderivative of A yields B,M,P,Fx-contributions which are added.
 The data A are overwritten by the zero polynomial }
var lF,e,ec: ^PT; k: integer; label 1,9;
begin
1: lF:=LargestPT(A,1); if lF=NIL then goto 9;
   e:=CopyPT(lF,A.level); ec:=CopyPT(e,A.level);
   DeletePolynomialTerm_direct(A,lF);
   k:=e^.degree[1];
   if k>1 then
   begin TermOp(P, ec,e,-1,0,1/(k-1),0,0); { cf. (5.12) }
      if ix=0 then
         TermOp(Fx,ec,e,-2,0,-1,     1,0)  { symbolic variable X } else
      if x<>0 then
         TermOp(Fx,ec,e,-2,0,-x,     0,0); { X evaluated by x }
         TermOp(A, ec,e,-1,1, 2,     0,0);
         TermOp(A, ec,e,-2,2,-1,     0,0);
      if iz=0 then
         TermOp(A, ec,e,-2,0,-1,     0,2)  { symbolic variable Z } else
      if z<>0 then
         TermOp(A, ec,e,-2,0,-sqr(z),0,0)  { Z evaluated by z }
   end else
   if k=0 then AddTerm(B,e) else           {cf. (5.10)}
   if k=1 then 
   begin TermOp(A,ec,e,-1,1,1,0,0);
         TermOp(M,ec,e,-1,0,1,0,0)         {cf. (5.11)}
   end;
   DeletePT(e,A.level); DeletePT(ec,A.level);
   goto 1;
9:
end; { Ax_Fx }

procedure AZx_Fx(var BZ,M,P,Fx,AZ: sparsePolynomial; ix,iz: integer; x,z: real);
{input are A and initial values of BZ,M,P,Fx.
 ix=0: x is the value of the parameter X, otherwise X is symbolic.
 iz=0: z is the value of the parameter Z, otherwise Z is symbolic.
 The x-antiderivative of A yields BZ,M,P,Fx-contributions which are added.
 The data AZ are overwritten by the zero polynomial }
var lF,e,ec: ^PT; k: integer; f: real; label 1,9;
begin
1: lF:=LargestPT(AZ,1); if lF=NIL then goto 9;
   e:=CopyPT(lF,AZ.level); ec:=CopyPT(e,AZ.level);
   DeletePolynomialTerm_direct(AZ,lF);
   k:=e^.degree[1];
   if k>1 then
   begin f:=1/(k-1); if iz=1 then f:=f*sqr(z); if f<>0 then 
         TermOp(P, ec,e,-1,0,f,0,2*(1-iz));
      f:=-1; if ix=1 then f:=f*x; if iz=1 then f:=f*sqr(z); if f<>0 then
         TermOp(Fx,ec,e,-2,0,f, 1-ix,2*(1-iz));
         TermOp(AZ,ec,e,-1,1, 2,0,0);
         TermOp(AZ,ec,e,-2,2,-1,0,0);
      f:=-1; if iz=1 then f:=f*sqr(z); if f<>0 then
         TermOp(AZ,ec,e,-2,0,f,0,2*(1-iz))
   end else
   if k=1 then 
   begin TermOp(AZ,ec,e,-1,1,1,0,0);
         f:=1; if iz=1 then f:=sqr(z); if f<>0 then
         TermOp(M ,ec,e,-1,0,f,0,2*(1-iz))  {cf. (5.11)}
   end else {k=0}
         TermOp(BZ, ec,e,0,0,1,0,0);         {cf. (5.10)}
   DeletePT(e,AZ.level); DeletePT(ec,AZ.level); goto 1;
9:
end; { AZx_Fx }

procedure AZy_Fy(var BZ,M,P,Fy,AZ: sparsePolynomial; ix,iz: integer; x,z: real);
{input are A and initial values of BZ,M,P,Fy.
 ix=0: x is the value of the parameter X, otherwise X is symbolic.
 iz=0: z is the value of the parameter Z, otherwise Z is symbolic.
 The y-antiderivative of A yields BZ,M,P,Fy-contributions which are added.
 The data AZ are overwritten by the zero polynomial }
var lF,e,ec: ^PT; k: integer; f: real; label 1,9;
begin
1: lF:=LargestPT(AZ,2); if lF=NIL then goto 9;
   e:=CopyPT(lF,AZ.level); ec:=CopyPT(e,AZ.level);
   DeletePolynomialTerm_direct(AZ,lF);
   k:=e^.degree[2];
   if k>1 then
   begin f:=1/(k-1); if iz=1 then f:=f*sqr(z); if f<>0 then 
         TermOp(P, ec,e,0,-1,f, 0,2*(1-iz));
      f:=-1; if ix=1 then f:=f*x; if iz=1 then f:=f*sqr(z); if f<>0 then
         TermOp(Fy,ec,e,0,-2,f, 1-ix,2*(1-iz));
         TermOp(AZ,ec,e,1,-1, 2,0,0);
         TermOp(AZ,ec,e,2,-2,-1,0,0);
      f:=-1; if iz=1 then f:=f*sqr(z); if f<>0 then
         TermOp(AZ,ec,e,0,-2,f, 0,2*(1-iz))
   end else
   if k=1 then 
   begin TermOp(AZ,ec,e,1,-1,1,0,0);
         f:=1; if iz=1 then f:=sqr(z); if f<>0 then
         TermOp(M,ec,e,0,-1,f,0,2*(1-iz))  {cf. (5.14)}
   end else {k=0}
         TermOp(BZ,ec,e,0,0,-1,0,0);        {cf. (5.13)} 
   DeletePT(e,AZ.level); DeletePT(ec,AZ.level); goto 1;
9:
end; { AZy_Fy }

procedure DAx_Fx(var M,P,Fx,AZ,DA: sparsePolynomial; ix: integer; x: real);
{input are A and initial values of M,P,Fx.
 ix=0: x is the value of the parameter X, otherwise X is symbolic.
 The x-antiderivative of DA yields BZ,M,P,Fx,AZ-contributions which are added.
 The data A are overwritten by the zero polynomial }
var lF,e,ec: ^PT; k: integer; label 1,9;
begin
1: lF:=LargestPT(DA,1); if lF=NIL then goto 9;
   e:=CopyPT(lF,DA.level); ec:=CopyPT(e,DA.level);
   DeletePolynomialTerm_direct(DA,lF);
   k:=e^.degree[1];
   if k=0 then AddTerm(M,e) {cf. (5.11)}
   else {k>=1}
   begin TermOp(P, ec,e,0,0,1/k,0,0); { cf. (5.12) }
      if ix=0 then
         TermOp(Fx,ec,e,-1,0,-1,1,0)  { symbolic variable X } else
      if x<>0 then
         TermOp(Fx,ec,e,-1,0,-x,0,0); { X evaluated by x }
         TermOp(DA,ec,e,-1,1, 1,0,0);
         TermOp(AZ,ec,e,-1,0,-1,0,0)
   end;
   DeletePT(e,DA.level); DeletePT(ec,DA.level); goto 1;
9:
end; { DAx_Fx }

procedure DAy_Fy(var M,P,Fy,AZ,DA: sparsePolynomial; ix: integer; x: real);
{input are A and initial values of M,P,Fy,AZ
 ix=0: x is the value of the parameter X, otherwise X is symbolic.
 The y-antiderivative of DA yields BZ,M,P,Fy,AZ-contributions which are added.
 The data A are overwritten by the zero polynomial }
var lF,e,ec: ^PT; k: integer; label 1,9;
begin
1: lF:=LargestPT(DA,2); if lF=NIL then goto 9;
   e:=CopyPT(lF,DA.level); ec:=CopyPT(e,DA.level);
   DeletePolynomialTerm_direct(DA,lF);
   k:=e^.degree[2];
   if k=0 then SubtractTerm(M,e)      {cf. (5.14)}
   else {k>=1}
   begin TermOp(P, ec,e,0,0,-1/k,0,0); { cf. (5.15) }
      if ix=0 then
         TermOp(Fy,ec,e,0,-1,1,1,0)  { symbolic variable X } else
      if x<>0 then
         TermOp(Fy,ec,e,0,-1,x,0,0); { X evaluated by x }
         TermOp(DA,ec,e,1,-1,1,0,0);
         TermOp(AZ,ec,e,0,-1,1,0,0)
   end;
   DeletePT(e,DA.level); DeletePT(ec,DA.level); goto 1;
9:
end; { DAy_Fy }

procedure Ay_Fy(var B,M,P,Fy,A: sparsePolynomial; ix,iz: integer; x,z: real);
{input are A and initial values of B,M,P,Fy.
 The y-antiderivative of A yields L,Fy-contributions which are added.
 The data DA are overwritten by the zero polynomial }
var lF,e,ec: ^PT; ll: integer; label 1,9;
begin
1: lF:=LargestPT(A,2); if lF=NIL then goto 9;
   e:=CopyPT(lF,A.level); ec:=CopyPT(e,A.level);
   DeletePolynomialTerm_direct(A,lF);
   ll:=e^.degree[2];
   if ll>1 then
   begin TermOp(P, ec,e,0,-1,1/(ll-1),0,0); { cf. (5.15) }
      if ix=0 then
         TermOp(Fy,ec,e,0,-2,-1,      1,0)  { symbolic variable X } else
      if x<>0 then
         TermOp(Fy,ec,e,0,-2,-x,      0,0); { X evaluated by x }
         TermOp(A, ec,e,1,-1, 2,      0,0);
         TermOp(A, ec,e,2,-2,-1,      0,0);
      if iz=0 then
         TermOp(A, ec,e,0,-2,-1,      0,2)  { symbolic variable Z } else
      if z<>0 then
         TermOp(A, ec,e,0,-2,-sqr(z), 0,0)  { Z evaluated by z }
   end else
   if ll=0 then SubtractTerm(B,e) else      {cf. (5.13)}
   if ll=1 then 
   begin TermOp(A,ec,e,1,-1,1,0,0);
         TermOp(M,ec,e,0,-1,1,0,0)          {cf. (5.14)}
   end;
   DeletePT(e,A.level); DeletePT(ec,A.level);
   goto 1;
9:
end; { Ay_Fy }

procedure Mx_Ax(var M,Ax,MM: sparsePolynomial);
{input are MM and initial values of M,Ax.
 The x-antiderivative of MM yields M,Ax-contributions which are added.
 The data MM are overwritten by the zero polynomial }
begin Gx_Fx(M,Ax,MM)                                   { cf. (5.16) }
end; { Mx_Ax }

procedure My_Ay(var M,Ay,MM: sparsePolynomial);
{input are MM and initial values of M,Ay.
 The y-antiderivative of MM yields M,Ay-contributions which are added.
 The data MM are overwritten by the zero polynomial }
begin Gy_Fy(M,Ay,MM)                                   { cf. (5.17) }
end; { My_Ay }

procedure Bx_Ax(var B,Ax,BB: sparsePolynomial);
{input are BB and initial values of B,Ax.
 The x-antiderivative of BB yields B,Ax-contributions which are added.
 The data BB are overwritten by the zero polynomial }
begin Lx_Fx(B,Ax,BB) end; { Bx_Ax }

procedure By_Ay(var B,Ay,BB: sparsePolynomial);
{input are BB and initial values of B,Ay.
 The y-antiderivative of BB yields B,Ay-contributions which are added.
 The data BB are overwritten by the zero polynomial }
begin Ly_Fy(B,Ay,BB)                                   { cf. (5.18) }
end; { By_Ay }

procedure Pxy(var P,PP: sparsePolynomial; ixy: integer);
{input are PP and initial values of P.
 The x (ixy=1) or y (ixy=2) antiderivative of PP yields P-contributions 
 which are added. The data PP are overwritten by the zero polynomial }
var lF,e,ec: ^PT;
begin while PP.sp<>NIL do
   begin lF:=PP.sp;
      e:=CopyPT(lF,PP.level); ec:=CopyPT(e,PP.level);
      DeletePolynomialTerm_direct(PP,lF);
      TermOp(P,ec,e,2-ixy,ixy-1,1/(1+e^.degree[ixy]),0,0);
      DeletePT(e,PP.level); DeletePT(ec,PP.level);
   end
end; { Pxy }

procedure Px(var P,PP: sparsePolynomial);
begin Pxy(P,PP,1) end;

procedure Py(var P,PP: sparsePolynomial);
begin Pxy(P,PP,2) end;

procedure Psx(var Ps,PP: sparsePolynomial);
{The functions are of the form polynomial(x,y)*sign(x-y).
 Input are PP and initial values of Ps.
 The x antiderivative of PP yields Ps-contributions 
 which are added. The data PP are overwritten by the zero polynomial }
var lF,e,ec: ^PT;
begin while PP.sp<>NIL do
   begin lF:=PP.sp;
      e:=CopyPT(lF,PP.level); ec:=CopyPT(e,PP.level);
      DeletePolynomialTerm_direct(PP,lF);
      TermOp(Ps,ec,e,1,0,1/(1+e^.degree[1]),0,0);
      TermOp(Ps,ec,e,-e^.degree[1],1+e^.degree[1],-1/(1+e^.degree[1]),0,0);
      DeletePT(e,PP.level); DeletePT(ec,PP.level);
   end
end; { Psx }

procedure Psy(var Ps,PP: sparsePolynomial);
{The functions are of the form polynomial(x,y)*sign(x-y).
 Input are PP and initial values of Ps.
 The y antiderivative of PP yields Ps-contributions 
 which are added. The data PP are overwritten by the zero polynomial }
var lF,e,ec: ^PT;
begin while PP.sp<>NIL do
   begin lF:=PP.sp;
      e:=CopyPT(lF,PP.level); ec:=CopyPT(e,PP.level);
      DeletePolynomialTerm_direct(PP,lF);
      TermOp(Ps,ec,e,0,1,1/(1+e^.degree[2]),0,0);
      TermOp(Ps,ec,e,1+e^.degree[2],-e^.degree[2],-1/(1+e^.degree[2]),0,0);
      DeletePT(e,PP.level); DeletePT(ec,PP.level);
   end
end; { Psx }

procedure C0x_Qx(var Qx,C0: sparsePolynomial);
{input are C0 and initial values of Qx.
 The x-antiderivative of C0 yields Qx-contributions which are added.
 The data C0 are overwritten by the zero polynomial }
var lF,e,ec: ^PT;
begin while C0.sp<>NIL do
   begin lF:=C0.sp;
      e:=CopyPT(lF,C0.level); ec:=CopyPT(e,C0.level);
      DeletePolynomialTerm_direct(C0,lF);
      TermOp(Qx,ec,e,1,0, 1,0,0);
      TermOp(Qx,ec,e,0,1,-1,0,0);			{cf. (5.28)}
      DeletePT(e,C0.level); DeletePT(ec,C0.level);
   end
end;

procedure C0y_Qy(var Qy,C0: sparsePolynomial);
begin C0x_Qx(Qy,C0) end;

procedure Cpx_Rpx(var Rpx,Cp: sparsePolynomial);
begin C0x_Qx(Rpx,Cp) end;

procedure Cmx_Rmx(var Rmx,Cm: sparsePolynomial);
begin C0x_Qx(Rmx,Cm) end;

procedure Cpy_Rpy(var Rpy,Cp: sparsePolynomial);
begin C0x_Qx(Rpy,Cp) end;

procedure Cmy_Rmy(var Rmy,Cm: sparsePolynomial);
begin C0x_Qx(Rmy,Cm) end;

procedure Qx_Dx(var Q,Dx,QQ: sparsePolynomial);
{input are QQ and initial values of Q,Dx.
 The x-antiderivative of QQ yields Q,Dx-contributions which are added.
 The data QQ are overwritten by the zero polynomial }
begin Lx_Fx(Q,Dx,QQ) end;

procedure Qy_Dy(var Q,Dy,QQ: sparsePolynomial);
{input are QQ and initial values of Q,Dy.
 The y-antiderivative of QQ yields Q,Dy-contributions which are added.
 The data QQ are overwritten by the zero polynomial }
begin Ly_Fy(Q,Dy,QQ) end;

procedure Dx(var Q,M,P,DD: sparsePolynomial; iy: integer; y: real);
{input are DD and initial values of Q,M,P.
 iy=0: y is the value of the parameter Y, otherwise Y is symbolic.
 The x-antiderivative of DD yields Q,M,P-contributions which are added.
 The data DD are overwritten by the zero polynomial. }
var lF,e,ec: ^PT; k: integer; label 1,9;
begin
1: lF:=LargestPT(DD,1); if lF=NIL then goto 9;
   e:=CopyPT(lF,DD.level); ec:=CopyPT(e,DD.level);
   DeletePolynomialTerm_direct(DD,lF);
   k:=e^.degree[1];
   if k>1 then
   begin if iy=0 then
         TermOp(P, ec,e,-1,0,1/(k-1),0,1) { symbolic variable Y } else if y<>0 then
         TermOp(P, ec,e,-1,0,y/(k-1),0,0);
         TermOp(DD,ec,e,-1,1,2,      0,0);
         TermOp(DD,ec,e,-2,2,-1,     0,0);
      if iy=0 then
         TermOp(DD,ec,e,-2,0,-1,     0,2)  { symbolic variable Y } else
      if y<>0 then
         TermOp(DD,ec,e,-2,0,-sqr(y),0,0)  { Y evaluated by y }
   end {cf. (5.21)} else
   if k=0 then AddTerm(Q,e) else         {cf. (5.19)}
   if k=1 then 
   begin if iy=0 then
         TermOp(M, ec,e,-1,0,1,0,1)  { symbolic variable Y } else if y<>0 then
         TermOp(M, ec,e,-1,0,y,0,0);
         TermOp(DD,ec,e,-1,1,1,0,0)
   end;                                  {cf. (5.20)}
   DeletePT(e,DD.level); DeletePT(ec,DD.level);
   goto 1;
9:
end; { Dx }

procedure Dy(var Q,M,P,DD: sparsePolynomial; iy: integer; y: real);
{input are DD and initial values of Q,M,P.
 iy=0: y is the value of the parameter Y, otherwise Y is symbolic.
 The y-antiderivative of DD yields Q,M,P-contributions which are added.
 The data DD are overwritten by the zero polynomial. }
var lF,e,ec: ^PT; k: integer; label 1,9;
begin
1: lF:=LargestPT(DD,2); if lF=NIL then goto 9;
   e:=CopyPT(lF,DD.level); ec:=CopyPT(e,DD.level);
   DeletePolynomialTerm_direct(DD,lF);
   k:=e^.degree[2];
   if k>1 then
   begin if iy=0 then
         TermOp(P, ec,e,0,-1,1/(k-1),0,1) { symbolic variable Y } else if y<>0 then
         TermOp(P, ec,e,0,-1,y/(k-1),0,0);
         TermOp(DD,ec,e,1,-1,2,      0,0);
         TermOp(DD,ec,e,2,-2,-1,     0,0);
      if iy=0 then
         TermOp(DD,ec,e,0,-2,-1,     0,2)  { symbolic variable Y } else
      if y<>0 then
         TermOp(DD,ec,e,0,-2,-sqr(y),0,0)  { Y evaluated by y }
   end else
   if k=0 then SubtractTerm(Q,e) else
   if k=1 then
   begin if iy=0 then
         TermOp(M, ec,e,0,-1,-1,0,1)  { symbolic variable Y } else if y<>0 then
         TermOp(M, ec,e,0,-1,-y,0,0);
         TermOp(DD,ec,e,1,-1, 1,0,0)
   end;                                  {cf. (5.21)}
   DeletePT(e,DD.level); DeletePT(ec,DD.level);
   goto 1;
9:
end; { Dy }

procedure RxEx_Fx(var R,E,RR,EE,Fx: sparsePolynomial;
                  ix,iy,ipm: integer; x,y: real);
{input are RR,EE and initial values of R,E,Fx.
 ix=0: x is the value of the parameter X, otherwise X is symbolic.
 iy=0: y is the value of the parameter Y, otherwise Y is symbolic.
 ipm=1: Rp,Ep / ipm=-1: Rm,Em.
 The x-antiderivatives of RR and EE yield R,E,Fx-contributions which are added.
 The data RR,EE are overwritten by the zero polynomial. }
var lF,t,tc: ^PT; k,kk: integer; f,ff: real; label 1,9;
begin if ix<>0 then ix:=1; if iy<>0 then iy:=1;
   if ipm<0 then begin k:=ix; ix:=iy; iy:=k; f:=x; x:=y; y:=f end;
1: t :=LargestPT(RR,1); if t =NIL then k :=-1 else k :=t^.degree[1];
   tc:=LargestPT(EE,1); if tc=NIL then kk:=-1 else kk:=tc^.degree[1];
   lF:=t; if kk>k then begin lF:=tc; k:=kk; kk:=2 end else kk:=1;
   if lf=NIL then goto 9;
   t:=CopyPT(lF,RR.level); tc:=CopyPT(t,RR.level);
   if kk=1 then DeletePolynomialTerm_direct(RR,lF)
   else DeletePolynomialTerm_direct(EE,lF);
   k:=t^.degree[1]; f:=1/(k+1);
   if kk=1 then { R-recursion }
   begin TermOp(R, tc,t,1,0,f,   0,0);
         TermOp(R, tc,t,0,1,-f,  0,0);
         TermOp(E, tc,t,0,0,-f,  0,0);
      if k>0 then
      begin
         TermOp(RR,tc,t,-1,1,k*f,0,0);
         TermOp(EE,tc,t,-1,0,k*f,0,0)
      end
   end else { E-recursion }
   begin TermOp(E, tc,t,1,0,f,       0,0);
         TermOp(E, tc,t,0,1,-f,      0,0);
      ff:=f; if (ix<>0) and (ipm=1) then ff:=ff*sqr(x) else 
      if (iy<>0) and (ipm=-1) then ff:=ff*sqr(y);
         TermOp(R, tc,t,0,0,ff,      (1+ipm)*(1-ix),(1-ipm)*(1-iy));
      ff:=-f; if ix<>0 then ff:=ff*x; if iy<>0 then ff:=ff*y;
         TermOp(Fx,tc,t,0,0,ff,      1-ix,1-iy);
      if k>0 then
      begin ff:=-k*f; if (ix<>0) and (ipm=1) then ff:=ff*sqr(x) else 
         if (iy<>0) and (ipm=-1) then ff:=ff*sqr(y);
         TermOp(RR,tc,t,-1,0,-k*f,   (1+ipm)*(1-ix),(1-ipm)*(1-iy));
         TermOp(EE,tc,t,-1,1,k*f,    0,0)
      end
   end;
   DeletePT(t,RR.level); DeletePT(tc,RR.level);
   goto 1;
9:
end;

procedure RyEy_Fy(var R,E,RR,EE,Fy: sparsePolynomial;
                  ix,iy,ipm: integer; x,y: real);
{input are RR,EE and initial values of R,E,Fy.
 ix=0: x is the value of the parameter X, otherwise X is symbolic.
 iy=0: y is the value of the parameter Y, otherwise Y is symbolic.
 ipm=1: Rp,Ep / ipm=-1: Rm,Em.
 The y-antiderivatives of RR and EE yield R,E,Fy-contributions which are added.
 The data RR,EE are overwritten by the zero polynomial. }
var lF,t,tc: ^PT; k,kk: integer; f,ff: real; label 1,9;
begin if ix<>0 then ix:=1; if iy<>0 then iy:=1;
   if ipm<0 then begin k:=ix; ix:=iy; iy:=k; f:=x; x:=y; y:=f end;
1: t :=LargestPT(RR,2); if t =NIL then k :=-1 else k :=t^.degree[2];
   tc:=LargestPT(EE,2); if tc=NIL then kk:=-1 else kk:=tc^.degree[2];
   lF:=t; if kk>k then begin lF:=tc; k:=kk; kk:=2 end else kk:=1;
   if lf=NIL then goto 9;
   t:=CopyPT(lF,RR.level); tc:=CopyPT(t,RR.level);
   if kk=1 then DeletePolynomialTerm_direct(RR,lF)
           else DeletePolynomialTerm_direct(EE,lF);
   k:=t^.degree[2]; f:=1/(k+1);
   if kk=1 then { R-recursion }
   begin TermOp(R, tc,t,1,0,-f,   0,0);
         TermOp(R, tc,t,0,1,f,    0,0);
         TermOp(E, tc,t,0,0,f,    0,0);
      if k>0 then
      begin
         TermOp(RR,tc,t,1,-1, k*f,0,0);
         TermOp(EE,tc,t,0,-1,-k*f,0,0)
      end
   end else { E-recursion }
   begin TermOp(E, tc,t,1,0,-f,   0,0);
         TermOp(E, tc,t,0,1, f,   0,0);
      ff:=-f; if (ix<>0) and (ipm=1) then ff:=ff*sqr(x) else 
      if (iy<>0) and (ipm=-1) then ff:=ff*sqr(y);
         TermOp(R, tc,t,0,0,ff,   (1+ipm)*(1-ix),(1-ipm)*(1-iy));
      ff:=-f; if ix<>0 then ff:=ff*x; if iy<>0 then ff:=ff*y;
         TermOp(Fy,tc,t,0,0,ff,   1-ix,1-iy);
      if k>0 then
      begin ff:=k*f; if (ix<>0) and (ipm=1) then ff:=ff*sqr(x) else
         if (iy<>0) and (ipm=-1) then ff:=ff*sqr(y);
         TermOp(RR,tc,t,0,-1,ff,  (1+ipm)*(1-ix),(1-ipm)*(1-iy));
         TermOp(EE,tc,t,1,-1,k*f, 0,0)
      end
   end;
   DeletePT(t,RR.level); DeletePT(tc,RR.level);
   goto 1;
9:
end;

function TeX_Brackets(ein: string): string; var s: string;
{WERT ist "ein" in Brackets}
begin writev(s,auf,ein,zu); TeX_Brackets:=s
end; { TeX_Brackets }

procedure TeX_VarName(var t: TeXExpr_T; s: string);
begin with t do begin
   p0:=s; writev(m0,'-',s); p1:=p0; m1:=m0;
   p2:=p0; m2:=TeX_Brackets(m0);
   writev(pv,'+',p0); mv:=m0;
   Vorz:=3
end end; { TeX_VarName }

procedure TeX_pi(var t: TeXExpr_T);
begin TeX_VarName(t,'\pi ') end; { TeX_pi }

procedure TeX_NegativeExpr(var tNegativ: TeXExpr_T; t: TeXExpr_T);
begin with tNegativ do begin
   p0:=t.m0; p1:=t.m1; p2:=t.m2; pv:=t.mv;
   m0:=t.p0; m1:=t.p1; m2:=t.p2; mv:=t.pv;
   if abs(t.Vorz)<=2 then vorz:=-t.vorz else vorz:=3
end end; { TeX_NegativeExpr }

procedure TeX_Integer(var t: TeXExpr_T; i: integer); var s: string;
begin writev(s,abs(i):1); TeX_VarName(t,s);
   if i<0 then TeX_NegativeExpr(t,t) else
      if i=0 then with t do begin m0:=p0; m1:=p0; m2:=p0 end;
   if abs(i)<=1 then t.Vorz:=i else if i>0 then t.Vorz:=2 else if i<0 then t.Vorz:=-2;
end; { TeX_Integer }

procedure TeX_Real(var t: TeXExpr_T; x: real); var s: string;
begin writev(s,abs(x):25:22); TeX_VarName(t,s);
   if x<0 then TeX_NegativeExpr(t,t) else
      if x=0 then with t do begin m0:=p0; m1:=p0; m2:=p0 end;
   if x=1 then t.Vorz:=1 else if x=-1 then t.Vorz:=-1 else
      if x=0 then t.Vorz:=0 else if x>0 then t.Vorz:=2 else if x<0 then t.Vorz:=-2
end; { TeX_Real }

procedure TeX_Sum(var s: TeXExpr_T; t1,t2: TeXExpr_T; Simplify: integer);
{ s := t1 + t2; falls Simplify=1: t1=0 bzw t2=0 entfaellt }
label 9;
begin with s do begin
   if Simplify=1 then
   begin if t1.vorz=0 then begin s:=t2; goto 9 end;
         if t2.vorz=0 then begin s:=t1; goto 9 end;
         if t1.p0=t2.m0 then begin TeX_Integer(s,0); goto 9 end
   end;
   if (abs(t1.vorz)<=1) and (abs(t2.vorz)<=1) then Vorz:=t1.vorz+t2.vorz
     else if (t1.Vorz=3) or (t2.Vorz=3) then Vorz:=3
       else if t1.vorz=0 then Vorz:=t2.vorz
         else if t2.vorz=0 then Vorz:=t1.vorz
           else if t1.vorz*t2.vorz>0 then
             begin if t1.vorz>0 then Vorz:=2 else if t1.vorz<0 then Vorz:=-2 else Vorz:=0
             end else Vorz:=3;
   writev(p0,t1.p0,t2.pv); writev(m0,t1.m0,t2.mv);
   p1:=TeX_Brackets(p0); m1:=TeX_Brackets(m0);
   p2:=p1; m2:=m1;
   writev(pv,t1.pv,t2.pv); writev(mv,t1.mv,t2.mv)
   end;
9:end; { TeX_Sum }

procedure TeX_Product(var p: TeXExpr_T; t1,t2: TeXExpr_T; Simplify: integer);
{ s := t1 * t2; if Simplify=1: t1=+/-1 or t2=+/-1 reduced }
var s1,s2: string; label 9;
begin with p do begin
   if Simplify=1 then
   begin if t1.vorz*t2.vorz=0 then begin TeX_Integer(p,0); goto 9 end;
      if t1.vorz=1  then begin p:=t2; goto 9 end;
      if t1.vorz=-1 then begin TeX_NegativeExpr(p,t2); goto 9 end;
      if t2.vorz=1  then begin p:=t1; goto 9 end;
      if t2.vorz=-1 then begin TeX_NegativeExpr(p,t1); goto 9 end
   end;
   Vorz:=t1.vorz*t2.vorz;
   if (t1.Vorz=3) or (t2.Vorz=3) then Vorz:=3
   else if (abs(t1.vorz)<=1) and (abs(t2.vorz)<=1) then 
   else if t1.vorz*t2.vorz=0 then
   else if Vorz>0 then Vorz:=2 else Vorz:=-2;
   writev(s1,t1.p1,'*',t2.p2); writev(s2,t1.m1,'*',t2.m2);
   if length(s1)<=length(s2) then p0:=s1 else p0:=s2;
   writev(s1,t1.m1,'*',t2.p2); writev(s2,t1.p1,'*',t2.m2);
   if length(s1)<=length(s2) then m0:=s1 else m0:=s2;
   p1:=p0; m1:=m0;
   writev(s1,t1.p2,'*',t2.p2); writev(s2,t1.m2,'*',t2.m2);
   if length(s1)<=length(s2) then p2:=s1 else p2:=s2;
   writev(s1,t1.m2,'*',t2.p2); writev(s2,t1.p2,'*',t2.m2);
   if length(s1)<=length(s2) then m2:=s1 else m2:=s2;
   writev(s1,'+',p2); writev(s2,'-',m2);
   if length(s1)<=length(s2) then pv:=s1 else pv:=s2;
   writev(s1,'-',p2); writev(s2,'+',m2);
   if length(s1)<=length(s2) then mv:=s1 else mv:=s2;
end;
9:end; { TeX_Product }

procedure TeX_Frac(var f: TeXExpr_T; t1,t2: TeXExpr_T; Simplify: integer);
{ f := t1 / t2; falls Simplify=1: t1=+/-1 bzw t2=+/-1 entfaellt }
var s1,s2: string; label 9;
begin with f do begin
   if t2.vorz=0 then writeln('#### WARNING: DIVISION BY ZERO IN TeX_Frac') else
   if Simplify=1 then
   begin if t1.vorz=0   then begin TeX_Integer(f,0); goto 9 end;
         if t2.vorz=1   then begin f:=t1; goto 9 end;
         if t2.vorz=-1  then begin TeX_NegativeExpr(f,t1); goto 9 end;
         if t1.p0=t2.p0 then begin TeX_Integer(f,1); goto 9 end;
         if t1.p0=t2.m0 then begin TeX_Integer(f,-1); goto 9 end
   end;
   Vorz:=t1.vorz*t2.vorz;
   if (t1.Vorz=3) or (t2.Vorz=3) then Vorz:=3
   else if (abs(t1.vorz)<=1) and (abs(t2.vorz)<=1) then 
   else if t1.vorz=0 then
   else if t2.vorz=0 then Vorz:=3
   else if Vorz>0 then Vorz:=2 else Vorz:=-2;
   writev(s1,fa,t1.p1,fm,t2.p2,fe);
   writev(s2,fa,t1.m1,fm,t2.m2,fe);
   if length(s1)<=length(s2) then p0:=s1 else p0:=s2;
   writev(s1,fa,t1.m1,fm,t2.p2,fe); writev(s2,fa,t1.p1,fm,t2.m2,fe);
   if length(s1)<=length(s2) then m0:=s1 else m0:=s2;
   p1:=p0; m1:=m0; p2:=p0; m2:=m0;
   writev(s1,'+',p2); writev(s2,'-',m2);
   if length(s1)<=length(s2) then pv:=s1 else pv:=s2;
   writev(s1,'-',p2); writev(s2,'+',m2);
   if length(s1)<=length(s2) then mv:=s1 else mv:=s2;
end;
9:end; { TeX_Frac }

procedure TeX_Root(var r: TeXExpr_T; t: TeXExpr_T; Simplify: integer);
{ r := Root von t1; falls Simplify=1: speziell fuer t1=1 und t1=0 }
label 9;
begin with r do begin
   if t.vorz<0 then writeln('#### WARNING: ROOT WITH NEGATIVE ARGUMENT IN TeX_Root');
   if Simplify=1 then
   begin if t.vorz=0 then begin TeX_Integer(r,0); goto 9 end;
         if t.vorz=1 then begin TeX_Integer(r,1); goto 9 end
   end;
   writev(p0,ra,t.p0,fe); TeX_VarName(r,p0);
   Vorz:=t.vorz; if t.Vorz<0 then Vorz:=3
   end;
9:end; { TeX_Root }

procedure TeX_RationalNumber(var t: TeXExpr_T; var r: RationalNumber_T);
var t1,t2: TeXExpr_T;
begin TeX_Integer(t1,r.Z); TeX_Integer(t2,r.N); TeX_Frac(t,t1,t2,1)
end; { TeX_RationalNumber }

procedure TeX_XAsRationalNumber(var t: TeXExpr_T; x: real); var r: RationalNumber_T;
var d: real;
begin d:=RationalFromReal0(r,x);
   if d<=TolRational then TeX_RationalNumber(t,r) else TeX_Real(t,x)
end; { TeX_XAsRationalNumber }

procedure TeX_Function(var t: TeXExpr_T; argument: TeXExpr_T; FctName: string);
var s: string;
begin writev(s,FctName,TeX_Brackets(argument.p0)); TeX_VarName(t,s)
end; { TeX_Function }

procedure TeX_OddFunction(var t: TeXExpr_T; argument: TeXExpr_T;
                          FctName: string; Simplify: integer);
var t1,t2: TeXExpr_T; var s: string; label 9;
begin if Simplify=1 then
   if argument.Vorz=0 then begin TeX_Integer(t,0); goto 9 end;
   with t do begin
      TeX_Function(t1,argument,FctName);
      TeX_NegativeExpr(t2,argument); TeX_Function(t2,t2,FctName);
      TeX_NegativeExpr(t2,t2);
      if length(t1.p0)<=length(t2.p0) then p0:=t1.p0 else p0:=t2.p0;
      if length(t1.p1)<=length(t2.p1) then p1:=t1.p1 else p1:=t2.p1;
      if length(t1.p2)<=length(t2.p2) then p2:=t1.p2 else p2:=t2.p2;
      if length(t1.pv)<=length(t2.pv) then pv:=t1.pv else pv:=t2.pv;
      if length(t1.m0)<=length(t2.m0) then m0:=t1.m0 else m0:=t2.m0;
      if length(t1.m1)<=length(t2.m1) then m1:=t1.m1 else m1:=t2.m1;
      if length(t1.m2)<=length(t2.m2) then m2:=t1.m2 else m2:=t2.m2;
      if length(t1.mv)<=length(t2.mv) then mv:=t1.mv else mv:=t2.mv;
      if argument.Vorz=0 then Vorz:=0 else Vorz:=3
   end;
9:end; { TeX_OddFunction }

procedure TeX_arctan(var t: TeXExpr_T; argument: TeXExpr_T; Simplify: integer);
begin TeX_OddFunction(t,argument,'\arctan ',Simplify);
   t.Vorz:=argument.Vorz; if abs(t.Vorz)=1 then t.Vorz:=2*t.Vorz
end; { TeX_arctan }

procedure TeX_ln(var t: TeXExpr_T; argument: TeXExpr_T; Simplify: integer); label 9;
begin if argument.Vorz<=0 then writeln('#### WARNING: LOGARITHMS WITH ARGUMENT <=0 IN TeX_ln');
   if Simplify=1 then if argument.Vorz=1 then begin TeX_Integer(t,0); goto 9 end;
   TeX_Function(t,argument,'\ln '); if abs(argument.Vorz)=1 then t.Vorz:=0 else t.Vorz:=3;
9:end; { TeX_ln }

procedure TeX_G(var s: string; var p: sparsePolynomial; xy,xx,yy: real);
var ff,w: real; t,t2: TeXExpr_T; label 9;
begin s:=''; w:=sqr(xy)+sqr(xx)+sqr(yy); if w=0 then goto 9;
   ff:=constant_from_Polynomial(p); if ff=0 then goto 9;
   TeX_XAsRationalNumber(t,ff); TeX_XAsRationalNumber(t2,w);
   TeX_Root(t2,t2,1); TeX_Product(t,t,t2,1);
   s:=t.pv;
9:
end; { TeX_G }

procedure TeX_F(var s: string; var p: sparsePolynomial; xy,xx,yy: real);
var ff,w: real; t,t2: TeXExpr_T; label 9;
begin s:=''; w:=sqr(xy)+sqr(xx)+sqr(yy); if w=0 then goto 9;
   ff:=constant_from_Polynomial(p); if ff=0 then goto 9;
   TeX_XAsRationalNumber(t,ff); TeX_XAsRationalNumber(t2,w);
   TeX_Root(t2,t2,1); TeX_Frac(t,t,t2,1);
   s:=t.pv;
9:
end; { TeX_F }

procedure TeX_L(var s: string; var p: sparsePolynomial; xy,xx,yy: real);
var t,tz,t2: TeXExpr_T; ff: real; label 9;
begin s:=''; ff:=constant_from_Polynomial(p); if ff=0 then goto 9;
   xx:=sqr(xy)+sqr(xx)+sqr(yy); if xy+sqrt(xx)=0 then goto 9;
   TeX_XAsRationalNumber(t,ff);
   TeX_XAsRationalNumber(tz,xy); TeX_XAsRationalNumber(t2,xx); TeX_Root(t2,t2,1);
   TeX_Sum(t2,tz,t2,1); TeX_ln(t2,t2,1); TeX_Product(t,t,t2,1);
   if t.Vorz<>0 then s:=t.pv;
9:
end; { TeX_L }

procedure TeX_M(var s: string; var p: sparsePolynomial; xy,xx,yy: real);
var t,tz,t2: TeXExpr_T; ff: real; label 9;
begin s:=''; ff:=constant_from_Polynomial(p); if ff=0 then goto 9;
   yy:=sqr(xy)+sqr(xx)+sqr(yy); if xx+sqrt(yy)=0 then goto 9;
   TeX_XAsRationalNumber(t,ff);
   TeX_XAsRationalNumber(tz,xx); TeX_XAsRationalNumber(t2,yy); TeX_Root(t2,t2,1);
   TeX_Sum(t2,tz,t2,1); TeX_ln(t2,t2,1); TeX_Product(t,t,t2,1);
   if t.Vorz<>0 then s:=t.pv;
9:
end; { TeX_M }

procedure TeX_Rp(var s: string; var p: sparsePolynomial; xy,xx,yy: real);
var t,tw,t2,tf: TeXExpr_T; ff: real;  label 9;
begin s:=''; if xy*xx*yy=0 then begin writeln('##### 0/0 in TeX_Rp'); goto 9 end;
   if yy*xy=0 then goto 9;
   ff:=constant_from_Polynomial(p); if ff=0 then goto 9;
   if yy<0 then begin ff:=-ff; yy:=-yy end;
   if xx<0 then begin ff:=-ff; xx:=-xx end;
   if xy<0 then begin ff:=-ff; xy:=-xy end;
   if xx=0 then ff:=ff/2; TeX_XAsRationalNumber(tf,ff);
   if xx=0 then begin TeX_pi(t2); TeX_Product(t,tf,t2,1); goto 9 end;
   TeX_XAsRationalNumber(t2,yy*xy/xx);
   TeX_XAsRationalNumber(tw,sqr(xy)+sqr(xx)+sqr(yy)); TeX_Root(tw,tw,1);
   TeX_Integer(t,1); TeX_Frac(tw,t,tw,1); TeX_Product(tw,t2,tw,1);
   TeX_arctan(tw,tw,1); TeX_XAsRationalNumber(tf,ff); TeX_Product(t,tf,tw,1);
   if t.Vorz<>0 then s:=t.pv;
9:
end; { TeX_Rp }

procedure TeX_Rm(var s: string; var p: sparsePolynomial; xy,xx,yy: real);
var t,tw,t2,tf: TeXExpr_T; ff: real;  label 9;
begin s:=''; if xy*xx*yy=0 then begin writeln('##### 0/0 in TeX_Rm'); goto 9 end;
   TeX_Rp(s,p,xy,yy,xx);
9:
end; { TeX_Rp }

procedure TeX_Q(var s: string; var p: sparsePolynomial; xy,yy: real);
var t,t2: TeXExpr_T; ff: real; label 9;
begin s:=''; if xy*yy=0 then begin writeln('##### 0/0 in TeX_Q'); goto 9 end;
   if xy=0 then goto 9;
   ff:=constant_from_Polynomial(p); if ff=0 then goto 9;
   if ((yy=0) and (xy<0)) or (yy*xy<0) then begin ff:=-ff; xy:=-xy end;
   if yy=0 then ff:=ff/2; TeX_XAsRationalNumber(t,ff);
   if yy=0 then begin TeX_pi(t2); TeX_Product(t,t,t2,1); goto 9 end;
   TeX_XAsRationalNumber(t2,xy/yy);
   TeX_arctan(t2,t2,1); TeX_Product(t,t,t2,1);
   if t.Vorz<>0 then s:=t.pv;
9:
end; { TeX_Q }

procedure TeX_PiE(var s: string; var p: sparsePolynomial);
var t,t2: TeXExpr_T; ff: real; label 9;
begin s:=''; ff:=constant_from_Polynomial(p); if ff=0 then goto 9;
   TeX_XAsRationalNumber(t,ff); TeX_pi(t2); TeX_Product(t,t,t2,1);
   s:=t.pv;
9:
end; { TeX_PiE }

procedure TeX_Ps(var s: string; var p: sparsePolynomial; xy: real);
var t,t2: TeXExpr_T; ff: real; label 9;
begin s:=''; ff:=constant_from_Polynomial(p); if ff=0 then goto 9;
   if xy<0 then ff:=-ff else if xy=0 then goto 9;
   TeX_XAsRationalNumber(t,ff/2); TeX_pi(t2); TeX_Product(t,t,t2,1);
   s:=t.pv;
9:
end; { TeX_Ps }

procedure TeX_P(var s: string; var p: sparsePolynomial);
var t,t2: TeXExpr_T; ff: real; label 9;
begin s:=''; ff:=constant_from_Polynomial(p); if ff=0 then goto 9;
   TeX_XAsRationalNumber(t,ff); s:=t.pv;
9:
end; { TeX_P }

procedure StartAT(var A: AT; fct,sx,sy,sxx,syy: string; x,y,xx,yy: real);
begin if sx='*' then sx:='x'; if sy='*' then sy:='y';
   if sxx='*' then sxx:='X'; if syy='*' then syy:='Y';
   A.fct:=fct; A.x:=x; A.y:=y; A.sx:=sx; A.sy:=sy; A.xx:=xx; A.yy:=yy; 
   A.sxx:=sxx; A.syy:=syy; A.post:=NIL; StartP(A.sp,1)
end; { StartAT }

function NewAT(fct,sx,sy,sxx,syy: string; x,y,xx,yy: real): ^AT;
var A: ^AT;
begin new(A); StartAT(A^,fct,sx,sy,sxx,syy,x,y,xx,yy); NewAT:=A end;

function CompareAT(fct1,sx1,sy1,sxx1,syy1: string; x1,y1,xx1,yy1: real;
                   fct2,sx2,sy2,sxx2,syy2: string; x2,y2,xx2,yy2: real): integer;
{ Result=0: components fct,,...,yy have equal meaning,
        <0: Data1 ordered before Data2,
        >0: Data1 ordered after Data2. }
var c,cxy: integer; d: real; label 0,1,2,9;
begin c:=0; cxy:=0;

   { difference in function name: }
   if fct1<fct2 then c:=-1 else if fct1>fct2 then c:=1; if c<>0 then goto 9;

   { difference in x-,y-evaluation status or symbolic names: }
   if (sx1<>'$') or (sx2<>'$') or (sy1<>'$') or (sy2<>'$') then
   begin if sx1<sx2 then c:=-1 else if sx1>sx2 then c:=1; if c<>0 then goto 9;
         if sy1<sy2 then c:=-1 else if sy1>sy2 then c:=1; if c<>0 then goto 9
   end; { now sx1=sx2, sy1=sy2 holds }

   { difference in xx-,yy-evaluation status or symbolic names: }
   if (sxx1<>'$') or (sxx2<>'$') or (syy1<>'$') or (syy2<>'$') then
   begin if sxx1<sxx2 then c:=-1 else if sxx1>sxx2 then c:=1; if c<>0 then goto 9;
         if syy1<syy2 then c:=-1 else if syy1>syy2 then c:=1; if c<>0 then goto 9
   end; { now sxx1=sxx2, syy1=syy2 holds }

   { case when yy is evaluated -> value yy decides }
   if syy1='$' then if (fct1<>'P') and (fct1<>'Ps') and (fct1<>'Pi') then
   begin d:=yy1-yy2; if d<-Tolerance then c:=-1 else if d>Tolerance then c:=1 end;

   { case when xx is evaluated -> value xx decides }
   if sxx1='$' then
   if (fct1<>'P') and (fct1<>'Ps') and (fct1<>'Pi')
       and (fct1<>'D') and (fct1<>'Q') and (fct1<>'C0') then
   begin d:=xx1-xx2; if d<-Tolerance then c:=-1 else if d>Tolerance then c:=1 end;

   { case when xx and yy are evaluated -> xx*xx+yy*yy may decide }
   if g34=1 then if c<>0 then if (sxx1='$') and (syy1='$') then
      if (fct1='F') or (fct1='G') or (fct1='L') then
      begin d:=sqr(xx1)+sqr(yy1)-(sqr(xx2)+sqr(yy2));
	 if abs(d)<=Tolerance then c:=0 else
            if d<-Tolerance then c:=-1 else c:=1
         end;

   { case when y is evaluated -> value y decides }
   if sy1='$' then
   begin d:=y1-y2; if d<-Tolerance then cxy:=-1 else if d>Tolerance then cxy:=1 end;

   { case when x is evaluated -> value x decides }
   if sx1='$' then
   begin d:=x1-x2; if d<-Tolerance then cxy:=-1 else if d>Tolerance then cxy:=1 end;

   if (c=0) and (cxy=0) then goto 9;
   if (sx1<>'$') or (sy1<>'$') then
   begin if cxy<>0 then c:=cxy; goto 9 end; { now sx1=sy1='$' }

   { if all parameters are evaluated, their values do not matter for P,Pi }
   if (fct1='P') or (fct1='Pi') then
      if (sxx1='$') and (syy1='$') then begin c:=0; goto 9 end;

   { case when x and y are evaluated -> difference x-y decides }
   if cxy<>0 then 
   begin if abs(x1-y1-x2+y2)<=Tolerance then cxy:=0 else
      if x1-y1-x2+y2<Tolerance then cxy:=-1 else cxy:=1
   end;
   if cxy<>0 then c:=cxy;

   if c=0 then goto 9;
   if ((fct1<>'F') and (fct1<>'G')) or (sxx1<>'$') then goto 9;

   { case of F,G -> value sqr(x-y)+sqr(xx) decides }
   if g13=1 then
   begin if syy1<>'$' then
   begin d:=sqr(x1-y1)+sqr(xx1)-(sqr(x2-y2)+sqr(xx2));
         if abs(d)<=Tolerance then c:=0 else if d<-Tolerance then c:=-1 else c:=1;
         if c=0 then goto 9
      end else
      begin { case of F,G -> value sqr(x-y)+sqr(xx)+sqr(yy) decides }
	 d:=sqr(x1-y1)+sqr(xx1)+sqr(yy1)-(sqr(x2-y2)+sqr(xx2)+sqr(yy2));
	 if abs(d)<=Tolerance then c:=0 else if d<-Tolerance then c:=-1 else c:=1
      end
   end;

9: CompareAT:=c
{   ;writeln('------------ >>>>>>>: ',c:1); if c=0 then
   if (fct1<>fct2) or (sx1<>sx2) or (sy1<>sy2) or (sxx1<>sxx2) or (syy1<>syy2)
      or (x1<>x1) or (y1<>y2) or (xx1<>xx2) or (yy1<>yy2) then readln}
end;

function CompareFct(A1,A2: ^AT): integer;
{ Result=0: components x,...,syy have equal meaning,
        <0: A1 before A2 in lexicographical ordering,
        >0: A1 after A2. }
begin CompareFct:=
   CompareAT(A1^.fct,A1^.sx,A1^.sy,A1^.sxx,A1^.syy,A1^.x,A1^.y,A1^.xx,A1^.yy,
             A2^.fct,A2^.sx,A2^.sy,A2^.sxx,A2^.syy,A2^.x,A2^.y,A2^.xx,A2^.yy)
end; { CompareFct }

function FindAT(list: ^AT; fct,sx,sy,sxx,syy: string; x,y,xx,yy: real): ^AT;
{ find the first element in the "list" with the characteristics fct,...,yy.
  Strings '*' allows any value.
  Result=NIL in the negative case }
var l: ^AT; Cfct,Csx,Csy,Csxx,Csyy: string; Cx,Cy,Cxx,Cyy: real; label 9;
begin l:=list; while l<>NIL do
   begin Cfct:=fct; if fct='*' then Cfct:=l^.fct;
      Csx:=sx;   if sx ='*' then Csx :=l^.sx;
      Csy:=sy;   if sy ='*' then Csy :=l^.sy;
      Csxx:=sxx; if sxx='*' then Csxx:=l^.sxx;
      Csyy:=syy; if syy='*' then Csyy:=l^.syy;
      Cx:=x; if (sx<>'$') or (Csx<>'$') then Cx:=l^.x;
      Cy:=y; if (sy<>'$') or (Csy<>'$') then Cy:=l^.y;
      Cxx:=xx; if (sxx<>'$') or (Csxx<>'$') then Cxx:=l^.xx;
      Cyy:=yy; if (syy<>'$') or (Csyy<>'$') then Cyy:=l^.yy;
      if CompareAT(l^.fct,l^.sx,l^.sy,l^.sxx,l^.syy,l^.x,l^.y,l^.xx,l^.yy,
                   Cfct,Csx,Csy,Csxx,Csyy,Cx,Cy,Cxx,Cyy)=0 then goto 9;
      l:=l^.post
   end; l:=NIL;
9: FindAT:=l
end;

procedure writeFct(var fct,sx,sy,sxx,syy: string; x,y,xx,yy: real);
var s: string; label 6,8,9;
begin if fct='P' then begin s:='1'; write_P(s); goto 9 end;
   if fct='Pi' then begin s:='pi'; write_P(s); goto 9 end;
   if fct='Ps' then begin s:='pi/2 * '; write_P(s); s:='sign' end else s:=fct;
   write_P(s); s:='('; write_P(s);
   if (sx='$') and (sy='$') then WriteReal(x-y) else
      if sx='$' then
      begin WriteReal(x); writev(s,' - ',sy); write_P(s) end else
         if sy='$' then
         begin write_P(sx); WriteReal(-y) end else
            begin writev(s,sx,'-',sy); write_P(s) end;
   if fct='Ps' then goto 8;
   s:=';'; write_P(s);
   if (fct='D') or (fct='Q') or (fct='C0') then goto 6;  { only 2nd parameter }
   if sxx<>'$' then write_P(sxx) else WriteReal(xx);
   s:=','; write_P(s);
6: if syy<>'$' then write_P(syy) else WriteReal(yy);
8: s:=')'; write_P(s);
9: writeln_P
end;

procedure WriteFctList(title: string; List: ^AT);
var l: ^AT; s: string;
begin l:=List; write_P(title); writeln_P;
   if l=NIL then begin s:='ZERO'; write_P(s); writeln_P end else
   while l<>NIL do with l^ do
   begin if sp.sp<>NIL then
      begin WritePolynomial(' + ',sp,sx,sy,sxx,syy);
            writeFct(Fct,sx,sy,sxx,syy,x,y,xx,yy)
      end;
      l:=post
   end
end;

procedure InsertAT(var list: ^AT; e: ^AT);
{ e is to be inserted into the list }
var l,lp: ^AT; label 9;
begin if e=NIL then goto 9; l:=list;
   if l=NIL then begin list:=e; e^.post:=NIL; goto 9 end;
   if CompareFct(l,e)>=0 then
   begin list:=e; e^.post:=l; goto 9 end;  { e inserted as first element }
   lp:=NIL; {previous element}
   repeat if CompareFct(l,e)>=0 then
      begin lp^.post:=e; e^.post:=l; goto 9 end; { e inserted between lp and l }
      lp:=l; l:=l^.post
   until l=NIL;
   lp^.post:=e; e^.post:=NIL; { e inserted after lp as last element }
9:
end;

function FindOrCreateAT(var list: ^AT; fct,sx,sy,sxx,syy: string; x,y,xx,yy: real): ^AT;
{ as FindAT. In the negative case create and insert an element }
var l: ^AT; sx_,sy_,sxx_,syy_: string;
begin sx_:=sx; sy_:=sy; sxx_:=sxx; syy_:=syy;
   l:=FindAT(list,fct,sx_,sy_,sxx_,syy_,x,y,xx,yy);
   if sx_ ='*' then sx_ :='$'; if sy_ ='*' then sy_ :='$';
   if sxx_='*' then sxx_:='$'; if syy_='*' then syy_:='$';
   if l=NIL then
   begin l:=NewAT(fct,sx_,sy_,sxx_,syy_,x,y,xx,yy); InsertAT(list,l) end;
   FindOrCreateAT:=l
end;

procedure ExcludeAT(var list: ^AT; a: ^AT);
{ remove "a" from the list without deleting "a" }
var l,lp: ^AT; label 9;
begin if a=NIL then goto 9; if list=NIL then goto 9;
   if list=a then begin list:=a^.post; goto 9 end;
   l:=list; while l<>NIL do
   begin if l^.post=a then begin l^.post:=a^.post; goto 9 end;
         l:=l^.post
   end;
9: 
end;

function DeleteAT(a: ^AT): ^AT;
begin if a<>NIL then begin
   DeletePolynomial(a^.sp);
   dispose(a)
   end;
   DeleteAT:=NIL
end;

function DeleteAFromList(list,a: ^AT): ^AT;
begin ExcludeAT(list,a); DeleteAT(a); DeleteAFromList:=list end;

function DeleteAList(list: ^AT): ^AT;
begin while list<>NIL do list:=DeleteAFromList(list,list); DeleteAList:=NIL end;

procedure DefineFct(var list: ^AT; sxx,syy: string);
var a: ^AT; fct,sx,sy: string;
begin
   writeln('##### input of a function ...');
   writeln('      ** The name should one of the following:');
   writeln('         AZ,BZ,C0,Cm,Cp,D,DA,Em,Ep,F,G,L,M,P,Pi,Ps,Q,Rm,Rp');
   write('      -> name = '); readln(fct);
   sx:='x'; sy:='y'; if sxx='$' then sxx:='X'; if syy='$' then syy:='Y';
   a:=NewAT(fct,sx,sy,sxx,syy,0,0,0,0); Polynomial1_input(a^.sp);
   InsertAT(list,a)
end;


function CopyAT(A: ^AT): ^AT;
{ creates a new copy of A }
var Acopy: ^AT;
begin with A^ do Acopy:=NewAT(A^.fct,sx,sy,sxx,syy,x,y,xx,yy);
      CopyPolynomialsp(Acopy^.sp,A^.sp);
      CopyAT:=Acopy
end;

procedure CopyAList(var listcopy: ^AT; listoriginal: ^AT);
var a: ^AT;
begin listcopy:=NIL; a:=listoriginal;
   while a<>NIL do begin InsertAT(listcopy,CopyAT(a)); a:=a^.post end
end;

procedure multiplyATbyPolynomial(list: ^AT; var p: sparsePolynomial);
{ list := list * p }
var l: ^AT; pp: ^sparsePolynomial;
begin l:=list; while l<>NIL do
   begin pp:=CopyPolynomial(l^.sp);
      DeletePolynomial(l^.sp); 
      multiply(l^.sp,pp^,p); { product restored on l^.sp }
      l:=l^.post
   end
end;

function reorganiseAT(list: ^AT): ^AT;
var l,lnew,a: ^AT; label 5;
begin lnew:=NIL; l:=list;
   while l<>NIL do with l^ do
   begin SparsifyPolynomial(sp);
      if sp.sp=NIL then goto 5; { zero function ignored }
      { add l-term to lnew }
      a:=FindOrCreateAT(lnew,fct,sx,sy,sxx,syy,x,y,xx,yy);
      P_plus_Q(a^.sp,sp);
5:    l:=post
   end;
   DeleteAList(list); reorganiseAT:=lnew
end;

function AddFct(l1,l2: ^AT): ^AT;
{ l1 := l1 + l2 }
var l,lnew,a: ^AT; label 5;
begin l:=l2;
   while l<>NIL do with l^ do
   begin
      a:=FindOrCreateAT(l1,fct,sx,sy,sxx,syy,x,y,xx,yy);
      P_plus_Q(a^.sp,sp);
5:    l:=post
   end;
   AddFct:=l1
end;

function xIntegral(list: ^AT): ^AT;
{ all functions from the list are integrated w.r.t. x.
  "list" is deleted after the computation.
  The resulting list is the output value}
var aa,ab,ad,ae,al,ag,af,am,ap,aq,ar,ax,listx: ^AT; s: string; ix,iy: integer; label 8,9;
begin xintegration:=xintegration+1; listx:=NIL; if list=NIL then goto 8;
   if list^.sx='$' then
   begin writeln('### functions already evaluated w.r.t. x'); goto 9 end;

					{ P -x-> P }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if (fct='P') or (fct='Pi') then 
      begin ap:=FindOrCreateAT(listx,fct,sx,sy,sxx,syy,x,y,xx,yy);
            Px(ap^.sp, sp)
      end;
      ax:=post
   end;
					{ Ps -x-> Ps }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='Ps' then 
      begin ap:=FindOrCreateAT(listx,fct,sx,sy,sxx,syy,x,y,xx,yy);
            Psx(ap^.sp, sp)
      end;
      ax:=post
   end;
					{ BZ -x-> BZ (AZ) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='BZ' then 
      begin aa:=FindOrCreateAT(list ,'AZ',sx,sy,sxx,syy,x,y,xx,yy);
            ab:=FindOrCreateAT(listx,'BZ',sx,sy,sxx,syy,x,y,xx,yy);
            Bx_Ax(ab^.sp, aa^.sp, sp)
      end;
      ax:=post
   end;
{	---------------------------------------     M -x-> M (A)
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='M' then 
      begin aa:=FindOrCreateAT(list ,'A',sx,sy,sxx,syy,x,y,xx,yy);
            am:=FindOrCreateAT(listx,'M',sx,sy,sxx,syy,x,y,xx,yy);
            Mx_Ax(am^.sp, aa^.sp, sp)
      end;
      ax:=post
   end;                                  ------     not in use }
					{ M -x-> M (DA) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='M' then 
      begin aa:=FindOrCreateAT(list,'DA',sx,sy,sxx,syy,x,y,xx,yy);
            am:=FindOrCreateAT(listx,'M',sx,sy,sxx,syy,x,y,xx,yy);
            Lx_Fx(am^.sp, aa^.sp, sp)
      end;
      ax:=post
   end;

					{ DA -x-> M,P (F,AZ) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='DA' then
      begin if sxx='$' then ix:=1 else ix:=0;
            aa:=FindOrCreateAT(list,'AZ',sx,sy,sxx,syy,x,y,xx,yy);
            am:=FindOrCreateAT(listx,'M',sx,sy,sxx,syy,x,y,xx,yy);
            ap:=FindOrCreateAT(listx,'P',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list ,'F',sx,sy,sxx,syy,x,y,xx,yy);
            DAx_Fx(am^.sp, ap^.sp, af^.sp, aa^.sp, sp, ix, xx)
      end;
      ax:=post
   end;

					{ AZ -x-> B,M,P (F) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='AZ' then
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            ab:=FindOrCreateAT(listx,'BZ',sx,sy,sxx,syy,x,y,xx,yy);
            am:=FindOrCreateAT(listx,'M',sx,sy,sxx,syy,x,y,xx,yy);
            ap:=FindOrCreateAT(listx,'P',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list ,'F',sx,sy,sxx,syy,x,y,xx,yy);
            AZx_Fx(ab^.sp, am^.sp, ap^.sp, af^.sp, sp, ix, iy, xx, yy)
      end;
      ax:=post
   end;
{------------------------------------------------------------  A -x-> B,M,P (F)
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='A' then
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            ab:=FindOrCreateAT(listx,'B',sx,sy,sxx,syy,x,y,xx,yy);
            am:=FindOrCreateAT(listx,'M',sx,sy,sxx,syy,x,y,xx,yy);
            ap:=FindOrCreateAT(listx,'P',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list ,'F',sx,sy,sxx,syy,x,y,xx,yy);
            Ax_Fx(ab^.sp, am^.sp, ap^.sp, af^.sp, sp, ix, iy, xx, yy)
      end;
      ax:=post
   end;
----------------------------------------------------------------- not in use }
					{ G -x-> G (F) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='G' then 
      begin af:=FindOrCreateAT(list ,'F',sx,sy,sxx,syy,x,y,xx,yy);
            ag:=FindOrCreateAT(listx,'G',sx,sy,sxx,syy,x,y,xx,yy);
            Gx_Fx(ag^.sp, af^.sp, sp)
      end;
      ax:=post
   end;
   					{ L -x-> L (F) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='L' then 
      begin af:=FindOrCreateAT(list ,'F',sx,sy,sxx,syy,x,y,xx,yy);
            al:=FindOrCreateAT(listx,'L',sx,sy,sxx,syy,x,y,xx,yy);
            Lx_Fx(al^.sp, af^.sp, sp)
      end;
      ax:=post
   end;
					{ C0 -x-> (Q) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='C0' then 
      begin aq:=FindOrCreateAT(list,'Q',sx,sy,sxx,syy,x,y,xx,yy);
            C0x_Qx(aq^.sp, sp)
      end;
      ax:=post
   end;
					{ Cp -x-> (Rp) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='Cp' then 
      begin ar:=FindOrCreateAT(list,'Rp',sx,sy,sxx,syy,x,y,xx,yy);
            Cpx_Rpx(ar^.sp, sp)
      end;
      ax:=post
   end;
					{ Cm -x-> (Rm) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='Cm' then 
      begin ar:=FindOrCreateAT(list,'Rm',sx,sy,sxx,syy,x,y,xx,yy);
            Cmx_Rmx(ar^.sp, sp)
      end;
      ax:=post
   end;
					{ Rp -x-> Rp,Ep (F) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='Rp' then 
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            aa:=FindOrCreateAT(list,'Ep',sx,sy,sxx,syy,x,y,xx,yy);
            ar:=FindOrCreateAT(listx,'Rp',sx,sy,sxx,syy,x,y,xx,yy);
            ae:=FindOrCreateAT(listx,'Ep',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list,'F',sx,sy,sxx,syy,x,y,xx,yy);
            RxEx_Fx(ar^.sp,ae^.sp,sp,aa^.sp,af^.sp,ix,iy,1,xx,yy);
      end;
      ax:=post
   end;
					{ Ep -x-> Rp,Ep (F) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='Ep' then 
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            aa:=FindOrCreateAT(list,'Rp',sx,sy,sxx,syy,x,y,xx,yy);
            ar:=FindOrCreateAT(listx,'Rp',sx,sy,sxx,syy,x,y,xx,yy);
            ae:=FindOrCreateAT(listx,'Ep',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list,'F',sx,sy,sxx,syy,x,y,xx,yy);
            RxEx_Fx(ar^.sp,ae^.sp,aa^.sp,sp,af^.sp,ix,iy,1,xx,yy);
      end;
      ax:=post
   end;
					{ Rm -x-> Rm,Em (F) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='Rm' then 
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            aa:=FindOrCreateAT(list,'Em',sx,sy,sxx,syy,x,y,xx,yy);
            ar:=FindOrCreateAT(listx,'Rm',sx,sy,sxx,syy,x,y,xx,yy);
            ae:=FindOrCreateAT(listx,'Em',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list,'F',sx,sy,sxx,syy,x,y,xx,yy);
            RxEx_Fx(ar^.sp,ae^.sp,sp,aa^.sp,af^.sp,ix,iy,-1,xx,yy);
      end;
      ax:=post
   end;
					{ Em -x-> Rm,Em (F) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='Em' then 
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            aa:=FindOrCreateAT(list,'Rm',sx,sy,sxx,syy,x,y,xx,yy);
            ar:=FindOrCreateAT(listx,'Rm',sx,sy,sxx,syy,x,y,xx,yy);
            ae:=FindOrCreateAT(listx,'Em',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list,'F',sx,sy,sxx,syy,x,y,xx,yy);
            RxEx_Fx(ar^.sp,ae^.sp,aa^.sp,sp,af^.sp,ix,iy,-1,xx,yy);
      end;
      ax:=post
   end;
					{ F -x-> L,G }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='F' then 
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            ag:=FindOrCreateAT(listx,'G',sx,sy,sxx,syy,x,y,xx,yy);
            al:=FindOrCreateAT(listx,'L',sx,sy,sxx,syy,x,y,xx,yy);
            Fx(ag^.sp, al^.sp, sp, ix, iy, xx, yy)
      end;
      ax:=post
   end;
   					{ Q -x-> Q (D) }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='Q' then 
      begin ad:=FindOrCreateAT(list ,'D',sx,sy,sxx,syy,x,y,xx,yy);
            aq:=FindOrCreateAT(listx,'Q',sx,sy,sxx,syy,x,y,xx,yy);
            Qx_Dx(aq^.sp, ad^.sp, sp)
      end;
      ax:=post
   end;
					{ D -x-> Q,M,P }
   ax:=list; while ax<>NIL do with ax^ do
   begin if sp.sp<>NIL then if fct='D' then 
      begin if syy='$' then iy:=1 else iy:=0;
            ap:=FindOrCreateAT(listx,'P',sx,sy,sxx,syy,x,y,xx,yy);
            aq:=FindOrCreateAT(listx,'Q',sx,sy,sxx,syy,x,y,xx,yy);
            am:=FindOrCreateAT(listx,'M',sx,sy,sxx,syy,x,y,xx,yy);
            Dx(aq^.sp, am^.sp, ap^.sp, sp, iy, yy)
      end;
      ax:=post
   end;

8: s:='### xIntegration performed'; write_P(s); writeln_P;
9: list:=DeleteAT(list); xIntegral:=listx
end; { xIntegral }

function yIntegral(list: ^AT): ^AT;
{ all functions from the list are integrated w.r.t. y.
  "list" is deleted after the computation.
  The resulting list is the output value}
var aa,ab,ad,ae,al,ag,af,am,ap,aq,ar,ay,listy: ^AT; s: string; ix,iy: integer; label 8,9;
begin yintegration:=yintegration+1; listy:=NIL; if list=NIL then goto 8;
   if list^.sy='$' then
   begin writeln('### functions already evaluated w.r.t. y'); goto 9 end;

					{ P -y-> P }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if (fct='P') or (fct='Pi') then 
      begin ap:=FindOrCreateAT(listy,fct,sx,sy,sxx,syy,x,y,xx,yy);
            Py(ap^.sp, sp)
      end;
      ay:=post
   end;
					{ Ps -y-> Ps }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='Ps' then 
      begin ap:=FindOrCreateAT(listy,fct,sx,sy,sxx,syy,x,y,xx,yy);
            Psy(ap^.sp, sp)
      end;
      ay:=post
   end;

					{ BZ -y-> BZ (AZ) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='BZ' then 
      begin aa:=FindOrCreateAT(list ,'AZ',sx,sy,sxx,syy,x,y,xx,yy);
            ab:=FindOrCreateAT(listy,'BZ',sx,sy,sxx,syy,x,y,xx,yy);
            By_Ay(ab^.sp, aa^.sp, sp)
      end;
      ay:=post
   end;
{	---------------------------------------    M -y-> M (DA)
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='M' then 
      begin aa:=FindOrCreateAT(list ,'A',sx,sy,sxx,syy,x,y,xx,yy);
            am:=FindOrCreateAT(listy,'M',sx,sy,sxx,syy,x,y,xx,yy);
            My_Ay(am^.sp, aa^.sp, sp)
      end;
      ay:=post
   end;                                 ------     not in use }
					{ M -y-> M (DA) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='M' then 
      begin aa:=FindOrCreateAT(list,'DA',sx,sy,sxx,syy,x,y,xx,yy);
            am:=FindOrCreateAT(listy,'M',sx,sy,sxx,syy,x,y,xx,yy);
            Ly_Fy(am^.sp, aa^.sp, sp)
      end;
      ay:=post
   end;
					{ DA -y-> M,P (F,AZ) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='DA' then
      begin if sxx='$' then ix:=1 else ix:=0;
            aa:=FindOrCreateAT(list,'AZ',sx,sy,sxx,syy,x,y,xx,yy);
            am:=FindOrCreateAT(listy,'M',sx,sy,sxx,syy,x,y,xx,yy);
            ap:=FindOrCreateAT(listy,'P',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list ,'F',sx,sy,sxx,syy,x,y,xx,yy);
            DAy_Fy(am^.sp, ap^.sp, af^.sp, aa^.sp, sp, ix, xx)
      end;
      ay:=post
   end;
					{ AZ -y-> BZ,M,P (F) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='AZ' then
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            ab:=FindOrCreateAT(listy,'BZ',sx,sy,sxx,syy,x,y,xx,yy);
            am:=FindOrCreateAT(listy,'M',sx,sy,sxx,syy,x,y,xx,yy);
            ap:=FindOrCreateAT(listy,'P',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list ,'F',sx,sy,sxx,syy,x,y,xx,yy);
            AZy_Fy(ab^.sp, am^.sp, ap^.sp, af^.sp, sp, ix, iy, xx, yy)
      end;
      ay:=post
   end;
{------------------------------------------------------------ A -y-> B,M,P (F)
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='A' then
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            ab:=FindOrCreateAT(listy,'B',sx,sy,sxx,syy,x,y,xx,yy);
            am:=FindOrCreateAT(listy,'M',sx,sy,sxx,syy,x,y,xx,yy);
            ap:=FindOrCreateAT(listy,'P',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list ,'F',sx,sy,sxx,syy,x,y,xx,yy);
            Ay_Fy(ab^.sp, am^.sp, ap^.sp, af^.sp, sp, ix, iy, xx, yy)
      end;
      ay:=post
   end;  --------------------------------------------------------- not in use }
					{ G -y-> G (F) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='G' then 
      begin af:=FindOrCreateAT(list ,'F',sx,sy,sxx,syy,x,y,xx,yy);
            ag:=FindOrCreateAT(listy,'G',sx,sy,sxx,syy,x,y,xx,yy);
            Gy_Fy(ag^.sp, af^.sp, sp)
      end;
      ay:=post
   end;

   					{ L -y-> L (F) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if ay^.fct='L' then 
      begin af:=FindOrCreateAT(list ,'F',sx,sy,sxx,syy,x,y,xx,yy);
            al:=FindOrCreateAT(listy,'L',sx,sy,sxx,syy,x,y,xx,yy);
            Ly_Fy(al^.sp, af^.sp, ay^.sp)
      end;
      ay:=post
   end;
					{ C0 -y-> (Q) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='C0' then 
      begin aq:=FindOrCreateAT(list,'Q',sx,sy,sxx,syy,x,y,xx,yy);
            C0y_Qy(aq^.sp, sp)
      end;
      ay:=post
   end;
					{ Cp -y-> (Rp) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='Cp' then 
      begin ar:=FindOrCreateAT(list,'Rp',sx,sy,sxx,syy,x,y,xx,yy);
            Cpy_Rpy(ar^.sp, sp)
      end;
      ay:=post
   end;
					{ Cm -y-> (Rm) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='Cm' then 
      begin ar:=FindOrCreateAT(list,'Rm',sx,sy,sxx,syy,x,y,xx,yy);
            Cmy_Rmy(ar^.sp, sp)
      end;
      ay:=post
   end;
					{ Rp -y-> Rp,Ep (F) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='Rp' then 
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            aa:=FindOrCreateAT(list,'Ep',sx,sy,sxx,syy,x,y,xx,yy);
            ar:=FindOrCreateAT(listy,'Rp',sx,sy,sxx,syy,x,y,xx,yy);
            ae:=FindOrCreateAT(listy,'Ep',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list,'F',sx,sy,sxx,syy,x,y,xx,yy);
            RyEy_Fy(ar^.sp,ae^.sp,sp,aa^.sp,af^.sp,ix,iy,1,xx,yy);
      end;
      ay:=post
   end;
					{ Ep -y-> Rp,Ep (F) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='Ep' then 
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            aa:=FindOrCreateAT(list,'Rp',sx,sy,sxx,syy,x,y,xx,yy);
            ar:=FindOrCreateAT(listy,'Rp',sx,sy,sxx,syy,x,y,xx,yy);
            ae:=FindOrCreateAT(listy,'Ep',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list,'F',sx,sy,sxx,syy,x,y,xx,yy);
            RyEy_Fy(ar^.sp,ae^.sp,aa^.sp,sp,af^.sp,ix,iy,1,xx,yy);
      end;
      ay:=post
   end;
					{ Rm -y-> Rm,Em (F) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='Rm' then 
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            aa:=FindOrCreateAT(list,'Em',sx,sy,sxx,syy,x,y,xx,yy);
            ar:=FindOrCreateAT(listy,'Rm',sx,sy,sxx,syy,x,y,xx,yy);
            ae:=FindOrCreateAT(listy,'Em',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list,'F',sx,sy,sxx,syy,x,y,xx,yy);
            RyEy_Fy(ar^.sp,ae^.sp,sp,aa^.sp,af^.sp,ix,iy,-1,xx,yy);
      end;
      ay:=post
   end;
					{ Em -y-> Rm,Em (F) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='Em' then 
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            aa:=FindOrCreateAT(list,'Rm',sx,sy,sxx,syy,x,y,xx,yy);
            ar:=FindOrCreateAT(listy,'Rm',sx,sy,sxx,syy,x,y,xx,yy);
            ae:=FindOrCreateAT(listy,'Em',sx,sy,sxx,syy,x,y,xx,yy);
            af:=FindOrCreateAT(list,'F',sx,sy,sxx,syy,x,y,xx,yy);
            RyEy_Fy(ar^.sp,ae^.sp,aa^.sp,sp,af^.sp,ix,iy,-1,xx,yy);
      end;
      ay:=post
   end;
					{ F -y-> L,G }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if ay^.fct='F' then 
      begin if sxx='$' then ix:=1 else ix:=0; if syy='$' then iy:=1 else iy:=0;
            ag:=FindOrCreateAT(listy,'G',sx,sy,sxx,syy,x,y,xx,yy);
            al:=FindOrCreateAT(listy,'L',sx,sy,sxx,syy,x,y,xx,yy);
            Fy(ag^.sp, al^.sp, sp, ix, iy, xx, yy)
      end;
      ay:=post
   end;
   					{ Q -y-> Q (D) }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='Q' then 
      begin ad:=FindOrCreateAT(list ,'D',sx,sy,sxx,syy,x,y,xx,yy);
            aq:=FindOrCreateAT(listy,'Q',sx,sy,sxx,syy,x,y,xx,yy);
            Qy_Dy(aq^.sp, ad^.sp, sp)
      end;
      ay:=post
   end;
					{ D -y-> Q,M,P }
   ay:=list; while ay<>NIL do with ay^ do
   begin if sp.sp<>NIL then if fct='D' then 
      begin if syy='$' then iy:=1 else iy:=0;
            ap:=FindOrCreateAT(listy,'P',sx,sy,sxx,syy,x,y,xx,yy);
            aq:=FindOrCreateAT(listy,'Q',sx,sy,sxx,syy,x,y,xx,yy);
            am:=FindOrCreateAT(listy,'M',sx,sy,'$',syy,x,y,0,yy);
            Dy(aq^.sp, am^.sp, ap^.sp, sp, iy, yy)
      end;
      ay:=post
   end;

8: list:=DeleteAT(list); yIntegral:=listy;
   s:='### yIntegration performed'; write_P(s); writeln_P;
9:
end;

procedure NameVar(var s: string; VariableNr: integer; l: ^AT);
begin if l<>NIL then with l^ do case VariableNr of
1:	s:=sx;
2:	s:=sy;
3:	s:=sxx;
4:	s:=syy
end end;

procedure replaceAT(list: ^AT; VariableNr: integer; a,b,c: real);
var l: ^AT; s,sxy,syx: string; v: integer; label 8,9;
begin l:=list; if l=NIL then goto 9;
   NameVar(sxy,VariableNr,list);
   v:=VariableNr; if v<=2 then v:=3-v else v:=7-v;
   NameVar(syx,v,list);
   while l<>NIL do
   begin NameVar(s,VariableNr,l); if s='$' then goto 8;
      ReplaceByQuadratic(l^.sp,VariableNr,a,b,c);
      case VariableNr of
1:	begin l^.sx :='$'; l^.x :=c end;
2:	begin l^.sy :='$'; l^.y :=c end;
3:	begin l^.sxx:='$'; l^.xx:=c end;
4:	begin l^.syy:='$'; l^.yy:=c end;
      end;
8:    l:=l^.post
   end;
   writev(s,'### replacement of ',sxy,' by ',c); write_P(s);
   if a<>0 then
   begin
      writev(s,' + ',a,' * ',syx); write_P(s)
   end;
   writeln_P;
9:
end;

procedure replaceATconstant(list: ^AT; VariableNr,par: integer; c: real);
var l: ^AT;
begin if par=0 then
   begin writeln('### input for replacement ...');
         write('    -> constant = '); readln(c)
   end;
   replaceAT(list,VariableNr,0,0,c)
end;

procedure replaceATlinear(list: ^AT; VariableNr: integer);
var l: ^AT; a,c: real;
begin
   writeln('### input of the linear polynomial a*z+c for replacement ...');
   write('    -> a = '); readln(a);
   write('    -> c = '); readln(c);
   replaceAT(list,VariableNr,0,a,c)
end;

procedure DifferenceAT(var list: ^AT; l2: ^AT);
{ list := list - l2 }
var l,ll: ^AT;
begin l:=l2; while l<>NIL do with l^ do
    begin ll:=FindOrCreateAT(list,fct,sx,sy,sxx,syy,x,y,xx,yy);
          P_minus_Q(ll^.sp,sp);
          l:=post
end end;

procedure evaluateATatLinearBounds(var list: ^AT; VarNr,par: integer; cl,cr: real);
var lc: ^AT;
begin CopyAList(lc,list);
   writeln('#### evaluation at right side ...');
   replaceATconstant(list,VarNr,par,cr);
   writeln('#### evaluation at left side ...');
   replaceATconstant(lc,VarNr,par,cl);
   DifferenceAT(list,lc); DeleteAT(lc)
end;

procedure evaluateBothBounds(var list: ^AT; par: integer; xl,xr,yl,yr: real);
begin
   writeln('###### evaluation for x-variable ...');
   evaluateATatLinearBounds(list,1,par,xl,xr);
   writeln('###### evaluation for y-variable ...');
   evaluateATatLinearBounds(list,2,par,yl,yr)
end;

procedure InterchangeXYFct(a: ^AT);
var z: real; s: string;
begin if a<>NIL then with a^ do
    begin interchangeXYPolynomial(sp);
       z:= xx;  xx:= yy;  yy:=z;
       s:=sxx; sxx:=syy; syy:=s
end end;

procedure stabilise(var list: ^AT);
var a,add,l,lnew: ^AT; z: real; p: sparsePolynomial; pc: ^sparsePolynomial;
    n: integer; s: string; t: IT; label 5;
begin t[1]:=0; t[1]:=1;
   repeat lnew:=NIL; add:=NIL; l:=list; n:=0;
   s:='### stabilise: started'; write_P(s); writeln_P; s:='stabilise: ';
   while l<>NIL do with l^ do
   begin SparsifyPolynomial(sp); if sp.sp=NIL then goto 5; { term = zero }

      { in the antisymmetric case, x=y leads to a zero term }
      if (sx='$') and (sy='$') then if abs(x-y)<=Tolerance then
        if (fct='BZ') or (fct='C') or (fct='Q')
           or (fct='Rm') or (fct='Rp') or (fct='Ps') then 
           begin writeln(s,fct,'-term omitted since x=y'); goto 5 end;
      
      { expressing BZ by Y*(Q-Rm) }
      if (fct='BZ') and (syy='$') then
      begin fct:='Rm'; sxP(sp,-yy); n:=n+1;
            add:=CopyAT(l); sxP(add^.sp,yy); add^.fct:='Q';
            writeln(s,'BZ expressed by Y*(Q-Rm) #### -> additional Q-term')
      end;

      { in the antisymmetric case, the sign is changed when x,y are interchanged }
      if (sx='$') and (sy='$') then if x<y then
        if (fct='BZ') or (fct='C') or (fct='Q')
           or (fct='Rm') or (fct='Rp') or (fct='Ps') then
           begin z:=x; x:=y; y:=z; sxP(sp,-1); n:=n+1;
                 writeln(s,'x,y interchanged and sign reversed in ',fct)
           end;

      { P, Ps, C0, D and Q do not depend on X }
      if sxx='$' then if xx<>0 then
        if (fct='C0') or (fct='Q') or (fct='P') or (fct='Ps') or (fct='Pi') or (fct='D') then
        begin xx:=0; writeln(s,'xx:=0 in ',fct) end;

      { P,Ps does not depend on Y }
      if syy='$' then if (fct='P') or (fct='Ps') or (fct='Pi') then if yy<>0 then
        begin yy:=0; writeln(s,'yy:=0 in ',fct) end;

      { Ep,Em,Rm,Cm=0 if xx=0 }
      if (sxx='$') and (xx=0) then
        if (fct='Ep') or (fct='Em') or (fct='Rm') or (fct='Cm') then 
        begin writeln(s,fct,'-term omitted since xx=0'); goto 5 end;

      { Ep,Em,Rp,Cp=0 if yy=0 }
      if (syy='$') and (yy=0) then
        if (fct='Ep') or (fct='Em') or (fct='Rp') or (fct='Cp') then
        begin writeln(s,fct,'-term omitted since yy=0'); goto 5 end;

      { ordering of the values xx and yy }
      if g34=1 then if (sxx='$') and (syy='$') then if xx>yy then
        if (fct='F') or (fct='G') or (fct='L') then
           begin z:=xx; xx:=yy; yy:=z; n:=n+1;
                 writeln(s,'xx,yy interchanged for ',fct)
           end;

      { expressing Ep,Em by M due to (2.12),(2.13) }
      if (sxx='$') and (syy='$') then if (fct='Ep') or (fct='Em') then
      begin if fct='Ep' then InterchangeXYFct(l);
         writeln(s,fct,' expressed by M  #### -> additional M-term');
         fct:='M'; n:=n+1; add:=CopyAT(l); sxP(sp,-yy); sxP(add^.sp,yy); add^.xx:=0
      end;

      { expressing M by L }
      if m12=1 then if fct='M' then if (sx='$') and (sy='$') and (sxx='$') then
      if x-y>=0 then
      begin writeln(s,fct,' expressed by L');
         fct:='L'; n:=n+1; z:=x-y; x:=xx; y:=0; xx:=z
      end;

      { if the x,y variables are evaluated, certain functions are symmetric
        w.r.t. the difference x-y. In these cases a negative value x-y is changed
        into a positive one by interchanging the x and y values. }
      if (sx='$') and (sy='$') then if x<y then
        if (fct='F') or (fct='G') or (fct='A') or (fct='D') 
           or (fct='Em') or (fct='Ep') or (fct='M') or (fct='C') 
           or (fct='Cp') or (fct='Cm') or (fct='C0') then
           begin z:=x; x:=y; y:=z; n:=n+1; writeln(s,'x,y interchanged in ',fct)
           end;

      { the limits X->0 for Rp }
      if sxx='$' then if abs(xx)<=Tolerance then
      begin if fct='Rp' then
          begin fct:='Ps'; n:=n+1; writeln(s,'Rp yields Ps since xx=0') end 
          else if fct='Cp' then
          begin StartP(p,1); n:=n+1;
             DefineCoefficient1_real(p,1,0,1); DefineCoefficient1_real(p,0,1,-1);
             pc:=CopyPolynomial(sp); DeletePolynomial(sp);
             multiply(sp,pc^,p); fct:='Ps';
             DeletePolynomialPointer(pc); DeletePolynomial(p);
             writeln(s,'Cp yields Ps since xx=0')
      end end;

      { Ps expressed by Pi }
      if fct='Ps' then if (sx='$') and (sy='$') then
      begin z:=1/2; n:=n+1;
	 if x<y then z:=-z else if x=y then z:=0; fct:='Pi'; sxP(sp,z);
         writeln(s,'Ps expressed by Pi')
      end;

      { Cm,Em expressed by Cp,Ep }
      if (fct='Cm') or (fct='Em') or ( fct='Rm') then
      if ((cem=1) and (sxx='$')  and (syy='$')) or (cXY=1) then
      begin InterchangeXYFct(l); n:=n+1;
         if fct='Rm' then
         begin fct:='Rp'; writeln(s,'Rm expressed by Rp') end else
         if fct='Cm' then
         begin fct:='Cp'; writeln(s,'Cm expressed by Cp') end else
         if fct='Em' then
         begin fct:='Ep'; writeln(s,'Em expressed by Ep') end
      end;

      { if x<y, function L is stabilised by L(-z;...) = - L(z;...) + 2*L(0;...) }
      if (sx='$') and (sy='$') then if x<y then if fct='L' then
        begin add:=CopyAT(l); { becomes additional term 2*L(0;...) }
           sxP(add^.sp,2); add^.x:=0; add^.y:=0;
           z:=x; x:=y; y:=z; sxP(sp,-1); n:=n+1;
           writeln(s,'L(-z;...) replaced by L(z;...) and L(0;...)   ### -> additional L-term')
        end;

      { if X<0, function M is stabilised by M(.;X,Y) = - M(.;-X,Y) + 2*M(.;0,Y) }
      if (sxx='$') and (xx<0) then if fct='M' then
        begin add:=CopyAT(l); { becomes additional term 2*M(.;0,Y) }
           sxP(add^.sp,2); add^.xx:=0;
           sxP(sp,-1); xx:=-xx; n:=n+1;
           writeln(s,'M(.;X,Y) = - M(.;-X,Y) + 2*M(.;0,Y)   ### -> additional M-term')
        end;

      { add possibly unchanged l-term to lnew }
      a:=FindOrCreateAT(lnew,fct,sx,sy,sxx,syy,x,y,xx,yy);
      P_plus_Q(a^.sp,l^.sp);

5:    if add<>NIL then
      begin a:=FindOrCreateAT(lnew,add^.fct,add^.sx,add^.sy,add^.sxx,add^.syy,
                              add^.x,add^.y,add^.xx,add^.yy);
            P_plus_Q(a^.sp,add^.sp); add:=DeleteAT(add);
            writeln(s,'additional term added')
      end;

      l:=post
   end;
   DeleteAList(list); list:=lnew
   until n=0
end;

procedure interpret12(var list: ^AT);
var a,l,lnew: ^AT; s: string; label 5;
begin lnew:=NIL; l:=list;
   s:='### interpret12 performed'; write_P(s); writeln_P; s:='### ERROR in Interpret12: ';
   while l<>NIL do with l^ do
   begin
      if sxx='$' then
      begin writeln(s,'Y variable (=x2-y2) already evaluated');
            goto 5
      end;
      if sx<>'$' then
      begin writeln(s,'x variable not yet evaluated'); goto 5 end;
      if sy<>'$' then
      begin writeln(s,'y variable not yet evaluated'); goto 5 end;

{G}   if fct='G' then
      begin sxx:='$'; xx:=x-y; writeln('interpret12: G -> G') end else

{L}   if fct='L' then
      begin fct:='M'; sxx:='$'; xx:=x-y; writeln('interpret12: L -> M') end else

      begin writeln(s,'The function ',fct,' should not appear');
            goto 5
      end;

      reinterpretationPolynomial(sp,3,4);
      sx:='x'; sy:='y'; x:=0; y:=0; { primary variables are again free variables }
      a:=FindOrCreateAT(lnew,fct,sx,sy,sxx,syy,x,y,xx,yy);
      P_plus_Q(a^.sp,sp);

5:    l:=post
   end;
   DeleteAList(list); list:=lnew; stabilise(list)
end;

procedure interpret23(var list: ^AT);
var a,l,lnew: ^AT; s: string; label 5;
begin lnew:=NIL; l:=list;
   s:='### interpret23 performed'; write_P(s); writeln_P; s:='### ERROR in Interpret23: ';
   while l<>NIL do with l^ do
   begin
      if sxx<>'$' then
      begin writeln(s,'X variable (=x1-y1) not yet evaluated');
            goto 5
      end;
      if syy='$' then
      begin writeln(s,'Z variable (=x3-y3) already evaluated');
            goto 5
      end;
      if sx<>'$' then
      begin writeln(s,'x variable not yet evaluated'); goto 5 end;
      if sy<>'$' then
      begin writeln(s,'y variable not yet evaluated'); goto 5 end;

{G,Cp,Cm}
      if (fct='G') or (fct='Cp') or (fct='Cm') then
      begin syy:='$'; yy:=x-y; if fct='G' then writeln('interpret23: G -> G') end else

{L}   if fct='L' then
      begin fct:='M'; syy:=sxx; yy:=xx; sxx:='$'; xx:=x-y; writeln('interpret23: L -> M') end
      else

{BZ}  if fct='BZ' then  {BZ -> Cp+Cm-C0}
      begin a:=CopyAT(l); a^.fct:='Cp'; InsertAT(list,a);
         a:=CopyAT(l); a^.fct:='Cm'; InsertAT(list,a);
         sxP(sp,-1); fct:='C0'; syy:='$'; sxx:='$'; xx:=0; yy:=x-y;
	 writeln('interpret23: BZ -> Cp,Cm,C0')
      end else

{M}   if fct='M' then
      begin fct:='M'; syy:='$'; yy:=x-y; writeln('interpret23: M -> M') end else

{P}   if fct='P' then writeln('interpret23: P -> P') else

      begin writeln(s,'The function ',fct,' should not appear');
            goto 5
      end;

      reinterpretationPolynomial(sp,4,3);
      sx:='x'; sy:='y'; x:=0; y:=0; { primary variables are again free variables }
      a:=FindOrCreateAT(lnew,fct,sx,sy,sxx,syy,x,y,xx,yy);
      P_plus_Q(a^.sp,sp);

5:    l:=post
   end;
   DeleteAList(list); list:=lnew; stabilise(list)
end;

procedure WriteExFct(l: ^AT);
begin writeln(ES,phase,xintegration,yintegration);
   while l<>NIL do with l^ do
   begin writeln(ES,'Function');
         writeln(ES,fct);
         writeln(ES,sx);
         writeln(ES,sy);
         writeln(ES,sxx);
         writeln(ES,syy);
         writeln(ES,x:25:22,' ',y:25:22,' ',xx:25:22,' ',yy:25:22);
         WriteExPolynomial(sp);
         l:=post
   end;
   writeln(ES,'NIL');
end;

procedure WriteTEX(l: ^AT);
var s: string; label 1,5;
begin writeln(TX,'% begin of expression'); writeln(TX,'$');
   while l<>NIL do with l^ do
   begin if fct='Ps' then 
         begin TeX_Ps(s,sp,x-y); 
1:             if s<>'' then writeln(TX,s,' \allowbreak'); goto 5
         end;
         if fct='P'  then begin TeX_P(s,sp); goto 1 end;
         if fct='Pi' then begin TeX_PiE(s,sp); goto 1 end;
         if fct='Ps' then begin TeX_Ps(s,sp,x-y); goto 1 end;
         if fct='Q'  then begin TeX_Q(s,sp,x-y,yy); goto 1 end;
         if fct='F'  then begin TeX_F(s,sp,x-y,xx,yy); goto 1 end;
         if fct='G'  then begin TeX_G(s,sp,x-y,xx,yy); goto 1 end;
         if fct='L'  then begin TeX_L(s,sp,x-y,xx,yy); goto 1 end;
         if fct='M'  then begin TeX_M(s,sp,x-y,xx,yy); goto 1 end;
         if fct='Rp' then begin TeX_Rp(s,sp,x-y,xx,yy); goto 1 end;
         if fct='Rm' then begin TeX_Rm(s,sp,x-y,xx,yy); goto 1 end;
         writeln('#### WriteTEX - output of ',fct,' not yet implemented');
5:       l:=post
   end;
   writeln(TX,'$'); writeln(TX,'% end of expression');
end;

function ReadExFct: ^AT;
var list,v,l: ^AT; s: string; p: ^sparsePolynomial; label 1,9;
begin list:=NIL; l:=NIL; readln(ES,phase,xintegration,yintegration);
1: readln(ES,s); if s='NIL' then goto 9; v:=l;
   new(l); if v=NIL then list:=l else v^.post:=l;
   with l^ do
   begin readln(ES,fct);
         readln(ES,sx);
         readln(ES,sy);
         readln(ES,sxx);
         readln(ES,syy);
         readln(ES,x,y,xx,yy);
         ReadExPolynomial(sp);
         post:=NIL
   end;
   goto 1;
9: ReadExFct:=list
end;

function arctan2arg(x,y: real): real; {arctan(x/y}
var a: real;
begin a:=0; if (x=0) and (y=0) then writeln('ERROR: arctan(x/y) undefined for x=y=0') else
   if x<>0 then 
   begin if y<>0 then a:=arctan(x/y) else
       begin writeln('WARNING: arctan(x/0) = sign(x)*pi/2 used');
           if x>0 then a:=1.570796326794897 else a:=-1.570796326794897
   end end;
   arctan2arg:=a
end;

function evaluateElementaryFct(var f: string; xy,xx,yy: real): real;
var z,w: real; e: string; label 9;
begin z:=0; w:=sqrt(sqr(xy)+sqr(xx)+sqr(yy)); e:='ERROR in evaluateElementaryFct: ';

{P} if f='P' then begin z:=1; goto 9 end;

{Pi}if f='Pi' then begin z:=3.141592653589793; goto 9 end;

{Ps}if f='Ps' then 
    begin z:=1.570796326794897;
          if xy<0 then z:=-z else if xy=0 then z:=0; goto 9
    end;

{G} if f='G' then begin z:=w; goto 9 end;

{L} if f='L' then
    begin z:=xy+w; if z<=0 then writeln(e,'ln(0) in L') else z:=ln(z);
          goto 9
    end;

{M} if f='M' then
    begin z:=xx+w; if z<=0 then writeln(e,'ln(0) in M') else z:=ln(z);
          goto 9
    end;

{Q} if f='Q' then
    begin if (yy=0) and (xy=0) then writeln(e,'0/0 in Q')
          else z:=arctan2arg(xy,yy);
      goto 9
    end;

{D} if f='D' then
    begin if (yy=0) and (xy=0) then writeln(e,'0/0 in D')
          else z:=yy/(sqr(xy)+sqr(yy));
      goto 9
    end;

{C0}if f='C0' then begin if xy<>0 then z:=xy*arctan2arg(xy,yy); goto 9 end;

{Rp}if f='Rp' then
    begin z:=yy*xy;
      if (z=0) and (xx=0) then writeln(e,'arctan(0/0) in Rp')
      else if z<>0 then z:=arctan2arg(z,xx*w);
      goto 9
    end;

{F} if f='F' then
    begin z:=sqr(xy)+sqr(xx)+sqr(yy);
      if z=0 then writeln(e,'Division by zero in F') else z:=1/sqrt(z);
      goto 9
    end;

{B} if f='BZ' then
    begin z:=xy/yy; if z<>0 then z:=yy*(arctan(z)-arctan2arg(z*xx,w));
      goto 9
    end;

{A} if f='A' then
    begin z:=w*(xx+w); if z=0 then writeln(e,'Division by zero in A') else z:=1/z;
      goto 9
    end;

{Ep}if f='Ep' then begin if xx*yy<>0 then z:=xx*ln((w-yy)/(w+yy))/2; goto 9 end;

{Rm}if f='Rm' then
    begin z:=xx*xy;
      if (z=0) and (yy=0) then writeln(e,'0/0 in Rm')
      else if z<>0 then z:=arctan2arg(z,yy*w);
      goto 9
    end;

{Em}if f='Em' then begin if xx*yy<>0 then z:=yy*ln((w-xx)/(w+xx))/2; goto 9 end;

writeln('ERROR in evaluateElementaryFct: unknown function ',f);

9:  evaluateElementaryFct:=z;
end;

function evaluateFct(list: ^AT; var absvalue: real): real;
var c,z,a,d: real; l: ^AT; s: string; label 5,9;
begin z:=0; a:=0; l:=list;
   while l<>NIL do with l^ do
   begin
      s:=sx; if s<>'$' then
      begin
5:       writeln('ERROR in evaluateFct: ',s,' in ',fct,' not yet evaluated');
         goto 9
      end;
      s:=sy;  if s<>'$' then goto 5;
      s:=sxx; if s<>'$' then goto 5;
      s:=syy; if s<>'$' then goto 5;
      c:=evaluatePolynomial(sp,x,y,xx,yy);
      if c<>0 then
      begin d:=c*evaluateElementaryFct(fct,x-y,xx,yy);
	 z:=z+d; absvalue:=absvalue+abs(d)
      end;
9:    l:=post
   end;
   evaluateFct:=z;
end;

procedure menuinternal;
var choice: char; s: string; label 1,9;
begin
1: writeln; writeln('###################################');
   writeln('### MENU for internal variables ###');
   writeln('> l: list status of internal variables');
   writeln('     representation for output ...');
   if RationalRepresentation=1 then
   writeln('> f: use floating point numbers') else
   writeln('> r: use rational numbers');
   writeln('> T: change value of Tolerance');
   writeln('> q: change value of TolRational');
   writeln('> a: change value of MaxDenominator');
   writeln('     files ...');
   writeln('> t: change name of external TEX file');
   writeln('> x: change name of external storage file');
   writeln('> p: change name of protocol file');
   writeln('     for further choices see >l');
   writeln('### return to main menu with " "');
   write('  -> selected character = '); readln(choice);
   writelnP; writev(s,'>>> internal menu item "',choice,'"'); writeP(s); writelnP;
   case choice of
'r':	RationalRepresentation:=1;
'f':	RationalRepresentation:=0;
't':	begin write('new name for external TEX file = '); readln(fileTex);
              if fileTex='' then fileTex:='coulomb.tex';
              writev(s,'Redefinition...  fileTex := ',fileTex); writeP(s); writelnP
        end;
'p':	begin close(protocol); write('new name for protocol file = ');
              readln(fileLog); if fileLog='' then fileLog:='coulomb.log';
              rewrite(protocol,fileLog)
	end;
'x':	begin write('new name for external storage file = '); readln(fileExt);
              if fileExt='' then fileExt:='coulomb.ext';
              writev(s,'Redefinition...  fileExt := ',fileExt); writeP(s); writelnP
        end;
'T':	begin write('-> Tolerance = '); readln(Tolerance);
              if Tolerance<0 then Tolerance:=100*MachineEps;
              writev(s,'Redefinition...  Tolerance := ',Tolerance); writeP(s); writelnP
        end;
'q':	begin write('-> TolRational = '); readln(TolRational);
              if TolRational<0 then TolRational:=1000*MachineEps;
              writev(s,'Redefinition...  TolRational := ',TolRational); writeP(s); writelnP
        end;
'a':	begin write('-> MaxDenominator = '); readln(MaxDenominator);
              if MaxDenominator<=10 then MaxDenominator:=10000;
              TreatAsRationalNumber:=1/(MaxDenominator+1);
              writev(s,'Redefinition...  MaxDenominator := ',MaxDenominator:1); writeP(s); writelnP
        end;
'm':	m12:=1-m12;
'c':	cem:=1-cem;
'C':	cXY:=1-cXY;
'g':	g13:=1-g13;
'G':	g34:=1-g34;
'l':	begin writeln('--- internal status:');
          write('    external storage file is named  "'); write(fileExt); writeln('"');
          write('    external TEX file is named      "'); write(fileTex); writeln('"');
          write('    file for protocol of actions is "'); write(fileLog); writeln('"');
	  if xintegration>0 then writeln('    x-integration already performed in phase ',phase:1);
	  if yintegration>0 then writeln('    y-integration already performed in phase ',phase:1);
          writeln('    x with |x|<=Tolerance is treated as zero for polynomial coefficients x, where Tolerance = ',Tolerance);
	  write('    number representation in output by ');
	     if RationalRepresentation=1 then write('rational')
	     else write('floating point'); writeln(' numbers');
          writeln('    rational numbers accepted if |x-r|<=TolRational, where TolRational = ',TolRational);
          writeln('    rational are m/n, where n is bounded by MaxDenominator = ',MaxDenominator:1);
	  writeln('    --- In the treatment of the expressions the following');
	  writeln('        replacements may be allowed or not. The variables');
	  writeln('        x-y, X, Y used below indicate evaluated ones. The');
	  writeln('        characters in brackets allow to revert the state ...');
	  write('    <g> In F,G:   F(x-y;X,.) <-> F(X;x-y,.) is ');
          if g13=0 then write('not '); writeln('allowed');
	  write('    <G> In F,G,L: F(.;X,Y)   <-> F(.;Y,X)   is ');
          if g34=0 then write('not '); writeln('allowed');
	  write('    <m> M(x-y;X,.) -> L(X;x-y,.) is ');
          if m12=0 then write('not '); writeln('allowed');
	  write('    <c> In Cm,Em,Rm: Cm(.;X,Y) -> Cp,Ep,Rp(.;Y,X) is ');
          if cem=0 then write('not '); writeln('allowed');
          write('    <X> ... same operation when X,Y are not both evaluated is ');
          if cXY=0 then write('not '); writeln('allowed');
	end;
otherwise      goto 9;
   end;
   goto 1;
9: writeln('###################################')
end; {menuinternal}

procedure menuout(var fl,fs: ^AT);
var choice: char; s: string; v,w: real; label 1,2,9;
begin
1: writeln; writeln('#####################################################################');
   writeln('### MENU for output and storage of the actual function expression ###');
   writeln('     screen ...');
   writeln('> l: list the actual value of the expression');
   if fl<>NIL then begin
	if (fl^.sx='$') and (fl^.sy='$') and (fl^.sxx='$') and (fl^.syy='$') then
	writeln('> n: compute the numerical value of the expression');
	writeln('     internal ...');
	writeln('> s: write into internal storage');
	writeln('     external ...');
	if (fl^.sx='$') and (fl^.sy='$') and (fl^.sxx='$') and (fl^.syy='$') then
	begin write('> t: write into TEX file "'); write(fileTex); writeln('"') end;
	write('> x: write into external storage "'); write(fileExt); writeln('"')
   end;
   writeln('return to main menu with " "');
   write('### -> selected character = '); readln(choice);
   writelnP; writev(s,'>>> output menu item "',choice,'"'); writeP(s); writelnP;
   case choice of
'l':	begin WriteFctList('### output of the function array ...',fl);
2:            write('### give RETURN to return to the menu'); readln
	end;
'n':    begin v:=evaluateFct(fl,w);
              writev(s,'value = ',v:25:22,' *** sum of moduli of terms = ',w);
              write_P(s); writeln_P; goto 2
	end;
's':	begin fs:=DeleteAList(fs); CopyAList(fs,fl) end;
't':	begin rewrite(TX,fileTex); WriteTex(fl); close(TX) end;
'x':	begin rewrite(ES,fileExt); WriteExFct(fl); close(ES) end;
otherwise      goto 9;
   end;
   goto 1;
9: writeln('#####################################################################')
end; {menuout}

procedure menufct(var fl,fs: ^AT);
var choice: char; p: sparsePolynomial; s: string; label 1,2,9;
begin
   case phase of
   1:		begin msxx:='Y'; msyy:='Z' end;
   2:		begin msxx:='X'; msyy:='Z' end;
   otherwise	begin msxx:='X'; msyy:='Y' end;
   end;
1: writeln; writeln('########################################################');
   writeln('### MENU for defining the actual function expression ###');
   writeln('     standard start ...');
   write('> f: define expression by F');
   if phase<>1 then write(' and return to phase:=1'); writeln;
   writeln('     special ...');
   if fl<>NIL then writeln('> 0: set expression to zero');
   writeln('> a: add a term "polynomial*function"');
   if fl<>NIL then writeln('  m: multiply by a polynomial (subject of input)');
   writeln('     from storage ...');
   if fs<>NIL then begin
      writeln('> s: take from internal storage (expression overwritten by storage content)');
      writeln('> +: add from internal storage (expression := expression + storage)')
   end;
   write('> x: read from external storage "'); write(fileExt); writeln('"');
   writeln('     other actions ...');
   writeln('> l: list the actual value of the expression');
   writeln('> b: stabilisation and further simplifications');
   writeln('> o: omit zero terms, combine suitable terms');
   writeln('> p: change phase number');
   writeln('### return to main menu with " "');
   write('  -> selected character = '); readln(choice);
   writelnP; writev(s,'>>> definition menu item "',choice,'"'); writeP(s); writelnP;
   case choice of
's':	begin fl:=DeleteAList(fl); CopyAList(fl,fs) end;
'+':	fl:=AddFct(fl,fs);
'a':	DefineFct(fl,msxx,msyy);
'f':    begin fl:=DeleteAList(fl); fl:=NewAT('F','x','y',msxx,msyy,0,0,0,0);
              DefineCoefficient1_real(fl^.sp,0,0,1); phase:=1;
2:            xintegration:=0; yintegration:=0
        end;
'm':	begin StartP(p,1); Polynomial1_input(p);
              multiplyATbyPolynomial(fl,p); DeletePolynomial(p)
	end;
'l':	begin WriteFctList('### output of the function array ...',fl);
              write('### give RETURN to return to the menu'); readln
	end;
'x':	begin fl:=DeleteAList(fl); reset(ES,fileExt); fl:=ReadExFct; close(ES) end;
'p':	begin write('-> Phase = '); readln(phase); goto 2 end;
'0':    fl:=DeleteAList(fl);
'o':	fl:=reorganiseAT(fl);
'b':	stabilise(fl);
otherwise      goto 9;
   end;
   if choice<>'f' then goto 1;
9: writeln('########################################################')
end; {menufct}

procedure menuint(var fl: ^AT);
var choice: char; s: string; label 1,9;
begin
1: writeln; writeln('###########################');
   writeln('### MENU for integraton ###');
   writeln('     integration ...');
   if (fl^.sx<>'$') and (xintegration=0) then
   begin
   if (fl^.sy<>'$') and (yintegration=0) then
	writeln('> 2: xy- integration including stabilisation');
	writeln('> x: x - integration followed by stabilisation');
	writeln('> X: .  .  .  .  .   without stabilisation')
   end;
   if (fl^.sy<>'$') and (yintegration=0) then
   begin
	writeln('> y: y - integration followed by stabilisation');
	writeln('> Y: .  .  .  .  .   without stabilisation')
   end;
   writeln('     other actions ...');
   writeln('> l: list actual expression');
   writeln('> b: stabilisation and further simplifications');
   writeln('> o: omit zero terms, combine suitable terms');
   writeln('### return to main menu with " "');
   write('  -> selected character = '); readln(choice);
   writelnP; writev(s,'>>> integration menu item "',choice,'"'); writeP(s); writelnP;
   case choice of
'x','X': begin fl:=xIntegral(fl); if choice='x' then stabilise(fl) end;
'y','Y': begin fl:=yIntegral(fl); if choice='y' then stabilise(fl) end;
'2':	 begin fl:=xIntegral(fl); stabilise(fl); fl:=yIntegral(fl); stabilise(fl) end;
'l':	 begin WriteFctList('### output of the function array ...',fl);
	    write('### give RETURN to return to the menu'); readln
	 end;
'o':	fl:=reorganiseAT(fl);
'b':	stabilise(fl);
otherwise      goto 9;
   end;
   if xintegration*xintegration=0 then goto 1;
9: writeln('###########################')
end; {menuint}

procedure menueva(var fl: ^AT);
var choice: char; s: string; label 1,9;
begin
1: writeln; writeln('###########################');
   writeln('### MENU for evaluation ###');
{  writeln('< 132: .  .  . x by a linear polynomial a*y+c >');
   writeln('< 133: .  .  . y by a linear polynomial a*x+c >'); }
   if fl^.sx<>'$'  then
   begin
	writeln('     evaluation of x ...');
   	writeln('> x: replace x by a real value (value is subject of input)');
   	writeln('> i: .  .  . p(x,y) by p(x2,y)-p(x1,y) for real values x1,x2');
	writeln('             ... same without stabilisation ...');
   	writeln('> a: replace x by a real value');
   	writeln('> A: .  .  . p(x,y) by p(x2,y)-p(x1,y) for real values x1,x2')
   end;
   if fl^.sy<>'$'  then
   begin writeln('     evaluation of y ...');
   	writeln('> y: replace y by a real value (value is subject of input)');
   	writeln('> j: .  .  . p(x,y) by p(x,y2)-p(x,y1) for real values y1,y2');
	writeln('             ... same without stabilisation ...');
   	writeln('> b: replace y by a real value');
   	writeln('> B: .  .  . p(x,y) by p(x,y2)-p(x,y1) for real values y1,y2')
   end;
   if (fl^.sx<>'$') and (fl^.sy<>'$') then
   begin writeln('     evaluation of x and y ...');
   	writeln('> c: .  .  . combine actions "x" and "y"');
   	writeln('> 1: .  .  . x=0,1 and y=0,1');
   	writeln('> 2: .  .  . x=0,1 and y=1,2');
   	writeln('> 3: .  .  . x=1,2 and y=0,1');
   	writeln('> 4: .  .  . x=1,2 and y=1,2')
   end;
   if (fl^.sxx<>'$') or (fl^.syy<>'$') then writeln('     evaluation of parameters X or Y ...');
   if fl^.sxx<>'$' then writeln('> X: replace ',fl^.sxx,' by a real value');
   if fl^.syy<>'$' then writeln('> Y: replace ',fl^.syy,' by a real value');
   writeln('     other actions ...');
   writeln('> l: list status of internal variables');
   writeln('### return to main menu with " "');
   write('  -> selected character = '); readln(choice);
   writelnP; writev(s,'>>> evaluation menu item "',choice,'"'); writeP(s); writelnP;
   case choice of
'l':	begin WriteFctList('### output of the function array ...',fl);
	    write('### give RETURN to return to the menu'); readln
	end;
'g':	begin replaceATlinear(fl,1); stabilise(fl) end;
'h':	begin replaceATlinear(fl,2); stabilise(fl) end;
'x':	begin replaceATconstant(fl,1,0,0); stabilise(fl) end;
'i':	begin evaluateATatLinearBounds(fl,1,0,0,0); stabilise(fl) end;
'a':	replaceATconstant(fl,1,0,0);
'A':	evaluateATatLinearBounds(fl,1,0,0,0);
'y':	begin replaceATconstant(fl,2,0,0); stabilise(fl) end;
'j':	begin evaluateATatLinearBounds(fl,2,0,0,0); stabilise(fl) end;
'b':	replaceATconstant(fl,2,0,0);
'B':	evaluateATatLinearBounds(fl,2,0,0,0);
'X':    replaceATconstant(fl,3,0,0);
'Y':    replaceATconstant(fl,4,0,0);
'c':    begin evaluateBothBounds(fl,0,0,0,0,0); stabilise(fl) end;
'1':    begin evaluateBothBounds(fl,1,0,1,0,1); stabilise(fl) end;
'2':    begin evaluateBothBounds(fl,1,0,1,1,2); stabilise(fl) end;
'3':    begin evaluateBothBounds(fl,1,1,2,0,1); stabilise(fl) end;
'4':    begin evaluateBothBounds(fl,1,1,2,1,2); stabilise(fl) end;
otherwise      goto 9;
   end;
   case choice of 'l','g','h','x','i','a','A','y','j','b','B','X','Y': goto 1;
	otherwise      goto 9;
 end;
9: writeln('###########################')
end; {menueva}

procedure menu;
var fl,fs: ^AT; s: string; choice: char; label 1,9;
begin fl:=NIL; fs:=NIL;
1: writeln; writeln('#######################');
   writeln('###### PHASE = ',phase:1,' ######');
   writeln('###### MAIN MENU ######');
   writeln('1:   define actual function expression');
   if fl<>NIL then
   begin if ((fl^.sx<>'$') and (xintegration=0)) or ((fl^.sy<>'$') and (yintegration=0)) then
	writeln('2:   integration');
	if (fl^.sx<>'$') or (fl^.sy<>'$') or (fl^.sxx<>'$') or (fl^.syy<>'$') then
	writeln('3:   evaluation');
	if (fl^.sx='$') and (fl^.sy='$') then if phase<3 then
	writeln('4:   interpretation of the function w.r.t. x',(phase+1):1,', y',(phase+1):1)
   end;
   writeln('5:   output and storage of actual function expression');
   writeln('6:   internal status');
   writeln('7,l: list the actual function expression');
   writeln('end of program with " "');
   write('  -> selected character = '); readln(choice);
   writelnP; writev(s,'>>>>>> main menu item "',choice,'"'); writeP(s); writelnP;
   case choice of
'1':	menufct(fl,fs);
'2':	if fl<>NIL then menuint(fl);
'3':	if fl<>NIL then menueva(fl);
'4':	begin if phase=1 then begin interpret12(fl); phase:=2 end else
		if phase=2 then begin interpret23(fl); phase:=3 end;
		xintegration:=0; yintegration:=0
	end;
'5':	menuout(fl,fs);
'6':	menuinternal;
'7','l':begin WriteFctList('### output of the function array ...',fl);
              write('### give RETURN to return to the menu'); readln
	end;
   otherwise goto 9
   end; {case}
   goto 1;
9: writeln('######## END ##########'); s:='END OF PROGRAM'; writeP(s); writelnP
end;

begin fileExt:='coulomb.ext'; fileTex:='coulomb.tex'; fileLog:='coulomb.log';
   rewrite(protocol,fileLog);
   RationalRepresentation:=1; phase:=1; xintegration:=0; yintegration:=0;
   m12:=1; g13:=1; g34:=1; cem:=1; cxy:=0; MaxDenominator:=10000;
   Tolerance:=100*MachineEps; TolRational:=1000*MachineEps;
   TreatAsRationalNumber:=1/(MaxDenominator+1);

   menu;
   
   close(protocol);
end.