// Software    : Scilab
// Function    : polyroot()
// License     : CeCILL-B : http://www.cecill.info/licences
// Author (C)  : Samuel GOUGEON - France
// Version     : 1.0
// Date        : 2015-06-16
// Release     : http://fileexchange.scilab.org/toolboxes/362000
// History     :
//   2010-05-23 : Creation
//   2015-06-16 : First release on FileExchange
//
// DISCLAIMER: This macro is provided as is, with NO WARRANTY of any 
//   kind. Use it to your own risks, awareness, and responsability.

function RootsMult = polyroots(rp, rounded?)
// 
// Gathers raw roots of a polynomial into main ones and their multiplicities
//
// CALLING SEQUENCE
//             polyroots   // displays this help
// RootsMult = polyroots(rp, rounded?) 
//
// PARAMETERS
// rp       : Polynomial (degree N) that's roots are searched for. Or,
//            Vector or matrix of N complex roots of a polynomial 
//             of degree N (as coming from roots()).
// rounded? : Scalar boolean. if True, main roots (in RootsMult(:,1)) 
//             are rounded to the %eps^(1/N)/10 relative accuracy. 
// RootsMult: Nx2 matrix of real or complex numbers. 
//            Column #1: Average gathered roots ("main roots")
//            Column #2: Related multiplicities
//
// DESCRIPTION
// When numerically computing the N roots of a polynomial of degree N
//  with a floating point accuracy of %eps (epsilon-machine), the 
//  relative accuracy of each computed root r can't be better than 
//  %eps^(1/N), that is to say: the true root can be anywhere in a
//  disk centered on r and of radius r*%eps^(1/N).
// Hence, for a root r of multiplicity m, the m computed raw roots 
//  coming out of roots() will lay in a disk centered on r and of
//  radius r*%eps^(1/N).
// After eventually calling roots(rp) whether rp is a polynomial,
//  polyroots() looks for clusters of raw roots within the r*%eps^(1/N)
//  vicinity criterion. Then, for each identified cluster, the actual
//  main root is computed as the average of raw roots covered by the disk,
//  and the number of raw roots covered sets main root's multiplicity.
//
// WARNING: When true roots become too near from each others, they 
//  might turn numerically undistinguishable. Returned main roots 
//  may then merge several true ones with cumulated multiplicities.
//  All from roots() as well as from polyroots() are numerically correct.
//
// EXAMPLES
// // Example #1:  roots = 0, (0.5, x14), (%pi, x15)
// p = (%s - %pi)^15 * (%s - 0.5)^14 * %s;
// r = roots(p)       // raw roots from roots()
// RM = polyroots(p)  // main roots from polyroots()
//
// // Example #2:  roots = (-%i, x7), (%i, x6)
// p = (%z -%i)^7 * (%z+%i)^6;
// r = roots(p)       // raw roots from roots()
// RM = polyroots(p)  // from polyroots()
//
// // Example #3:  randomized roots
// R = unique(grand(1,8,"uin",1,9)); 
// M = grand(R,"uin",1,4);
// gsort([R' M'],"lr")  // true roots and multiplicities
// p = prod((%s - R).^M);
// r = roots(p)       // raw roots from roots()
// RM = polyroots(p)  // main roots from polyroots()
//
// // Example #4:  root = (1-%i, x6), (2-3*%i, x4)
// p = %s * (%s - 1+%i)^6 * (%s - 2+3*%i)^4;
// r = roots(p)       // raw roots from roots()
// RM = polyroots(p)  // main roots from polyroots()
//
// // Example #5:  root = (0.789, x5),  (1.23456, x6)
// p = (%s - 0.789)^5 * (%s - 1.23456)^6;
// r = roots(p)       // raw roots from roots()
// RM = polyroots(p)  // main roots from polyroots()
//
  
  fname = "polyroots"
  
  if argn(2)==0
    head_comments(fname)
    RootsMult = []
    return
  end
  
  if type(rp)==2 then
      r = roots(rp)
  else
      r = rp(:)
  end
  nri = length(r)
  rea = real(r) 
  ima = imag(r)
  dr = rea*ones(1,nri) - ones(nri,1)*rea.'
  di = ima*ones(1,nri) - ones(nri,1)*ima.'
  d  = sqrt(dr.^2 + di.^2)

  NfreeRoots = nri
  ibr = 1                 // index of next big root

  // Maximal relative radius of a group
  rrball = %eps^(1/nri)  // theoretical value
  rrball = 1.5*rrball    // empirical one
  
  m = abs(r)*ones(1,nri) // Matrix of roots moduli
  rball = rrball.*m  // rball(i,:) : maximal radius of a group around the root r(i)
  smallest = 2*number_properties("tiniest")
  k = find(rball<smallest)
  rball(k) = smallest
  
  // "Annealing" => Going from large groups to isolated roots:
  k = ones(1:nri)
  while length(k)>0 & NfreeRoots>0
    k = sum(d<rball,"r")
    if k ~= [],
      [n,col] = max(k)
      p = find(d(:,col) < rball(:,col))
      RootsMult(ibr) = mean(r(p))
      mult(ibr)  = n
      NfreeRoots = NfreeRoots - n
      ibr = ibr+1
      d(p,:) = %inf  // Tagging already processed raw roots
    end 
  end
  
// Rounding bigRoots to their box resolution
  if isdef("rounded?","local") & rounded?
    pow = round(10^(1-floor(log10(%eps^(1/nri)))))
    RootsMult = round(RootsMult * pow) / pow
  end
  RootsMult = clean(RootsMult)
  
  // Sorting by decreasing moduli
  [tmp, k] = gsort(abs(RootsMult),"g")
  RootsMult = RootsMult(k)
  RootsMult(:,2) = mult(k)
endfunction
