From 85dbc2d3065e50f7abeaee733ef309124f927688 Mon Sep 17 00:00:00 2001 From: Robert Persson Date: Tue, 2 Jul 2013 16:01:41 +0200 Subject: Added poly and numeric modules --- lib/pure/numeric.nim | 113 ++++++++++++++ lib/pure/poly.nim | 411 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 524 insertions(+) create mode 100644 lib/pure/numeric.nim create mode 100644 lib/pure/poly.nim (limited to 'lib') diff --git a/lib/pure/numeric.nim b/lib/pure/numeric.nim new file mode 100644 index 000000000..902dadea3 --- /dev/null +++ b/lib/pure/numeric.nim @@ -0,0 +1,113 @@ +# +# +# Nimrod's Runtime Library +# (c) Copyright 2013 Robert Persson +# +# See the file "copying.txt", included in this +# distribution, for details about the copyright. +# + + +type TOneVarFunction* =proc (x:float):float + +proc brent*(xmin,xmax:float ,function:TOneVarFunction, rootx,rooty:var float,tol:float,maxiter=1000):bool= + ## 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 false is returned. + ## Otherwise there exists at least one root and true is always returned. + ## This root is searched for at most `maxiter` iterations. + ## If `tol` tolerance is reached within `maxiter` iterations + ## the root refinement stops and true is returned. + + # see http://en.wikipedia.org/wiki/Brent%27s_method + var + a=xmin + b=xmax + c=a + d=1.0e308 + fa=function(a) + fb=function(b) + fc=fa + s=0.0 + fs=0.0 + mflag:bool + i=0 + tmp2:float + + if fa*fb>=0: + if abs(fa)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 + (not mflag and (abs(s - b) >= (abs(c - d) / 2.0))): + s=(a+b)/2.0 + mflag=true + else: + if ((mflag and (abs(b - c) < tol)) or (not mflag and (abs(c - d) < tol))): + s=(a+b)/2.0 + mflag=true + else: + mflag=false + fs = function(s) + d = c + c = b + fc = fb + if fa * fs<0.0: + b=s + fb=fs + else: + a=s + fa=fs + if abs(fa)maxiter: + break + + rootx=b + rooty=fb + return true + + + +when isMainModule: + var rootx=0.0 + var rooty=0.0 + var err=0.0 + + var cnt=0 + proc myf(x:float):float= + inc cnt + echo ($cnt & " : " & $x) + return x*x-200 + + + var suc=brent(-10000.0,0.2,myf,rootx,rooty,1.0e-3) + + + echo suc + echo rootx + echo rooty + + + + discard(readline(stdin)) diff --git a/lib/pure/poly.nim b/lib/pure/poly.nim new file mode 100644 index 000000000..8333398e5 --- /dev/null +++ b/lib/pure/poly.nim @@ -0,0 +1,411 @@ +# +# +# Nimrod's Runtime Library +# (c) Copyright 2013 Robert Persson +# +# See the file "copying.txt", included in this +# distribution, for details about the copyright. +# + +import math +import strutils +import numeric + + +type + TPoly* = object + cofs:seq[float] + + +proc initPolyFromDegree(n:int):TPoly= + ## internal usage only + ## caller must initialize coefficients of poly + ## and possibly `clean` away zero exponents + var numcof=n+1 #num. coefficients is one more than degree + result.cofs=newSeq[float](numcof) + +proc degree*(p:TPoly):int= + ## Returns the degree of the polynomial, + ## that is the number of coefficients-1 + return p.cofs.len-1 + + +proc eval*(p:TPoly,x:float):float= + ## Evaluates a polynomial function value for `x` + ## quickly using Horners method + var n=p.degree + result=p.cofs[n] + dec n + while n>=0: + result = result*x+p.cofs[n] + dec n + +proc `[]` *(p:TPoly;idx:int):float= + ## Gets a coefficient of the polynomial. + ## p[2] will returns the quadric term, p[3] the cubic etc. + ## Out of bounds index will return 0.0. + if idx<0 or idx>p.degree: + return 0.0 + return p.cofs[idx] + +proc `[]=` *(p:var TPoly;idx:int,v:float)= + ## Sets an coefficient of the polynomial by index. + ## p[2] set the quadric term, p[3] the cubic etc. + ## If index is out of range for the coefficients, + ## the polynomial grows to the smallest needed degree. + if idx<0: + return + + if idx>p.degree: #polynomial must grow + echo("GROW!") + var oldlen=p.cofs.len + p.cofs.setLen(idx+1) + for q in oldlen.. =0: + yield p[i] + dec i + +proc clean*(p:var TPoly;zerotol=0.0)= + ## Removes leading zero coefficients of the polynomial. + ## An optional tolerance can be given for what's considered zero. + var n=p.degree + var relen=false + + while n>0 and abs(p[n])<=zerotol: # >0 => keep at least one coefficient + dec n + relen=true + + if relen: p.cofs.setLen(n+1) + + +proc `$` *(p:TPoly):string = + ## Gets a somewhat reasonable string representation of the polynomial + ## The format should be compatible with most online function plotters, + ## for example directly in google search + result="" + var first=true #might skip + sign if first coefficient + + for idx in countdown(p.degree,0): + var a=p[idx] + + if a==0.0: + continue + + if a>= 0.0 and not first: + result.add('+') + first=false + + if a!=1.0 or idx==0: + result=result & formatFloat(a,ffDefault,0) + if idx>=2: + result.add("x^" & $idx) + elif idx==1: + result.add("x") + + if result=="": + result="0" + + +proc derivative*(p:TPoly):TPoly= + ## Returns a new polynomial, which is the derivative of `p` + newSeq[float](result.cofs,p.degree) + for idx in 0..high(result.cofs): + result.cofs[idx]=p.cofs[idx+1]*float(idx+1) + +proc diff*(p:TPoly,x:float):float= + ## Evaluates the differentiation of a polynomial with + ## respect to `x` quickly using a modifed Horners method + var n=p.degree + result=p[n]*float(n) + dec n + while n>=1: + result = result*x+p[n]*float(n) + dec n + +proc integral*(p:TPoly):TPoly= + ## Returns a new polynomial which is the indefinite + ## integral of `p`. The constant term is set to 0.0 + result=initPolyFromDegree(p.degree+1) + result.cofs[0]=0.0 #constant arbitrary term, use 0.0 + for i in 1..high(result.cofs): + result.cofs[i]=p.cofs[i-1]/float(i) + + +proc integrate*(p:TPoly;xmin,xmax:float):float= + ## Computes the definite integral of `p` between `xmin` and `xmax` + # TODO: this can be done faster using a modified horners method, + # see 'diff' function above. + var igr=p.integral + result=igr.eval(xmax)-igr.eval(xmin) + +proc initPoly*(cofs:varargs[float]):TPoly= + ## Initializes a polynomial with given coefficients. + ## The most significant coefficient is first, so to create x^2-2x+3: + ## intiPoly(1.0,-2.0,3.0) + if len(cofs)<=0: + result.cofs= @[0.0] #need at least one coefficient + else: + # reverse order of coefficients so indexing matches degree of + # coefficient... + result.cofs= @[] + for idx in countdown(cofs.len-1,0): + result.cofs.add(cofs[idx]) + + result.clean #remove leading zero terms + + +proc divMod*(p,d:TPoly;q,r:var TPoly)= + ## Divides `p` with `d`, and stores the quotinent in `q` and + ## the remainder in `d` + var + pdeg=p.degree + ddeg=d.degree + power=p.degree-d.degree + ratio:float + + if power<0: #denominator is larger than numerator + q=initPoly(0.0) #division result is 0 + r=p #remainder is numerator + return + + q=initPolyFromDegree(power) + r=p + + for i in countdown(pdeg,ddeg): + ratio=r[i]/d[ddeg] + + q[i-ddeg]=ratio + r[i]=0.0 + + for j in countup(0,=res[high(res)]+mergetol: res.add(rootx) #dont add equal roots. + else: + #this might be a 'touching' case, check function value against + #zero tolerance + if abs(rooty)<=zerotol: + if res==nil: res= @[rootx] + elif rootx>=res[high(res)]+mergetol: res.add(rootx) #dont add equal roots. + + +proc roots*(p:TPoly,tol=1.0e-9,zerotol=1.0e-6,mergetol=1.0e-12):seq[float]= + ## Computes the real roots of the polynomial `p` + ## `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 + ## 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. + ## Returns a sequence with the solutions, or nil in case of no solutions. + var deg=p.degree + var res:seq[float]=nil + if deg<=0: + return nil + elif p.degree==1: + var linrt= -p.cofs[0]/p.cofs[1] + if linrt==inf or linrt==neginf: + return nil #constant only => no roots + return @[linrt] + elif p.degree==2: + return solveQuadric(p.cofs[2],p.cofs[1],p.cofs[0],zerotol) + else: + # degree >=3 , find min/max points of polynomial with recursive + # derivative and do a numerical search for root between each min/max + var x0,x1:float + p.getRangeForRoots(x0,x1) + var minmax=p.derivative.roots(tol,zerotol,mergetol) + if minmax!=nil: #ie. we have minimas/maximas in this function + for x in minmax.items: + addRoot(p,res,x0,x,tol,zerotol,mergetol) + x0=x + addRoot(p,res,x0,x1,tol,zerotol,mergetol) + + return res + +when isMainModule: + var ply=initPoly(1.0,-6.0,5.0,2.0) + var ply2 =initPoly(4.0,5.0,6.0) + + echo ply + echo ply2 + echo ply2-ply + + + + + var rts=ply.roots + if rts!=nil: + for i in rts: + echo formatFloat(i,ffDefault,0) + + + discard readLine(stdin) + + + + + + + + + + + + + + + + + + + \ No newline at end of file -- cgit 1.4.1-2-gfad0 From d1a90c6ec6c64f3fcdf39da2e54ffd67b70339c4 Mon Sep 17 00:00:00 2001 From: Robert Persson Date: Tue, 2 Jul 2013 20:55:27 +0200 Subject: Cleanup of poly an numeric modules Removed some test code --- lib/pure/numeric.nim | 25 ------------------------- lib/pure/poly.nim | 38 -------------------------------------- 2 files changed, 63 deletions(-) (limited to 'lib') diff --git a/lib/pure/numeric.nim b/lib/pure/numeric.nim index 902dadea3..593370fbb 100644 --- a/lib/pure/numeric.nim +++ b/lib/pure/numeric.nim @@ -86,28 +86,3 @@ proc brent*(xmin,xmax:float ,function:TOneVarFunction, rootx,rooty:var float,tol rootx=b rooty=fb return true - - - -when isMainModule: - var rootx=0.0 - var rooty=0.0 - var err=0.0 - - var cnt=0 - proc myf(x:float):float= - inc cnt - echo ($cnt & " : " & $x) - return x*x-200 - - - var suc=brent(-10000.0,0.2,myf,rootx,rooty,1.0e-3) - - - echo suc - echo rootx - echo rooty - - - - discard(readline(stdin)) diff --git a/lib/pure/poly.nim b/lib/pure/poly.nim index 8333398e5..e90bc96e8 100644 --- a/lib/pure/poly.nim +++ b/lib/pure/poly.nim @@ -371,41 +371,3 @@ proc roots*(p:TPoly,tol=1.0e-9,zerotol=1.0e-6,mergetol=1.0e-12):seq[float]= addRoot(p,res,x0,x1,tol,zerotol,mergetol) return res - -when isMainModule: - var ply=initPoly(1.0,-6.0,5.0,2.0) - var ply2 =initPoly(4.0,5.0,6.0) - - echo ply - echo ply2 - echo ply2-ply - - - - - var rts=ply.roots - if rts!=nil: - for i in rts: - echo formatFloat(i,ffDefault,0) - - - discard readLine(stdin) - - - - - - - - - - - - - - - - - - - \ No newline at end of file -- cgit 1.4.1-2-gfad0 From 2cae55ae042c17c66ca1831092a3b0929d5566b7 Mon Sep 17 00:00:00 2001 From: Robert Persson Date: Wed, 3 Jul 2013 01:15:31 +0200 Subject: Fixed a mixed bag of stuff poly and numeric --- lib/pure/numeric.nim | 29 ++++++++++-------------- lib/pure/poly.nim | 64 ++++++++++++++++++++++++---------------------------- 2 files changed, 42 insertions(+), 51 deletions(-) (limited to 'lib') diff --git a/lib/pure/numeric.nim b/lib/pure/numeric.nim index 593370fbb..8ef5fabda 100644 --- a/lib/pure/numeric.nim +++ b/lib/pure/numeric.nim @@ -10,15 +10,16 @@ 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` ## 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 false is returned. - ## Otherwise there exists at least one root and true is always returned. + ## 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 - ## the root refinement stops and true is returned. + ## the root refinement stops and success=true. # see http://en.wikipedia.org/wiki/Brent%27s_method var @@ -37,13 +38,9 @@ proc brent*(xmin,xmax:float ,function:TOneVarFunction, rootx,rooty:var float,tol if fa*fb>=0: if abs(fa) 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))): + 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 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 mflag=true else: @@ -83,6 +80,4 @@ proc brent*(xmin,xmax:float ,function:TOneVarFunction, rootx,rooty:var float,tol if i>maxiter: break - rootx=b - rooty=fb - return true + return (b,fb,true) diff --git a/lib/pure/poly.nim b/lib/pure/poly.nim index e90bc96e8..cc7720565 100644 --- a/lib/pure/poly.nim +++ b/lib/pure/poly.nim @@ -53,11 +53,9 @@ proc `[]=` *(p:var TPoly;idx:int,v:float)= ## p[2] set the quadric term, p[3] the cubic etc. ## If index is out of range for the coefficients, ## the polynomial grows to the smallest needed degree. - if idx<0: - return + assert(idx>=0) if idx>p.degree: #polynomial must grow - echo("GROW!") var oldlen=p.cofs.len p.cofs.setLen(idx+1) for q in oldlen.. =0: @@ -94,7 +92,7 @@ proc `$` *(p:TPoly):string = var first=true #might skip + sign if first coefficient for idx in countdown(p.degree,0): - var a=p[idx] + let a=p[idx] if a==0.0: continue @@ -104,7 +102,7 @@ proc `$` *(p:TPoly):string = first=false if a!=1.0 or idx==0: - result=result & formatFloat(a,ffDefault,0) + result.add(formatFloat(a,ffDefault,0)) if idx>=2: result.add("x^" & $idx) elif idx==1: @@ -248,7 +246,7 @@ proc `-` *(p1:TPoly,p2:TPoly):TPoly= result.clean # drop zero coefficients in remainder 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) for i in 0..high(p.cofs): result[i]=p.cofs[i]/f @@ -260,7 +258,7 @@ proc `/` *(p,q:TPoly):TPoly= p.divMod(q,result,dummy) proc `mod` *(p,q:TPoly):TPoly= - ## computes the polynomial modulo operation, + ## Computes the polynomial modulo operation, ## that is the remainder op `p`/`q` var dummy:TPoly 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]= ## Solves the quadric equation `ax^2+bx+c`, with a possible ## tolerance `zerotol` to find roots of curves just 'touching' - ## the x axis. Returns sequence with 1 or 2 solutions, or nil - ## in case of no real solution. - result=nil + ## the x axis. Returns sequence with 0,1 or 2 solutions. 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.. var linrt= -c/b if linrt==inf or linrt==neginf: #constant only - return nil + return @[] return @[linrt] 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 var err=a*p*p-b*p+c #evaluate error at parabola center axis if(err<=zerotol): return @[-p] - return nil + return @[] else: var sr=sqrt(d) 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 ## quickly computes a range, guaranteed to contain ## all the real roots of the polynomial @@ -319,36 +315,36 @@ proc getRangeForRoots(p:TPoly;xmin,xmax:var float)= bound2=bound2+c bound2=max(1.0,bound2) - xmax=min(bound1,bound2) - xmin= -xmax + result.xmax=min(bound1,bound2) + 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 ## try to do a numeric search for a single root in range xp0-xp1, ## adding it to `res` (allocating `res` if nil) - var rootx,rooty:float - - if brent(xp0,xp1, proc(x:float):float=p.eval(x),rootx,rooty,tol): - if res==nil: res= @[rootx] - elif rootx>=res[high(res)]+mergetol: res.add(rootx) #dont add equal roots. + var br=brent(xp0,xp1, proc(x:float):float=p.eval(x),tol) + if br.success: + if res.len==0 or br.rootx>=res[high(res)]+mergetol: #dont add equal roots. + res.add(br.rootx) else: #this might be a 'touching' case, check function value against #zero tolerance - if abs(rooty)<=zerotol: - if res==nil: res= @[rootx] - elif rootx>=res[high(res)]+mergetol: res.add(rootx) #dont add equal roots. + if abs(br.rooty)<=zerotol: + if res.len==0 or br.rootx>=res[high(res)]+mergetol: #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` ## `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 ## 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. - ## 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 res:seq[float]=nil + result= @[] if deg<=0: return nil 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: # degree >=3 , find min/max points of polynomial with recursive # derivative and do a numerical search for root between each min/max - var x0,x1:float - p.getRangeForRoots(x0,x1) + var range=p.getRangeForRoots() var minmax=p.derivative.roots(tol,zerotol,mergetol) if minmax!=nil: #ie. we have minimas/maximas in this function for x in minmax.items: - addRoot(p,res,x0,x,tol,zerotol,mergetol) - x0=x - addRoot(p,res,x0,x1,tol,zerotol,mergetol) + addRoot(p,result,range.xmin,x,tol,zerotol,mergetol,maxiter) + range.xmin=x + addRoot(p,result,range.xmin,range.xmax,tol,zerotol,mergetol,maxiter) - return res + + \ No newline at end of file -- cgit 1.4.1-2-gfad0 From ce1399db2d60c085bfd97d64dc97f5c5754e60b7 Mon Sep 17 00:00:00 2001 From: Robert Persson Date: Wed, 3 Jul 2013 12:44:26 +0200 Subject: Optimized divMod in poly Made a huge speed improvement in debug mode. Only a few % in release, but the memory load should be somewhat lower. --- lib/pure/poly.nim | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) (limited to 'lib') diff --git a/lib/pure/poly.nim b/lib/pure/poly.nim index cc7720565..d5f59cfa1 100644 --- a/lib/pure/poly.nim +++ b/lib/pure/poly.nim @@ -12,6 +12,8 @@ import strutils import numeric +import times #todo: remove + type TPoly* = object cofs:seq[float] @@ -169,28 +171,25 @@ proc divMod*(p,d:TPoly;q,r:var TPoly)= power=p.degree-d.degree ratio:float + r.cofs = p.cofs #initial remainder=numerator if power<0: #denominator is larger than numerator - q=initPoly(0.0) #division result is 0 - r=p #remainder is numerator - return + q.cofs= @ [0.0] #quotinent is 0.0 + return # keep remainder as numerator - q=initPolyFromDegree(power) - r=p + q.cofs=newSeq[float](power+1) for i in countdown(pdeg,ddeg): - ratio=r[i]/d[ddeg] + ratio=r.cofs[i]/d.cofs[ddeg] - q[i-ddeg]=ratio - r[i]=0.0 + q.cofs[i-ddeg]=ratio + r.cofs[i]=0.0 for j in countup(0, Date: Sun, 7 Jul 2013 01:26:00 +0200 Subject: Optimized integrate function in module poly Rewrote the integrate function since the old one was quite hacky. The new one is about 7 times faster in release and 3 times faster in debug --- lib/pure/poly.nim | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) (limited to 'lib') diff --git a/lib/pure/poly.nim b/lib/pure/poly.nim index d5f59cfa1..609e58bdc 100644 --- a/lib/pure/poly.nim +++ b/lib/pure/poly.nim @@ -11,9 +11,6 @@ import math import strutils import numeric - -import times #todo: remove - type TPoly* = object cofs:seq[float] @@ -141,11 +138,22 @@ proc integral*(p:TPoly):TPoly= proc integrate*(p:TPoly;xmin,xmax:float):float= ## Computes the definite integral of `p` between `xmin` and `xmax` - # TODO: this can be done faster using a modified horners method, - # see 'diff' function above. - var igr=p.integral - result=igr.eval(xmax)-igr.eval(xmin) + ## quickly using a modified version of Horners method + var + n=p.degree + s1=p[n]/float(n+1) + s2=s1 + fac:float + dec n + while n>=0: + fac=p[n]/float(n+1) + s1 = s1*xmin+fac + s2 = s2*xmax+fac + dec n + + result=s2*xmax-s1*xmin + proc initPoly*(cofs:varargs[float]):TPoly= ## Initializes a polynomial with given coefficients. ## The most significant coefficient is first, so to create x^2-2x+3: @@ -258,7 +266,7 @@ proc `/` *(p,q:TPoly):TPoly= proc `mod` *(p,q:TPoly):TPoly= ## Computes the polynomial modulo operation, - ## that is the remainder op `p`/`q` + ## that is the remainder of `p`/`q` var dummy:TPoly p.divMod(q,dummy,result) @@ -363,5 +371,3 @@ proc roots*(p:TPoly,tol=1.0e-9,zerotol=1.0e-6,mergetol=1.0e-12,maxiter=1000):seq addRoot(p,result,range.xmin,x,tol,zerotol,mergetol,maxiter) range.xmin=x addRoot(p,result,range.xmin,range.xmax,tol,zerotol,mergetol,maxiter) - - \ No newline at end of file -- cgit 1.4.1-2-gfad0 From 5a1810081606cb687ec15a5e2efa69cfb74fd460 Mon Sep 17 00:00:00 2001 From: Robert Persson Date: Sun, 7 Jul 2013 01:50:10 +0200 Subject: Fixed some minor stuff in module poly Removed the stupid initPolyFromDegree which only served ro re-allocate results. Also fixed some minor stuff with nil return values in roots. --- lib/pure/poly.nim | 24 +++++++++--------------- 1 file changed, 9 insertions(+), 15 deletions(-) (limited to 'lib') diff --git a/lib/pure/poly.nim b/lib/pure/poly.nim index 609e58bdc..45e528604 100644 --- a/lib/pure/poly.nim +++ b/lib/pure/poly.nim @@ -16,13 +16,6 @@ type cofs:seq[float] -proc initPolyFromDegree(n:int):TPoly= - ## internal usage only - ## caller must initialize coefficients of poly - ## and possibly `clean` away zero exponents - var numcof=n+1 #num. coefficients is one more than degree - result.cofs=newSeq[float](numcof) - proc degree*(p:TPoly):int= ## Returns the degree of the polynomial, ## that is the number of coefficients-1 @@ -130,7 +123,7 @@ proc diff*(p:TPoly,x:float):float= proc integral*(p:TPoly):TPoly= ## Returns a new polynomial which is the indefinite ## integral of `p`. The constant term is set to 0.0 - result=initPolyFromDegree(p.degree+1) + newSeq(result.cofs,p.cofs.len+1) result.cofs[0]=0.0 #constant arbitrary term, use 0.0 for i in 1..high(result.cofs): result.cofs[i]=p.cofs[i-1]/float(i) @@ -227,7 +220,7 @@ proc `*` *(p1:TPoly,p2:TPoly):TPoly= proc `*` *(p:TPoly,f:float):TPoly= ## Multiplies the polynomial `p` with a real number - result=initPolyFromDegree(p.degree) + newSeq(result.cofs,p.cofs.len) for i in 0..high(p.cofs): result[i]=p.cofs[i]*f result.clean @@ -254,7 +247,7 @@ proc `-` *(p1:TPoly,p2:TPoly):TPoly= proc `/`*(p:TPoly,f:float):TPoly= ## Divides polynomial `p` with a real number `f` - result=initPolyFromDegree(p.degree) + newSeq(result.cofs,p.cofs.len) for i in 0..high(p.cofs): result[i]=p.cofs[i]/f result.clean @@ -351,13 +344,12 @@ proc roots*(p:TPoly,tol=1.0e-9,zerotol=1.0e-6,mergetol=1.0e-12,maxiter=1000):seq ## `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 - result= @[] - if deg<=0: - return nil - elif p.degree==1: + if deg<=0: #constant only => no roots + return @[] + elif p.degree==1: #linear var linrt= -p.cofs[0]/p.cofs[1] if linrt==inf or linrt==neginf: - return nil #constant only => no roots + return @[] #constant only => no roots return @[linrt] elif p.degree==2: return solveQuadric(p.cofs[2],p.cofs[1],p.cofs[0],zerotol) @@ -366,8 +358,10 @@ proc roots*(p:TPoly,tol=1.0e-9,zerotol=1.0e-6,mergetol=1.0e-12,maxiter=1000):seq # derivative and do a numerical search for root between each min/max var range=p.getRangeForRoots() var minmax=p.derivative.roots(tol,zerotol,mergetol) + result= @[] if minmax!=nil: #ie. we have minimas/maximas in this function for x in minmax.items: addRoot(p,result,range.xmin,x,tol,zerotol,mergetol,maxiter) range.xmin=x addRoot(p,result,range.xmin,range.xmax,tol,zerotol,mergetol,maxiter) + -- cgit 1.4.1-2-gfad0