with Strings, IPowers, QPowers, Fun_Series, Newton, Roots, Globals;
use Strings, Globals;

pragma Elaborate_All (Strings,IPowers,QPowers,Fun_Series,Newton,Roots);

package body Chebs is

   use SM,SV,SV_More,Flt_Vectors,SV_Ops;

   procedure AddProd(S: in Scalar; C1: in CPoly1; C2: in out CPoly1; Tmp: in out Scalar) is
   begin
      if not IsZero(S) then
         for I in Freq loop AddProd(S,C1(I),C2(I),Tmp); end loop;
      end if;
   end AddProd;

   procedure SubProd(S: in Scalar; C1: in CPoly1; C2: in out CPoly1; Tmp: in out Scalar) is
   begin
      if not IsZero(S) then
         for I in Freq loop SubProd(S,C1(I),C2(I),Tmp); end loop;
      end if;
   end SubProd;

   --- basic

   function Info(Dummy: Cheb) return Scalar_Info is
      pragma Unreferenced(Dummy);
   begin
      return FInfo;
   end Info;

   function IsSharp(P: Cheb) return Boolean is
      C: CPoly1 renames P.C;
   begin
      return STrunc or else (IsSharp(P.C) and then IsZero(P.E));
   end IsSharp;

   function IsZero(P: Cheb) return Boolean is
   begin
      return IsZero(P.C) and then (STrunc or else IsZero(P.E));
   end IsZero;

   function "="(P1,P2: Cheb) return Boolean is
      C1: CPoly1 renames P1.C;
      C2: CPoly1 renames P2.C;
   begin
      if (C1 /= C2) then return False; end if;
      if STrunc then return True; end if;
      if (P1.E /= P2.E) then return False; end if;
      if IsZero(P1.E) then return True; end if;
      return (P1.R = P2.R);
   end "=";

   procedure SetZero(P: in out Cheb) is
   begin
      SetZero(P.C);
      if not STrunc then SetZero(P.E); end if;
      P.R := PolyFac*RadC;
   end SetZero;

   procedure Copy(P1: in Cheb; P2: in out Cheb) is
   begin
      Copy(P1.C,P2.C);
      if not STrunc then Copy(P1.E,P2.E); end if;
      P2.R := P1.R;
   end Copy;

   procedure Swap(P1,P2: in out Cheb) is
      R1: constant Flt := P1.R;
   begin
      Swap(P1.C,P2.C);
      if not STrunc then Swap(P1.E,P2.E); end if;
      P1.R := P2.R; P2.R := R1;
   end Swap;

   --- sets

   function Center0(P: Cheb) return Boolean is
   begin
      return Center0(P.C);
   end Center0;

   function Contains0(P: Cheb) return Logical is
      --- very crude
      L: constant Logical := Contains0(P.C);
   begin
      if L=True or else STrunc then return L; end if;
      return Uncertain;
   end Contains0;

   function Contains(P1,P2: Cheb) return Logical is
      --- very crude
      C1: CPoly1 renames P1.C;
      C2: CPoly1 renames P2.C;
      E1: Errs1 renames P1.E;
      E2: Errs1 renames P2.E;
      L: Logical := Contains(C1,C2);
   begin
      if STrunc then return L; end if;
      for I in Freq loop
         L := L and Value(E1(I) >= E2(I));
      end loop;
      if L=True then
         if IsPoly(P2) then return True; end if;
         if (P1.R>P2.R) then return False; end if;
      end if;
      return Uncertain;
   end Contains;

   procedure BallAt0(R: in Flt; P: in out Cheb) is
   begin
      SetZero(P.C);
      if not STrunc then
         SetZero(P.E);
         P.E(0) := abs(R);
      end if;
      P.R := RadC;
   end BallAt0;

   function BallAt0(R: Flt) return Cheb is
      P: Cheb;
   begin
      BallAt0(R,P);
      return P;
   end BallAt0;

   procedure ToErr(P: in out Cheb) is
      C: CPoly1 renames P.C;
   begin
      for I in Freq loop ToErr(C(I)); end loop;
   end ToErr;

   procedure ToErr(P1: in Cheb; P2: in out Cheb) is
      C1: CPoly1 renames P1.C;
      E2: Errs1 renames P2.E;
   begin
      SetZero(P2.C);
      if not STrunc then
         for I in Freq loop E2(I) := MaxNorm(C1(I)); end loop;
      end if;
      P2.R := P1.R;
   end ToErr;

   function ToErr(P: Cheb) return Cheb is
      Q: Cheb;
   begin
      ToErr(P,Q);
      return Q;
   end ToErr;

   procedure Center(P: in out Cheb) is
   begin
      if not STrunc then
         Center(P.C);
         SetZero(P.E);
      end if;
   end Center;

   procedure Center(P1: in Cheb; P2: in out Cheb) is
   begin
      Center(P1.C,P2.C);
      if not STrunc then SetZero(P2.E); end if;
      P2.R := PolyFac*RadC;
   end Center;

   function Center(P: Cheb) return Cheb is
      Q: Cheb;
   begin
      Center(P,Q);
      return Q;
   end Center;

   procedure ModCenter(P: in out Cheb) is
   begin
      ModCenter(P.C);
   end ModCenter;

   procedure ModCenter(P1: in Cheb; P2: in out Cheb) is
   begin
      ModCenter(P1.C,P2.C);
      if not STrunc then Copy(P1.E,P2.E); end if;
      P2.R := P1.R;
   end ModCenter;

   function ModCenter(P: Cheb) return Cheb is
      Q: Cheb;
   begin
      ModCenter(P,Q);
      return Q;
   end ModCenter;

   procedure ErrMult(R: in Radius; P: in out Cheb) is
      C: CPoly1 renames P.C;
      E: Errs1 renames P.E;
   begin
      if not STrunc then
         ErrMult(R,P.C);
         for I in EFreq loop E(I) := R*E(I); end loop;
      end if;
   end ErrMult;

   procedure Union(P1: in Cheb; P2: in out Cheb) is
      E1: Errs1 renames P1.E;
      E2: Errs1 renames P2.E;
   begin
      Union(P1.C,P2.C);
      if not STrunc then
         for I in EFreq loop E2(I) := RMax(E1(I),E2(I)); end loop;
      end if;
      if P2.R>P1.R then P2.R := P1.R; end if;
   end Union;

   function Union(P1,P2: Cheb) return Cheb is
      Q: Cheb := P2;
   begin
      Union(P1,Q);
      return Q;
   end Union;

   function QuasiDeg(C: CPoly1) return Natural is
   begin
      for D in reverse 1 .. CDeg loop
         if not IsZero(C(D)) then return D; end if;
      end loop;
      return 0;
   end QuasiDeg;

   function QuasiDeg(E: Errs1) return Natural is
   begin
      if not STrunc then
         for I in reverse 1 .. CDeg loop
            if E(I) /= Zero then return I; end if;
         end loop;
      end if;
      return 0;
   end QuasiDeg;

   procedure Intersection(P1: in Cheb; P2: in out Cheb; Empty: out Logical) is
      --- far from optimal
      C1: CPoly1 renames P1.C;
      C2: CPoly1 renames P2.C;
      F1: constant Integer := TrueF(P1.E);
      F2: constant Integer := TrueF(P2.E);
      E: Logical;
      Tmp: LF_Pointer;
   begin
      if F1>F2 then
         Pool.Allocate(Tmp);
         Copy(P2,Tmp.Data);
         Copy(P1,P2);
         Intersection(Tmp.Data,P2,Empty);
         Pool.Recycle(Tmp);
         return;
      end if;
      if P2.R<P1.R then P2.R := P1.R; end if;
      Intersection(C1(0),C2(0),Empty);
      for I in 1 .. IMin(F1-1,CDeg) loop
         Intersection(C1( I),C2( I),E);
         Empty := Empty or E;
      end loop;
   end Intersection;

   --- order

   function Sign(P: Cheb) return Integer is
   begin
      if not Frequency0(P) then
         raise Not_Certain;
      end if;
      return Sign(P.C(0));
   end Sign;

   function Compare(P1,P2: Cheb) return Integer is
   begin
      if not (Frequency0(P1) and then Frequency0(P2)) then
         raise Not_Certain;
      end if;
      return Compare(P1.C(0),P2.C(0));
   end Compare;

   function "<"(P1,P2: Cheb) return Boolean is
   begin
      return (Compare(P1,P2) < 0);
   end "<";

   function "<="(P1,P2: Cheb) return Boolean is
   begin
      return (Compare(P1,P2) <= 0);
   end "<=";

   function ">="(P1,P2: Cheb) return Boolean is
   begin
      return (Compare(P1,P2) >= 0);
   end ">=";

   function ">"(P1,P2: Cheb) return Boolean is
   begin
      return (Compare(P1,P2) > 0);
   end ">";

   procedure Min(P1: in Cheb; P2: in out Cheb) is
   begin
      if P1<P2 then Copy(P1,P2); end if;
   end Min;

   procedure Min(P1,P2: in Cheb; P3: in out Cheb) is
   begin
      if P1<P2 then Copy(P1,P3); else Copy(P2,P3); end if;
   end Min;

   function Min(P1,P2: Cheb) return Cheb is
   begin
      if P1<P2 then return P1; else return P2; end if;
   end Min;

   procedure Max(P1: in Cheb; P2: in out Cheb) is
   begin
      if P1>P2 then Copy(P1,P2); end if;
   end Max;

   procedure Max(P1,P2: in Cheb; P3: in out Cheb) is
   begin
      if P1>P2 then Copy(P1,P3); else Copy(P2,P3); end if;
   end Max;

   function Max(P1,P2: Cheb) return Cheb is
   begin
      if P1>P2 then return P1; else return P2; end if;
   end Max;

   function Sup(P: Cheb) return Flt is
   begin
      if not Frequency0(P) then
         raise Not_Implemented with "Fouriers1.Sup: P not constant";
      end if;
      return Sup(P.C(0));
   end Sup;

   function Inf(P: Cheb) return Flt is
   begin
      if not Frequency0(P) then
         raise Not_Implemented with "Fouriers1.Inf: P not constant";
      end if;
      return Inf(P.C(0));
   end Inf;

   --- addition and multiplication etc.

   procedure Neg(P: in out Cheb) is
   begin
      Neg(P.C);
   end Neg;

   procedure Neg(P1: in Cheb; P2: in out Cheb) is
   begin
      Neg(P1.C,P2.C);
      if not STrunc then Copy(P1.E,P2.E); end if;
      P2.R := P1.R;
   end Neg;

   function "-"(P: Cheb) return Cheb is
      Q: Cheb;
   begin
      Neg(P.C,Q.C);
      if not STrunc then Copy(P.E,Q.E); end if;
      Q.R := P.R;
      return Q;
   end "-";

   procedure Add(I: in Integer; P: in out Cheb) is
   begin
      Add(I,P.C(0));
   end Add;

   procedure Add(P1: in Cheb; P2: in out Cheb) is
   begin
      Add(P1.C,P2.C);
      if not STrunc then Add(P1.E,P2.E); end if;
      if P2.R>P1.R then P2.R := P1.R; end if;
   end Add;

   procedure Sum(P1,P2: in Cheb; P3: in out Cheb) is
   begin
      Sum(P1.C,P2.C,P3.C);
      if not STrunc then Sum(P1.E,P2.E,P3.E); end if;
      P3.R := RMin(P1.R,P2.R);
   end Sum;

   function "+"(P1,P2: Cheb) return Cheb is
      P3: Cheb;
   begin
      Sum(P1.C,P2.C,P3.C);
      if not STrunc then Sum(P1.E,P2.E,P3.E); end if;
      P3.R := RMin(P1.R,P2.R);
      return P3;
   end "+";

   procedure Sub(P1: in Cheb; P2: in out Cheb) is
   begin
      Sub(P1.C,P2.C);
      if not STrunc then Add(P1.E,P2.E); end if;
      if P2.R>P1.R then P2.R := P1.R; end if;
   end Sub;

   procedure Diff(P1,P2: in Cheb; P3: in out Cheb) is
   begin
      Diff(P1.C,P2.C,P3.C);
      if not STrunc then Sum(P1.E,P2.E,P3.E); end if;
      P3.R := RMin(P1.R,P2.R);
   end Diff;

   function "-"(P1,P2: Cheb) return Cheb is
      P3: Cheb;
   begin
      Diff(P1.C,P2.C,P3.C);
      if not STrunc then Sum(P1.E,P2.E,P3.E); end if;
      P3.R := RMin(P1.R,P2.R);
      return P3;
   end "-";

   procedure Mult(R: in Flt; V: in out Flt_Vector) is
   begin
      for J in V'Range loop
         V(J) := R*V(J);
      end loop;
   end Mult;

   procedure Prod(R: in Flt; V1: in Flt_Vector; V2: in out Flt_Vector) is
   begin
      for J in V1'Range loop
         V2(J) := R*V1(J);
      end loop;
   end Prod;

   procedure AddProd(R: in Flt; V1: in Flt_Vector; V2: in out Flt_Vector) is
   begin
      for J in V1'Range loop
         V2(J) := V2(J)+R*V1(J);
      end loop;
   end AddProd;

   procedure Mult(R: in Flt; P,Tmp: in out Cheb) is
      STmp: Scalar renames Tmp.C(0);
   begin
      Mult(R,P.C,STmp);
      if not STrunc then Mult(abs(R),P.E); end if;
   end Mult;

   procedure Mult(R: in Flt; P: in out Cheb) is
      STmp: Scalar;
   begin
      Mult(R,P.C,STmp);
      if not STrunc then Mult(abs(R),P.E); end if;
   end Mult;

   procedure Prod(R: in Flt; P1: in Cheb; P2: in out Cheb) is
   begin
      if R=Zero then
         SetZero(P2);
      else
         Prod(R,P1.C,P2.C);
         if not STrunc then Prod(abs(R),P1.E,P2.E); end if;
         P2.R := P1.R;
      end if;
   end Prod;

   function "*"(R: Flt; P: Cheb) return Cheb is
      Q: Cheb;
   begin
      Prod(R,P,Q);
      return Q;
   end "*";

   procedure AddProd(R: in Flt; P1: in Cheb; P2,Tmp: in out Cheb) is
      STmp: Scalar renames Tmp.C(0);
   begin
      if R /= Zero then
         AddProd(R,P1.C,P2.C,STmp);
         if not STrunc then AddProd(abs(R),P1.E,P2.E); end if;
         if P2.R>P1.R then P2.R := P1.R; end if;
      end if;
   end AddProd;

   procedure AddProd(R: in Flt; P1: in Cheb; P2: in out Cheb) is
   begin
      if R /= Zero then
         declare
            STmp: Scalar;
         begin
            AddProd(R,P1.C,P2.C,STmp);
         end;
         if not STrunc then AddProd(abs(R),P1.E,P2.E); end if;
         if P2.R>P1.R then P2.R := P1.R; end if;
      end if;
   end AddProd;

   procedure Mult(Q: in Rational; P: in out Cheb) is
   begin
      if LNum(Q)=0 then
         SetZero(P);
      else
         declare
            Tmp: Scalar;
         begin
            Mult(Flt(LNum(Q)),P.C,Tmp);
            Div(Flt(LDen(Q)),P.C,Tmp);
         end;
      end if;
   end Mult;

   function UpSumEven(E: Errs1) return Radius is
      R: Flt := Zero;
   begin
      for I in 0 .. EDeg loop R := R+E(I); end loop;
      return R;
   end UpSumEven;

   procedure Mult(P1: in Cheb; P2: in out Cheb) is
      Tmp: LF_Pointer;
   begin
      Pool.Allocate(Tmp);
      Prod(P1,P2,Tmp.Data);
      Copy(Tmp.Data,P2);
      Pool.Recycle(Tmp);
   end Mult;

   procedure Mult1(P1: in Cheb; P2: in out Cheb) is
      Tmp: LF_Pointer;
   begin
      Pool.Allocate(Tmp);
      Prod(P1,P2,Tmp.Data);
      Copy(Tmp.Data,P2);
      Pool.Recycle(Tmp);
   end Mult1;

   --  2*cos(m)*cos(n) =  cos(m-n) + cos(m+n)
   --  2*sin(m)*sin(n) =  cos(m-n) - cos(m+n)
   --  2*cos(m)*sin(n) = -sin(m-n) + sin(m+n)
   --  2*sin(m)*cos(n) =  sin(m-n) + sin(m+n)

   procedure Prod(P1,P2: in Cheb; P3: in out Cheb) is
      C1: CPoly1 renames P1.C;
      C2: CPoly1 renames P2.C;
      C3: CPoly1 renames P3.C;
      NotCos1: constant Boolean := False;
      NotSin1: constant Boolean := True;
      NotCos2: constant Boolean := False;
      NotSin2: constant Boolean := True;
      NotCos3: constant Boolean := False;
      NotSin3: constant Boolean := True;
      Tmp0: Scalar;

      procedure EvenPlus(K: in Natural; S: in out Scalar) is
         N: Integer;
         ---   cos(K.)  K=M+N
         ---   cos(m)*cos(n) --> +cos(m+n)
      begin
         if K <= CDeg then
            AddProd(C1(0),C2(K),S,Tmp0);             -- M=0
            if K=0 then return; end if;
            AddProd(C1(K),C2(0),S,Tmp0);             -- N=0
                                                     --        if K=1 then return; end if;
         end if;
         N := K/2;
         if 2*N=K then                             -- M=N
            AddProd(C1( N),C2( N),S,Tmp0);           --  +cos(m)*cos(n)
         end if;
         --      if K=2 then return; end if;
         for M in IMax(1,K-CDeg) .. (K-1)/2 loop -- 0<M<N<K
            N := K-M;
            AddProd(C1( M),C2( N),S,Tmp0);         --  +cos(m)*cos(n)
            AddProd(C1( N),C2( M),S,Tmp0);         --  +cos(n)*cos(m)
         end loop;
      end EvenPlus;

      procedure EvenMinus(K: in Natural; S: in out Scalar) is
         N: Integer;
         ---   cos(K.)  K=|M-N|
         ---   cos(m)*cos(n) --> +cos(m-n)
      begin
         if K=0 then                          -- M=N
            for M in 0 .. CDeg loop
               AddProd(C1( M),C2( M),S,Tmp0); --  +cos(m)*cos(n)
            end loop;
         else                                 -- K>0
            for M in K .. CDeg loop          -- M>N>=0
               N := M-K;
               AddProd(C1( M),C2( N),S,Tmp0); --  +cos(m)*cos(n)
               AddProd(C1( N),C2( M),S,Tmp0); --  +cos(n)*cos(m)
            end loop;
         end if;
      end EvenMinus;

      procedure EvenPlusMinus is
      begin
         if NotSin3 then
            for K in 0 .. CDeg loop                     -- coeff * coeff --> coeff (even)
               EvenPlus(K,C3(K));
               EvenMinus(K,C3(K));
            end loop;
         end if;
      end EvenPlusMinus;

      task type EvenPlusMinus_Task is end EvenPlusMinus_Task;
      task body EvenPlusMinus_Task is begin EvenPlusMinus; end EvenPlusMinus_Task;

      R3: Flt renames P3.R;
      E3: Errs1 renames P3.E;
   begin
      SetZero(C3);
      R3 := RMin(P1.R,P2.R);
      if FProd_Parallel then
         declare
            Run_EvenPlusMinus: EvenPlusMinus_Task;
         begin
            if not STrunc then
               ErrProd(P1,P2,R3,E3);
            end if;                                     -- all * err --> err
         end;
      else
         EvenPlusMinus;
         if not STrunc then
            ErrProd(P1,P2,R3,E3);
         end if;                                       -- all * err --> err
      end if;

      if not STrunc then
         declare
            Q3: constant Flt := One/R3;
            PR3,PQ3: Flt := Half;
            E0: Scalar;

            procedure EvenPlusE is
               PowR3,PowQ3: Flt;
            begin
               if NotSin3 then
                  PowR3 := PR3;
                  PowQ3 := PQ3;
                  for K in CDeg+1 .. EDeg loop            -- coeff * coeff --> err (even, sum)
                     PowR3 := PowR3*R3;
                     PowQ3 := PowQ3*Q3;
                     SetZero(E0);
                     EvenPlus(K,E0);
                     E3(K) := E3(K)+(PowR3+PowQ3)*MaxNorm(E0);
                  end loop;
               end if;
            end EvenPlusE;

         begin
            for K in 1 .. CDeg loop
               PR3 := R3*PR3;
               PQ3 := Q3*PQ3;
            end loop;
            EvenPlusE;
         end;
      end if;
      Mult(Half,P3);
   end Prod;

   procedure ErrProd(P1,P2: in Cheb; R3: in Radius; E3: out Errs1) is
      --- twice the product involving errors of P1 and P2
      C1: CPoly1 renames P1.C;
      C2: CPoly1 renames P2.C;
      E1: Errs1 renames P1.E;
      E2: Errs1 renames P2.E;
      Q3: constant Flt := One/R3;
      ETmp1,ETmp2: Flt := Zero;
      ETmp,EE3,NE1,NE2,Fac: Flt;
      N1,N2: Errs1;
   begin
      NE1 := MaxNorm(C1(0));
      NE2 := MaxNorm(C2(0));
      for K in EFreq loop                    --- coeff0 * err --> err  (freq diff)
         E3(K) := NE1*E2(K)+E1(K)*NE2;
         N1(K) := Zero;
         N2(K) := Zero;
      end loop;
      N1(0) := NE1; -- used later
      N2(0) := NE2; -- used later

      NE1 := UpSumEven(E1);
      NE2 := UpSumEven(E2);
      E3(0)  := E3( 0)+NE1*NE2;      --- err * err --> err (freq diff, even)

      Fac := One;
      for I in 1 .. CDeg loop                --- coeff * err --> err  (freq diff)
         Fac := Fac*R3;
         NE1 := Fac*MaxNorm(C1( I));
         NE2 := Fac*MaxNorm(C2( I));
         ETmp1 := ETmp1+E1( I);
         ETmp2 := ETmp2+E2( I);
         E3( 0) := E3( 0)+ETmp1*NE2+NE1*ETmp2;
         for K in I+1 .. EDeg loop
            E3(K-I) := E3(K-I)+E1( K)*NE2+NE1*E2( K);
         end loop;
         N1( I) := NE1; -- used below
         N2( I) := NE2; -- used below
      end loop;

      ETmp1 := N1(0)+E1(0);
      ETmp2 := N2(0)+E2(0);
      E3(0) := E3(0)+ETmp1*E2(0)+E1(0)*N2(0);
      for D in 1 .. 2*EDeg loop                                    --- all * err --> err (freq sum)
         if D <= EDeg then
            EE3 := ETmp1*E2( D)+E1(0)*N2( D)+E1( D)*ETmp2+N1( D)*E2(0);
         else
            EE3 := Zero;
         end if;
         for K in IMax(1,D-CDeg) .. IMin(D-1,CDeg) loop
            ETmp := N1( K)+E1( K);
            EE3 := EE3+ETmp*E2(D-K)+E1(K)*N2(D-K);
         end loop;
         if D <= EDeg then
            E3( D) := E3( D)+EE3;
         else
            E3( EDeg) := E3( EDeg)+EE3;
         end if;
      end loop;
   end ErrProd;

   function "*"(P1,P2: Cheb) return Cheb is
      P3: Cheb;
   begin
      Prod(P1,P2,P3);
      return P3;
   end "*";

   procedure AddProd(P1,P2: in Cheb; P3,Tmp: in out Cheb) is
   begin
      Prod(P1,P2,Tmp);
      Add(Tmp,P3);
   end AddProd;

   procedure AddProd(P1,P2: in Cheb; P3: in out Cheb) is
      Tmp: LF_Pointer;
   begin
      Pool.Allocate(Tmp);
      Prod(P1,P2,Tmp.Data);
      Add(Tmp.Data,P3);
      Pool.Recycle(Tmp);
   end AddProd;

   procedure SubProd(P1,P2: in Cheb; P3,Tmp: in out Cheb) is
   begin
      Prod(P1,P2,Tmp);
      Sub(Tmp,P3);
   end SubProd;

   procedure Div(R: in Flt; P,Tmp: in out Cheb) is
      E: Errs1 renames P.E;
      A: Flt renames Tmp.E(0);
      STmp: Scalar renames Tmp.C(0);
   begin
      Div(R,P.C,STmp);
      if not STrunc then
         A := abs(R);
         for I in EFreq loop E(I) := E(I)/A; end loop;
      end if;
   end Div;

   procedure Div(R: in Flt; P: in out Cheb) is
      E: Errs1 renames P.E;
      A: Flt;
      STmp: Scalar;
   begin
      Div(R,P.C,STmp);
      if not STrunc then
         A := abs(R);
         for I in EFreq loop E(I) := E(I)/A; end loop;
      end if;
   end Div;

   procedure Quot(P1: in Cheb; R: in Flt; P2: in out Cheb) is
      E1: Errs1 renames P1.E;
      E2: Errs1 renames P2.E;
      A: Flt;
   begin
      Quot(P1.C,R,P2.C);
      if not STrunc then
         A := abs(R);
         for I in EFreq loop E2(I) := E1(I)/A; end loop;
      end if;
      P2.R := P1.R;
   end Quot;

   function "/"(P: Cheb; R: Flt) return Cheb is
      Q: Cheb;
   begin
      Quot(P,R,Q);
      return Q;
   end "/";

   procedure Div(P1: in Cheb; P2: in out Cheb) is
      Tmp1,Tmp2: LF_Pointer;
   begin
      Pool.Allocate(Tmp1);
      Inv(P1,Tmp1.Data);
      Pool.Allocate(Tmp2);
      Prod(Tmp1.Data,P2,Tmp2.Data);
      Pool.Recycle(Tmp1);
      Copy(Tmp2.Data,P2);
      Pool.Recycle(Tmp2);
   end Div;

   procedure Quot(P1,P2: in Cheb; P3: in out Cheb) is
      Tmp: LF_Pointer;
   begin
      Pool.Allocate(Tmp);
      Inv(P2,Tmp.Data);
      Prod(P1,Tmp.Data,P3);
      Pool.Recycle(Tmp);
   end Quot;

   function "/"(P1,P2: Cheb) return Cheb is
      P,P3: Cheb;
   begin
      Inv(P1,P);
      Prod(P,P2,P3);
      return P3;
   end "/";

   procedure Mult(P1: in Cheb; P2,Tmp: in out Cheb) is
   begin
      Prod(P1,P2,Tmp);
      Copy(Tmp,P2);
   end Mult;

   function Simple(P: Cheb) return Boolean is
      --- only P.C(0) nonzero
      C: CPoly1 renames P.C;
   begin
      for I in 1 .. CDeg loop
         if not IsZero(C(I)) then return False; end if;
      end loop;
      return True;
   end Simple;

   procedure Inv(P: in out Cheb) is
      Tmp: LF_Pointer;
   begin
      Pool.Allocate(Tmp);
      Copy(P,Tmp.Data);
      Inv(Tmp.Data,P);
      Pool.Recycle(Tmp);
   end Inv;

   function Inv(P: Cheb) return Cheb is
      Q: Cheb;
   begin
      Inv(P,Q);
      return Q;
   end Inv;

   --- functions

   function IsReal(P: Cheb) return Boolean is
      C: CPoly1 renames P.C;
   begin
      if FInfo.IsComplex then
         if not IsReal(P.C) then return False; end if;
         if not IsPoly(P) then raise Not_Certain; end if;
      end if;
      return True;
   end IsReal;

   procedure Adjoint(P1: in Cheb; P2: in out Cheb) is
   begin
      Adjoint(P1.C,P2.C);
      if not STrunc then Copy(P1.E,P2.E); end if;
      P2.R := P1.R;
   end Adjoint;

   procedure Adjoint(P: in out Cheb) is
      C: CPoly1;
   begin
      Adjoint(P.C,C);
      Copy(C,P.C);
   end Adjoint;

   function Adjoint(P: Cheb) return Cheb is
   begin
      if FInfo.IsComplex then
         declare
            Q: Cheb;
         begin
            Adjoint(P,Q);
            return Q;
         end;
      end if;
      return P;
   end Adjoint;

   procedure Real_Part(P: in out Cheb) is
   begin
      if FInfo.IsComplex then
         Real_Part(P.C);
      end if;
   end Real_Part;

   procedure Imag_Part(P: in out Cheb) is
   begin
      Imag_Part(P.C);
      if not (STrunc or else FInfo.IsComplex) then
         SetZero(P.E);
      end if;
   end Imag_Part;

   procedure Eval0(P: in out Cheb) is
      C: CPoly1 renames P.C;
      E: Errs1 renames P.E;
      CSum: Scalar renames C(CDeg);
      ESum: Flt renames E(EDeg);
   begin
      Eval0(CSum);
      for I in reverse 0 .. CDeg-1 loop
         Eval0(C(I));
         Add(C(I),CSum);
      end loop;
      for I in -CDeg .. CDeg-1 loop SetZero(C(I)); end loop;
      Copy(CSum,C(0));
      if not STrunc then
         for I in reverse 0 .. EDeg-1 loop ESum := ESum+E(I); end loop;
         for I in -EDeg .. EDeg-1 loop E(I) := Zero; end loop;
         BallAt0(ESum,CSum);
         Eval0(CSum);
         Add(CSum,C(0));
         ESum := Zero;
      end if;
      SetZero(CSum);
   end Eval0;

   procedure Norm(P: in Cheb; Q: in out Cheb) is
   begin
      SetZero(Q);
      Norm1(P,Q.C(0));
   end Norm;

   function Norm(P: Cheb) return Cheb is
      Q: Cheb;
   begin
      SetZero(Q);
      Norm1(P,Q.C(0));
      return Q;
   end Norm;

   function MaxNorm(P: Cheb) return Radius is
      C: CPoly1 renames P.C;
      R: Flt renames P.R;
      PowR: Flt := One;
      N: Flt := MaxNorm(C(0));
   begin
      for I in 1 .. QuasiDeg(C) loop
         PowR := PowR*R;
         N := N+PowR*MaxNorm(C(I));
      end loop;
      if not STrunc then N := N+UpSum(P.E); end if;
      return N;
   end MaxNorm;

   procedure Sqr(P: in out Cheb) is
      Tmp: LF_Pointer;
   begin
      Pool.Allocate(Tmp);
      Prod(P,P,Tmp.Data);
      Copy(Tmp.Data,P);
      Pool.Recycle(Tmp);
   end Sqr;

   function Sqr(P: Cheb) return Cheb is
      Q: Cheb;
   begin
      Prod(P,P,Q);
      return Q;
   end Sqr;

   function Sqrt(P: Cheb) return Cheb is
      Q: Cheb;
   begin
      Sqrt(P,Q);
      return Q;
   end Sqrt;

   function Root(K: Positive; P: Cheb) return Cheb is
      Q: Cheb;
   begin
      Root(K,P,Q);
      return Q;
   end Root;

   function Exp(P: Cheb) return Cheb is
      Q: Cheb;
   begin
      Exp(P,Q);
      return Q;
   end Exp;

   function ArcCos(P: Cheb) return Cheb is
      Q: Cheb;
   begin
      ArcCos(P,Q);
      return Q;
   end ArcCos;

   function ArcSin(P: Cheb) return Cheb is
      Q: Cheb;
   begin
      ArcSin(P,Q);
      return Q;
   end ArcSin;

   procedure Cos(P1: in Cheb; P2: in out Cheb) is
      Tmp: LF_Pointer;
   begin
      Pool.Allocate(Tmp);
      CosSin(P1,P2,Tmp.Data);
      Pool.Recycle(Tmp);
   end Cos;

   function Cos(P: Cheb) return Cheb is
      PC,PS: Cheb;
   begin
      CosSin(P,PC,PS);
      return PC;
   end Cos;

   procedure Sin(P1: in Cheb; P2: in out Cheb) is
      Tmp: LF_Pointer;
   begin
      Pool.Allocate(Tmp);
      CosSin(P1,Tmp.Data,P2);
      Pool.Recycle(Tmp);
   end Sin;

   function Sin(P: Cheb) return Cheb is
      PC,PS: Cheb;
   begin
      CosSin(P,PC,PS);
      return PS;
   end Sin;

   procedure Simple_Random(P: in out Cheb) is
      R: constant Flt := P.R*Flt(9)/Flt(8); -- ad hoc
      Q: constant Flt := One/R;
      C: CPoly1 renames P.C;
      PowR,PowQ,Fac: Flt := Half;
   begin
      Simple_Random(C(0));
      for I in 1 .. CDeg loop
         PowR := PowR*R;
         PowQ := PowQ*Q;
         Fac := One/(PowR+PowQ);
         Simple_Random(C(I));
         Mult(Fac,C(I));
      end loop;
      if not STrunc then SetZero(P.E); end if;
   end Simple_Random;

   --- conversion and i/o

   function Approx(P: Cheb) return Flt is
   begin
      return Approx(P.C(0));
   end Approx;

   procedure Assign(I: in Integer; P: in out Cheb) is
   begin
      SetZero(P);
      Assign(I,P.C(0));
   end Assign;

   procedure Assign(Q: in Rational; P: in out Cheb) is
   begin
      SetZero(P);
      Assign(Q,P.C(0));
   end Assign;

   procedure Assign(R: in Flt; P: in out Cheb) is
   begin
      SetZero(P);
      Assign(R,P.C(0));
   end Assign;

   function Scal(I: Integer) return Cheb is
      P: Cheb;
   begin
      SetZero(P);
      Assign(I,P.C(0));
      return P;
   end Scal;

   function Scal(Q: Rational) return Cheb is
      P: Cheb;
   begin
      SetZero(P);
      Assign(Q,P.C(0));
      return P;
   end Scal;

   function Scal(R: Flt) return Cheb is
      P: Cheb;
   begin
      SetZero(P);
      Assign(R,P.C(0));
      return P;
   end Scal;

   procedure Enclose(R1,R2: in Flt; P: in out Cheb) is
   begin
      SetZero(P);
      Enclose(R1,R2,P.C(0));
   end Enclose;

   function Enclose(R1,R2: Flt) return Cheb is
      P: Cheb;
   begin
      SetZero(P);
      Enclose(R1,R2,P.C(0));
      return P;
   end Enclose;

   procedure Show1(N: in String; P: in Cheb; Hide0: in Boolean := True) is
   begin
      if not IsZero(P) then
         Show1(N & "R ",P.R);
         Show1(N & "C ",P.C);
         if not STrunc then
            Show1(N & "E ",P.E);
         end if;
      elsif not Hide0 then
         Show0(N & "0");
      end if;
   end Show1;

   procedure Show2(N: in String; P1,P2: in Cheb; Hide0: in Boolean := True) is
   begin
      if not (IsZero(P1) and IsZero(P2)) then
         Show2(N & "R ",P1.R,P2.R);
         Show2(N & "C ",P1.C,P2.C);
         if not STrunc then
            Show2(N & "E ",P1.E,P2.E);
         end if;
      elsif not Hide0 then
         Show0(N & "0 0");
      end if;
   end Show2;

   procedure Put(F: in File_Type; P: in Cheb; Decimal: in Boolean := False) is
      C: CPoly1 renames P.C;
      E: Errs1 renames P.E;
      DC: constant Integer := QuasiDeg(C);
      DE: constant Integer := QuasiDeg(E);
   begin
      Put(F,P.R,Decimal);
      Put(F,DC);
      for I in 0 .. DC loop
         Put(F,C(I),Decimal);
      end loop;
      if STrunc then
         Put(F,0);
         Put(F,Zero,Decimal);
      else
         Put(F,DE);
         for I in 0 .. DE loop
            Put(F,E(I),Decimal);
         end loop;
      end if;
   end Put;

   procedure Get(F: in File_Type; P: in out Cheb; Decimal: in Boolean := False) is
      E: Errs1 renames P.E;
      DC,DE: Integer := 0;
      R: Flt := Zero;
      S: Scalar;
   begin
      SetZero(P);
      Get(F,P.R,Decimal);
      Get(F,DC);
      for I in 0 .. DC loop
         Get(F,S,Decimal);
         AddCoeff(I,S,P);
      end loop;
      Get(F,DE);
      for I in 0 .. DE loop
         Get(F,R,Decimal);
         AddErr(I,R,E);
      end loop;
   end Get;

   procedure Write(FileName: in String; P: in Cheb; Decimal: in Boolean := False) is
      F: File_Type;
   begin
      if Verbosity>0 then Show0("Writing " & FileName); end if;
      Create(F,Out_File,FileName);
      Put(F,P,Decimal);
      Close(F);
   end Write;

   procedure Read(FileName: in String; P: in out Cheb; Decimal: in Boolean := False) is
      F: File_Type;
   begin
      if Verbosity>1 then Show0("Reading " & FileName); end if;
      Open(F,In_File,FileName);
      Get(F,P,Decimal);
      Close(F);
   end Read;

   --- misc

   function Get_Precision(P: Cheb) return Positive is
   begin
      return Get_Precision(P.C(0));
   end Get_Precision;

   --------------- other standard procedures

   function IsPoly(P: Cheb) return Boolean is
   begin
      return STrunc or else IsZero(P.E);
   end IsPoly;

   function Frequency0(P: Cheb) return Boolean is
      C: CPoly1 renames P.C;
   begin
      for I in 1 .. CDeg loop
         if not IsZero(C(I)) then return False; end if;
      end loop;
      return STrunc or else IsZero(P.E);
   end Frequency0;

   function Rho(P: Cheb) return Radius is
   begin
      return P.R;
   end Rho;

   procedure SetRho(R: in Radius; P: in out Cheb; Check: in Boolean := True) is
   begin
      if (R>P.R) and then Check and then TrueF(P.E) <= EDeg then
         raise Domain_Error with Show0("Fouriers1.SetRho: cannot increase domain");
      end if;
      P.R := R;
   end SetRho;

   procedure SetZero(R: in Radius; P: in out Cheb) is
   begin
      SetZero(P.C);
      if not STrunc then SetZero(P.E); end if;
      P.R := R;
   end SetZero;

   function TrueF(E: Errs1) return Natural is
      --- smallest |K| for which E(K) is nonzero, or EDeg+1
   begin
      if STrunc then return EDeg+1; end if;
      if E(0)=Zero then
         for K in 0 .. EDeg loop
            if E(K) /= Zero then return K; end if;
         end loop;
         return EDeg+1;
      end if;
      return 0;
   end TrueF;

   procedure LiftErrs(P: in out Cheb; F: in Natural := 0) is
      C: CPoly1 renames P.C;
      E: Errs1 renames P.E;
      F0: constant Integer := IMin(F,CDeg+1);
      PF: constant Integer := TrueF(E);
      B: Scalar;
   begin
      if PF >= F0 then return; end if;
      BallAt0(E(0),B);
      Add(B,C(0));
      E(1) := E(1)+E(0);
      E(0) := Zero;
      if F0=1 then return; end if;
      declare
         R: Flt renames P.R;
         Fac: Flt := One/(-R);
         Tmp: Scalar;
      begin
         for I in 1 .. F0-1 loop --- distribute E(I)
            if I >= PF then
               Fac := One/(-Fac);
               BallAt0(E(I),B);
               AddProd(Fac,B,C(I),Tmp);
               E(I+1) := E(I+1)+E(I);
            end if;
         end loop;
      end;
   end LiftErrs;

   procedure Assign(S: in Scalar; P: in out Cheb) is
   begin
      SetZero(P);
      Copy(S,P.C(0));
   end Assign;

   function Scal(S: Scalar) return Cheb is
      P: Cheb;
   begin
      SetZero(P);
      Copy(S,P.C(0));
      return P;
   end Scal;

   function UpInvNormFac(R: Radius; I: Integer) return Radius is
      PowR,PowQ: Flt := NegHalf;
   begin
      for K in 1 .. abs(I) loop
         PowR := PowR*R;
         PowQ := PowQ/R;
      end loop;
      return One/(-(PowR+PowQ));
   end UpInvNormFac;

   procedure EvenCoeff(I: in Natural; P: in Cheb; S: in out Scalar) is
      E: Errs1 renames P.E;
      ESum: Flt := Zero;
   begin
      if not STrunc then
         for K in reverse 0 .. IMin(I,EDeg) loop
            ESum := ESum+E(K);
         end loop;
      end if;
      if ESum=Zero then
         SetZero(S);
      else
         BallAt0(UpInvNormFac(P.R,I)*ESum,S);
      end if;
      if I <= CDeg then Add(P.C(I),S); end if;
   end EvenCoeff;

   procedure Coeff(I: in Integer; P: in Cheb; S: in out Scalar) is
      C: CPoly1 renames P.C;
      E: Errs1 renames P.E;
   begin
      if STrunc then
         if abs(I) <= CDeg then Copy(C(I),S); else SetZero(S); end if;
      elsif I=0 then
         BallAt0(E(0),S);
         Add(P.C(0),S);
      else
         EvenCoeff(I,P,S);
      end if;
   end Coeff;

   function Coeff(I: Integer; P: Cheb) return Scalar is
      S: Scalar;
   begin
      Coeff(I,P,S);
      return S;
   end Coeff;

   function UpNormFac(R: Radius; I: Integer) return Radius is
      Q: constant Flt := One/R;
      PowR,PowQ: Flt := Half;
   begin
      for K in 1 .. abs(I) loop
         PowR := PowR*R;
         PowQ := PowQ*Q;
      end loop;
      return PowR+PowQ;
   end UpNormFac;

   procedure AddErr(I: in Integer; R: in Flt; E: in out Errs1) is
   begin
      if STrunc then
         null;
      elsif abs(I) <= EDeg then
         E(I) := E(I)+abs(R);
      else
         E( EDeg) := E( EDeg)+abs(R);
      end if;
   end AddErr;

   procedure AddCoeff(I: in Integer; S: in Scalar; P: in out Cheb) is
   begin
      if abs(I) <= CDeg then
         Add(S,P.C(I));
      elsif not STrunc then
         AddErr(I,UpNormFac(P.R,I)*MaxNorm(S),P.E);
      end if;
   end AddCoeff;

   procedure Add(S: in Scalar; P: in out Cheb) is
   begin
      Add(S,P.C(0));
   end Add;

   procedure Mult(S: in Scalar; P: in out Cheb) is
      R: Flt;
   begin
      if IsZero(S) then
         SetZero(P);
      else
         Mult(S,P.C);
         R := MaxNorm(S);
         if not STrunc then Mult(R,P.E); end if;
      end if;
   end Mult;

   procedure Prod(S: in Scalar; P1: in Cheb; P2: in out Cheb) is
   begin
      if IsZero(S) then
         SetZero(P2);
      else
         Prod(S,P1.C,P2.C);
         P2.R := P1.R;
         if not STrunc then Prod(MaxNorm(S),P1.E,P2.E); end if;
      end if;
   end Prod;

   function "*"(S: Scalar; P: Cheb) return Cheb is
      Q: Cheb;
   begin
      Prod(S,P,Q);
      return Q;
   end "*";

   procedure AddProd(S: in Scalar; P1: in Cheb; P2: in out Cheb) is
   begin
      if not IsZero(S) then
         declare
            Tmp: Scalar;
         begin
            AddProd(S,P1.C,P2.C,Tmp);
         end;
         if not STrunc then
            AddProd(MaxNorm(S),P1.E,P2.E);
         end if;
         if P2.R>P1.R then P2.R := P1.R; end if;
      end if;
   end AddProd;

   procedure Val(P: in Cheb; S1: in Scalar; S2: in out Scalar) is
      C: CPoly1 renames P.C;
      E: Errs1  renames P.E;
      D: constant Integer := QuasiDeg(C);
      R,Err: Radius;
