mirrors/Nim
1
0
Fork 0
mirror of https://github.com/nim-lang/Nim.git synced 2026-01-01 02:42:05 +00:00
Code Issues Packages Projects Releases Wiki Activity

Removed deprecated numeric and poly module from the stdlib

Browse Source
This commit is contained in:
lcrees
2017-11-16 12:34:13 -07:00
parent a4d40d137e
commit 6de04da3a5
3 changed files with 2 additions and 458 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
changelog.md
View File

@@ -82,6 +82,8 @@ This now needs to be written as:
Nimble package.
- Removed deprecated gentabs module from the stdlib and published it as separate
Nimble package.
- Removed deprecated poly and numeric modules from the stdlib and published them
as one separate Nimble package.
- The ``nim doc`` command is now an alias for ``nim doc2``, the second version of
the documentation generator. The old version 1 can still be accessed
via the new ``nim doc0`` command.

87
lib/pure/numeric.nim
View File

@@ -1,87 +0,0 @@
#
#
# Nim's Runtime Library
# (c) Copyright 2013 Robert Persson
#
# See the file "copying.txt", included in this
# distribution, for details about the copyright.
#
## **Warning:** This module will be moved out of the stdlib and into a
## Nimble package, don't use it.
type OneVarFunction* = proc (x: float): float
{.deprecated: [TOneVarFunction: OneVarFunction].}
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`
## 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
## the root refinement stops and success=true.
# 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)<abs(fb):
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
(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)<abs(fb):
swap(a,b)
swap(fa,fb)
inc i
if i>maxiter:
break
return (b,fb,true)

371
lib/pure/poly.nim
View File

@@ -1,371 +0,0 @@
#
#
# Nim's Runtime Library
# (c) Copyright 2013 Robert Persson
#
# See the file "copying.txt", included in this
# distribution, for details about the copyright.
#
## **Warning:** This module will be moved out of the stdlib and into a
## Nimble package, don't use it.
import math
import strutils
import numeric
type
Poly* = object
cofs:seq[float]
{.deprecated: [TPoly: Poly].}
proc degree*(p:Poly):int=
## Returns the degree of the polynomial,
## that is the number of coefficients-1
return p.cofs.len-1
proc eval*(p:Poly,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:Poly;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 Poly;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.
assert(idx>=0)
if idx>p.degree: #polynomial must grow
var oldlen=p.cofs.len
p.cofs.setLen(idx+1)
for q in oldlen.. <high(p.cofs):
p.cofs[q]=0.0 #new-grown coefficients set to zero
p.cofs[idx]=v
iterator items*(p:Poly):float=
## Iterates through the coefficients of the polynomial.
var i=p.degree
while i>=0:
yield p[i]
dec i
proc clean*(p:var Poly;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:Poly):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):
let 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.add(formatFloat(a,ffDefault,0))
if idx>=2:
result.add("x^" & $idx)
elif idx==1:
result.add("x")
if result=="":
result="0"
proc derivative*(p: Poly): Poly=
## 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:Poly,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:Poly):Poly=
## Returns a new polynomial which is the indefinite
## integral of `p`. The constant term is set to 0.0
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)
proc integrate*(p:Poly;xmin,xmax:float):float=
## Computes the definite integral of `p` between `xmin` and `xmax`
## 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]):Poly=
## 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:Poly;q,r:var Poly)=
## 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
r.cofs = p.cofs #initial remainder=numerator
if power<0: #denominator is larger than numerator
q.cofs= @ [0.0] #quotinent is 0.0
return # keep remainder as numerator
q.cofs=newSeq[float](power+1)
for i in countdown(pdeg,ddeg):
ratio=r.cofs[i]/d.cofs[ddeg]
q.cofs[i-ddeg]=ratio
r.cofs[i]=0.0
for j in countup(0,<ddeg):
var idx=i-ddeg+j
r.cofs[idx] = r.cofs[idx] - d.cofs[j]*ratio
r.clean # drop zero coefficients in remainder
proc `+` *(p1:Poly,p2:Poly):Poly=
## Adds two polynomials
var n=max(p1.cofs.len,p2.cofs.len)
newSeq(result.cofs,n)
for idx in countup(0,n-1):
result[idx]=p1[idx]+p2[idx]
result.clean # drop zero coefficients in remainder
proc `*` *(p1:Poly,p2:Poly):Poly=
## Multiplies the polynomial `p1` with `p2`
var
d1=p1.degree
d2=p2.degree
n=d1+d2
idx:int
newSeq(result.cofs,n)
for i1 in countup(0,d1):
for i2 in countup(0,d2):
idx=i1+i2
result[idx]=result[idx]+p1[i1]*p2[i2]
result.clean
proc `*` *(p:Poly,f:float):Poly=
## Multiplies the polynomial `p` with a real number
newSeq(result.cofs,p.cofs.len)
for i in 0..high(p.cofs):
result[i]=p.cofs[i]*f
result.clean
proc `*` *(f:float,p:Poly):Poly=
## Multiplies a real number with a polynomial
return p*f
proc `-`*(p:Poly):Poly=
## Negates a polynomial
result=p
for i in countup(0,<result.cofs.len):
result.cofs[i]= -result.cofs[i]
proc `-` *(p1:Poly,p2:Poly):Poly=
## Subtract `p1` with `p2`
var n=max(p1.cofs.len,p2.cofs.len)
newSeq(result.cofs,n)
for idx in countup(0,n-1):
result[idx]=p1[idx]-p2[idx]
result.clean # drop zero coefficients in remainder
proc `/`*(p:Poly,f:float):Poly=
## Divides polynomial `p` with a real number `f`
newSeq(result.cofs,p.cofs.len)
for i in 0..high(p.cofs):
result[i]=p.cofs[i]/f
result.clean
proc `/` *(p,q:Poly):Poly=
## Divides polynomial `p` with polynomial `q`
var dummy:Poly
p.divMod(q,result,dummy)
proc `mod` *(p,q:Poly):Poly=
## Computes the polynomial modulo operation,
## that is the remainder of `p`/`q`
var dummy:Poly
p.divMod(q,dummy,result)
proc normalize*(p:var Poly)=
## Multiplies the polynomial inplace by a term so that
## the leading term is 1.0.
## This might lead to an unstable polynomial
## if the leading term is zero.
p=p/p[p.degree]
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 0,1 or 2 solutions.
var p,q,d:float
p=b/(2.0*a)
if p==Inf or p==NegInf: #linear equation..
var linrt= -c/b
if linrt==Inf or linrt==NegInf: #constant only
return @[]
return @[linrt]
q=c/a
d=p*p-q
if d<0.0:
#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 @[]
else:
var sr=sqrt(d)
result= @[-sr-p,sr-p]
proc getRangeForRoots(p:Poly):tuple[xmin,xmax:float]=
## helper function for `roots` function
## quickly computes a range, guaranteed to contain
## all the real roots of the polynomial
# see http://www.mathsisfun.com/algebra/polynomials-bounds-zeros.html
var deg=p.degree
var d=p[deg]
var bound1,bound2:float
for i in countup(0,deg):
var c=abs(p.cofs[i]/d)
bound1=max(bound1,c+1.0)
bound2=bound2+c
bound2=max(1.0,bound2)
result.xmax=min(bound1,bound2)
result.xmin= -result.xmax
proc addRoot(p:Poly,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 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(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:Poly,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 used 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 being the same root.
## `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
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 @[] #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 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)