// =============================================================================
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) 2008 - 2009 - INRIA - Michael Baudin
// Copyright (C) 2011 - DIGITEO - Michael Baudin
// Copyright (C) 2012 - INRIA - Serge Steer
//
//  This file is distributed under the same license as the Scilab package.
// =============================================================================
// <-- CLI SHELL MODE -->
function sortedRoots=sortRoots(rootsToSort)
    //Sort roots using rounded values to avoid rounding errors
    //Here 10000 is ok due to roots values
    [tmp,kRoots]=gsort(round(10000*[real(rootsToSort) imag(rootsToSort)]),"lr","i");
    sortedRoots = rootsToSort(kRoots);
endfunction
function checkroots(p,expectedroots,varargin)
     // Checks the roots function against given roots.
     //
     // 1. Check default algorithm
     myroots=roots(p);
     computedroots = sortRoots(myroots);
     expectedroots  = sortRoots(expectedroots);
     assert_checkalmostequal(computedroots,expectedroots,varargin(:));
     //
     // 2. Check "e" algorithm
     myroots=roots(p,"e");
     computedroots = sortRoots(myroots);
     expectedroots  = sortRoots(expectedroots);
     assert_checkalmostequal(computedroots,expectedroots,varargin(:));
     //
     // 3. Check "f" algorithm
     if ( isreal(p) ) then
         myroots=roots(p,"f");
         computedroots = sortRoots(myroots);
         expectedroots  = sortRoots(expectedroots);
         assert_checkalmostequal(computedroots,expectedroots,varargin(:));
     end