--      TA,TB,TC: Scalar;
      AS: Scalar := ArcCos(S1);
   begin
      if MaxNorm(S1)>One then raise Undefined with "Chebs.Val error: argument too large"; end if;
      if IsPoly(P) then
         SetZero(S2);
      else
         Err := Zero;
         R := One;
         for I in EFreq loop
            Err := Err+R*E(I);
            R := R/P.R;
         end loop;
         BallAt0(Err,S2);
      end if;
      for I in 0..D loop
         Add(Cos(Flt(I)*AS)*C(I),S2);
      end loop;
      --  Add(C(0)+C(1)*S1,S2);
      --  TA := Scal(1);
      --  TB := S1;
      --  for I in 2 .. D loop
      --     TC := Two*S1*TB-TA;
      --     if MaxNorm(TC)>One then
      --        Add(C(I)*BallAt0(One),S2);
      --     else
      --        Add(TC*C(I),S2);
      --     end if;
      --     Copy(TB,TA);
      --     Copy(TC,TB);
      --  end loop;
   end Val;

   procedure SetIV(P: in out Cheb; T: in Scalar) is
      C: CPoly1 renames P.C;
      E: Errs1  renames P.E;
      D: constant Integer := QuasiDeg(C);
      R,Err: Radius;
      S: Scalar;
   begin
      if IsPoly(P) then
         SetZero(S);
      else
         Err := Zero;
         R := One;
         for I in EFreq loop
            Err := Err+R*E(I);
            R := R/P.R;
         end loop;
         BallAt0(Err,S);
      end if;
      for I in 1 .. D loop
         if I mod 2 = 0 then
            Add(C(I),S);
         else
            Sub(C(I),S);
         end if;
      end loop;
      C(0) := T-S;
      E(0) := Zero;
   end SetIV;

   procedure SumEven(P: in Cheb; S: in out Scalar) is
      C: CPoly1 renames P.C;
      E: Errs1  renames P.E;
      D: constant Integer := QuasiDeg(C);
      R,Err: Radius;
   begin
      if IsPoly(P) then
         SetZero(S);
      else
         Err := Zero;
         R := One;
         for I in EFreq loop
            Err := Err+R*E(I);
            R := R/P.R;
         end loop;
         BallAt0(Err,S);
      end if;
      for I in 2 .. D loop
         if I mod 2 = 0 then Add(C(I),S); end if;
      end loop;
   end SumEven;

   procedure SumOdd(P: in Cheb; S: in out Scalar) is
      C: CPoly1 renames P.C;
      E: Errs1  renames P.E;
      D: constant Integer := QuasiDeg(C);
      R,Err: Radius;
   begin
      if IsPoly(P) then
         SetZero(S);
      else
         Err := Zero;
         R := One;
         for I in EFreq loop
            Err := Err+R*E(I);
            R := R/P.R;
         end loop;
         BallAt0(Err,S);
      end if;
      for I in 3 .. D loop
         if I mod 2 = 1 then Add(C(I),S); end if;
      end loop;
   end SumOdd;

   function Val(P: Cheb; S: Scalar) return Scalar is
      V: Scalar;
   begin
      Val(P,S,V);
      return V;
   end Val;

   function SumEven(P: Cheb) return Scalar is
      V: Scalar;
   begin
      SumEven(P,V);
      return V;
   end SumEven;

   function SumOdd(P: Cheb) return Scalar is
      V: Scalar;
   begin
      SumOdd(P,V);
      return V;
   end SumOdd;

   procedure Val(P1: in Cheb; S: in Scalar; P2: in out Cheb) is
   begin
      SetZero(P2);
      Val(P1,S,P2.C(0));
   end Val;

   procedure Norm1(R: in Radius; C: in CPoly1; S: in out Scalar) is
      N,Fac,Tmp: Scalar;
   begin
      SetZero(S);
      Fac := SOne;
      for I in 0 .. CDeg loop
         Norm(C(I),N);
         AddProd(Fac,N,S,Tmp);
         Mult(R,Fac,Tmp);
      end loop;
   end Norm1;

   function UpSum(E: Errs1) return Radius is
      R: Flt := Zero;
   begin
      for I in 0 .. EDeg loop R := R+E(I); end loop;
      return R;
   end UpSum;

   procedure Norm1(P: in Cheb; S: in out Scalar) is
   begin
      Norm1(P.R,P.C,S);
      if not STrunc then
         declare
            E: Scalar;
         begin
            Enclose(Zero,UpSum(P.E),E);
            Add(E,S);
         end;
      end if;
   end Norm1;

   function Norm1(P: Cheb) return Scalar is
      S: Scalar;
   begin
      Norm1(P,S);
      return S;
   end Norm1;

   procedure Monom(I: in Freq; P: in out Cheb) is
   begin
      SetZero(P);
      Assign(1,P.C(I));
   end Monom;

   procedure HiUnitBall(I: in Integer; P: in out Cheb; R: in Radius := RadC) is
   begin
      SetZero(P);
      AddErr(I,One,P.E);
      P.R := R;
   end HiUnitBall;

   procedure Compose(C: in SV.Polynom1; P1: in Cheb; P2: in out Cheb) is
      D: constant Integer := QuasiDeg(C);
      Tmp: LF_Pointer;
   begin
      SetZero(P2);
      Add(C(D),P2.C(0));
      Pool.Allocate(Tmp);
      for I in reverse 0 .. D-1 loop
         Mult(P1,P2,Tmp.Data);
         Add(C(I),P2.C(0));
      end loop;
      Pool.Recycle(Tmp);
   end Compose;

   procedure ScalProd(P,Q: in Cheb; S: in out Scalar) is
      Err,R: Radius;
   begin
      if IsPoly(P) and then IsPoly(Q) then
         SetZero(S);
      else
         Err := Zero;
         R := One;
         for I in EFreq loop
            Err := Err+R*P.E(I)*Q.E(I);
            R := R/(P.R*Q.R);
         end loop;
         BallAt0(Err,S);
      end if;
      for I in Freq loop
         S := S+P.C(I)*Q.C(I);
      end loop;
   end ScalProd;

   procedure DerAtOne(P: in Cheb; S: in out Scalar) is
   begin
      SetZero(S);
      if MaxNorm(P)>Zero then
         for I in 1..QuasiDeg(P.C) loop
            Add(Flt(I*I)*P.C(I),S);
         end loop;
         if not STrunc then
            raise Undefined;
         end if;
      end if;
   end DerAtOne;

   procedure AntiDer(C: in out CPoly1) is
      D: CPoly1;
      QD: Natural := IMin(QuasiDeg(C)+1,CDeg-1);
   begin
      SetZero(D);
      D(1) := C(0)-Half*C(2);
      for I in 2 .. QD  loop
         D(I) := (C(I-1)-C(I+1))/Flt(2*I);
      end loop;
      Swap(C,D);
   end AntiDer;

   procedure AntiDer(P: in out Cheb) is
      E: Errs1 renames P.E;
      QD: Natural := IMin(QuasiDeg(E),EDeg-1);
      EE: Errs1;
      Tmp: Flt renames E(0);
   begin
      if not STrunc then
         SetZero(EE);
         EE(1) := E(0)+Half*E(2);
         for I in 2 .. QD loop
            EE(I) := (E(I-1)+E(I+1))/Flt(2*I);
         end loop;
         Swap(E,EE);
      end if;
      AntiDer(P.C);
   end AntiDer;

   procedure AntiDer(P1: in Cheb; P2: in out Cheb) is
   begin
      Copy(P1,P2);
      AntiDer(P2);
   end AntiDer;

   procedure AntiDerP(P: in out Cheb) is
      V1,V2: Scalar;
   begin
      AntiDer(P);
      Val(P,Scal(-1),V1);
      Val(P,Scal(1),V2);
      Add(Half*(V1-V2),P.C(1));   ----- the average of P' is 0
   end AntiDerP;

   procedure AntiDer2D(P: in out Cheb) is
      Even,Odd: Scalar;
   begin
      AntiDer(P);
      AntiDer(P);
      SumEven(P,Even);
      SumOdd(P,Odd);
      P.C(0) := -Even;            ----  P( 1)=0
      P.C(1) := -Odd;             ----  P(-1)=0
      if Not_STrunc then
         P.E(0) := Zero;
         P.E(1) := Zero;
      end if;
   end AntiDer2D;

   function Integral(P: Cheb; S1,S2: Scalar) return Scalar is
      Tmp: LF_Pointer;
      V1,V2: Scalar;
   begin
      Pool.Allocate(Tmp);
      AntiDer(P,Tmp.Data);
      Val(Tmp.Data,S1,V1);
      Val(Tmp.Data,S2,V2);
      Pool.Recycle(Tmp);
      Sub(V1,V2);
      return V2;
   end Integral;

   procedure NewQuot(S0,S1: in Scalar; Q: in out Flt; N: in out Integer) is
   begin
      if not (IsZero(S0) or else IsZero(S1)) then
         Q := Q*MaxNorm(S0)/MaxNorm(S1);
         N := N+1;
      end if;
   end NewQuot;

   function GuessRho(P: Cheb) return Radius is
      use Flt_EF;
      DMin: constant Integer := CDeg/5;
      DMax: constant Integer := (4*CDeg)/5;
      C: CPoly1 renames P.C;
      N1: Integer := 0;
      Q1: Flt := One;
   begin
      for I in DMin .. DMax loop
         NewQuot(C( I),C( I+1),Q1,N1);
      end loop;
      return Exp(Log(Q1)/Flt(N1));
   end GuessRho;

   procedure Test_Some_Ops is
      R: constant Flt := Flt(5)/Flt(4);
      Fac: constant Flt := 1.0/16.0;
      X: constant Scalar := Scal(1.0/5.0);
      Q: constant Rational := 2/3;
      S0,S: Scalar;
      F0,F: Cheb;
   begin
      SetZero(R,F0);
      Simple_Random(F0);
      Mult(Fac/MaxNorm(F0),F0);
      Add(1,F0);
      S0 := Val(F0,X);              -- S0 := F0(X);
      S := Inv(S0);
      F := Inv(F0);
      Show2("Inv(F(X)) ",Val(F,X),S);
      S := Sqrt(S0);
      F := Sqrt(F0);
      Show2("Sqrt(F(X)) ",Val(F,X),S);
      S := Root(7,S0);
      F := Root(7,F0);
      Show2("Root(7,F(X)) ",Val(F,X),S);
      S := Exp(S0);
      F := Exp(F0);
      Show2("Exp(F(X)) ",Val(F,X),S);
      S := Cos(S0);
      F := Cos(F0);
      Show2("Cos(F(X)) ",Val(F,X),S);
      S := Sin(S0);
      F := Sin(F0);
      Show2("Sin(F(X)) ",Val(F,X),S);
      S := Log(S0);
      F := Log(F0);
      Show2("Log(F(X)) ",Val(F,X),S);
      QPower(Q,S0,S);
      QPower(Q,F0,F);
      Show2("(F(X)^Q) ",Val(F,X),S);
      S := ArcCos(S0/Four);
      F := ArcCos(F0/Four);
      Show2("ArcCos(F(X)/4) ",Val(F,X),S);
   end Test_Some_Ops;

   -----------------------------------------------------------------------
   -----------------------------------------------------------------------
   package FIP is new IPowers (Scalar => Cheb);
   package FQP is new QPowers (Scalar => Cheb, IPower => FIP.IPower);
   -----------------------------------------------------------------------
   -----------------------------------------------------------------------

   function "**"(P: Cheb; I: Integer) return Cheb is
      Q: Cheb;
   begin
      FIP.IPower(I,P,Q);
      return Q;
   end "**";

   procedure QPower(Q: in Rational; P1: in Cheb; P2: in out Cheb) renames FQP.QPower;

   -----------------------------------------------
   --- stuff for Fun_Series ----------------------
   -----------------------------------------------

   procedure Exp_Split(P1: in Cheb; P2: in out Cheb; S: in out Scalar) is
   begin
      Copy(P1,P2);
      LiftErrs(P2,1);
      Exp(P2.C(0),S);
      SetZero(P2.C(0));
   end Exp_Split;

   procedure CosSin_Split(P1: in Cheb; P2: in out Cheb; C,S: in out Scalar) is
   begin
      Copy(P1,P2);
      LiftErrs(P2,1);
      Cos(P2.C(0),C);
      Sin(P2.C(0),S);
      SetZero(P2.C(0));
   end CosSin_Split;

  procedure Log_Split(P1: in Cheb; P2: in out Cheb; S: in out Scalar) is
  begin
    Copy(P1,P2);
    LiftErrs(P2,1);
    Inv(P2.C(0),S);
    Mult(S,P2);
    SetZero(P2.C(0));
    Log(P1.C(0),S);
  end Log_Split;

  procedure Inv_Split(P1: in Cheb; P2: in out Cheb; S: in out Scalar) is
  begin
    Inv(P1.C(0),S);
    Prod(S,P1,P2);
    SetZero(P2.C(0));
  end Inv_Split;

  procedure QPower_Split(Q: in Rational; P1: in Cheb; P2: in out Cheb; S: in out Scalar) is
  begin
    Copy(P1,P2);
    LiftErrs(P2,1);
    Inv(P2.C(0),S);
    Mult(S,P2);
    SetZero(P2.C(0));
    QPower(Q,P1.C(0),S);
  end QPower_Split;

   --------------------------------------------------------------------------------------
   --------------------------------------------------------------------------------------
   package FFS is new Fun_Series (Numeric => STrunc, Scalar => Scalar, Fun => Cheb,
                                  Inv_Split => Inv_Split, Log_Split => Log_Split,
                                  QPower_Split => QPower_Split);
   --------------------------------------------------------------------------------------
   --------------------------------------------------------------------------------------

   function ExpTerms(R: Radius; Prec: Positive) return Natural is
      Eps: constant Flt := Half**Prec;
      A: Flt := One;
   begin
      for N in 1 .. 4096 loop
         A := A*R/Flt(N);
         if A<Eps then return N; end if;
      end loop;
      raise Sorry with "Chebs.ExpTerms: should not get here";
   end ExpTerms;

   procedure Exp(P1: in Cheb; P2: in out Cheb) is
      Iter: Integer; --- see below
      K: Integer := 0;
      R: Flt := MaxNorm(P1);
      Tmp: LF_Pointer;
   begin
      if R=Zero then Assign(1,P2); return; end if;
      while R>One loop
         K := K+1;
         R := Half*R;
      end loop;
      Iter := ExpTerms(R,Get_Precision(P1));
      if not STrunc then Iter := Iter+4; end if;
      if K=0 then
         FFS.Series_Exp(P1,P2,Iter);
      else
         Pool.Allocate(Tmp);
         Prod(Half**K,P1,Tmp.Data);
         FFS.Series_Exp(Tmp.Data,P2,Iter);
         for J in 1 .. K loop
            Prod(P2,P2,Tmp.Data);
            Copy(Tmp.Data,P2);
         end loop;
         Pool.Recycle(Tmp);
      end if;
   end Exp;

   procedure Log(P1: in Cheb; P2: in out Cheb) is
      Iter: constant Integer := Choose(STrunc,CDeg,IMax(CDeg+8,32));
   begin
      FFS.Series_Log(P1,P2,Iter);
   end Log;

   function Log(P: Cheb) return Cheb is
      Q: Cheb;
   begin
      Log(P,Q);
      return Q;
   end Log;

   procedure CosSin(P: in Cheb; PC,PS: in out Cheb) is
      Iter: Integer; --- see below
      K: Integer := 0;
      R: Flt := MaxNorm(P);
      Tmp: LF_Pointer;
   begin
      if R=Zero then Assign(1,PC); SetZero(PS); return; end if;
      while R>One loop
         K := K+1;
         R := Half*R;
      end loop;
      Iter := ExpTerms(R*R,Get_Precision(P));
      if not STrunc then Iter := Iter+4; end if;
      if K=0 then
         FFS.Series_CosSin(P,PC,PS,Iter);
      else
         Pool.Allocate(Tmp);
         Prod(Half**K,P,Tmp.Data);
         FFS.Series_CosSin(Tmp.Data,PC,PS,Iter);
         for J in 1 .. K loop
            Prod(PC,PS,Tmp.Data);
            Prod(Two,Tmp.Data,PS);  -- PS := 2*PC*PS
            Prod(PC,PC,Tmp.Data);
            Prod(Two,Tmp.Data,PC);  -- PC := 2*PC*PC-1
            Add(-1,PC);
         end loop;
         Pool.Recycle(Tmp);
      end if;
   end CosSin;

   -------------------------------------------------------------------------------------------
   -------------------------------------------------------------------------------------------
   package GDFT is new Generic_DFT (L => CDeg, Scalar => Scalar, Scalar_Array => SV.Vector);
   -------------------------------------------------------------------------------------------
   -------------------------------------------------------------------------------------------

   procedure Approx_Inv(P: in out Cheb) is
      procedure F1 is new GDFT.F1 (F => Inv);
   begin
      Center(P);
      F1(P.C);
      if not STrunc then SetZero(P.E); end if;
   end Approx_Inv;

   procedure Approx_Sqrt(P: in out Cheb) is
      procedure F1 is new GDFT.F1 (F => Sqrt);
   begin
      Center(P);
      F1(P.C);
      if not STrunc then SetZero(P.E); end if;
   end Approx_Sqrt;

   procedure Approx_Root(K: in Positive; P: in out Cheb) is
      procedure RootK(S1: in Scalar; S2: in out Scalar) is
      begin
         Root(K,S1,S2);
      end RootK;
      procedure F1 is new GDFT.F1 (F => RootK);
   begin
      Center(P);
      F1(P.C);
      if not STrunc then SetZero(P.E); end if;
   end Approx_Root;

   procedure Approx_ArcCos(P: in out Cheb) is
      procedure F1 is new GDFT.F1 (F => ArcCos);
   begin
      Center(P);
      F1(P.C);
      if not STrunc then SetZero(P.E); end if;
   end Approx_ArcCos;

   procedure Approx_ArcSin(P: in out Cheb) is
      procedure F1 is new GDFT.F1 (F => ArcSin);
   begin
      Center(P);
      F1(P.C);
      if not STrunc then SetZero(P.E); end if;
   end Approx_ArcSin;

   ------------------------------------------------------------------------------------------
   ------------------------------------------------------------------------------------------
   package Newt is new Newton(Numeric => STrunc, Scalar => Cheb, Operator => Cheb);
   ------------------------------------------------------------------------------------------
   ------------------------------------------------------------------------------------------

   procedure Inv(P1: in Cheb; P2: in out Cheb) is
      Steps: constant Integer := 64;  -- ad hoc
   begin
      Copy(P1,P2);
      if Frequency0(P1) then
         Inv(P1.C(0),P2.C(0));
      else
         Approx_Inv(P2);
         Newt.InvNewton(P1,P2,Steps);
      end if;
   end Inv;

   procedure Sqrt(P1: in Cheb; P2: in out Cheb) is
   begin
      Copy(P1,P2);
      if Frequency0(P1) then
         Sqrt(P1.C(0),P2.C(0));
      else
         Approx_Sqrt(P2);
         declare
            Steps: constant Integer := IMax(CDeg,64);  --- need to be chosen better
            procedure F(X: in Cheb; Y: in out Cheb) is
            begin
               Prod(X,X,Y);
               Sub(P1,Y);
            end F;
            procedure DF(X: in Cheb; Y: in out Cheb) is
            begin
               Prod(Two,X,Y);
            end DF;
            procedure Improve is new Newt.FindZero(F => F, DF => DF);
         begin
            Improve(P2,Steps);
         end;
      end if;
   end Sqrt;

   procedure Root(K: in Positive; P1: in Cheb; P2: in out Cheb) is
   begin
      if K=1 then Copy(P1,P2); return; end if;
      if K=2 then Sqrt(P1,P2); return; end if;
      Copy(P1,P2);
      if Frequency0(P1) then
         Root(K,P1.C(0),P2.C(0));
      else
         Approx_Root(K,P2);
         declare
            Steps: constant Integer := IMax(CDeg,64);  --- need to be chosen better
            procedure F(X: in Cheb; Y: in out Cheb) is
            begin
               FIP.IPower(K,X,Y);
               Sub(P1,Y);
            end F;
            procedure DF(X: in Cheb; Y: in out Cheb) is
            begin
               FIP.IPower(K-1,X,Y);
               Mult(Flt(K),Y);
            end DF;
            procedure Improve is new Newt.FindZero(F => F, DF => DF);
         begin
            Improve(P2,Steps);
         end;
      end if;
   end Root;

   procedure ArcCos(P1: in Cheb; P2: in out Cheb) is
   begin
      Copy(P1,P2);
      if Frequency0(P1) then
         ArcCos(P1.C(0),P2.C(0));
      else
         Approx_ArcCos(P2);
         declare
            Steps: constant Integer := IMax(CDeg,64);  --- need to be chosen better
            procedure F(X: in Cheb; Y: in out Cheb) is
            begin
               Cos(X,Y);
               Sub(P1,Y);
            end F;
            procedure DF(X: in Cheb; Y: in out Cheb) is
            begin
               Sin(X,Y);
               Neg(Y);
            end DF;
            procedure Improve is new Newt.FindZero(F => F, DF => DF);
         begin
            Improve(P2,Steps);
         end;
      end if;
   end ArcCos;

   procedure ArcSin(P1: in Cheb; P2: in out Cheb) is
   begin
      Copy(P1,P2);
      if Frequency0(P1) then
         ArcSin(P1.C(0),P2.C(0));
      else
         Approx_ArcSin(P2);
         declare
            Steps: constant Integer := IMax(CDeg,64);  --- need to be chosen better
            procedure F(X: in Cheb; Y: in out Cheb) is
            begin
               Sin(X,Y);
               Sub(P1,Y);
            end F;
            procedure DF(X: in Cheb; Y: in out Cheb) is
            begin
               Cos(X,Y);
            end DF;
            procedure Improve is new Newt.FindZero(F => F, DF => DF);
         begin
            Improve(P2,Steps);
         end;
      end if;
   end ArcSin;

   ----------------------------------------------
   ----------------------------------------------
   package FR is new Roots (Scalar => Cheb);
   ----------------------------------------------
   ----------------------------------------------

   procedure Roots2(B,C: in Cheb; U1,U2,V: in out Cheb) is
      Steps: constant Integer := 64;
   begin
      if Frequency0(B) and then Frequency0(C) then
         SetZero(U1);
         SetZero(U2);
         SetZero(V);
         Roots2(B.C(0),C.C(0),U1.C(0),U2.C(0),V.C(0));
      else
         FR.NewtonRoots2(B,C,U1,U2,V,Steps);
      end if;
   end Roots2;

   procedure Roots3(B,C,D: in Cheb; U0,U1,U2,V: in out Cheb) is
      Steps: constant Integer := IMax(CDeg,64);
   begin
      if Frequency0(B) and then Frequency0(C) and then Frequency0(D) then
         SetZero(U0);
         SetZero(U1);
         SetZero(U2);
         SetZero(V);
         Roots3(B.C(0),C.C(0),D.C(0),U0.C(0),U1.C(0),U2.C(0),V.C(0));
      else
         FR.NewtonRoots3(B,C,D,U0,U1,U2,V,Steps);
      end if;
   end Roots3;

end Chebs;
