with text_io,integer_io;                 use text_io,integer_io;
with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Standard_Floating_Numbers_io;       use Standard_Floating_Numbers_io;
with Standard_Random_Numbers;            use Standard_Random_Numbers;
with Standard_Mathematical_Functions;    use Standard_Mathematical_Functions;
with Standard_Integer_Vectors;
with Standard_Integer_Vectors_io;        use Standard_Integer_Vectors_io;
with Standard_Floating_Vectors;
with Standard_Floating_Matrices;         use Standard_Floating_Matrices;
with vLpRs_Algorithm;                    use vLpRs_Algorithm;

procedure ts_vlprs is

-- DESCRIPTION :
--   Test on the extrapolation algorithm on power series.

-- GENERATION AND EVALUATION OF A POWER SERIES :

  function Power_Series_Exponents
               ( n : natural ) return Standard_Integer_Vectors.Vector is

  -- DESCRIPTION :
  --   Asks for a vector of n exponents.

    res : Standard_Integer_Vectors.Vector(1..n) := (1..n => 0);

  begin
    put("Give "); put(n,1); put_line(" integer numbers : ");
    for i in res'range loop
      get(res(i));
    end loop;
    return res;
  end Power_Series_Exponents;

  function Power_Series_Coefficients
               ( n : natural ) return Standard_Floating_Vectors.Vector is

  -- DESCRIPTION :
  --   Returns a vector of range 0..n of randomly generated floating-point
  --   numbers, which represent the coefficients of a power series.

    res : Standard_Floating_Vectors.Vector(1..n);

  begin
    for i in res'range loop
      res(i) := Random;
    end loop;
    return res;
  end Power_Series_Coefficients;

  function Eval_Power_Series
                ( s : double_float; a : Standard_Floating_Vectors.Vector;
                  m : Standard_Integer_Vectors.Vector ) return double_float is

  -- DESCRIPTION :
  --   Returns the sum of a_i*s**m(i).

    res : double_float := 0.0;

  begin
    for i in m'range loop
      res := res + a(i)*s**m(i);
    end loop;
    return res;
  end Eval_Power_Series;

  procedure Eval_Log_Power_Series
               ( s : in double_float; a : in Standard_Floating_Vectors.Vector;
                 m : in Standard_Integer_Vectors.Vector;
                 logs,logx : out double_float ) is

  -- DESCRIPTION :
  --   On return are the logarithms of s and the result of the power series.

  begin
    logs := LOG10(s);
    logx := LOG10(Eval_Power_Series(s,a,m));
  end Eval_Log_Power_Series;

  procedure Eval_Log_Power_Series
               ( s,a : in Standard_Floating_Vectors.Vector;
                 m : in Standard_Integer_Vectors.Vector;
                 logs,logx : out Standard_Floating_Vectors.Vector ) is

  -- DESCRIPTION :
  --   On return are the logarithms of s and the result of the power series.

  begin
    for i in s'range loop
      logs(i) := LOG10(s(i));
      logx(i) := LOG10(abs(Eval_Power_Series(s(i),a,m)));
    end loop;
  end Eval_Log_Power_Series;

-- AUXILIARY :

  procedure Write ( l,v : in Standard_Floating_Vectors.Vector;
                    m : in integer ) is

    w,err : double_float;

  begin
    for i in v'range loop
      put(v(i),3,3,3); put(l(i),3,3,3); 
      w := v(i)/l(i); put(" "); put(w);
      err := abs(w - double_float(m));
      put(err,3,3,3); new_line;
    end loop;
  end Write;

-- CALLING the EXTRAPOLATOR :

  procedure Extrapolate ( n,m,r : in natural;
                          a : in Standard_Floating_Vectors.Vector;
                          exp : in Standard_Integer_Vectors.Vector ) is

    s,logs,logx : Standard_Floating_Vectors.Vector(0..m-1);
    srp,dsp : Standard_Floating_Vectors.Vector(1..r-1) := (1..r-1 => 0.0);
    p : Standard_Floating_Vectors.Vector(0..r-1) := (0..r-1 => 0.0);
    l,v : Standard_Floating_Vectors.Vector(0..r) := (0..r => 0.0);
    rt1,rt2 : Matrix(1..r-1,1..r-1);

  begin
    put("Give an decreasing sequence of "); put(m,1); 
    put_line(" floating-point numbers :");
    for i in s'range loop
      get(s(i));
    end loop;
    Eval_Log_Power_Series(s,a,exp,logs,logx);
    put_line("The logs of the s-values : ");
    for i in logs'range loop
      put(logs(i)); new_line;
    end loop;
    put_line("The logs of the x-values : ");
    for i in logx'range loop
      put(logx(i)); new_line;
    end loop;
    rt1(1,1) := 0.0; rt2(1,1) := 0.0;
    vlprs_full(r,s,logs,logx,srp,dsp,p,l,v,rt1,rt2);
    put_line("The extrapolated values : ");
    Write(l,v,exp(exp'first));
    put_line("The errors of the extrapolated values : ");
    vlprs_pipe(Standard_Output,r,s,logs,logx,srp,dsp,p,l,v,rt1,rt2);
    put_line("The extrapolated values : ");
    Write(l,v,exp(exp'first));
  end Extrapolate;

  procedure Test_vlprs is

    n,m,r : natural;
 
  begin
    put("Give the order of the extrapolator : "); get(r);
    loop
      put("Give the number of points : "); get(m);
      exit when m > r;
      put("should be larger than "); put(r,1); put_line(".  Please retry.");
    end loop;
    put("Give the length of the power series : "); get(n);
    declare
      exp : Standard_Integer_Vectors.Vector(1..n) := Power_Series_Exponents(n);
      pcf : Standard_Floating_Vectors.Vector(1..n)
          := Power_Series_Coefficients(n);
    begin
      put("the exponents : "); put(exp); new_line;
      put_line("the coefficients : ");
      for i in pcf'range loop
        put(pcf(i)); new_line;
      end loop;
      Extrapolate(n,m,r,pcf,exp);
    end;
  end Test_vlprs;

begin
  Test_vlprs;
end ts_vlprs;