# USAGE: leastroot( expr, x = a .. b) # finds floating-point approximation to least root of the expression # or equation expr in the closed finite interval [a,b] # greatestroot( expr, x = a .. b) # does the same for the greatest root. # allsolve(expr, x = a .. b) # finds all the roots. # # 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". For example, # 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; # Maple forgot this leastroot:= proc(eq::{equation,algebraic}, r::(name = realcons .. realcons)) options `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved.`; allsolve(eq, r, -1); end; greatestroot:= proc(eq::{equation,algebraic}, r::(name = realcons .. realcons)) options `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved.`; allsolve(eq, r, 1); end; 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, flag; global `allsolve/eps1`, `allsolve/eps2`, `allsolve/eps3`; options `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 1997 by Robert B. Israel. All rights reserved.`; # third argument -1 for least, 1 for greatest, 0 for all # -2 or 2 for least or greatest not including 0 if nargs = 2 then flag:= 0 else flag:= args[3] fi; if not member(flag,{-2,-1,0,1,2}) then ERROR("Third argument of allsolve must be integer -2 to 2, not ",flag); fi; 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; if flag = 0 then 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,0)); od; RETURN(sort([res])); elif flag <= -1 then 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:= allsolve(fl[2*i-1],x=ai .. bi,flag); if res <> NULL then RETURN(res) fi; od; RETURN(res) else for i from n to 1 by -1 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:= allsolve(fl[2*i-1],x=ai .. bi, flag); if res <> NULL then RETURN(res) fi; od; RETURN(res) fi; fi; if has(r, infinity) then Digits:= Digits-5; fn:= proc(t) local a,b; if t = 0 then NULL elif type(t,specfunc(anything,INTERVAL)) then a:= op([1,1],t); if a = 0 then a:= infinity else a:= 1/a fi; b:= op([1,2],t); if b = 0 then b:= -infinity else b:= 1/b fi; INTERVAL(b..a) else 1/t fi end; if a = Float(-infinity) then if b = Float(infinity) then # two-sided infinite # on (-infinity, -10] transform v = 1/x, v = -.1 .. 0 # on [10, infinity) transform v = 1/x, v = 0 .. .1 if flag = 0 then res:= map(fn, { op(allsolve(subs(x=1/v, f), v = -0.1 .. 0, 0)), op(allsolve(subs(x=1/v, f), v = 0 .. 0.1, 0)) }); res:= { op(res), op(allsolve(f, x = -10 .. 10, 0))}; elif flag <= -1 then res:= allsolve(subs(x=1/v, f), v = -0.1 .. 0, 2); if res <> NULL then RETURN(fn(res)) fi; res:= allsolve(f, x = -10 .. 10, -1); if res <> NULL then RETURN(res) fi; res:= allsolve(subs(x=1/v, f), v = 0 .. 0.1, 2); if res <> NULL then RETURN(fn(res)) fi; RETURN(NULL); else res:= allsolve(subs(x=1/v,f), v = 0 .. 0.1, -2); if res <> NULL then RETURN(fn(res)) fi; res:= allsolve(f, x = -10 .. 10, 1); if res <> NULL then RETURN(res) fi; res:= allsolve(subs(x=1/v,f), v = -0.1 .. 0, -2); if res <> NULL then RETURN(fn(res)) fi; RETURN(NULL); fi; else # one-sided left c:= -10 - abs(b); if flag = 0 then res:= map(fn, { op(allsolve(subs(x=1/v, f), v = 1/c .. 0, 0)) }); res:= { op(res), op(allsolve(f, x = c .. b,0))}; elif flag <= -1 then res:= allsolve(subs(x=1/v, f), v = 1/c .. 0, 2); if res <> NULL then RETURN(fn(res)) fi; RETURN(allsolve(f, x = c .. b, -1)); else res:= allsolve(f, x = c .. b, 1); if res <> NULL then RETURN(res) fi; res:= allsolve(subs(x=1/v, f), v=1/c .. 0, -2); if res <> NULL then RETURN(fn(res)) fi; RETURN(NULL); fi fi else # one-sided right c:= 10 + abs(a); if flag = 0 then res:= map(fn, { op(allsolve(subs(x=1/v, f), v = 0 .. 1/c, 0 )) }); res:= { op(res), op(allsolve(f, x = a .. c, 0))}; elif flag <= -1 then res:= allsolve(f, x = a .. c, -1); if res <> NULL then RETURN(res) fi; res:= allsolve(subs(x=1/v, f), v = 0 .. 1/c, 2); if res <> NULL then RETURN(fn(res)) fi; RETURN(NULL) else res:= allsolve(subs(x=1/v, f), v = 0 .. 1/c, -2); if res <> NULL then RETURN(fn(res)) fi; RETURN(allsolve(f, x = a .. c, 1)); fi; fi; RETURN(sort([op(res)])); fi; fn := unapply(f, x); fa := traperror(evalf(fn(a))); if (fa = lasterror) or type(fa,undefined) 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) or type(fb,undefined) 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, flag); Digits:= Digits-5; if flag = 0 then sort([op(evalf({res}))]) else evalf(res) fi; end; `allsolve/allsolve` := proc(fn, fnp, fnpp, a, b, fa, fb, flag) local fr, fpp1, fpp2, fp1, fp2, f1, f2, m, fm, v, w, fpa, m1, fm1, fbb, fppb, fpb, m2, fm2, g1, q, res; options `Maple Advisor Database 1.01 for Maple 6`, `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(evalr(q))); if (fr = lasterror) or hastype(fr, undefined) then fr:= traperror(evalf(evalr(normal(q)))) fi; if (fr = lasterror) or hastype(fr, undefined) 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 fi; q:= fnp(INTERVAL(a .. b)); fr := traperror(evalf(evalr(q))); if (fr = lasterror) or hastype(fr,undefined) then fr:= traperror(evalf(evalr(normal(q)))) fi; if (fr = lasterror) or hastype(fr,undefined) 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 fa = 0 then RETURN(INTERVAL(a..b)) else RETURN(NULL) fi elif fa = 0 and (flag <> -2 or a <> 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, flag)) 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) or hastype(fr,undefined) then fr:= traperror(evalf(evalr(normal(q)))) fi; if (fr = lasterror) or hastype(fr,undefined) 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) or type(fpa,undefined) 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, flag)) else RETURN(NULL) fi fi; m1 := a - fpa/fpp2; fm1 := evalf(fn(m1)); if fm1 <= 0 then if flag = 0 then RETURN(`allsolve/monotone`(fn, fnp, a, m1, fa, fm1, flag), `allsolve/allsolve`(-fn, -fnp, -fnpp, m1, b, -fm1, -fb, flag)) elif flag <= -1 then res:= `allsolve/monotone`(fn, fnp, a, m1, fa, fm1, flag); if res <> NULL then RETURN(res) fi; RETURN(`allsolve/allsolve`(-fn, -fnp, -fnpp, m1, b, -fm1, -fb, flag)); else res:= `allsolve/allsolve`(-fn, -fnp, -fnpp, m1, b, -fm1, -fb, flag); if res <> NULL then RETURN(res) fi; RETURN(`allsolve/monotone`(fn, fnp, a, m1, fa, fm1, flag)); fi; fi; else m1 := a; fm1 := fa fi; fpb:= traperror(evalf(fnp(b))); if (fpb = lasterror) or type(fpb,undefined) 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,flag)) 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 flag <= -1 then g1:= `allsolve/allsolve`(fn, fnp, fnpp, m1, m2, fm1, fm2, flag); if g1 <> NULL then RETURN(g1) fi; if 0 <= fm2 then RETURN(`allsolve/monotone`(fn, fnp, m2, b, fm2, fb, flag)) else RETURN(`allsolve/monotone`(-fn, -fnp, m2, b, -fm2, -fb,flag)) fi; else if 0 <= fm2 then g1 := `allsolve/monotone`(fn, fnp, m2, b, fm2, fb, flag) else g1 := `allsolve/monotone`(-fn, -fnp, m2, b, -fm2, -fb, flag) fi; if flag >= 1 then if g1 <> NULL then RETURN(g1) else RETURN(`allsolve/allsolve`(fn, fnp, fnpp, m1, m2, fm1, fm2, flag)) fi; fi; RETURN( `allsolve/allsolve`(fn, fnp, fnpp, m1, m2, fm1, fm2, flag), g1) fi; fi; else m2 := b; fm2 := fb fi; if -`allsolve/eps2`/(b-a)^2 < fpp1 then RETURN(`allsolve/convex`(fn, fnp, m1, m2, fm1, fm2, flag)) elif fpp2 < `allsolve/eps2`/(b-a)^2 then if 0 < fm2 then if fm1 < `allsolve/eps2` then RETURN(m1) else RETURN(NULL) fi elif flag = 0 then RETURN(`allsolve/reverse`([`allsolve/convex`( unapply(-fn(-v),v),unapply(fnp(-v),v),-m2,-m1,-fm2,-fm1,0)])); else g1:= `allsolve/convex`( unapply(-fn(-v),v),unapply(fnp(-v),v),-m2,-m1,-fm2,-fm1,-flag); if g1 <> NULL then RETURN(-g1) else RETURN(NULL) fi; 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)); if flag = 0 then g1 := `allsolve/allsolve`(fn, fnp, fnpp, m1, m, fm1, fm,flag); if 0 <= fm then g1, `allsolve/allsolve`(fn, fnp, fnpp, m, m2, fm, fm2, flag) else g1, `allsolve/allsolve`(-fn, -fnp, -fnpp, m, m2, -fm, -fm2, flag) fi elif flag <= -1 then g1 := `allsolve/allsolve`(fn, fnp, fnpp, m1, m, fm1, fm,flag); if g1 <> NULL then RETURN(g1) fi; if 0 <= fm then `allsolve/allsolve`(fn, fnp, fnpp, m, m2, fm, fm2, flag) else `allsolve/allsolve`(-fn, -fnp, -fnpp, m, m2, -fm, -fm2, flag) fi else if 0 <= fm then g1:=`allsolve/allsolve`(fn, fnp, fnpp, m, m2, fm, fm2, flag) else g1:=`allsolve/allsolve`(-fn, -fnp, -fnpp, m, m2, -fm, -fm2, flag) fi; if g1 <> NULL then RETURN(g1) fi; `allsolve/allsolve`(fn, fnp, fnpp, m1, m, fm1, fm,flag); fi; end; `allsolve/reverse`:= proc(L) local n; options `Maple Advisor Database 1.01 for Maple 6`, `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, flag) local fpa, fpb, v, m, fm, fmp, ap, bp, fap, fbp, res; options `Maple Advisor Database 1.01 for Maple 6`, `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 if flag <> -2 or a <> 0 then RETURN(a) else RETURN(NULL) fi 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, flag)) 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, flag)) 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` and (flag <> -2 or ap <> 0) then RETURN(ap) else RETURN(NULL) fi fi; if fpb <= 0 then if fbp < `allsolve/eps2` and (flag <> 2 or bp <> 0) 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 if flag = 0 then `allsolve/convex`(fn, fnp, ap, m, fap, fm, 0), `allsolve/reverse`([`allsolve/convex`( unapply(fn(-v),v),unapply(-fnp(-v),v),-bp,-m,fbp,fm, 0)]) elif flag <= -1 then res:= `allsolve/convex`(fn, fnp, ap, m, fap, fm, flag); if res <> NULL then RETURN(res) fi; res:= `allsolve/convex`( unapply(fn(-v),v),unapply(-fnp(-v),v),-bp,-m,fbp,fm, -flag); if res <> NULL then -res else NULL fi else res:= `allsolve/convex`( unapply(fn(-v),v),unapply(-fnp(-v),v),-bp,-m,fbp,fm, -flag); if res <> NULL then RETURN(-res) fi; `allsolve/convex`(fn, fnp, ap, m, fap, fm, flag); fi 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, flag) else `allsolve/convex`(fn, fnp, m, bp, fm, fbp, flag) fi fi end; `allsolve/monotone` := proc(fn, fnp, a, b, fa, fb, flag) local d, fpm, m1, fm1, m2, fm2, m, fm, res; options `Maple Advisor Database 1.01 for Maple 6`, `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 res:= m2 else res:= m1 fi; if abs(flag) <> 2 or res <> 0 then res else NULL fi; end;