Add bisection and secant algorithms for root finding.

.gitignore

@@ -1,5 +1,7 @@
 # Object file
 *.o
+obj/*
+lib/*
 
 # Ada Library Information
 *.ali

README.md

@@ -17,6 +17,13 @@ for some of the operations, if that's important to you.
 ## Differential Equations
 
 ## Root Finding
+These are used to find roots (zero crossings) of functions. There are a number
+of different methods with diffrent pros and cons. It is important to understand
+the nature of the function before selecting the root finding algorithm.
+
+The implemented algorithms are:
+* Bisection
+* Secant
 
 ## Vectors
 

src/complex/BBS-complex_real.adb

@@ -108,8 +108,10 @@ package body BBS.complex_real is
  -- effectively if value is within a circle centered on self.
  --
  function near_cir(self : in complex; value : complex; tol : f) return Boolean is
+ dr : f'Base := self.r - value.r;
+ di : f'Base := self.i - value.i;
  begin
- return ;
+ return (dr*dr + di*dr) < tol;
  end;
  --
  -- Additional functions

src/complex/bbs-roots_real.adb

@@ -0,0 +1,66 @@
+with Ada.Text_IO;
+with Ada.Float_Text_IO;
+package body BBS.roots_real is
+ --
+ -- Bisection algorithm for finding a root.
+ --
+ function bisection(test : test_func; lower, upper : in out f; limit : Positive; err : out errors) return f is
+ val_high : f;
+ val_low : f;
+ val_mid : f;
+ mid_limit : f;
+ begin
+ val_high := test(upper);
+ val_low := test(lower);
+ if (val_high * val_low) > 0.0 then
+ err := bad_args;
+ return 0.0;
+ end if;
+ err := none;
+ for i in 0 .. limit loop
+ mid_limit := (lower + upper)/2.0;
+ val_mid := test(mid_limit);
+ exit when val_mid = 0.0;
+ if (val_high * val_mid) > 0.0 then
+ upper := mid_limit;
+ val_high := val_mid;
+ else
+ lower := mid_limit;
+ val_low := val_mid;
+ end if;
+ end loop;
+ return (lower + upper)/2.0;
+ end;
+ --
+ function seacant(test : test_func; lower, upper : in out f; limit : Positive; err : out errors) return f is
+ val_high : f;
+ val_low : f;
+ val_mid : f;
+ mid_limit : f;
+ begin
+ val_high := test(upper);
+ val_low := test(lower);
+ if (val_high * val_low) > 0.0 then
+ err := bad_args;
+ return 0.0;
+ end if;
+ err := none;
+ for i in 0 .. limit loop
+ mid_limit := upper - val_high*(upper - lower)/(val_high - val_low);
+ val_mid := test(mid_limit);
+ exit when val_mid = 0.0;
+ if (val_high * val_mid) > 0.0 then
+ upper := mid_limit;
+ val_high := val_mid;
+ else
+ lower := mid_limit;
+ val_low := val_mid;
+ end if;
+ end loop;
+ if val_mid = 0.0 then
+ upper := mid_limit;
+ lower := mid_limit;
+ end if;
+ return upper - val_high*(upper - lower)/(val_high - val_low);
+ end;
+end;

src/complex/bbs-roots_real.ads

@@ -0,0 +1,33 @@
+generic
+ type F is digits <>;
+package BBS.roots_real is
+ type errors is (none, bad_args, no_solution);
+ type test_func is access function (x : f) return f;
+
+ --
+ -- If there is an odd number of roots between the lower and upper limits,
+ -- the bisection algorithm will always converge. Each iteration simply halves
+ -- the interval between the limits, keeping a root between them. The result
+ -- is the midpoint of the final interval after the specified number of iterations.
+ --
+ -- The lower and upper values are updated during the iterations to provide an
+ -- interval containing the root.
+ --
+ function bisection(test : test_func; lower, upper : in out f; limit : Positive; err : out errors) return f;
+ --
+ -- If there is an odd number of roots between the lower and upper limits,
+ -- the secant method will always converge. Each stage of the iteration
+ -- identifies a point where a line between the lower and upper values crosses
+ -- the axis. This point is used instead of the midpoint in the bisection
+ -- algorithm.
+ --
+ -- Depending on the function, this can converge to a root much faster than
+ -- the bisection algorithm. On the other had, it can also converge much
+ -- slower.
+ --
+ -- The lower and upper values are updated during the iterations to provide an
+ -- interval containing the root. If the root is exact, the lower and upper
+ -- values are equal to the returned value.
+ --
+ function seacant(test : test_func; lower, upper : in out f; limit : Positive; err : out errors) return f;
+end;