#    USAGE:   allsolve( expr, x = a .. b)  
#  finds list of floating-point approximations to all roots of the expression
#  or equation expr in the closed finite interval [a,b]
#
#  expr, and all of its subexpressions, should be real-valued and twice 
#  continuously differentiable on [a,b].
#  In particular, you must make sure that expr does not go to infinity on
#  this interval, nor does it involve functions that do so.  
#
# Note:  This code is experimental, and is not guaranteed to work.  
# In particular, "allsolve" uses "evalr" to do interval arithmetic, and is
# therefore subject to the weaknesses of "evalr".  Among the bugs and
# weaknesses I know about:
# 1) It doesn't work with "GAMMA".
# 2) It doesn't work with the two-variable version of "arctan".
#
# It also suffers from the usual weaknesses of rounding error in floating-point 
# computations.
# To some extent this can be helped by increasing Digits.  Execution
# can be quite slow, particularly for functions that are complicated or
# change direction rapidly on the given interval. 
#
# An idea for further investigation: use the interval arithmetic
# package of Connell and Corless in the share library instead of "evalr".
#
# I would appreciate your sending comments and bug reports to
# israel@math.ubc.ca.

if kernelopts(version)='kernelopts(version)' then 
  ERROR(`This package requires at least Release 4`)
fi;

`evalf/INTERVAL`:= proc() map(evalf,INTERVAL(args)) end;
# Release 4 forgot this


allsolve := proc(eq::{equation,algebraic}, r::(name = realcons .. realcons))
local x, a, b, fn, fnp, fnpp, fa, fb, f, fl, i, n, ai, bi, c, res;
global `allsolve/eps1`, `allsolve/eps2`, `allsolve/eps3`;
options 
`Maple Advisor Database 1.00 for Maple V Release 4 and 5`,
`Copyright (c) 1997 by Robert B. Israel.  All rights reserved.`;
    Digits:= Digits+5;
    x := op(1, r);
    a := evalf(op(1, op(2, r)));
    b := evalf(op(2, op(2, r)));
    if a > b then
        ERROR(`invalid interval `, [a, b])
    fi;
    if type(eq,`=`) then f:= lhs(eq)-rhs(eq)
    else f:= eq
    fi;
    if has(f,{Heaviside,abs,signum,piecewise,min,max}) then
      fl:= convert(f,pwlist,x);
      if has(fl,{Heaviside,abs,signum,piecewise,min,max}) then
        ERROR(`could not convert to piecewise list`)
      fi;
      n:= (nops(fl)+1)/2;
      Digits:= Digits-5;
      res:= NULL;
      for i from 1 to n do
        if i=1 then ai:= a else ai:= max(a,evalf(fl[2*i-2],Digits+5)) fi;
        if i=n then bi:= b else bi:= min(b,evalf(fl[2*i],Digits+5)) fi;
        if ai > bi then next fi;
        res:= res, op(allsolve(fl[2*i-1],x=ai .. bi));
      od;
      RETURN(sort([op({res})]));
    fi;
    if has(r, infinity) then
      Digits:= Digits-5;
      if a = -infinity then
        if b = infinity then # two-sided infinite
