#
#
#            Nim'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 
  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 polynomipx; padding-right: 5px; }
td.linenos .special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; }
span.linenos.special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; }
.highlight .hll { background-color: #ffffcc }
.highlight .c { color: #888888 } /* Comment */
.highlight .err { color: #a61717; background-color: #e3d2d2 } /* Error */
.highlight .k { color: #008800; font-weight: bold } /* Keyword */
.highlight .ch { color: #888888 } /* Comment.Hashbang */
.highlight .cm { color: #888888 } /* Comment.Multiline */
.highlight .cp { color: #cc0000; font-weight: bold } /* Comment.Preproc */
.highlight .cpf { color: #888888 } /* Comment.PreprocFile */
.highlight .c1 { color: #888888 } /* Comment.Single */
.highlight .cs { color: #cc0000; font-weight: bold; background-color: #fff0f0 } /* Comment.Special */
.highlight .gd { color: #000000; background-color: #ffdddd } /* Generic.Deleted */
.highlight .ge { font-style: italic } /* Generic.Emph */
.highlight .ges { font-weight: bold; font-style: italic } /* Generic.EmphStrong */
.highlight .gr { color: #aa0000 } /* Generic.Error */
.highlight .gh { color: #333333 } /* Generic.Heading */
.highlight .gi { color: #000000; background-color: #ddffdd } /* Generic.Inserted */
.highlight .go { color: #888888 } /* Generic.Output */
.highlight .gp { color: #555555 } /* Generic.Prompt */
.highlight .gs { font-weight: bold } /* Generic.Strong */
.highlight .gu { color: #666666 } /* Generic.Subheading */
.highlight .gt { color: #aa0000 } /* Generic.Traceback */
.highlight .kc { color: #008800; font-weight: bold } /* Keyword.Constant */
.highlight .kd { color: #008800; font-weight: bold } /* Keyword.Declaration */
.highlight .kn { color: #008800; font-weight: bold } /* Keyword.Namespace */
.highlight .kp { color: #008800 } /* Keyword.Pseudo */
.highlight .kr { color: #008800; font-weight: bold } /* Keyword.Reserved */
.highlight .kt { color: #888888; font-weight: bold } /* Keyword.Type */
.highlight .m { color: #0000DD; font-weight: bold } /* Literal.Number */
.highlight .s { color: #dd2200; background-color: #fff0f0 } /* Literal.String */
.highlight .na { color: #336699 } /* Name.Attribute */
.highlight .nb { color: #003388 } /* Name.Builtin */
.highlight .nc { color: #bb0066; font-weight: bold } /* Name.Class */
.highlight .no { color: #003366; font-weight: bold } /* Name.Constant */
.highlight .nd { color: #555555 } /* Name.Decorator */
.highlight .ne { color: #bb0066; font-weight: bold } /* Name.Exception */
.highlight .nf { color: #0066bb; font-weight: bold } /* Name.Function */
.highlight .nl { color: #336699; font-style: italic } /* Name.Label */
.highlight .nn { color: #bb0066; font-weight: bold } /* Name.Namespace */
.highlight .py { color: #336699; font-weight: bold } /* Name.Property */
.highlight .nt { color: #bb0066; font-weight: bold } /* Name.Tag */
.highlight .nv { color: #336699 } /* Name.Variable */
.highlight .ow { color: #008800 } /* Operator.Word */
.highlight .w { color: #bbbbbb } /* Text.Whitespace */
.highlight .mb { color: #0000DD; font-weight: bold } /* Literal.Number.Bin */
.highlight .mf { color: #0000DD; font-weight: bold } /* Literal.Number.Float */
.highlight .mh { color: #0000DD; font-weight: bold } /* Literal.Number.Hex */
.highlight .mi { color: #0000DD; font-weight: bold } /* Literal.Number.Integer */
.highlight .mo { color: #0000DD; font-weight: bold } /* Literal.Number.Oct */
.highlight .sa { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Affix */
.highlight .sb { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Backtick */
.highlight .sc { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Char */
.highlight .dl { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Delimiter */
.highlight .sd { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Doc */
.highlight .s2 { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Double */
.highlight .se { color: #0044dd; background-color: #fff0f0 } /* Literal.String.Escape */
.highlight .sh { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Heredoc */
.highlight .si { color: #3333bb; background-color: #fff0f0 } /* Literal.String.Interpol */
.highlight .sx { color: #22bb22; background-color: #f0fff0 } /* Literal.String.Other */
.highlight .sr { color: #008800; background-color: #fff0ff } /* Literal.String.Regex */
.highlight .s1 { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Single */
.highlight .ss { color: #aa6600; background-color: #fff0f0 } /* Literal.String.Symbol */
.highlight .bp { color: #003388 } /* Name.Builtin.Pseudo */
.highlight .fm { color: #0066bb; font-weight: bold } /* Name.Function.Magic */
.highlight .vc { color: #336699 } /* Name.Variable.Class */
.highlight .vg { color: #dd7700 } /* Name.Variable.Global */
.highlight .vi { color: #3333bb } /* Name.Variable.Instance */
.highlight .vm { color: #336699 } /* Name.Variable.Magic */
.highlight .il { color: #0000DD; font-weight: bold } /* Literal.Number.Integer.Long */</style><div class="highlight"><pre><span></span>*~
*.pyc
*.pyo
stuff/*
</pre></div>
en 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.
  ## `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)