Updates to polynomial division.

src/BBS-Numerical-polynomial_real.adb

@@ -95,21 +95,18 @@ package body BBS.Numerical.polynomial_real is
  -- u/v => q, r
  --
  procedure divide(u, v : poly; q : out poly; r : out poly) is
+ temp : poly(u'Range) := u;
  begin
- for i in u'Range loop
- r(i) := u(i);
- end loop;
  q := (others => 0.0);
+ r := (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);
+ q(k) := temp(v'Last + k)/v(v'Last);
+ for j in k .. v'Last + k - 1 loop
+ temp(j) := temp(j) - q(k)*v(j-k);
  end loop;
  end loop;
- for j in v'Last .. u'Last loop
- r(j) := 0.0;
- end loop;
+ r(0 .. v'Last - 1) := temp(0 .. v'Last - 1);
  end;
  --
  -- Evaluate a polynomial at x.
@@ -141,6 +138,17 @@ package body BBS.Numerical.polynomial_real is
  end;
  end;
  --
+ -- Return the order of a polynomial
+ --
+ function order(p : poly) return Natural is
+ limit : Natural := p'Last;
+ begin
+ while (p(limit) = 0.0) and (limit > 0) loop
+ limit := limit - 1;
+ end loop;
+ return limit;
+ end;
+ --
  -- Basic calculus
  --
  function integrate(p : poly; c : f'Base) return poly is

src/BBS-Numerical-polynomial_real.ads

@@ -24,9 +24,15 @@ package BBS.Numerical.polynomial_real is
  --
  -- u/v => q, r
  --
+ -- The quotient and remainder may need to be trimmed to remove leading
+ -- zero coefficients.
+ --
  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);
+ (q'First = 0) and (r'First = 0) and
+ (q'Last >= (u'Last - v'Last)) and
+ (r'Last >= (v'Last - 1)) and
+ (u'Last >= v'Last);
  --
  -- Evaluate a polynomial at x.
  --
@@ -39,6 +45,11 @@ package BBS.Numerical.polynomial_real is
  with pre => (p'First = 0),
  post => (trim'Result'Last <= p'Last);
  --
+ -- Return the order of a polynomial
+ --
+ function order(p : poly) return Natural
+ with pre => (p'First = 0);
+ --
  -- Basic calculus
  --
  function integrate(p : poly; c : f'Base) return poly

test/test_poly.adb

@@ -29,8 +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);
+ p8 : poly.poly(0 .. 2);
+ p9 : poly.poly(0 .. 1);
  x : real;
  r : real;
  l : real;
@@ -40,7 +40,7 @@ procedure test_poly is
 
  function test(x : real) return real is
  begin
- return poly.evaluate(p1, x);
+ return poly.evaluate(p4, x);
  end;
 begin
  Ada.Text_IO.Put_Line("Testing some of the numerical routines.");
@@ -65,9 +65,9 @@ begin
  Ada.Text_IO.New_Line;
  poly.divide(p4, p1, p8, p9);
  Ada.Text_IO.Put(" p8 = p4/p1 = ");
- put_poly(poly.trim(p8));
+ put_poly(p8);
  Ada.Text_IO.Put(", remainder p9 = ");
- put_poly(poly.trim(p9));
+ put_poly(p9);
  Ada.Text_IO.New_Line;
  --
  Ada.Text_IO.Put_Line("Evaluations of polynomials");
@@ -88,7 +88,7 @@ begin
  Ada.Text_IO.New_Line;
  end loop;
  --
- Ada.Text_IO.Put_Line("Find a root of P1");
+ Ada.Text_IO.Put_Line("Find a root of P4");
  l := 0.0;
  u := 1.0;
  iter := 13;