with Ints,Strings,Globals;
use Ints,Strings,Globals;

pragma Elaborate_All (Globals);

package body Fun_Series is

   procedure Series(X: in Fun; Y: in out Fun; Decreasing: in Boolean; Iter: in Positive) is
      --- returns 1+F(1)*X*[1+F(2)*X*[1+F(3)*X*[1+...]]]
      Tmp: Fun;
   begin
      SetZero(Y);
      if IsZero(X) then Add(1,Y); return; end if;
      if not Numeric then
         declare
            N: Radius := MaxNorm(X);
            Q: Rational;
         begin
            if Decreasing then -- assume Flt operations are rounded up
               Q := Fac(Iter+1);
               N := N*Flt(abs(INum(Q)));
               N := N/Flt(IDen(Q));
            end if;
            if N >= One then
               raise Undefined with "Fun_Series.Series error: argument too large";
            end if;
            N := N*(NegOne/(N+NegOne));
            BallAt0(N,Y);
         end;
      end if;
      Add(1,Y);
      for I in reverse 1 .. Iter loop
         Mult(X,Y,Tmp);
         Mult(Fac(I),Y);
         Add(1,Y);
      end loop;
   end Series;

   procedure Series_Inv(X: in Fun; Y: in out Fun; Iter: in Positive) is
      S: Scalar;
   begin
      Inv_Split(X,Y,S); -- Y=S*X-1
      if IsZero(Y) then
         Add(S,Y);
         return;
      end if;
      declare
         function Fac(Dummy: Positive) return Rational is
            pragma Unreferenced(Dummy);
         begin
            return (-1)/1;
         end Fac;
         procedure InvSeries is new Series (Fac => Fac); --- Inv(1+Y)
         Z: Fun;
      begin
         Copy(Y,Z);
         InvSeries(Z,Y,False,Iter);
         Mult(S,Y);
      end;
   end Series_Inv;

   procedure Series_QPower(Q: in Rational; X: in Fun; Y: in out Fun; Iter: in Positive) is
      N: constant Integer := INum(Q);
      D: constant Integer := IDen(Q);
      S: Scalar;
   begin
      if abs(N)>D then
         raise Not_Implemented with "Fun_Series.Series_QPower: |Q|>1";
      end if;
      QPower_Split(Q,X,Y,S);  -- Y=X/(S**(1/Q))-1
      if IsZero(Y) then
         Add(S,Y);
         return;
      end if;
      declare
         function Fac(I: Positive) return Rational is
         begin
            return (N-D*(I-1))/(D*I);
         end Fac;
         procedure QPowerSeries is new Series (Fac => Fac); --- (1+Y)**Q
         Z: Fun;
      begin
         Copy(Y,Z);
         QPowerSeries(Z,Y,False,Iter);
         Mult(S,Y);
      end;
   end Series_QPower;

   procedure Series_Log(X: in Fun; Y: in out Fun; Iter: in Positive) is
      S: Scalar;
   begin
      Log_Split(X,Y,S);  -- Y=X/X0-1  S=Log(X0)
      if IsZero(Y) then
         Add(S,Y);
         return;
      end if;
      declare
         function Fac(I: Positive) return Rational is
         begin
            return (-I)/(I+1);
         end Fac;
         procedure LogSeries is new Series (Fac => Fac); --- Log(1+Y)/Y
         Z,Tmp: Fun;
      begin
         Copy(Y,Z);
         LogSeries(Z,Y,False,Iter);
         Mult(Z,Y,Tmp);
         Add(S,Y);
      end;
   end Series_Log;

   procedure Series_Exp(X: in Fun; Y: in out Fun; Iter: in Positive) is
      S: Scalar;
   begin
      Exp_Split(X,Y,S);  -- Y=X-X0  S=Exp(X0)
      if IsZero(Y) then
         Add(S,Y);
         return;
      end if;
      declare
         function Fac(I: Positive) return Rational is
         begin
            return 1/I;
         end Fac;
         procedure ExpSeries is new Series (Fac => Fac); --- Exp(Y)
         Z: Fun;
      begin
         Copy(Y,Z);
         ExpSeries(Z,Y,True,Iter);
         Mult(S,Y);
      end;
   end Series_Exp;

   procedure Series_CosSqrt(X: in Fun; Y: in out Fun; Iter: in Positive) is
      --- Cos(Sqrt(X))
      I2: Integer;
      function Fac(I: Positive) return Rational is
      begin
         I2 := 2*I;
         return 1/(I2*(1-I2));
      end Fac;
      procedure CosSqrtSeries is new Series (Fac => Fac);
   begin
      CosSqrtSeries(X,Y,True,Iter);
   end Series_CosSqrt;

   procedure Series_SinSqrt(X: in Fun; Y: in out Fun; Iter: in Positive) is
      --- Sin(Sqrt(X))/Sqrt(X)
      I2: Integer;
      function Fac(I: Positive) return Rational is
      begin
         I2 := 2*I;
         return 1/(I2*(-1-I2));
      end Fac;
      procedure SinSqrtSeries is new Series (Fac => Fac);
   begin
      SinSqrtSeries(X,Y,True,Iter);
   end Series_SinSqrt;

   procedure Series_CosSin(X: in Fun; Yc,Ys: in out Fun; Iter: in Positive) is
      C,S: Scalar;
      Z: Fun renames Yc;
   begin
      CosSin_Split(X,Z,C,S);  -- Z=X-X0  C=Cos(X0)  S=Sin(X0)
      if IsZero(Z) then
         SetZero(Ys);
         Add(C,Yc);
         Add(S,Ys);
         return;
      end if;
      declare
         Z1,Z2: Fun;
      begin
         Copy(Z,Z1);
         Prod(Z1,Z1,Z2);
         Series_CosSqrt(Z2,Yc,Iter);
         Series_SinSqrt(Z2,Ys,Iter);
         Mult(Z1,Ys,Z2); -- Z2 is Tmp
         Prod(C,Yc,Z1);
         Prod(S,Ys,Z2);
         Sub(Z2,Z1);     -- Z1 := C*Yc-S*Ys
         Prod(C,Ys,Z2);
         Prod(S,Yc,Ys);
         Add(Z2,Ys);     -- Ys := C*Ys+S*Yc
         Copy(Z1,Yc);    -- Yc := Z1
      end;
   end Series_CosSin;

   procedure Series_ArcSinSqrt(X: in Fun; Y: in out Fun; Iter: in Positive) is
      --- ArcSin(Sqrt(X))/Sqrt(X)
      I2: Integer;
      function Fac(I: Positive) return Rational is
      begin
         I2 := 2*I;
         return ((I2-1)*(I2-1))/(I2*(I2+1));
      end Fac;
      procedure ArcSinSqrtSeries is new Series (Fac => Fac);
   begin
      ArcSinSqrtSeries(X,Y,False,Iter);
   end Series_ArcSinSqrt;

   procedure Series_ArcSin(X: in Fun; Y: in out Fun; Iter: in Positive) is
      Z: Fun;
   begin
      Prod(X,X,Z);
      Series_ArcSinSqrt(Z,Y,Iter);
      Prod(X,Y,Z);
      Copy(Z,Y);
   end Series_ArcSin;

end Fun_Series;
