Fix polynomial trim(). Add polynomial division. Division probably nee…

…ds more work.

src/BBS-Numerical-polynomial_real.adb

@@ -88,6 +88,32 @@ package body BBS.Numerical.polynomial_real is
  return res;
  end;
  --
+ -- Division based on poldiv() in "Numerical Recipes in C", second
+ -- edition, 1992 by William H Press, Saul A. Tuekolsky, William T.
+ -- Vetterlink, and Brian P. Flannery, section 5.3.
+ --
+ -- u/v => q, r
+ --
+ procedure divide(u, v : poly; q : out poly; r : out poly) is
+ begin
+ for i in u'Range loop
+ r(i) := u(i);
+ end loop;
+ q := (others => 0.0);
+ --
+ for k in reverse 0 .. u'Last - v'Last loop
+ q(k) := r(v'Last + k)/v(v'Last);
+ for j in reverse k .. v'Last+k - 1 loop
+ r(j) := r(j) - q(k)*v(j-k);
+ end loop;
+ end loop;
+ for j in v'Last .. u'Last loop
+ r(j) := 0.0;
+ end loop;
+ end;
+ --
+ -- Evaluate a polynomial at x.
+ --
  function evaluate(p : poly; x : f'Base) return f'Base is
  accum : f'Base := 0.0;
  begin
@@ -97,6 +123,24 @@ package body BBS.Numerical.polynomial_real is
  return accum + p(0);
  end;
  --
+ -- Trims a polynomial by removing leading zero coefficients.
+ --
+ function trim(p : poly) return poly is
+ limit : Natural := p'Last;
+ begin
+ while (p(limit) = 0.0) and (limit > 0) loop
+ limit := limit - 1;
+ end loop;
+ declare
+ res : poly(0 .. limit);
+ begin
+ for i in 0 .. limit loop
+ res(i) := p(i);
+ end loop;
+ return res;
+ end;
+ end;
+ --
  -- Basic calculus
  --
  function integrate(p : poly; c : f'Base) return poly is

src/BBS-Numerical-polynomial_real.ads

@@ -18,13 +18,34 @@ package BBS.Numerical.polynomial_real is
  with pre => (Left'First = 0),
  post => ("*"'Result'Last = Left'Last);
  --
- function evaluate(p : poly; x : f'Base) return f'Base;
+ -- Division based on poldiv() in "Numerical Recipes in C", second
+ -- edition, 1992 by William H Press, Saul A. Tuekolsky, William T.
+ -- Vetterlink, and Brian P. Flannery, section 5.3.
+ --
+ -- u/v => q, r
+ --
+ procedure divide(u, v : poly; q : out poly; r : out poly)
+ with pre => (u'First = 0) and (v'First = 0) and
+ (q'First = 0) and (r'First = 0);
+ --
+ -- Evaluate a polynomial at x.
+ --
+ function evaluate(p : poly; x : f'Base) return f'Base
+ with pre => (p'First = 0);
+ --
+ -- Trims a polynomial by removing leading zero coefficients.
+ --
+ function trim(p : poly) return poly
+ with pre => (p'First = 0),
+ post => (trim'Result'Last <= p'Last);
  --
  -- Basic calculus
  --
  function integrate(p : poly; c : f'Base) return poly
- with post => (integrate'Result'Last = (p'Last + 1));
+ with pre => (p'First = 0),
+ post => (integrate'Result'Last = (p'Last + 1));
  function derivative(p : poly) return poly
- with post => (((derivative'Result'Last = (p'Last - 1)) and (p'Last > 0)) or
+ with pre => (p'First = 0),
+ post => (((derivative'Result'Last = (p'Last - 1)) and (p'Last > 0)) or
  ((derivative'Result'Last = 0) and (p'Last = 0)));
 end BBS.Numerical.polynomial_real;

test/test_poly.adb

@@ -29,6 +29,8 @@ procedure test_poly is
  p5 : poly.poly(0 .. 2);
  p6 : poly.poly(0 .. 3);
  p7 : poly.poly(0 .. 2);
+ p8 : poly.poly(0 .. 4);
+ p9 : poly.poly(0 .. 4);
  x : real;
  r : real;
  l : real;
@@ -61,10 +63,16 @@ begin
  Ada.Text_IO.Put(" p5 = -p2 = ");
  put_poly(p5);
  Ada.Text_IO.New_Line;
+ poly.divide(p4, p1, p8, p9);
+ Ada.Text_IO.Put(" p8 = p4/p1 = ");
+ put_poly(poly.trim(p8));
+ Ada.Text_IO.Put(", remainder p9 = ");
+ put_poly(poly.trim(p9));
+ Ada.Text_IO.New_Line;
  --
  Ada.Text_IO.Put_Line("Evaluations of polynomials");
  Ada.Text_IO.Put_Line(" x p1 p2 p3 p4 p5");
- for i in -20 .. 20 loop
+ for i in -15 .. 15 loop
  x := real(i)*0.1;
  float_io.put(x, 3, 2, 0);
  Ada.Text_IO.Put(" ");