File : taylors2-num.adb

with Reps.Ops, Messages;
use Reps.Ops, Messages;

package body Taylors2.Num is

  function NumInv(F: Taylor2) return Taylor2 is
--- inverse by power series
    C: Scalar;
    Ft,G: Taylor2;
  begin
    C := F.C(0,0);
    if C=ScalZero then
      Error("Taylors2.Num.NumInv error: division by 0");
    end if;
    C := Inv(C);
    Ft := (-C)*F;
    Ft.C(0,0) := ScalZero;
    G := Ft;
    G.C(0,0) := ScalOne;
    for I in 1 .. PDeg loop
      Mult(Ft,G);
      G.C(0,0) := ScalOne;
    end loop;
    return C*G;
  end NumInv;

  function GuessRho(D0,D1: Power; F: Taylor2) return RadPair is
    J,N1,N2: Integer := 0;
    Q1,Q2: Rep := One;
    R: RadPair := (9.9,9.9);
    procedure NewQuot(S0,S1: in Scalar; Q: in out Rep; N: in out Integer) is
    begin
      if (S0 /= ScalZero) and then (S1 /= ScalZero) then
        Q := Q*MaxNorm(S0)/MaxNorm(S1);
        N := N+1;
      end if;
    end NewQuot;
  begin
    if QuasiDeg(F) >= D1 then
      for D in D0 .. D1-1 loop
        for I in 0 .. D loop
          J := D-I;
          NewQuot(F.C(I,J),F.C(I+1,J),Q1,N1);
          NewQuot(F.C(I,J),F.C(I,J+1),Q2,N2);
        end loop;
      end loop;
      R(1) := Exp(Log(Q1)/Rep(N1));
      R(2) := Exp(Log(Q2)/Rep(N2));
    end if;
    return R;
  end GuessRho;

end Taylors2.Num;