# This procedure finds a zero of a function using bisection.  The parameters
# are an expression in an unknown defining the function, the name of the
# unknown, a range known to contain a zero, and a tolerance value for the
# accuracy to which the root is found.  The zero is found using floating
# point computation.
#
# The value of the function at one of the end-points of the range passed 
# must be less than or equal to zero, and the other must be greater than
# or equal to zero.
#
# If the expression defining the function contains a call of a procedure,
# it will probably be necessary to put it in single-quotes (') to prevent 
# immediate evaluation.
#
# If a very small tolerance value is passed as the last parameter, the
# zero will be computed to full floating-point accuracy, except if it
# very close to zero.
#
# Example:  bisection (x^2-2, x, 0..10, 1e-100)  should return sqrt(2).

bisection := proc(f,x,rng,tolerance)

    local lx, hx, mx, dx, lf, hf, mf;

    if tolerance<=0 then
        ERROR(`Fourth operand must be a positive tolerance value`);
    fi;
    
    if not type(rng,range) then 
        ERROR(`Third operand must be a range`);
    fi;

    lx := evalf(op(1,rng));
    hx := evalf(op(2,rng));
    
    if hx<lx then
        ERROR(`Low bound of range is greater than high bound`);
    fi;

    lf := evalf(eval(subs(x=lx,f)));
    hf := evalf(eval(subs(x=hx,f)));

    if lf=0 then RETURN(lx) fi;
    if hf=0 then RETURN(hx) fi;

    if lf>0 and hf>0 or lf<0 and hf<0 then
        ERROR(`Range does not straddle a zero of the function`);
    fi;

    do # loop until we find zero, or interval is smaller than tolerance

        dx := hx-lx;
        mx := lx + dx/2;

        if abs(mx-lx)<tolerance or abs(mx-hx)<tolerance then
            RETURN(mx);
        fi;

        mf := evalf(eval(subs(x=mx,f)));

        if mf=0 then
            RETURN(mx);
        fi;

        if mf>0 and lf>0 or mf<0 and lf<0 then
            lx := mx;
        else
            hx := mx;
        fi;

    od;

end: