File : funs2.adb

with Ints.IO, Reps.IO, Messages;
use Ints.IO, Reps.IO, Messages;

package body Funs2 is

  function QuasiDeg2(D,I: Integer; C: Components) return Integer is
  begin
    for J in reverse 0 .. D loop
      if C(I,J) /= ScalZero then return J; end if;
    end loop;
    return -1;
  end QuasiDeg2;

  procedure MultAdd(S: in Scalar; P: in Affine; F: in out Taylor) is
  begin
    MultAdd(S,P.A,F.C(1,0));
    MultAdd(S,P.B,F.C(0,1));
    MultAdd(S,P.C,F.C(0,0));
  end MultAdd;

  procedure Mult(P: in Affine; F: in out Taylor; Tmp: out Taylor) is
--- Tmp is scratch space
  begin
    Tmp := F;
    F := TayZero;
    MultAdd(P,Tmp,F);
  end Mult;

  procedure Mult(P: in Affine; F: in out Taylor) is
    Tmp: constant Taylor := F;
  begin
    F := TayZero;
    MultAdd(P,Tmp,F);
  end Mult;

  procedure MultAdd(P: in Affine; F: in Taylor; G: in out Taylor) is
    DF: constant Integer := QuasiDeg(F);
    D1: constant Integer := IMin(DF+1,PDeg);
    E: Scalar;
  begin
    if DF<0 then return; end if;
    MultAdd(P.C,F,G);
    if P.A /= ScalZero then
      for I in 1 .. D1 loop
        for J in 0 .. D1-I loop
          MultAdd(P.A,F.C(I-1,J),G.C(I,J));
        end loop;
      end loop;
    elsif P.B=ScalZero then return; end if;
    if P.B /= ScalZero then
      for J in 1 .. D1 loop
        for I in 0 .. D1-J loop
          MultAdd(P.B,F.C(I,J-1),G.C(I,J));
        end loop;
      end loop;
    end if;
    G.R := Min(F.R,G.R);
    if Trunc or else DF<PDeg then return; end if;
    E := P.A*BallAt0(F.R(1))+P.B*BallAt0(F.R(2));
    for I in Power loop
      MultAdd(E,F.C(I,PDeg-I),G.C(I,PDeg-I));
    end loop;
  end MultAdd;

  procedure Comp(U,V: in Affine; F: in Taylor; G: out Taylor) is
    D: constant Integer := QuasiDeg(F);
  begin
    G := TayZero;
    if D<0 then return; end if;
    declare
      D1: constant Integer := QuasiDeg1(D,F);
      R:  constant RadPair := Dom(U,V,F.R);
      D2: Integer;
      H,Tmp: Taylor2;
    begin
      for I in reverse 0 .. D1 loop
        SetZero(R,H);
        D2 := QuasiDeg2(D,I,F.C);
        if D2 >=0 then
          H.C(0,0) := F.C(I,D2);
          for J in reverse 0 .. D2-1 loop
            Mult(V,H,Tmp);
            Add(F.C(I,J),H.C(0,0));
          end loop;
        end if;
        Mult(U,G,Tmp);
        Add(H,G);
      end loop;
    end;
  end Comp;

  procedure FillZero(D: in Natural; G: out TaylorMat) is
  begin
    for N in 0 .. D loop
      for I in 0 .. N loop
        G(I,N-I) := TayZero;
      end loop;
    end loop;
  end FillZero;

-------------------------------

  procedure CheckArgs(F: in Fun; G: in out Fun; OK: out Boolean) is
  begin
    if F.A=G.A or else IsZero(F) then
      OK := True;
    else
      G.A := F.A;
      OK := IsZero(G);
    end if;
  end CheckArgs;

------------------------------

  function Rho1(F: Fun) return Radius is
  begin
    return RMin(F.P.R(1),F.Q.R(1));
  end Rho1;

  function Rho2(F: Fun) return Radius is
  begin
    return RMin(F.P.R(2),F.Q.R(2));
  end Rho2;

  function Rho(F: Fun) return RadPair is
  begin
    return Min(F.P.R,F.Q.R);
  end Rho;

  procedure SetRho(R: in RadPair; F: in out Fun) is
  begin
    SetRho(R,F.P);
    SetRho(R,F.Q);
  end SetRho;

  function Dom(F: Fun) return Domain is
    D: Domain;
  begin
    D.A := F.A;
    D.R := Rho(F);
    return D;
  end Dom;

  function IsZero(F: Fun) return Boolean is
  begin
    return IsZero(F.P) and then IsZero(F.Q);
  end IsZero;

  procedure FillZero(F: in out Fun) is
