summary refs log tree commit diff stats
path: root/lib/pure/poly.nim
diff options
context:
space:
mode:
authorlcrees <lcrees@gmail.com>2017-11-16 12:34:13 -0700
committerlcrees <lcrees@gmail.com>2017-11-16 12:40:56 -0700
commit6de04da3a54276fccac880dcec7b23bc4c98e27f (patch)
tree10ea2423afc3137ad01cedc714f687821ea728a5 /lib/pure/poly.nim
parenta4d40d137e402acb02d557baeba159b843ebda3f (diff)
downloadNim-6de04da3a54276fccac880dcec7b23bc4c98e27f.tar.gz
Removed deprecated numeric and poly module from the stdlib
Diffstat (limited to 'lib/pure/poly.nim')
-rw-r--r--lib/pure/poly.nim371
1 files changed, 0 insertions, 371 deletions
diff --git a/lib/pure/poly.nim b/lib/pure/poly.nim
deleted file mode 100644
index e286c5d17..000000000
--- a/lib/pure/poly.nim
+++ /dev/null
@@ -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)
-