# on (-infinity, -10] transform v = 1/x, v = -.1 .. 0
          res:= map(proc(t) if t=0 then NULL else 1/t fi end,
            { op(allsolve(subs(x=1/v, f), v = -0.1 .. 0)),
              op(allsolve(subs(x=1/v, f), v = 0 .. 0.1)) });
          res:= { op(res), op(allsolve(f, x = -10 .. 10))};
        else  # one-sided left
          c:= -10 - abs(b);
          res:= map(proc(t) if t=0 then NULL else 1/t fi end,
            { op(allsolve(subs(x=1/v, f), v = 1/c .. 0)) });
          res:= { op(res), op(allsolve(f, x = c .. b))};
        fi
      else # one-sided right
        c:= 10 + abs(a);
        res:= map(proc(t) if t=0 then NULL else 1/t fi end,
            { op(allsolve(subs(x=1/v, f), v = 0 .. 1/c )) });
        res:= { op(res), op(allsolve(f, x = a .. c))};
      fi;
      RETURN(sort([op(res)]));
    fi;
    fn := unapply(f, x);
    fa := traperror(evalf(fn(a)));
    if fa = lasterror then fa:= evalf(limit(f,x=a,right)) fi;
    if not type(fa,realcons) then
      ERROR(Limit(f,x=a,right),`does not evaluate to a real constant`)
    fi;
    fb := traperror(evalf(fn(b)));
    if fb = lasterror then fb:= evalf(limit(f,x=b,left)) fi;
    if not type(fb, realcons) then
        ERROR(Limit(f,x=a,right),`does not evaluate to a real constant`)
    fi;
    if fa < 0 then
        fa := -fa;
        fb := -fb;
        fn := unapply(-f, x);
        fnp := unapply(-normal(diff(f, x)), x);
        fnpp := unapply(-normal(diff(f, x $ 2)), x)
    else
        fnp := unapply(normal(diff(f, x)), x);
        fnpp := unapply(normal(diff(f, x $ 2)), x)
    fi;
    `allsolve/eps1`:= Float(1,5-Digits)*(abs(b)+abs(a));
    if has([fa,fb],infinity) then `allsolve/eps2`:= Float(1,5-Digits)
    else `allsolve/eps2` := Float(1, 5-Digits)*(abs(fa)+abs(fb)+.01)/2
    fi;
    if b = a then `allsolve/eps3`:= 1
    else `allsolve/eps3`:= `allsolve/eps2`/(b-a)
    fi;
    res:={`allsolve/allsolve`(fn, fnp, fnpp, a, b, fa, fb)};
    Digits:= Digits-5;
    sort([op(evalf(res))]);
end;


`allsolve/allsolve` := proc(fn, fnp, fnpp, a, b, fa, fb)
local fr, fpp1, fpp2, fp1, fp2, f1, f2, m, fm, v, w, fpa,
m1, fm1, fbb, fppb, fpb, m2, fm2, g1, q;
options 
`Maple Advisor Database 1.00 for Maple V Release 4 and 5`,
 `Copyright (c) 1997 by Robert B. Israel.  All rights reserved.`;

    userinfo(4, allsolve, `Checking interval `, a .. b, `end values`,
        fa, fb);
    if b - a < `allsolve/eps1` then
        if abs(fa + fb) < `allsolve/eps2` then RETURN(1/2*a + 1/2*b)
        else RETURN(NULL)
        fi
    fi;
    q:= fn(INTERVAL(a .. b));
    fr := traperror(evalf(readlib(evalr)(q)));
    if fr = lasterror then fr:= traperror(evalf(evalr(normal(q)))) fi;
    if fr = lasterror then 
      assume(w>=a,w<=b);
      if member(signum(0, fn(w), 0),{-1,1}) then # cut off
        RETURN(NULL)
      else fr:= INTERVAL(-infinity .. infinity)
      fi
    fi;
    if type(fr, {specfunc(range, INTERVAL), list(range)}) then
        f1 := op(1, op(1,fr)); f2 := op(2, op(-1,fr))
    elif type(fr, realcons) then f1 := fr; f2 := fr
    else ERROR(`evalr returned `, fr)
    fi;
    if has(f1,infinity) then f1:= Float(-1,Digits+10) fi;
    if has(f2,infinity) then f2:= Float(1,Digits+10) fi;
    if signum(0, f1, 0) = 1 or signum(0, f2, 0) = -1 then
        RETURN(NULL)
    fi;
    if abs(f1)+abs(f2) < .01*`allsolve/eps2` then RETURN(INTERVAL(a .. b)) fi;
    q:= fnp(INTERVAL(a .. b));
    fr := traperror(evalf(evalr(q)));
    if fr = lasterror then fr:= traperror(evalf(evalr(normal(q)))) fi;
    if fr = lasterror then 
      assume(w>=a,w<=b);
      q:= fnp(w);
      if signum(0,q,1) = 1 then fr:= INTERVAL(0 .. infinity)
      elif signum(0,q,-1) = -1 then fr:= INTERVAL(-infinity .. 0)
      else fr:= INTERVAL(-infinity .. infinity)
      fi
    fi;
    if type(fr, {specfunc(range, INTERVAL), list(range)}) then
        fp1 := op(1, op(1,fr)); fp2 := op(2, op(-1,fr))
    elif type(fr, realcons) then fp1 := fr; fp2 := fr
    else ERROR(`evalr returned `, fr)
    fi;
    if has(fp1,infinity) then fp1:= Float(-1,Digits+10) fi;
    if has(fp2,infinity) then fp2:= Float(1,Digits+10) fi;
    if fp1 >= -`allsolve/eps3` then ## nondecreasing; note fa >= 0
      if fp2 <= `allsolve/eps3` then ## function is flat
        if a = 0 then RETURN(INTERVAL(a..b))
        else RETURN(NULL)
        fi
      elif fa = 0 then RETURN(a)
      else RETURN(NULL)
      fi
    fi;
    if fp2 <= `allsolve/eps3` then ## nonincreasing with fa >= 0
        if fb <= 0 then
            RETURN(`allsolve/monotone`(fn, fnp, a, b, fa, fb))
        else RETURN(NULL)
        fi
    fi;