--- set F.A elsewhere
  begin
    F.E := True;
    F.P := TayZero;
    F.Q := TayZero;
  end FillZero;

  procedure SetZero(A: in Args; F: out Fun) is
  begin
    F.A := A;
    FillZero(F);
  end SetZero;

  procedure SetZero(Template: in Fun; F: out Fun) is
  begin
    F.A := Template.A;
    FillZero(F);
  end SetZero;

  procedure ToBalls(F: in out Fun) is
  begin
    ToBalls(F.P);
    ToBalls(F.Q);
  end ToBalls;

  procedure Symm(F: in out Fun) is
--- extract symmetric part
    R: constant RadPair := Rho(F);
  begin
    if not (F.E or else IsZero(F.Q)) then
      SetRho(R,F);
      MultAdd(TNorm(F.A,R),ModNumbers(F.Q),F.P);
    end if;
    F.Q := TayZero;
    F.E := True;
  end Symm;

  function Symm(F: Fun) return Fun is
    G: Fun := F;
  begin
    Symm(G);
    return G;
  end Symm;

  procedure AddConst(C: in Scalar; F: in out Fun) is
  begin
    AddConst(C,F.P);
  end AddConst;

  procedure AddBall(C: in Scalar; F: in out Fun) is
  begin
    if not Trunc then
      AddConst(BallAt0(C),F.P);
      F.E := F.E and then IsZero(F.Q);
    end if;
  end AddBall;

  procedure AddT1(C: in Scalar; F: in out Fun) is
  begin
    AddConst(C,F.Q);
  end AddT1;

  procedure AddS1(C: in Scalar; F: in out Fun) is
  begin
    MultAdd(C,CoordS1(F.A),F.P);
  end AddS1;

  procedure AddX(C: in Scalar; F: in out Fun) is
  begin
    AddT1(Half*C,F);
    AddS1(Half*C,F);
  end AddX;

  procedure AddY(C: in Scalar; F: in out Fun) is
  begin
    AddT1(Half*C,F);
    AddS1((-Half)*C,F);
  end AddY;

  function IsNumbers(F: Fun) return Boolean is
  begin
    return IsNumbers(F.P) and then IsNumbers(F.Q);
  end IsNumbers;

  function ModNumbers(F: Fun) return Fun is
    G: Fun := F;
  begin
    ModNumbers(G.P);
    ModNumbers(G.Q);
    return G;
  end ModNumbers;

  function Norm(R: RadPair; F: Fun) return Scalar is
  begin
    return Norm(R,F.P)+TNorm(F.A,R)*Norm(R,F.Q);
  end Norm;

  function Norm(F: Fun) return Scalar is
  begin
    return Norm(Rho(F),F);
  end Norm;

  procedure Neg(F: in out Fun) is
  begin
    Neg(F.P);
    Neg(F.Q);
  end Neg;

  procedure Neg(F: in Fun; G: out Fun) is
  begin
    G := F;
    Neg(G);
  end Neg;

  function "-"(F: Fun) return Fun is
    G: Fun := F;
  begin
    Neg(G);
    return G;
  end "-";

  procedure Add(F: in Fun; G: in out Fun) is
    OK: Boolean;
  begin
    CheckArgs(F,G,OK);
    if not OK then Error("Funs2.Add error"); end if;
    if G.E then G.E := F.E; end if;
    Add(F.P,G.P);
    Add(F.Q,G.Q);
  end Add;

  procedure Sum(F,G: in Fun; H: out Fun) is
  begin
    H := F;
    Add(G,H);
  end Sum;

  function "+"(F,G: Fun) return Fun is
    H: Fun := F;
  begin
    Add(G,H);
    return H;
  end "+";

  procedure Sub(F: in Fun; G: in out Fun) is
    OK: Boolean;
  begin
    CheckArgs(F,G,OK);
    if not OK then Error("Funs2.Sub error"); end if;
    if G.E then G.E := F.E; end if;
    Sub(F.P,G.P);
    Sub(F.Q,G.Q);
  end Sub;

  procedure Diff(F,G: in Fun; H: out Fun) is
  begin
    H := F;
    Sub(G,H);
  end Diff;

  function "-"(F,G: Fun) return Fun is
    H: Fun := F;
  begin
    Sub(G,H);
    return H;
  end "-";

  procedure Mult(S: in Scalar; F: in out Fun) is
  begin
    Mult(S,F.P);
    Mult(S,F.Q);
    if S=ScalZero then F.E := True; end if;
  end Mult;

  procedure Prod(S: in Scalar; F: in Fun; G: out Fun) is
  begin
    if S=ScalZero then
      SetZero(F.A,G);
    else
      G := F;
      Mult(S,G);
    end if;
  end Prod;

  function "*"(S: Scalar; F: Fun) return Fun is
    G: Fun;
  begin
    Prod(S,F,G);
    return G;
  end "*";

  procedure MultAdd(S: in Scalar; F: in Fun; G: in out Fun) is
    OK: Boolean;
  begin
    if S /= ScalZero then
      CheckArgs(F,G,OK);
      if not OK then Error("Funs2.MultAdd error"); end if;
      MultAdd(S,F.P,G.P);
      MultAdd(S,F.Q,G.Q);
    end if;
  end MultAdd;

  procedure Mult(F: in Fun; G: in out Fun; Tmp: out Fun) is
