(******************************************************************************)

(* :Title: Roots                                                              *)
(* :Author: Jeff Olson                                                        *)
(* :Summary:
	Alternative methods for finding roots of equations. Mathematica sometimes
	has trouble with the obvious.
                                                                              *)
(* :Package Version: 1.0                                                      *)
(* :Copyright:                                                                *)
(* :Context: JOlson`Roots`                                                    *)
(* :Category: General                                                         *)
(* :History:                                                                  *) 
(* :Keywords:                                                                 *)
(* :Sources:                                                                  *)
(* :Warnings:                                                                 *)
(* :Mathematica Version: 4.0                                                  *)
(* :Limitations:                                                              *)
(* :Discussion:                                                               *)
(* :Requirements:                                                             *)
(* :Examples:                                                                 *)

(******************************************************************************)
(******************************************************************************)
BeginPackage["JOlson`Roots`"]

Through[{Unprotect, ClearAll}[
	Bisect, FindRootOnInterval, FindRoots, RootIntervals, SamplePoints
]]

Bisect::usage = "Bisect[f, {a, b}, cond] can be used to find the roots of a
	continuous function when used in conjuction with FixedPoint or NestWhile.
	If f contains a root in the interval (a, b), then Bisect will return a
	subinterval containing the first root.
	
	Syntax: f should be a pure function with a root in (a, b). cond[f[x], f[y]]
	should yield true if there is at least one root between x and y. By default
	cond evaluates to true if its arguments have the opposition sign."

FindRootOnInterval::usage = "FindRootOnInterval[f[x], x, {a, b}, opts]"

FindRoots::usage = "FindRoots[f[x], x, {a, b}, opts]"

RootIntervals::usage = "RootIntervals[f[x], x, {a, b}, opts]"

SamplePoints::usage = "SamplePoints is an option to FindRoots. Default is 25."

FindRootOnInterval::nosol = "Could not find a root on the interval `1`";

(******************************************************************************)
(******************************************************************************)
Begin["`Private`"]

Needs["Utilities`FilterOptions`"]

$DefaultCond = Composition[Negative, N, Times]

Options[RootIntervals] = 
	Options[FindRoots] = {WorkingPrecision -> Automatic, SamplePoints -> 25};

FindRootOnInterval[f_, x_, int:{_, _}, opts___?OptionQ] := 
Module[
	{cand, findOpts = FilterOptions[FindRoot, opts]}, 
	
	cand = (x /. FindRoot[f == 0, {x, Random[Real, int]}, Evaluate[findOpts]]);
	If[Head[cand] === List, cand = First[cand]];
	If[IntervalMemberQ[Interval[int], cand], cand, 
		Message[FindRootOnInterval::nosol, int]; $Failed
	]
]

RootIntervals[f_, x_, int:{a_, b_}, opts___?OptionQ] := 
Module[
	{myopts, prec, points, dx, i, x1, x2, list = {}}, 
	
	myopts = Flatten[{opts, Options[RootIntervals]}];
	prec = WorkingPrecision /. myopts;
	points = SamplePoints /. myopts;
	If[prec === Automatic, prec = $MachinePrecision];
	dx = (b - a)/points;
	Do[
		x1 = i;
		x2 = Min[i + dx, b];
		If[$DefaultCond[f /. x -> x1, f /. x -> x2],
			AppendTo[list, {x1, x2}]
		],
		{i, a, b - dx, dx}
	];
	list
]

FindRoots[f_, x_, int:{_, _}, opts___?OptionQ] := 
	(FindRootOnInterval[f, x, #, opts]& /@ RootIntervals[f, x, int, opts])

Bisect[f_, {a_, b_}, cond_:$DefaultCond] :=
With[{m = (b + a)/2},
    If[cond[f[a], f[m]], {a, m}, {m, b}]
]

End[ ]
(******************************************************************************)
(******************************************************************************)

Protect[
	Bisect, FindRootOnInterval, FindRoots, RootIntervals, SamplePoints
]

EndPackage[ ]
(******************************************************************************)
(******************************************************************************)