## now not monotone, check for convex/concave
    q:= fnpp(INTERVAL(a .. b));
    fr := traperror(evalf(evalr(q)));
    if fr = lasterror then fr:= traperror(evalf(evalr(normal(q)))) fi;
    if fr = lasterror then      
      assume(w > a, w < b);
      q:= fnpp(w);
      if   signum(0,q,1) = 1 then fr:= INTERVAL(0 .. infinity)
      elif signum(0,q,-1) = -1 then fr:= INTERVAL(-infinity .. 0)
      else  fr:= INTERVAL(-infinity .. infinity) 
      fi
    fi;
    if type(fr, {specfunc(range, INTERVAL),list(range)}) then
            fpp1 := op(1, op(1,fr)); fpp2 := op(2, op(-1,fr))
    elif type(fr, realcons) then fpp1 := fr; fpp2 := fr
    fi;
    if has(fpp1,infinity) then fpp1:= Float(-1,Digits+10) fi;
    if has(fpp2,infinity) then fpp2:= Float(1,Digits+10) fi;
    fpa := traperror(evalf(fnp(a)));
    if fpa = lasterror then fpa:= evalf(limit(fnp(v),v=a,right)) fi;
    if not type(fpa,realcons) then
      ERROR(Limit(fnp(v),v=a,right),`does not evaluate to a real constant`)
    fi;
    if fpa < 0 then
        if fpa + fpp2*(b - a) <= 0 then
            if fb <= 0 then
                RETURN(`allsolve/monotone`(fn, fnp, a, b, fa, fb))
            else RETURN(NULL)
            fi
        fi;
        m1 := a - fpa/fpp2;
        fm1 := evalf(fn(m1));
        if fm1 <= 0 then
            RETURN(`allsolve/monotone`(fn, fnp, a, m1, fa, fm1),
                `allsolve/allsolve`(-fn, -fnp, -fnpp, m1, b, -fm1,
                -fb))
        fi;
    else m1 := a; fm1 := fa
    fi;
    fpb:= traperror(evalf(fnp(b)));
    if fpb = lasterror then fpa:= evalf(limit(fnp(x),x=b,left)) fi;
    if fb < 0 then fbb := -fb; fppb := -fpp1; fpb := -fpb
    else fbb := fb; fppb := fpp2;
    fi;
    if 0 < fpb then
        if 0 <= fpb + fppb*(m1 - b) then
            if fb <= 0 then
                RETURN(
                `allsolve/monotone`(fn, fnp, m1, b, fm1, fb))
            elif fm1 < `allsolve/eps2` then RETURN(m1)
            else RETURN(NULL)
            fi
        fi;
        m2 := b - fpb/fppb;
        fm2 := evalf(fn(m2));
        if signum(0, fm2, 0) <> signum(0, fb, 2) then
            if 0 <= fm2 then
                g1 :=
                `allsolve/monotone`(fn, fnp, m2, b, fm2, fb)
            else g1 :=
                `allsolve/monotone`(-fn, -fnp, m2, b, -fm2, -fb)
            fi;
            RETURN(
                `allsolve/allsolve`(fn, fnp, fnpp, m1, m2, fm1,
                fm2), g1)
        fi;
    else m2 := b; fm2 := fb
    fi;
    if -`allsolve/eps2`/(b-a)^2 < fpp1 then
        RETURN(`allsolve/convex`(fn, fnp, m1, m2, fm1, fm2))
    elif fpp2 < `allsolve/eps2`/(b-a)^2 then
        if 0 < fm2 then 
          if fm1 < `allsolve/eps2` then RETURN(m1)
          else RETURN(NULL)
          fi
        else 
          RETURN(`allsolve/reverse`([`allsolve/convex`(
         unapply(-fn(-v),v),unapply(fnp(-v),v),-m2,-m1,-fm2,-fm1)]));
        fi
    fi;
    if `allsolve/eps2` < fm1 then
        fpa := evalf(fnp(m1));
        if 0 <= fpa then
            m1 := m1 - (fpa + sqrt(fpa^2 - 2*fm1*fpp1))/fpp1
        else m1 := m1 - 2*fm1/(fpa - sqrt(fpa^2 - 2*fm1*fpp1))
        fi;
        fm1 := evalf(fn(m1));
    fi;
    if fm2 < 0 then fbb := -fm2; fppb := -fpp2; fpb := -evalf(fnp(m2))
    else fbb := fm2; fppb := fpp1; fpb := evalf(fnp(m2))
    fi;
    if  `allsolve/eps2` < fm2 then
        if fpb <= 0 then
            m2 := m2 - (fpb - sqrt(fpb^2 - 2*fbb*fppb))/fppb
        else m2 := m2 - 2*fbb/(fpb + sqrt(fpb^2 - 2*fbb*fppb))
        fi;
        if m2 < m1 then RETURN(NULL) fi
    fi;
    m := 1/2*m1 + 1/2*m2;
    fm := evalf(fn(m));
    g1 := `allsolve/allsolve`(fn, fnp, fnpp, m1, m, fm1, fm);
    if 0 <= fm then
        g1, `allsolve/allsolve`(fn, fnp, fnpp, m, m2, fm, fm2)
    else
        g1,
        `allsolve/allsolve`(-fn, -fnp, -fnpp, m, m2, -fm, -fm2)
    fi
