with Standard_Mathematical_Functions;    use Standard_Mathematical_Functions;

package body Chebychev_Polynomials is

  function Create ( k : natural ) return Vector is
  begin
    if k = 0
     then declare
            res : Vector(0..0);
          begin
            res(0) := 1.0;
            return res;
          end;
     elsif k = 1
         then declare
                res : Vector(0..1);
              begin
                res(0) := 0.0;
                res(1) := 1.0;
                return res;
              end;
         else declare
                res : Vector(0..k);
                tk1 : constant Vector := Create(k-1);
                tk2 : constant Vector := Create(k-2);
              begin
                for i in tk1'range loop
                  res(i+1) := 2.0*tk1(i);
                end loop;
                res(0) := 0.0;
                for i in tk2'range loop
                  res(i) := res(i) - tk2(i);
                end loop;
                return res;
              end;
    end if;
  end Create;

  function Eval ( k : natural; x : double_float ) return double_float is

    res : double_float := COS(double_float(k)*ARCCOS(x));

  begin
    return res;
  end Eval;

  function Eval ( p : Vector; x : double_float ) return double_float is

    res : double_float := p(0);
    pow : double_float := x;

  begin
    for i in 1..p'last loop
      res := res + p(i)*pow;
      pow := pow*x;
    end loop;
    return res;
  end Eval;

  function Diff ( p : Vector ) return Vector is

    res : Vector(0..p'last-1);

  begin
    for i in res'range loop
      res(i) := double_float(i+1)*p(i+1);
    end loop;
    return res;
  end Diff;

  function Diff ( p : Vector; k : natural ) return Vector is
  begin
    if k = 0
     then return p;
     elsif k = 1
         then return Diff(p);
         else return Diff(Diff(p),k-1);
	end if;
  end Diff;

  function Int ( p : Vector ) return Vector is

    res : Vector(0..p'last+1);

  begin
    res(0) := 0.0;
    for i in p'range loop
      res(i+1) := p(i)/double_float(i+1);
    end loop;
    return res;
  end Int;

  function Int ( p : Vector; k : natural ) return Vector is
  begin
    if k = 0
     then return p;
     elsif k = 1
         then return Int(p);
         else return Int(Int(p),k-1);
    end if;
  end Int;

end Chebychev_Polynomials;