--- Tmp is scratch space
  begin
    Tmp := G;
    FillZero(G);
    MultAdd(F,Tmp,G);
  end Mult;

  procedure Mult(F: in Fun; G: in out Fun) is
    Tmp: Fun;
  begin
    Mult(F,G,Tmp);
  end Mult;

  procedure Prod(F,G: in Fun; H: out Fun) is
  begin
    SetZero(F.A,H);
    MultAdd(F,G,H);
  end Prod;

  function "*"(F,G: Fun) return Fun is
    H: Fun;
  begin
    SetZero(F.A,H);
    MultAdd(F,G,H);
    return H;
  end "*";

  function Sqr(F: Fun) return Fun is
    H: Fun;
  begin
    SetZero(F.A,H);
    MultAdd(F,F,H);
    return H;
  end Sqr;

  procedure MultAdd(F,G: in Fun; H: in out Fun) is
    OK: Boolean;
  begin
    MultAdd(F.P,G.P,H.P);
    MultAdd(F.Q,G.P,H.Q);
    MultAdd(F.P,G.Q,H.Q);
    MultAdd(CoordT2(F.A),F.Q*G.Q,H.P);
    if H.E then H.E := (F.E and G.E); end if;
    CheckArgs(F,H,OK);
    if not OK then Error("Funs2.MultAdd error F"); end if;
    CheckArgs(G,H,OK);
    if not OK then Error("Funs2.MultAdd error G"); end if;
  end MultAdd;

  procedure ChangeVar(A: in Args; F: in out Fun) is
    G: constant Fun := F;
    B: Mat3;
  begin
    if F.A=A then return; end if;
    if IsZero(F) then F.A := A; return; end if;
    if not F.E then
      Error("Funs2.ChangeVar error: not implemented");
    end if;
    Coord(F.A,A,B);
    F.A := A;
    Comp(B.U,B.V,G.P,F.P);
    Comp(B.U,B.V,G.Q,F.Q);
  end ChangeVar;

  procedure Scale(S: in Scalar; A: in Args; F: in out Fun) is
    G: constant Fun := F;
    B: Mat3;
  begin
    if not F.E then
      Error("Funs2.Scale error: not implemented");
    end if;
    Scale(S,F.A,A,B);
    F.A := A;
    Comp(B.U,B.V,G.P,F.P);
    Comp(B.U,B.V,G.Q,F.Q);
  end Scale;

  procedure Scale(S: in Scalar; F: in out Fun) is
    A: Args;
    B: Mat3;
  begin
    Convert(F.A,B);
    Scale(S,B);
    Approx(B,A);
    Scale(S,A,F);
  end Scale;

  procedure DerTS(FA: in Args; DP: in TaylorMat; F1,F2: in out Fun) is
--- not using/setting Args or setting/checking E of F1,F2
  begin
