## Copyright (C) 2015 Doug Stewart <address@hidden>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} address@hidden =} newroots3 (@var{a})
## @deftypefnx {Function File} address@hidden =} newroots3 (@var{a}, "m")
## @deftypefnx {Function File} address@hidden =} newroots3 (@var{a}, "m", @var{Tol})
## Find the roots of Polynomials
##
##   res=newroots3([1 9 26 24])
##             the roots are -2 -3 -4
##
##   res=newroots3([1 6 12 8],'m')
##     use the multiple roots method 
##
##   res=newroots3([1 6 12 8],'m',1e-5)
##     use the multiple roots method and set the
##       tolerace for differentiating different sets of roots
##     use this Tolerance only when there are repeated roots
##       that are close to other roots
##
##
address@hidden
address@hidden
##    newroots3([1 6 12 8 ])
##  @result{} 
##  -1.99998027897029 + 0.00000000000000i
##  -2.00000986051485 + 0.00001707909463i
##  -2.00000986051485 - 0.00001707909463i
address@hidden group
address@hidden example
##
address@hidden
address@hidden
##    newroots3([1 6 12 8 ],'m')
##  @result{} 
##  -2.00000000000000
##  -2.00000000000000
##  -2.00000000000000
address@hidden group
address@hidden example
##
##
## @end deftypefn

## version 0.1 nov. 2015
function res=newroots3(a,multiple,Tol=[])
##
  if (nargin<1 || nargin>3)
    print_usage ();
  endif
  if (nargin==1)
    roots (a)
    break
  endif
  multiple
  
  if (multiple != 'm' && multiple != 'M')
      print_usage ();
  endif


## use old roots to find all roots
b=roots(a);
## get the multiplicity numbers
[mm ii]=mpoles(b,Tol);
idxcount=1;
good=true;
## find all the single roots first and factor them out
while (good==true)
  [r good ]=findm11(a); # Find a simgle root m=1
  if (good) # if there is one then save it and factor it out
    [a rr]=deconv(a,[1 -r]); # ToDo what it rr is big?
    res(idxcount)=r;
    idxcount ++;
  endif
endwhile

b=roots(a);
[mm ii]=mpoles(b,Tol);
## get the largest multiplicity number
[w iw]=max(mm);

while (w>1)  # loop till all are procesed 
  sub=b(iw-w+1:iw); # get 1 group of roots
  mea=mean(sub);
  for kk=1:w
     res(idxcount)=mea; # add them to the output list
    idxcount ++;
  endfor
  rett=convvn1([1 -mea ],w); # make a ploynomial from this group of roots
  [a rr]=deconv(a,rett); # factore them out of the original ploy.
  b=roots(a); # get a new set of roots
  [mm ii]=mpoles(b,1e-2);
  [w iw]=max(mm); # see if there are any more groups 
endwhile
b=roots(a);
for kk=1 :length(b)
  res(idxcount)=b(kk);
    idxcount ++;
  endfor
 res=sort(res(:)); 
 
  
  
  function ret=convvn1(inp,n)
    ret=inp;
    for k=1:n-1
      ret=conv(ret,inp);
    endfor
  endfunction
  
  function [ret good]=findm11(p,Tol=[])
  ret=0;
  good=false;
  b=roots(p);
  [mm ii]=mpoles(b,Tol);
  if (length(mm)>1)
    for k=1 : length(mm)-1
      if (mm(k)==1 && mm(k+1)==1)
        ret=b(ii(k));
        good=true;
      endif
    endfor
    
  elseif (length(b)==1)
    ret=b;
    good=true;
  else
    ret=0;
  endif
  endfunction

  
  
endfunction