end;

`allsolve/reverse`:= proc(L) local n;
options
`Maple Advisor Database 1.00 for Maple V Release 4 and 5`,
 `Copyright (c) 1997 by Robert B. Israel.  All rights reserved.`;
      seq(-L[nops(L)-n], n=0 .. nops(L)-1);
      end;

`allsolve/convex` := proc(fn, fnp, a, b, fa, fb)
local fpa, fpb, v, m, fm, fmp,  ap, bp, fap, fbp;
options `Maple Advisor Database 1.00 for Maple V Release 4 and 5`,
 `Copyright (c) 1997 by Robert B. Israel.  All rights reserved.`;
    userinfo(5, allsolve, `Checking interval `, a .. b, `end values`,
        fa, fb);
    fpa := evalf(fnp(a));
    fpb := evalf(fnp(b));
    if  fpa > -`allsolve/eps3` then
        if `allsolve/eps2` < fa then RETURN(NULL)
        elif `allsolve/eps2` < fb then RETURN(a)
        else RETURN(1/2*a + 1/2*b)
        fi
    fi;
    if b  <= a + `allsolve/eps1` then
        if abs(fa + fb) < `allsolve/eps2` then RETURN(1/2*a + 1/2*b)
        else RETURN(NULL)
        fi
    fi;
    if fpb <= `allsolve/eps3` then
        if fb <= 0 then
            RETURN(`allsolve/monotone`(fn, fnp, a, b, fa, fb))
        else RETURN(NULL)
        fi
    fi;
    if fb < 0 then
        if fpa + fpb < 0 then m := b - 2*fpb*(b - a)/(fpb - fpa)
        else m := (b*fa - a*fb)/(fa - fb)
        fi;
        fm := evalf(fn(m));
        if 0 <= fm then 
             ap:= m; fap:= fm; bp:= b; fbp:= fb;
             m:= (m+b)/2;
             fm:= evalf(fn(m));
             if 0 <= fm then ap:= m; fap:= fm 
             else bp:= m; fbp:= fm;
             fi
        else 
             ap:= a; fap:= fa; bp:= m; fbp:= fm;
             fmp:= evalf(fnp(m));            
             if 0 <= fmp then
                  m:= (a+m)/2;
                  fm:= evalf(fn(m));
                  if 0 <= fm then ap:= m; fap:= fm
                  else bp:= m; fbp:= fm
                  fi
             fi
        fi; 
        RETURN(`allsolve/convex`(fn, fnp, ap, bp, fap, fbp))
    fi;
    ap := a - fa/fpa;
    bp := b - fb/fpb;
    if bp < ap then
        RETURN(NULL)
    fi;
    fap := evalf(fn(ap));
    if fap < 0 then
        ap := .9*ap + .1*a; fap := evalf(fn(ap))
    fi;
    fbp := evalf(fn(bp));
    if fbp < 0 then
        bp := .9*bp + .1*b; fbp := evalf(fn(bp))
    fi;
    fpa := evalf(fnp(ap));
    fpb := evalf(fnp(bp));
    if 0 <= fpa then
        if fap < `allsolve/eps2` then RETURN(ap) else RETURN(NULL) fi
    fi;
    if fpb <= 0 then
        if fbp < `allsolve/eps2` then RETURN(bp) else RETURN(NULL) fi
    fi;
    m := (ap*fpb - bp*fpa)/(fpb - fpa);
    if m <= ap or bp <= m then m := 1/2*ap + 1/2*bp fi;
    fm := evalf(fn(m));
    if fm < 0 then
        `allsolve/convex`(fn, fnp, ap, m, fap, fm),
        `allsolve/reverse`([`allsolve/convex`(
         unapply(fn(-v),v),unapply(-fnp(-v),v),-bp,-m,fbp,fm)]) 
    else
        fmp := evalf(fnp(m));
        if abs(fmp) < `allsolve/eps3` and
          abs(fm) < `allsolve/eps2` then
            m
        elif 0 < fmp then
            `allsolve/convex`(fn, fnp, ap, m, fap, fm)
        else
            `allsolve/convex`(fn, fnp, m, bp, fm, fbp)
        fi
    fi
