-- deriv.adb  compute coefficients for numerical derivatives a(0)...
--            order is order of derivative, 1 = first derivative, 2 = second
--            npoints is number of points where value of function is known
--            f(x1), f(x2), f(x3) ... f(x npoints)
--            term is the point where derivative is computed
--             f'(xt) =  (aa(1) * f(x1) + aa(2) * f(x2)
--                        + ... + aa(points) * f(x points)) / h^order)
--                        note: caller does divide by h^order

with Ada.text_io;
use  Ada.text_io;
with Real_Arrays;
use  Real_Arrays;

procedure deriv(order   : integer;
                npoints : integer;
                point   : integer;
                aa      : out real_vector) is

  type integer_vector is array(integer range <>) of integer;
  type integer_matrix is array(integer range <>, integer range <>) of integer;
  a : integer_vector(1..50);         -- coefficients of f(x0), f(x1), ...
  h : integer_vector(1..99);           -- coefficients of h
  x : integer_vector(1..99);           -- coefficients of x,
                                       -- variable for differentiating
  numer : integer_matrix(1..99,1..99); -- numerator of a term numer(term,pos)
  denom : integer_vector(1..99);       -- denominator of coefficient
  j, k: integer;
  ipower, isum, iat: integer;
  debug : integer := 0;

begin
  if npoints<=order then
    put_line("ERROR in call to deriv, npoints=" & integer'image(npoints) &
             "< order=%d+1=" & integer'image(order));
    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; end if; -- no index j in jth term
      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_Line("denom("&Integer'Image(i)&")="&Integer'Image(denom(i)));
      for j in 1..npoints-1 loop
        Put_Line("numer("&Integer'Image(i)&","&Integer'Image(i)&")="&
                 Integer'Image(numer(i,j)));
      end loop;
    end loop;
  end if; -- debug

  -- substitute xj = x1 + j*h, just keep j value in numer(term,i)
  -- don't need to count x1 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(1) := 0;
    h(2) := 0;
    k := 1;
    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 by 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+1) := 0;
      h(k+1) := 0;
    end loop;
    for j in 1..k loop
      numer(term,j) := x(j);
      if debug>0 then
        Put_Line("p(x)= "&Integer'Image(numer(term,j))&
                 " x^"&Integer'Image(j-1)&" at term="&Integer'Image(term));

      end if; -- debug
    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_Line("p^("&Integer'Image(order)&")(x)= "&Integer'Image(numer(term,i))&
                 " x^"&Integer'Image(i-1)&" at term="&Integer'Image(term));
      end loop;
    end if; -- debug
  end loop; -- end 'term' loop  for one 'order', for one 'npoints'

  -- substitute each point for x, evaluate, get coefficients of f(x(j))
  k := npoints-order;
  if debug>0 then
    for i in 1..npoints loop
      for j in 1..npoints loop
        Put_Line("numer("&Integer'Image(i)&","&Integer'Image(j)&
                 ")="&Integer'Image(numer(i,j)));
      end loop;
    end loop;
    Put_Line("k="&Integer'Image(k));
  end if; -- debug
  iat := point-1; -- evaluate at x(iat), substitute
  for jterm in 1..npoints loop -- get each term
    ipower := iat;
    isum := numer(jterm,1);
    for j in 2..k loop
      isum := isum + ipower * numer(jterm,j);
      ipower := ipower * iat;
    end loop;
    a(jterm) := isum;
    if debug>0 then
      Put_Line("f^("&Integer'Image(order)&")(x("&Integer'Image(iat)&
               ")) = (1/h^"&Integer'Image(order)&") ("&
               Integer'Image(a(jterm))&"/"&Integer'Image(denom(jterm))&
               " f(x("&Integer'Image(jterm)&") + ");
    end if;
  end loop;
  for jterm in 1..npoints loop -- adjust for divisor
      aa(jterm) := Real(a(jterm)) / Real(denom(jterm));
  end loop; -- end computing terms of coefficients
  -- see rderiv that does the divide by h^order
end deriv;