--- P1 := Q+2*(T^2)*[DQ(1,0)+A*DQ(0,1)]]
    F1.P := TayZero;
--- Q1 := 2*[DP(1,0)+A*DP(0,1)]
    F1.Q := DP(1,0);
    MultAdd(Scal(FA.A),DP(0,1),F1.Q);
    Mult(ScalTwo,F1.Q);
--- P2 := B*DP(1,0)+DP(0,1)
    F2.P := DP(0,1);
    MultAdd(Scal(FA.B),DP(1,0),F2.P);
--- Q2 := B*DQ(1,0)+DQ(0,1)
    F2.Q := TayZero;
  end DerTS;

  procedure DerTS(FA: in Args; DP: in TaylorMat; F11,F12,F22: out Fun) is
--- not using/setting Args or setting/checking E of F11,F12,F22
  begin
--- P11 := 2*[DP(1,0)+A*DP(0,1)]+4*(T^2)[DP(2,0)+2*A*DP(1,1)+A*A*DP(0,2)]
    Prod(ScalTwo,DP(1,1),F11.P);
    MultAdd(Scal(FA.A),DP(0,2),F11.P);
    Mult(Scal(FA.A),F11.P);
    Add(DP(2,0),F11.P);
    Mult(CoordT2(FA),F11.P);
    Mult(ScalFour,F11.P);
    MultAdd(ScalTwo,DP(1,0),F11.P);
    MultAdd(FA.A*ScalTwo,DP(0,1),F11.P);
--- Q11 := 6*[DQ(1,0)+A*DQ(0,1)]+4*(T^2)[DQ(2,0)+2*A*DQ(1,1)+A*A*DQ(0,2)]
    F11.Q := TayZero;
--- P12 := [B*DQ(1,0)+DQ(0,1)]+2*(T^2)*[B*DQ(2,0)+(1+A*B)*DQ(1,1)+A*DQ(0,2)]
    F12.P := TayZero;
--- Q12 := 2*[B*DP(2,0)+(1+A*B)*DP(1,1)+A*DP(0,2)]
    Prod(ScalOne+FA.A*Scal(FA.B),DP(1,1),F12.Q);
    MultAdd(Scal(FA.A),DP(0,2),F12.Q);
    MultAdd(Scal(FA.B),DP(2,0),F12.Q);
    Mult(ScalTwo,F12.Q);
--- P22 := B*B*DP(2,0)+2*B*DP(1,1)+DP(0,2)
    F22.P := DP(0,2);
    MultAdd(FA.B*ScalTwo,DP(1,1),F22.P);
    MultAdd(FA.B*Scal(FA.B),DP(2,0),F22.P);
--- Q22 := B*B*DQ(2,0)+2*B*DQ(1,1)+DQ(0,2)
    F22.Q := TayZero;
  end DerTS;

  DP: TaylorMat(0..2,0..2);           -- used by derivative procedures below

  procedure MDer(R: in RadPair; F: in Fun; G: out Fun) is
--- generator of scaling
--- uses global DP, as storage only
    H: Taylor;
  begin
    if not F.E then
      Error("Funs2.MDer error: not implemented");
    end if;
    G := F;
    Ders(1,R,F.P,DP);
--- P2 := S*[B*DP(1,0)+DP(0,1)]
    G.P := DP(0,1);
    MultAdd(Scal(F.A.B),DP(1,0),G.P);
    Mult(CoordS1(F.A),G.P);
--- P1 := 2*(T^2)*[DP(1,0)+A*DP(0,1)]
    H := DP(1,0);
    MultAdd(Scal(F.A.A),DP(0,1),H);
    Mult(CoordT2(F.A),H);
    MultAdd(ScalTwo,H,G.P);
  end MDer;

  procedure DerTS(R: in RadPair; F: in Fun; F1,F2: out Fun) is
--- uses global DP, as storage only
  begin
    if not F.E then
      Error("Funs2.DerTS error: not implemented");
    end if;
    Ders(1,R,F.P,DP);
    DerTS(F.A,DP,F1,F2);
    F1.A := F.A;
    F1.E := True;
    F2.A := F.A;
    F2.E := True;
  end DerTS;

  procedure Der(R: in RadPair; F: in Fun; F1,F2: out Fun) is
    FT,FS: Fun;
  begin
    DerTS(R,F,FT,FS);
    Sum(FT,FS,F1);
    Diff(FT,FS,F2);
  end Der;

  procedure DerTS(R: in RadPair; F: in Fun; F1,F2,F11,F12,F22: out Fun) is