end;



`allsolve/monotone` := proc(fn, fnp, a, b, fa, fb)
local d, fpm, m1, fm1, m2, fm2, m, fm;
options `Maple Advisor Database 1.00 for Maple V Release 4 and 5`,
`Copyright (c) 1997 by Robert B. Israel.  All rights reserved.`;
    userinfo(5, allsolve, `Checking interval `, a .. b, `end values`,
        fa, fb);
    if fa = fb then RETURN(1/2*a + 1/2*b) fi;
    m1 := a;
    fm1 := fa;
    m2 := b;
    fm2 := fb;
    d := m2 - m1;
    while `allsolve/eps2` < fm1 and fm2 < -`allsolve/eps2` and
          `allsolve/eps1` < d do
        m := (fm1*m2 - fm2*m1)/(fm1 - fm2);
        if m1 < m and m < m2 then
            fm := evalf(fn(m));
            if 0 < fm then m1 := m; fm1 := fm
            else m2 := m; fm2 := fm
            fi;
        fi;
        fpm := evalf(fnp(m1));
        if fpm < 0 then
            m := m1 - fm1/fpm;
            if m < m2 then
                fm := evalf(fn(m));
                if 0 < fm then m1 := m; fm1 := fm
                else m2 := m; fm2 := fm
                fi;
            fi
        fi;
        fpm := evalf(fnp(m2));
        if fpm < 0 then
            m := m2 - fm2/fpm;
            if m1 < m then
                fm := evalf(fn(m));
                if 0 < fm then m1 := m; fm1 := fm
                else m2 := m; fm2 := fm
                fi;
            fi
        fi;
        if d < 2*m2 - 2*m1 then
            m := 1/2*m1 + 1/2*m2;
            fm := evalf(fn(m));
            if 0 < fm then m1 := m; fm1 := fm
            else m2 := m; fm2 := fm
            fi;
        fi;
        d := m2 - m1
    od;
    if -fm2 < fm1 then m2 else m1 fi
end;