Some refactoring and updating comments in root finding. Add trapezoid…

… rule for integration.

README.md

@@ -11,6 +11,8 @@ is that complex numbers are a tagged record so that object notation can be used
 for some of the operations, if that's important to you.
 
 ## Integration
+Implemented algorithms are:
+* Composite trapezoid
 
 ## Differentiation
 

src/bbs-integration_real.adb

@@ -0,0 +1,19 @@
+package body BBS.integration_real is
+ --
+ function trapezoid(test : test_func; lower, upper : f; steps : Positive) return f is
+ accumulator : f := 0.0;
+ last_func : f := test(lower);
+ next_func : f;
+ base : f := lower;
+ step_size : f := (upper - lower)/f(steps);
+ begin
+ for i in 0 .. steps - 1 loop
+ base := base + step_size;
+ next_func := test(base);
+ accumulator := accumulator + ((next_func + last_func)/2.0)*step_size;
+ last_func := next_func;
+ end loop;
+ return accumulator;
+ end;
+ --
+end;

src/bbs-integration_real.ads

@@ -0,0 +1,10 @@
+generic
+ type F is digits <>;
+package BBS.integration_real is
+ type test_func is access function (x : f) return f;
+ --
+ -- Integrate the provided function between the lower and upper limits using
+ -- the composite trapezoid method with the specified number of steps.
+ --
+ function trapezoid(test : test_func; lower, upper : f; steps : Positive) return f;
+end;

src/bbs-roots_real.adb

@@ -75,47 +75,48 @@ package body BBS.roots_real is
  end;
  --
  function mueller(test : test_func; x0, x2 : in out f; limit : Positive; err : out errors) return f is
- x1 : f := (x0 + x2)/2.0;
- h1 : f := x1 - x0;
- h2 : f := x2 - x1;
- d1 : f := (test(x1) - test(x0))/h1;
- d2 : f := (test(x2) - test(x1))/h2;
- d_small : f := (d2 - d1) / (h2 + h1);
- d_big : f;
+ x1 : f := (x0 + x2)/2.0;
+ step1 : f := x1 - x0;
+ step2 : f := x2 - x1;
+ delta1 : f := (test(x1) - test(x0))/step1;
+ delta2 : f := (test(x2) - test(x1))/step2;
+ d_small : f := (delta2 - delta1) / (step2 + step1);
+ d_big : f;
+ discriminant : f;
  b : f;
- h : f;
+ value : f;
  e : f;
  root : f;
  temp : f;
  begin
  err := none;
  for i in 0 .. limit loop
- b := d2 + (h2 * d_small);
- temp := b*b - (4.0*d_small*test(x2));
- if temp < 0.0 then
+ b := delta2 + (step2 * d_small);
+ discriminant := b*b - (4.0*d_small*test(x2));
+ if discriminant < 0.0 then
  err := no_solution;
  return 0.0;
  end if;
- d_big := elem.Sqrt(temp);
+ d_big := elem.Sqrt(discriminant);
  if abs(b - d_big) < abs(b + d_big) then
  e := b + d_big;
  else
  e := b - d_big;
  end if;
- h := -2.0*test(x2)/e;
- root := x2 + h;
- if h = 0.0 then
+ value := -2.0*test(x2)/e;
+ root := x2 + value;
+ if value = 0.0 then
  return root;
  end if;
  x0 := x1;
  x1 := x2;
  x2 := root;
- h1 := x1 - x0;
- h2 := x2 - x1;
+ step1 := x1 - x0;
+ step2 := x2 - x1;
  temp := test(x1);
- d1 := (temp - test(x0))/h1;
- d2 := (test(x2) - temp)/h2;
- d_small := (d2 - d1)/(h2 + h1);
+ delta1 := (temp - test(x0))/step1;
+ delta2 := (test(x2) - temp)/step2;
+ d_small := (delta2 - delta1)/(step2 + step1);
  end loop;
  return root;
  end;

src/bbs-roots_real.ads

@@ -36,8 +36,10 @@ package BBS.roots_real is
  function seacant(test : test_func; lower, upper : in out f; limit : Positive; err : out errors) return f;
  --
  -- The Mueller algorithm uses three points to model a quadratic curve and uses
- -- that to find a candidate root. This can potentially be used to find complex
- -- roots, however this implementation does not.
+ -- that to find a candidate root. Unlike the bisection method, Mueller's
+ -- method does not require a root to be located within the three points.
+ -- This can potentially be used to find complex roots, however this
+ -- implementation does not.
  --
  -- This method will fail if the function value at the three points is equal.
  --