--- uses global DP, as storage only
  begin
    if not F.E then
      Error("Funs2.DerTS error: not implemented");
    end if;
    Ders(2,R,F.P,DP);
    DerTS(F.A,DP,F1,F2);
    DerTS(F.A,DP,F11,F12,F22);
    F1.A := F.A;
    F1.E := True;
    F2.A := F.A;
    F2.E := True;
    F11.A := F.A;
    F11.E := True;
    F12.A := F.A;
    F12.E := True;
    F22.A := F.A;
    F22.E := True;
  end DerTS;

  procedure Der(R: in RadPair; F: in Fun; F1,F2,F11,F12,F22: out Fun) is
    FT,FS,FTT,FTS,FSS: Fun;
  begin
    DerTS(R,F,FT,FS,FTT,FTS,FSS);
    Sum(FT,FS,F1);
    Diff(FT,FS,F2);
    Sum(FTT,FSS,F11);
    F22 := F11;
    MultAdd(ScalTwo,FTS,F11);
    MultAdd(-ScalTwo,FTS,F22);
    Diff(FTT,FSS,F12);
  end Der;

  function ValTS(T,S: Scalar; F: Fun) return Scalar is
    U,V: Scalar;
  begin
    ValTS(T,S,F.A,U,V);
    return Val(U,V,F.P)+T*Val(U,V,F.Q);
  end ValTS;

  function Val(X,Y: Scalar; F: Fun) return Scalar is
  begin
    return ValTS(X+Y,X-Y,F);
  end Val;

  DF1,DF2,DF11,DF12,DF22: Fun;   --- used by ValDerTS

  procedure ValDerTS(T,S: in Scalar; F: in Fun; D1,D2: out Scalar) is
--- uses globals DF1,DF2
    R: RadPair;
    U,V: Scalar;
  begin
    ValTS(T,S,F.A,U,V);
    R := (MaxNorm(U),MaxNorm(V));
    DerTS(R,F,DF1,DF2);
    D1 := ValTS(T,S,DF1);
    D2 := ValTS(T,S,DF2);
  end ValDerTS;

  procedure ValDer(X,Y: in Scalar; F: in Fun; D1,D2: out Scalar) is
    DT,DS: Scalar;
  begin
    ValDerTS(X+Y,X-Y,F,DT,DS);
    D1 := DT+DS;
    D2 := DT-DS;
  end ValDer;

  procedure ValDerTS(T,S: in Scalar; F: in Fun; D1,D2,D11,D12,D22: out Scalar) is
--- uses globals DF1,DF2,DF11,DF12,DF22
    R: RadPair;
    U,V: Scalar;
  begin
    ValTS(T,S,F.A,U,V);
    R := (MaxNorm(U),MaxNorm(V));
    DerTS(R,F,DF1,DF2,DF11,DF12,DF22);
    D1 := ValTS(T,S,DF1);
    D2 := ValTS(T,S,DF2);
    D11 := ValTS(T,S,DF11);
    D12 := ValTS(T,S,DF12);
    D22 := ValTS(T,S,DF22);
  end ValDerTS;

  procedure ValDer(X,Y: in Scalar; F: in Fun; D1,D2,D11,D12,D22: out Scalar) is
    DT,DS,DTT,DTS,DSS: Scalar;
  begin
    ValDerTS(X+Y,X-Y,F,DT,DS,DTT,DTS,DSS);
    D1 := DT+DS;
    D2 := DT-DS;
    D11 := DTT+DSS+Two*DTS;
    D22 := DTT+DSS-Two*DTS;
    D12 := DTT-DSS;
  end ValDer;

  procedure SplitParity(F: in Fun; F0,F1: out Fun) is
  begin
    F0 := F;
    F1 := F;
    if F.E then        -- only applies when Trunc
      F0.Q := TayZero;
      F1.P := TayZero;
    else
      ModNumbers(F0.Q);
      ModNumbers(F1.P);
    end if;
  end SplitParity;

  procedure Reflect(F: in out Fun) is
  begin
    Neg(F.Q);
  end Reflect;

  function Reflect(F: Fun) return Fun is
    G: Fun := F;
  begin
    Reflect(G);
    return G;
  end Reflect;

  function "*"(L: Flop; F: Fun) return Fun is
