-- nderiv.adb compute formulas for numerical derivitives
--            order is order of derivitive, 1 := first derivitive, 2 := second
--            points is number of points where value of function is known
--            f(x1), f(x2), f(x3) ... f(x points)
--            term is the point where derivitive is computer
--            f'(x1) := (1/bh)*( a(1)*f(x1) + a(2)*f(x2)
--                     + ... + a(points)*f(x points)

-- algorithm: use divided differences to get polynomial p(x) that
--            approximates  f(x).  f(x):=p(x)+error term
--            f'(x) := p'(x) + error term'
--            substitute xj := x1 + (j-1)*h
--            substitute x := x1 to get  p'(x1) etc

with Integer_Arrays; use Integer_Arrays;
with Ada.Text_IO; use Ada.Text_Io;

procedure  nderiv(Order   : Integer; -- first derivitive is 1, second is 2, etc
                  Npoints : Integer; -- number of points used to compute
                  Point   : Integer; -- compute derivitive at this point 1, 2,
                  A       : out Integer_Vector; -- integer coefficients
                  BB      : out Integer) is    -- denominator

   package Int_IO is new Integer_IO(Integer);
   use Int_IO;

  function gcd(a : Integer; b : Integer) return Integer is
    a1, b1, q, r : integer;
  begin
    if a=0 or b=0 then return 1; end if;
    if abs(a)>abs(b) then
      a1 := abs(a);
      b1 := abs(b);
    else
       a1 := abs(b);
       b1 := abs(a);
    end if;
    r:=1;
    while r/=0 loop
      q := a1 / b1;
      r := a1 - q * b1;
      a1 := b1;
      b1 := r;
    end loop;
    return a1;
  end gcd;

  -- compute the exact coefficients to numerically compute a derivitive
  -- of order 'order' using 'npoints' at ordinate point,
  -- where order>=1, npoints>=order+1, 1 <= point <= npoints,
  -- 'a' array returned with numerator of coefficients,
  -- 'bb' returned with denominator of h^order

  h : Integer_Vector(1..100);     -- coefficients of h
  x : Integer_Vector(1..100);     -- coefficients of x, var for differentiating
  numer : Integer_Matrix(1..100,1..100); -- numerator of term numer(term)(pos)
  denom : Integer_Vector(1..100); -- denominator of coefficient
  j, k, b : Integer;
  ipower, isum, iat, r : Integer;
  debug : Integer := 0;
begin
  if Npoints <= Order then
    Put_Line("ERROR in call to deriv, npoints < order+1 ");
    return;
  end if;
  for term in 1..Npoints loop
    denom(term) := 1;
    j := 1;
    for i in 1..Npoints-1 loop
       if j=term then
          j:=j+1; -- no index j in jth term
       end if;
      numer(term,i) := -(J-1); -- from (x-x(j))
      denom(term) := denom(term)*(term-j);
      j:=j+1;
    end loop;
  end loop; -- end term setting up numerator and denominator
  if debug>0 then
    for i in 1..Npoints loop
      Put("denom(");
      Put(i, 0);
      Put("):=");
      Put(denom(i), 0);
      New_Line;
      for J in 1..Npoints-1 loop
         Put("numer(");
         Put(i, 0);
         Put(",");
         Put(j, 0);
         Put("):=");
         Put(numer(i,j), 0);
         New_Line;
      end loop;
    end loop;
  end if;

  -- substitute xj := x0 + j*h, just keep j value in numer(term)(i)
  -- don't need to count x0 because same as x's
  -- done by above setup

  -- multiply out polynomial in x with coefficients of h
  -- starting from (x+j1*h)*(x+j2*h) * ...
  -- then differentiate d/dx
  for term in 1..Npoints loop
    x(1) := 1; -- initialize
    x(2) := 0;
    h(2) := 0;
    k := 2;
    for i in 1..Npoints-1 loop
      for J in 1..k loop            -- multiply by (x + j h)
        h(j) := x(j)*numer(term,i); -- multiply be constant x(j)
      end loop;
      for j in 1..k loop         -- multiply by x
        h(j+1) := h(j+1)+x(j); -- multiply by x and add
      end loop;
      k:=k+1;
      for j in 1..k loop
        x(j) := h(j);
      end loop;
      x(k) := 0;
      h(k) := 0;
    end loop;
    k := k-1;
    for j in 1..k loop
      numer(term,j) := x(j);
      if debug>0 then
         Put("p(x):= ");
         Put(numer(term,j), 0);
         Put("  x^");
         Put((j-1), 0); -- adjust power
         Put(" at term:=");
         Put(term, 0);
         New_Line;
      end if;
    end loop; -- have p(x) for this 'term'

    -- differentiate 'order' number of times
    for jorder in 1..Order loop
      for i in 2..k loop -- differentiate p'(x)
        numer(term,i-1) := (i-1)*numer(term,i);
      end loop;
      k:=k-1;
    end loop; -- have p^(order) for this term
    if debug>0 then
      for i in 1..k loop
         Put("p^(");
         Put(order, 0);
         Put(")(x):= ");
         Put(numer(term,i), 0);
         Put(" x^");
         Put(i-1, 0); -- adjust power
         Put(" at term:=");
         Put(term, 0);
         New_Line;
      end loop;
    end if;
  end loop; -- end 'term' loop  for one 'order', for one 'npoints'

  -- substitute each point for x, evaluate, get coeficients of f(x(j))
  iat := point; -- evaluate at x(iat), substitute
  for jterm in 1..Npoints loop -- get each term
    ipower := iat-1;
    isum := numer(jterm,1);
    for j in 2..k loop
      if Debug>0 then
        Put("ipower=");
        Put(ipower);
        New_Line;
      end if;
      isum := isum + ipower * numer(jterm,j);
      ipower := ipower * (Iat-1);
    end loop;
    a(jterm) := isum;
    if debug>0 then
      Put("f^(");
      Put(order, 0);
      Put(")(x(");
      Put(iat, 0);
      Put(")) := (1/h^");
      Put(order, 0);
      Put(") (");
      Put(a(jterm), 0);
      Put("/");
      Put(denom(jterm), 0);
      Put(" f(x(");
      Put(jterm, 0);
      Put_Line(")) +");
    end if;
  end loop;
  if debug>0 then
    New_Line;
  end if;
  b := 0;
  for jterm in 1..Npoints loop -- clean up each term
    r := gcd(a(jterm),denom(jterm));
    a(jterm) := a(jterm) / r;
    denom(jterm) := denom(jterm) /r;
    if denom(jterm)<0 then
      denom(jterm) := -denom(jterm);
      a(jterm) := - a(jterm);
    end if;
    if denom(jterm)>b then
       b:=denom(jterm); -- largest denominator
    end if;
  end loop;
  for jterm in 1..Npoints loop -- check for divisor
    r := (a(jterm) * b) / denom(jterm);
    if r * denom(jterm) /= a(jterm) * b then
      r := gcd(b, denom(jterm));
      b := b * (denom(jterm) / r);
    end if;
  end loop;
  for jterm in 1..Npoints loop -- check for divisor
    a(jterm) := (a(jterm) * b) / denom(jterm);
    if debug>0 then
       Put("f^(");
       Put(order, 0);
       Put(")(x(");
       Put(iat, 0);
       Put(")):=(1/");
       Put(b, 0);
       Put("h^");
       Put(order, 0);
       Put(")(");
       Put(a(jterm), 0);
       Put(" f(x(");
       Put(jterm, 0);
       Put_Line(")) +");
    end if;
  end loop; -- end computing terms of coefficients
  bb := b;
end Nderiv;