endfunction
//   Check the computation of the roots of a polynomial
//   with different kinds of polynomials and different 
//   kinds of roots :
//   - real poly,
//   - complex poly,
//   - real roots,
//   - complex roots.
//roots : 3 real roots
p=-6+11*%s-6*%s^2+%s^3;
expectedroots  = [1; 2; 3];
checkroots(p,expectedroots,100*%eps);
//roots : 3 real roots + polynomials algebra
p=-6+11*%s-6*%s^2+%s^3;
q = p+0;
expectedroots  = [1; 2; 3];
checkroots(q,expectedroots,100*%eps);
//roots : 3 complex roots
p=-6-%i*6+(11+%i*5)*%s+(-6-%i)*%s^2+%s^3;
expectedroots  = [1+%i; 2 ; 3];
checkroots(p,expectedroots,1d-12,20*%eps);
//roots : 3 complex roots + polynomials algebra
p=-6-%i*6+(11+%i*5)*%s+(-6-%i)*%s^2+%s^3;
q = p+0;
expectedroots  = [1+%i; 2 ; 3];
checkroots(p,expectedroots,1d-12,20*%eps);
// roots : no root at all
p=1;
v=[];
checkroots(p,[]);
q = p+0;
checkroots(q,[]);
//roots : 2 complex roots
p=1+%s+%s^2;
expectedroots  = [-0.5 - sqrt(3.)/2.*%i; -0.5 + sqrt(3.)/2.*%i ];
checkroots(p,expectedroots,10*%eps);
//roots : 2 roots equals 0
p=%s^2;
expectedroots  = [0. ; 0. ];
checkroots(p,expectedroots,%eps);
// 2 real roots with a zero derivate at the root
p=(%s-%pi)^2;
expectedroots  = [%pi;%pi];
checkroots(p,expectedroots,10*%eps);
//
// Caution !
// The following are difficult root-finding problems
// with expected precision problems.
// See "Principles for testing polynomial 
// zerofinding programs"
// Jenkins, Traub
// 1975
// p.28
// "The accuracy which one may expect to achieve in calculating
// zeros is limited by the condition of these zeros. In particular,
// for multiple zeros perturbations of size epsilon in the 
// coefficients cause perturbations of size epsilon^(1/m)
// in the zeros."
// 
//
// 3 real roots with a zero derivate at the root
// Really difficult problem : only simple precision computed, instead of double precision ***
p=(%s-%pi)^3;
expectedroots  = [%pi;%pi;%pi];
checkroots(p,expectedroots,%eps^(1/3),5*%eps^(1/3));
// 4 real roots with a zero derivate at the root
// Really difficult problem : only simple precision
p=(%s-%pi)^4;
expectedroots  = [%pi;%pi;%pi;%pi];
checkroots(p,expectedroots,%eps^(1/4),5*%eps^(1/4))
// 10 real roots with a zero derivate at the root
// Really difficult problem : only one correct digit
p=(%s-%pi)^10;
expectedroots  = [%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi];
checkroots(p,expectedroots,%eps^(1/10),8*%eps^(1/10))
// "Numerical computing with Matlab", Cleve Moler.
A = diag(1:20);
p = poly(A,'x');
e = [1:20]';
checkroots(p,e,%eps,0.2);
// Tests from CPOLY
// M. A. Jenkins and J. F. Traub. 1972.
// Algorithm 419: zeros of a complex polynomial.
// Commun. ACM 15, 2 (February 1972), 97-99.
//
// EXAMPLE 1.  POLYNOMIAL WITH ZEROS 1,2,...,10.
P=[];
PI=[];
P(1)=1;
P(2)=-55;
P(3)=1320;
P(4)=-18150;
P(5)=157773;
P(6)=-902055;
P(7) = 3416930;
P(8)=-8409500;
P(9)=12753576;
P(10)=-10628640;
P(11)=3628800;
PI(1:11) = 0;
P = complex(P,PI);
E = (1:10)';
R = roots(P);
E = sortRoots(E);
R = sortRoots(R);
assert_checkalmostequal(R, E, 1.e-10);
// EXAMPLE 2. ZEROS ON IMAGINARY AXIS DEGREE 3.
// x^3-10001.0001*i*x^2-10001.0001*x+i
P = [];
PI=[];
P(1)=1;
P(2)=0;
P(3)=-10001.0001;
P(4)=0;
PI(1)=0;
PI(2)=-10001.0001;
PI(3)=0;
PI(4)=1;
P = complex(P,PI);
E = [
0.0001*%i
%i
10000*%i
];
R = roots(P);
E = sortRoots(E);
R = sortRoots(R);
assert_checkalmostequal(R, E, 1.e-15, 1.e-10);
// plot(real(R),imag(R),"bo")
// xtitle("Roots","Real","Imaginary")
// EXAMPLE 3. ZEROS AT 1+I,1/2*(1+I)....1/(2**-9)*(1+I)
P = [];
PI=[];
P(1)=1.0;
P(2)=-1.998046875;
P(3)=0.0;
P(4)=.7567065954208374D0;
P(5)=-.2002119533717632D0;
P(6)=1.271507365163416D-2;
P(7)=0;
P(8)=-1.154642632172909D-5;
P(9)=1.584803612786345D-7;
P(10)=-4.652065399568528D-10;
P(11)=0;
PI(1)=0;
PI(2)=P(2);
PI(3)=2.658859252929688D0;
PI(4)=-7.567065954208374D-1;
PI(5)=0;
PI(6)=P(6);
PI(7)=-7.820779428584501D-4;
PI(8)=-P(8);
PI(9)=0;
PI(10)=P(10);
PI(11)=9.094947017729282D-13;
P = complex(P,PI);
R = roots(P);
E = (1+%i)*2.^((0:-1:-9)');
E = sortRoots(E);
R = sortRoots(R);
assert_checkalmostequal(R, E, 1.e-13, 1.e-14);
// EXAMPLE 4. MULTIPLE ZEROS
// Real part:
// 288 - 1344*x + 2204*x^2 - 920*x^3 - 1587*x^4 + 2374*x^5 - 1293*x^6 + 284*x^7 + 3*x^8 - 10*x^9 + x^10
// Imaginary part:
// 504*x - 2352*x^2 + 4334*x^3 - 3836*x^4 + 1394*x^5 + 200*x^6 - 334*x^7 + 100*x^8 - 10*x^9
P = [];
PI=[];
P(1)=1;
P(2)=-10;
P(3)=3;
P(4)=284;
P(5)=-1293;
P(6)=2374;
P(7)=-1587;
P(8)=-920;
P(9)=2204;
P(10)=-1344;
P(11)=288;
PI(1)=0;
PI(2)=-10;
PI(3)=100;
PI(4)=-334;
PI(5)=200;
PI(6)=1394;
PI(7) =-3836;
PI(8)=4334;
PI(9)=-2352;
PI(10)=504;
PI(11)=0;
P = complex(P,PI);
R = roots(P);
E = [
1
1
1
1
2*%i
2*%i
2*%i
3
3
4*%i
];
E = sortRoots(E);
R = sortRoots(R);
assert_checkalmostequal(R, E, 1.e-3, 1.e-3);
// EXAMPLE 5. 12 ZEROS EVENLY DISTRIBUTE ON A CIRCLE OF 
// RADIUS 1. CENTERED AT 0+2I
// Real part:
// 4095 - 67584*x^2 + 126720*x^4 - 59136*x^6 + 7920*x^8 - 264*x^10 + x^12
// Imaginary part:
// 24576*x - 112640x^3 + 101376x^5 - 25344x^7 + 1760x^9 - 24x^11
P = [];
PI=[];
P(1)=1;
P(2)=0;
P(3)=-264;
P(4)=0;
P(5)=7920;
P(6)=0;
P(7)=-59136;
P(8)=0;
P(9)=126720;
P(10)=0;
P(11)=-67584;
P(12)=0;
P(13)=4095;
PI(1)=0;
PI(2)=-24;
PI(3)=0;
PI(4)=1760;
PI(5)=0;
PI(6)=-25344;
PI(7)=0;
PI(8)=101376;
PI(9)=0;
PI(10)=-112640;
PI(11)=0;
PI(12)=24576;
PI(13)=0;
P = complex(P,PI);
R = roots(P);
S3=sqrt(3);
E = [
-1 + 2*%i
%i
3*%i
1+2*%i
(1/2)*(-S3+3*%i)
(1/2)*(-S3+5*%i)
-(1/2)*%i*(S3+(-4-%i))
(1/2)*((1+4*%i)-%i*S3)
(1/2)*%i*(S3+(4+%i))
(1/2)*((1+4*%i)+%i*S3)
(1/2)*(S3+3*%i)
(1/2)*(S3+5*%i)
];
E = sortRoots(E);
R = sortRoots(R);
assert_checkalmostequal(R, E, 1.e-10, 1.e-8);
assert_checkequal(roots([4 3 2 1]), roots(poly([1 2 3 4], 'x', 'coeff')));
assert_checkequal(roots([4 3 2 1] + [1 2 3 4]*%i), roots(poly([1 2 3 4]+[4 3 2 1]*%i,'x','coeff')));