--- L(1)*F+L(2)*Reflect(F)
    G: Fun := F;
  begin
    Reflect(G);
    Mult(L(2),G);
    MultAdd(L(1),F,G);
    return G;
  end "*";

  function "*"(L1,L2: Flop) return Flop is
    L: Flop;
  begin
    L(1) := L1(1)*L2(1)+L1(2)*Reflect(L2(2));
    L(2) := L1(1)*L2(2)+L1(2)*Reflect(L2(1));
    return L;
  end "*";

  function "-"(L1,L2: Flop) return Flop is
    L: Flop;
  begin
    Diff(L1(1),L2(1),L(1));
    Diff(L1(2),L2(2),L(2));
    return L;
  end "-";

  function Norm(R: RadPair; L: Flop) return Scalar is
  begin
    return Norm(R,L(1))+Norm(R,L(2));
  end Norm;

  function CompTSRadii(FT,FS: Fun; A: Args) return RadPair is
    FU,FV: Fun;
  begin
    CompTS(FT,FS,A,FU,FV);
    return (Sup(Norm(FU)),Sup(Norm(FV)));
  end CompTSRadii;

  function CompVWRadii(V: Fun; A: Args) return RadPair is
--- domain radii needed to compose with (V,W)
    V0,V1: Fun;
  begin
    SplitParity(ScalTwo*V,V0,V1);
    return CompTSRadii(V1,V0,A);
  end CompVWRadii;

  function CompRadii(FX,FY: Fun; A: Args) return RadPair is
    FU,FV: Fun;
  begin
    Comp(FX,FY,A,FU,FV);
    return (Sup(Norm(FU)),Sup(Norm(FV)));
  end CompRadii;

  function CompXRadii(FX: Fun; A: Args) return RadPair is
    FY: Fun;
  begin
    SetZero(FX.A,FY);
    AddY(ScalOne,FY);
    return CompRadii(FX,FY,A);
  end CompXRadii;

  function CompYRadii(FY: Fun; A: Args) return RadPair is
    FX: Fun;
  begin
    SetZero(FY.A,FX);
    AddX(ScalOne,FX);
    return CompRadii(FX,FY,A);
  end CompYRadii;

  procedure CompTS(FT,FS: in Fun; A: in Args; FU,FV: out Fun) is
    TT0: constant Scalar := A.T0*Scal(A.T0);
  begin
    Prod(FT,FT,FU);
    FV := FS;
    MultAdd(Scal(A.A),FU,FV);
    MultAdd(Scal(A.B),FS,FU);
    AddConst(-(TT0+A.B*Scal(A.S0)),FU.P);
    AddConst(-(A.A*TT0+Scal(A.S0)),FV.P);
  end CompTS;

  procedure Comp(FX,FY: in Fun; A: in Args; FU,FV: out Fun) is
    FT,FS: Fun;
  begin
    Sum(FX,FY,FT);
    Diff(FX,FY,FS);
    CompTS(FT,FS,A,FU,FV);
  end Comp;

  procedure Eval(V1,V2: in Fun; F: in Taylor; V3: out Fun) is
    D:  constant Integer := QuasiDeg(F);
    RV: constant RadPair := Min(Rho(V1),Rho(V2));
    RF: constant RadPair := (Sup(Norm(RV,V1)),Sup(Norm(RV,V2)));
    V4,Tmp: Fun;
  begin
    SetZero(V1,V3);
    if D<0 then return; end if;
    if not InDomain(RF,F) then
      Error("Taylors2.Eval domain error");
    end if;
    if V1.A /= V2.A then
      Error("Funs2.Eval error: arguments differ");
    end if;
    for I in reverse 0 .. QuasiDeg1(D,F) loop
      SetZero(V2,V4);
      for J in reverse 0 .. QuasiDeg2(D,I,F.C) loop
        Mult(V2,V4,Tmp);
        AddConst(F.C(I,J),V4);
      end loop;
      Mult(V1,V3,Tmp);
      Add(V4,V3);
    end loop;
    SetRho(RV,V3);
    if V3.E then V3.E := IsNumbers(F); end if;
  end Eval;

  procedure CompTS(FT,FS,G: in Fun; H: out Fun) is
    FU,FV,HQ: Fun;
  begin
    CompTS(FT,FS,G.A,FU,FV);
    Eval(FU,FV,G.P,H);
    Eval(FU,FV,G.Q,HQ);
    MultAdd(FT,HQ,H);
    if (not G.E) and then Norm(Rho(H),FT)>TNorm(G.A,Rho(G)) then
      Error("Funs2.CompTS error: FT too large",T3);
    end if;
  end CompTS;

  procedure Comp(FX,FY,G: in Fun; H: out Fun) is
    FT,FS: Fun;
  begin
    Sum(FX,FY,FT);
    Diff(FX,FY,FS);
    CompTS(FT,FS,G,H);
  end Comp;

  procedure CompX(FX,G: in Fun; H: out Fun) is
    FY: Fun;
  begin
    SetZero(FX,FY);
    AddY(ScalOne,FY);
    Comp(FX,FY,G,H);
  end CompX;

  function CompX(FX,G: Fun) return Fun is
    H: Fun;
  begin
    CompX(FX,G,H);
    return H;
  end CompX;

  procedure CompY(FY,G: in Fun; H: out Fun) is
    FX: Fun;
  begin
    SetZero(FY,FX);
    AddX(ScalOne,FX);
    Comp(FX,FY,G,H);
  end CompY;

  function CompY(FY,G: Fun) return Fun is
    H: Fun;
  begin
    CompY(FY,G,H);
    return H;
  end CompY;

  function CompVW(V,F: Fun) return Fun is
