jyapayne/Nim
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Fixed a mixed bag of stuff poly and numeric

Browse source
This commit is contained in:
Robert Persson 2013-07-03 01:15:31 +02:00
commit 2cae55ae04
2 changed files with 42 additions and 51 deletions
Download patch file Download diff file Expand all files Collapse all files

29
lib/pure/numeric.nim
View file

@ -10,15 +10,16 @@
type TOneVarFunction* =proc (x:float):float type TOneVarFunction* =proc (x:float):float
proc brent*(xmin,xmax:float ,function:TOneVarFunction, rootx,rooty:var float,tol:float,maxiter=1000):bool= proc brent*(xmin,xmax:float ,function:TOneVarFunction, 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 ## 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 ## same sign, `rootx`/`rooty` is set too the extrema value closest to x-axis
## and false is returned. ## and succes is set to false.
## Otherwise there exists at least one root and true is always returned. ## Otherwise there exists at least one root and success is set to true.
## This root is searched for at most `maxiter` iterations. ## 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 true is returned. ## the root refinement stops and success=true.
# see http://en.wikipedia.org/wiki/Brent%27s_method # see http://en.wikipedia.org/wiki/Brent%27s_method
var var
@ -37,13 +38,9 @@ proc brent*(xmin,xmax:float ,function:TOneVarFunction, rootx,rooty:var float,tol
if fa*fb>=0: if fa*fb>=0:
if abs(fa)<abs(fb): if abs(fa)<abs(fb):
rootx=a return (a,fa,false)
rooty=fa
return false
else: else:
rootx=b return (b,fb,false)
rooty=fb
return false
if abs(fa)<abs(fb): if abs(fa)<abs(fb):
swap(fa,fb) swap(fa,fb)
@ -55,13 +52,13 @@ proc brent*(xmin,xmax:float ,function:TOneVarFunction, rootx,rooty:var float,tol
else: #secant rule else: #secant rule
s = b - fb * (b - a) / (fb - fa) s = b - fb * (b - a) / (fb - fa)
tmp2 = (3.0 * a + b) / 4.0 tmp2 = (3.0 * a + b) / 4.0
if (not(((s > tmp2) and (s < b)) or ((s < tmp2) and (s > b)))) 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 (mflag and abs(s - b) >= (abs(b - c) / 2.0)) or
(not mflag and (abs(s - b) >= (abs(c - d) / 2.0))): (not mflag and abs(s - b) >= abs(c - d) / 2.0):
s=(a+b)/2.0 s=(a+b)/2.0
mflag=true mflag=true
else: else:
if ((mflag and (abs(b - c) < tol)) or (not mflag and (abs(c - d) < tol))): if (mflag and (abs(b - c) < tol)) or (not mflag and (abs(c - d) < tol)):
s=(a+b)/2.0 s=(a+b)/2.0
mflag=true mflag=true
else: else:
@ -83,6 +80,4 @@ proc brent*(xmin,xmax:float ,function:TOneVarFunction, rootx,rooty:var float,tol
if i>maxiter: if i>maxiter:
break break
rootx=b return (b,fb,true)
rooty=fb
return true

64
lib/pure/poly.nim
View file

@ -53,11 +53,9 @@ proc `[]=` *(p:var TPoly;idx:int,v:float)=
## p[2] set the quadric term, p[3] the cubic etc. ## p[2] set the quadric term, p[3] the cubic etc.
## If index is out of range for the coefficients, ## If index is out of range for the coefficients,
## the polynomial grows to the smallest needed degree. ## the polynomial grows to the smallest needed degree.
if idx<0: assert(idx>=0)
return
if idx>p.degree: #polynomial must grow if idx>p.degree: #polynomial must grow
echo("GROW!")
var oldlen=p.cofs.len var oldlen=p.cofs.len
p.cofs.setLen(idx+1) p.cofs.setLen(idx+1)
for q in oldlen.. <high(p.cofs): for q in oldlen.. <high(p.cofs):
@ -66,7 +64,7 @@ proc `[]=` *(p:var TPoly;idx:int,v:float)=
p.cofs[idx]=v p.cofs[idx]=v
iterator coefficients*(p:TPoly):float= iterator items*(p:TPoly):float=
## Iterates through the corfficients of the polynomial. ## Iterates through the corfficients of the polynomial.
var i=p.degree var i=p.degree
while i>=0: while i>=0:
@ -94,7 +92,7 @@ proc `$` *(p:TPoly):string =
var first=true #might skip + sign if first coefficient var first=true #might skip + sign if first coefficient
for idx in countdown(p.degree,0): for idx in countdown(p.degree,0):
var a=p[idx] let a=p[idx]
if a==0.0: if a==0.0:
continue continue
@ -104,7 +102,7 @@ proc `$` *(p:TPoly):string =
first=false first=false
if a!=1.0 or idx==0: if a!=1.0 or idx==0:
result=result & formatFloat(a,ffDefault,0) result.add(formatFloat(a,ffDefault,0))
if idx>=2: if idx>=2:
result.add("x^" & $idx) result.add("x^" & $idx)
elif idx==1: elif idx==1:
@ -248,7 +246,7 @@ proc `-` *(p1:TPoly,p2:TPoly):TPoly=
result.clean # drop zero coefficients in remainder result.clean # drop zero coefficients in remainder
proc `/`*(p:TPoly,f:float):TPoly= proc `/`*(p:TPoly,f:float):TPoly=
## Divides polynomial `p`with real number `f` ## Divides polynomial `p` with a real number `f`
result=initPolyFromDegree(p.degree) result=initPolyFromDegree(p.degree)
for i in 0..high(p.cofs): for i in 0..high(p.cofs):
result[i]=p.cofs[i]/f result[i]=p.cofs[i]/f
@ -260,7 +258,7 @@ proc `/` *(p,q:TPoly):TPoly=
p.divMod(q,result,dummy) p.divMod(q,result,dummy)
proc `mod` *(p,q:TPoly):TPoly= proc `mod` *(p,q:TPoly):TPoly=
## computes the polynomial modulo operation, ## Computes the polynomial modulo operation,
## that is the remainder op `p`/`q` ## that is the remainder op `p`/`q`
var dummy:TPoly var dummy:TPoly
p.divMod(q,dummy,result) p.divMod(q,dummy,result)
@ -277,9 +275,7 @@ proc normalize*(p:var TPoly)=
proc solveQuadric*(a,b,c:float;zerotol=0.0):seq[float]= proc solveQuadric*(a,b,c:float;zerotol=0.0):seq[float]=
## Solves the quadric equation `ax^2+bx+c`, with a possible ## Solves the quadric equation `ax^2+bx+c`, with a possible
## tolerance `zerotol` to find roots of curves just 'touching' ## tolerance `zerotol` to find roots of curves just 'touching'
## the x axis. Returns sequence with 1 or 2 solutions, or nil ## the x axis. Returns sequence with 0,1 or 2 solutions.
## in case of no real solution.
result=nil
var p,q,d:float var p,q,d:float
@ -288,7 +284,7 @@ proc solveQuadric*(a,b,c:float;zerotol=0.0):seq[float]=
if p==inf or p==neginf: #linear equation.. if p==inf or p==neginf: #linear equation..
var linrt= -c/b var linrt= -c/b
if linrt==inf or linrt==neginf: #constant only if linrt==inf or linrt==neginf: #constant only
return nil return @[]
return @[linrt] return @[linrt]
q=c/a q=c/a
@ -298,12 +294,12 @@ proc solveQuadric*(a,b,c:float;zerotol=0.0):seq[float]=
#check for inside zerotol range for neg. roots #check for inside zerotol range for neg. roots
var err=a*p*p-b*p+c #evaluate error at parabola center axis var err=a*p*p-b*p+c #evaluate error at parabola center axis
if(err<=zerotol): return @[-p] if(err<=zerotol): return @[-p]
return nil return @[]
else: else:
var sr=sqrt(d) var sr=sqrt(d)
result= @[-sr-p,sr-p] result= @[-sr-p,sr-p]
proc getRangeForRoots(p:TPoly;xmin,xmax:var float)= proc getRangeForRoots(p:TPoly):tuple[xmin,xmax:float]=
## helper function for `roots` function ## helper function for `roots` function
## quickly computes a range, guaranteed to contain ## quickly computes a range, guaranteed to contain
## all the real roots of the polynomial ## all the real roots of the polynomial
@ -319,36 +315,36 @@ proc getRangeForRoots(p:TPoly;xmin,xmax:var float)=
bound2=bound2+c bound2=bound2+c
bound2=max(1.0,bound2) bound2=max(1.0,bound2)
xmax=min(bound1,bound2) result.xmax=min(bound1,bound2)
xmin= -xmax result.xmin= -result.xmax
proc addRoot(p:TPoly,res:var seq[float],xp0,xp1,tol,zerotol,mergetol:float)= proc addRoot(p:TPoly,res:var seq[float],xp0,xp1,tol,zerotol,mergetol:float,maxiter:int)=
## helper function for `roots` function ## helper function for `roots` function
## try to do a numeric search for a single root in range xp0-xp1, ## try to do a numeric search for a single root in range xp0-xp1,
## adding it to `res` (allocating `res` if nil) ## adding it to `res` (allocating `res` if nil)
var rootx,rooty:float var br=brent(xp0,xp1, proc(x:float):float=p.eval(x),tol)
if br.success:
if brent(xp0,xp1, proc(x:float):float=p.eval(x),rootx,rooty,tol): if res.len==0 or br.rootx>=res[high(res)]+mergetol: #dont add equal roots.
if res==nil: res= @[rootx] res.add(br.rootx)
elif rootx>=res[high(res)]+mergetol: res.add(rootx) #dont add equal roots.
else: else:
#this might be a 'touching' case, check function value against #this might be a 'touching' case, check function value against
#zero tolerance #zero tolerance
if abs(rooty)<=zerotol: if abs(br.rooty)<=zerotol:
if res==nil: res= @[rootx] if res.len==0 or br.rootx>=res[high(res)]+mergetol: #dont add equal roots.
elif rootx>=res[high(res)]+mergetol: res.add(rootx) #dont add equal roots. res.add(br.rootx)
proc roots*(p:TPoly,tol=1.0e-9,zerotol=1.0e-6,mergetol=1.0e-12):seq[float]= proc roots*(p:TPoly,tol=1.0e-9,zerotol=1.0e-6,mergetol=1.0e-12,maxiter=1000):seq[float]=
## Computes the real roots of the polynomial `p` ## Computes the real roots of the polynomial `p`
## `tol` is the tolerance use to break searching for each root when reached. ## `tol` is the tolerance use to break searching for each root when reached.
## `zerotol` is the tolerance, which is 'close enough' to zero to be considered a root ## `zerotol` is the tolerance, which is 'close enough' to zero to be considered a root
## and is used to find roots for curves that only 'touch' the x-axis. ## and is used to find roots for curves that only 'touch' the x-axis.
## `mergetol` is the tolerance, of which two x-values are considered beeing the same root. ## `mergetol` is the tolerance, of which two x-values are considered beeing the same root.
## Returns a sequence with the solutions, or nil in case of no solutions. ## `maxiter` can be used to limit the number of iterations for each root.
## Returns a (possibly empty) sorted sequence with the solutions.
var deg=p.degree var deg=p.degree
var res:seq[float]=nil result= @[]
if deg<=0: if deg<=0:
return nil return nil
elif p.degree==1: elif p.degree==1:
@ -361,13 +357,13 @@ proc roots*(p:TPoly,tol=1.0e-9,zerotol=1.0e-6,mergetol=1.0e-12):seq[float]=
else: else:
# degree >=3 , find min/max points of polynomial with recursive # degree >=3 , find min/max points of polynomial with recursive
# derivative and do a numerical search for root between each min/max # derivative and do a numerical search for root between each min/max
var x0,x1:float var range=p.getRangeForRoots()
p.getRangeForRoots(x0,x1)
var minmax=p.derivative.roots(tol,zerotol,mergetol) var minmax=p.derivative.roots(tol,zerotol,mergetol)
if minmax!=nil: #ie. we have minimas/maximas in this function if minmax!=nil: #ie. we have minimas/maximas in this function
for x in minmax.items: for x in minmax.items:
addRoot(p,res,x0,x,tol,zerotol,mergetol) addRoot(p,result,range.xmin,x,tol,zerotol,mergetol,maxiter)
x0=x range.xmin=x
addRoot(p,res,x0,x1,tol,zerotol,mergetol) addRoot(p,result,range.xmin,range.xmax,tol,zerotol,mergetol,maxiter)
return res