diff options
Diffstat (limited to 'lib/pure/numeric.nim')
| -rw-r--r-- | lib/pure/numeric.nim | 20 |
1 files changed, 10 insertions, 10 deletions
diff --git a/lib/pure/numeric.nim b/lib/pure/numeric.nim index 9b298c0a0..71adf19b3 100644 --- a/lib/pure/numeric.nim +++ b/lib/pure/numeric.nim @@ -11,23 +11,23 @@ type OneVarFunction* = proc (x: float): float {.deprecated: [TOneVarFunction: OneVarFunction].} -proc brent*(xmin,xmax:float, function:OneVarFunction, tol:float,maxiter=1000): +proc brent*(xmin,xmax:float, function:OneVarFunction, tol:float,maxiter=1000): tuple[rootx, rooty: float, success: bool]= - ## Searches `function` for a root between `xmin` and `xmax` + ## Searches `function` for a root between `xmin` and `xmax` ## using brents method. If the function value at `xmin`and `xmax` has the ## same sign, `rootx`/`rooty` is set too the extrema value closest to x-axis ## and succes is set to false. ## Otherwise there exists at least one root and success is set to true. ## This root is searched for at most `maxiter` iterations. - ## If `tol` tolerance is reached within `maxiter` iterations + ## If `tol` tolerance is reached within `maxiter` iterations ## the root refinement stops and success=true. # see http://en.wikipedia.org/wiki/Brent%27s_method - var + var a=xmin b=xmax c=a - d=1.0e308 + d=1.0e308 fa=function(a) fb=function(b) fc=fa @@ -42,19 +42,19 @@ proc brent*(xmin,xmax:float, function:OneVarFunction, tol:float,maxiter=1000): return (a,fa,false) else: return (b,fb,false) - + if abs(fa)<abs(fb): swap(fa,fb) swap(a,b) - + while fb!=0.0 and abs(a-b)>tol: if fa!=fc and fb!=fc: # inverse quadratic interpolation s = a * fb * fc / (fa - fb) / (fa - fc) + b * fa * fc / (fb - fa) / (fb - fc) + c * fa * fb / (fc - fa) / (fc - fb) else: #secant rule s = b - fb * (b - a) / (fb - fa) tmp2 = (3.0 * a + b) / 4.0 - if not((s > tmp2 and s < b) or (s < tmp2 and s > b)) or - (mflag and abs(s - b) >= (abs(b - c) / 2.0)) or + if not((s > tmp2 and s < b) or (s < tmp2 and s > b)) or + (mflag and abs(s - b) >= (abs(b - c) / 2.0)) or (not mflag and abs(s - b) >= abs(c - d) / 2.0): s=(a+b)/2.0 mflag=true @@ -80,5 +80,5 @@ proc brent*(xmin,xmax:float, function:OneVarFunction, tol:float,maxiter=1000): inc i if i>maxiter: break - + return (b,fb,true) |