--- compose F with (V,W) where W(x,y)=-V(-y,-x)
    H,V0,V1: Fun;
  begin
    Prod(ScalTwo,V,H);
    SplitParity(H,V0,V1);
    CompTS(V1,V0,F,H);
    return H;
  end CompVW;

  function FF(Fx,Gy,V: Fun) return Fun is
--- V -> Gy(x,V)+Fx(V,W), W=-SV
  begin
    return CompY(V,Gy)+CompVW(V,Fx);
  end FF;

  function DFF(Fxx,Fxy,Gyy,V: Fun) return Flop is
    L: Flop;
  begin
    L(1) := CompY(V,Gyy)+CompVW(V,Fxx);
    L(2) := -CompVW(V,Fxy);
    return L;
  end DFF;

  procedure Show1(N: in String; H: in Fun) is
  begin
    Show1(N,H.A);
    New_Line;
    Show1(N & " type",Int(not H.E));
    New_Line;
    Show1(N & " even",H.P);
    New_Line;
    Show1(N & "  odd",H.Q);
    New_Line;
    New_Line;
  end Show1;

  procedure Show2(N: in String; H1,H2: in Fun) is
  begin
    Show2(N,H1.A,H2.A);
    New_Line;
    Show2(N & " types",Int(not H1.E),Int(not H2.E));
    New_Line;
    Show2(N & "  even",H1.P,H2.P);
    New_Line;
    Show2(N & "   odd",H1.Q,H2.Q);
    New_Line;
    New_Line;
  end Show2;

  procedure Get(F: in File_Type; H: out Fun; Decimal: in Boolean := False) is
  begin
    Get(F,H.A.T0,Decimal);
    Get(F,H.A.S0,Decimal);
    Get(F,H.A.B,Decimal);
    Get(F,H.A.A,Decimal);
    Get(F,H.P,Decimal);
    Get(F,H.Q,Decimal);
    H.E := IsNumbers(H) or else IsZero(H.Q);
  end Get;

  procedure Put(F: in File_Type; H: in Fun; Decimal: in Boolean := False) is
--- H.E is not saved
  begin
    Put(F,H.A.T0,Decimal);
    Put(F,H.A.S0,Decimal);
    Put(F,H.A.B,Decimal);
    Put(F,H.A.A,Decimal);
    Put(F,H.P,Decimal);
    Put(F,H.Q,Decimal);
    if H.E /= (IsNumbers(H) or else IsZero(H.Q)) then
      Error("Funs2.Put error: not implemented");
    end if;
  end Put;

  procedure Read(Name: in String; H: out Fun; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    Show("Reading " & Name,1);
    Open(F,In_File,Name);
    Get(F,H,Decimal);
    Close(F);
  end Read;

  procedure Write(Name: in String; H: in Fun; Decimal: in Boolean := False) is
    F: File_Type;
  begin
    Show("Writing " & Name,1);
    Create(F,Out_File,Name);
    Put(F,H,Decimal);
    Close(F);
  end Write;

------------------------------------ some implicit equations

  procedure RVal(G: in Fun; X,C: in Scalar; S: in out Scalar) is
--- solves G(x,s)=c
    R0: constant Scalar := -ScalTwo;  -- guess
    R1: constant Scalar :=  ScalTwo;  -- guess
    GX,GY: Scalar;

    function H(Y: Scalar) return Scalar is
    begin
      return Val(X,Y,G)-C;
    end H;

    function DH(Y: Scalar) return Scalar is
    begin
      ValDer(X,Y,G,GX,GY);
      return GY;
    end DH;

    function Find is new FindZero(F => H, DF => DH);

  begin
    S := Find(S,R0,R1,Eps1,Eps2);
  end RVal;

  procedure LVal(F: in Fun; Y,C: in Scalar; S: in out Scalar) is
--- solves F(s,y)=c
    R0: constant Scalar := -ScalOne;  -- guess
    R1: constant Scalar :=  ScalOne;  -- guess
    FX,FY: Scalar;

    function H(X: Scalar) return Scalar is
    begin
      return Val(X,Y,F)-C;
    end H;

    function DH(X: Scalar) return Scalar is
    begin
      ValDer(X,Y,F,FX,FY);
      return FX;
    end DH;

    function Find is new FindZero(F => H, DF => DH);

  begin
     S := Find(S,R0,R1,Eps1,Eps2);
  end LVal;

  procedure ZVal(F,G: in Fun; X,Y: Scalar; S: in out Scalar) is
--- solve F(x,s)+G(s,y)=0
    R0: constant Scalar := -ScalOne;  -- guess
    R1: constant Scalar :=  ScalOne;  -- guess
    FX,FY,GX,GY: Scalar;

    function H(Z: Scalar) return Scalar is
    begin
      return Val(X,Z,F)+Val(Z,Y,G);
    end H;

    function DH(Z: Scalar) return Scalar is
    begin
      ValDer(X,Z,F,FX,FY);
      ValDer(Z,Y,G,GX,GY);
      return FY+GX;
    end DH;

    function Find is new FindZero(F => H, DF => DH);

  begin
    S := Find(S,R0,R1,Eps1,Eps2);
  end ZVal;

  procedure VVal(F,G: in Fun; X: Scalar; S: in out Scalar) is
--- solve G(x,s)+F(s,-s)=0
    R0: constant Scalar := -ScalOne;  -- guess
    R1: constant Scalar :=  ScalOne;  -- guess
    GX,GY,FX,FY: Scalar;

    function H(V: Scalar) return Scalar is
    begin
      return Val(X,V,G)+Val(V,-V,F);
    end H;

    function DH(V: Scalar) return Scalar is
    begin
      ValDer(X,V,G,GX,GY);
      ValDer(V,-V,F,FX,FY);
      return GY+FX-FY;
    end DH;

    function Find is new FindZero(F => H, DF => DH);

  begin
     S := Find(S,R0,R1,Eps1,Eps2);
  end VVal;

  procedure YVal(G: in Fun; S: in out Scalar) is
--- estimate solution of Gx(s^2,s)=[Gy(s^2,s)]^2
    R0: constant Scalar := Scal(-1454)/Rep(2048);  -- guess
    R1: constant Scalar := Scal(-1441)/Rep(2048);  -- guess
    GX,GY,GXX,GXY,GYY: Scalar;

    function H(Y: Scalar) return Scalar is
    begin
      ValDer(Y*Y,Y,G,GX,GY);
      return GX-Sqr(GY);
    end H;

    function DH(Y: Scalar) return Scalar is
    begin
      ValDer(Y*Y,Y,G,GX,GY,GXX,GXY,GYY);
      return Two*Y*GXX+GXY-Two*GY*(Two*Y*GXY+GYY);
    end DH;

    function Find is new FindZero(F => H, DF => DH);

  begin
    S := Find(S,R0,R1,Eps1,Eps2);
  end YVal;

end Funs2;