summaryrefslogtreecommitdiffstats
path: root/scilab/modules
diff options
context:
space:
mode:
authorMichael Baudin <michael.baudin@scilab.org>2011-05-02 11:13:46 +0200
committerVincent COUVERT <vincent.couvert@scilab.org>2011-05-06 11:06:40 +0200
commita441ccc87520c803497d91b70b193e90ec2c350e (patch)
tree842170b7e5ff420df6d62a68f28e226739fa2dfc /scilab/modules
parenta5fb769b63e2aab8a70c9d41a9aed1efd332c50c (diff)
downloadscilab-a441ccc87520c803497d91b70b193e90ec2c350e.zip
scilab-a441ccc87520c803497d91b70b193e90ec2c350e.tar.gz
polynomials: bug 7101 fixed: the roots had convergence problems for some polynomials.
Change-Id: I66e50b40a2b68dbef35bf0bd5b076190fd6f742b
Diffstat (limited to 'scilab/modules')
-rw-r--r--scilab/modules/polynomials/help/en_US/roots.xml127
-rw-r--r--scilab/modules/polynomials/help/fr_FR/roots.xml140
-rw-r--r--scilab/modules/polynomials/sci_gateway/fortran/sci_f_roots.f136
-rw-r--r--scilab/modules/polynomials/tests/nonreg_tests/bug_415.tst4
-rw-r--r--scilab/modules/polynomials/tests/nonreg_tests/bug_7101.dia.ref272
-rw-r--r--scilab/modules/polynomials/tests/nonreg_tests/bug_7101.tst271
-rw-r--r--scilab/modules/polynomials/tests/unit_tests/poly_roots.dia.ref177
-rw-r--r--scilab/modules/polynomials/tests/unit_tests/poly_roots.tst145
-rw-r--r--scilab/modules/polynomials/tests/unit_tests/roots.dia.ref624
-rw-r--r--scilab/modules/polynomials/tests/unit_tests/roots.tst609
10 files changed, 2090 insertions, 415 deletions
diff --git a/scilab/modules/polynomials/help/en_US/roots.xml b/scilab/modules/polynomials/help/en_US/roots.xml
index 59078f0..bfe72c8 100644
--- a/scilab/modules/polynomials/help/en_US/roots.xml
+++ b/scilab/modules/polynomials/help/en_US/roots.xml
@@ -2,6 +2,7 @@
2<!-- 2<!--
3 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 3 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4 * Copyright (C) 2006-2008 - INRIA 4 * Copyright (C) 2006-2008 - INRIA
5 * Copyright (C) 2011 - DIGITEO - Michael Baudin
5 * 6 *
6 * This file must be used under the terms of the CeCILL. 7 * This file must be used under the terms of the CeCILL.
7 * This source file is licensed as described in the file COPYING, which 8 * This source file is licensed as described in the file COPYING, which
@@ -17,10 +18,6 @@
17 xmlns:ns5="http://www.w3.org/1999/xhtml" 18 xmlns:ns5="http://www.w3.org/1999/xhtml"
18 xmlns:mml="http://www.w3.org/1998/Math/MathML" 19 xmlns:mml="http://www.w3.org/1998/Math/MathML"
19 xmlns:db="http://docbook.org/ns/docbook"> 20 xmlns:db="http://docbook.org/ns/docbook">
20 <info>
21 <pubdate>$LastChangedDate: 2008-07-11 10:31:18 +0200 (ven., 11 juil. 2008)
22 $</pubdate>
23 </info>
24 21
25 <refnamediv> 22 <refnamediv>
26 <refname>roots</refname> 23 <refname>roots</refname>
@@ -31,8 +28,10 @@
31 <refsynopsisdiv> 28 <refsynopsisdiv>
32 <title>Calling Sequence</title> 29 <title>Calling Sequence</title>
33 30
34 <synopsis>[x]=roots(p) 31 <synopsis>
35[x]=roots(p,'e')</synopsis> 32 x=roots(p)
33 x=roots(p,algo)
34 </synopsis>
36 </refsynopsisdiv> 35 </refsynopsisdiv>
37 36
38 <refsection> 37 <refsection>
@@ -43,35 +42,114 @@
43 <term>p</term> 42 <term>p</term>
44 43
45 <listitem> 44 <listitem>
46 <para>polynomial with real or complex coefficients or vector of the 45 <para>
47 polynomial coefficients in decreasing degree order (Matlab 46 a polynomial with real or complex coefficients, or
48 compatibility).</para> 47 a m-by-1 or 1-by-m matrix of doubles, the polynomial coefficients in decreasing degree order.
48 </para>
49 </listitem> 49 </listitem>
50 </varlistentry> 50 </varlistentry>
51
52 <varlistentry>
53 <term>algo</term>
54
55 <listitem>
56 <para>
57 a string, the algorithm to be used (default algo="f").
58 If algo="e", then the eigenvalues of the companion matrix are returned.
59 If algo="f", then the Jenkins-Traub method is used (if the polynomial is real and
60 has degree lower than 100).
61 If algo="f" and the polynomial is complex, then an error is generated.
62 If algo="f" and the polynomial has degree greater than 100, then an error is generated.
63 </para>
64 </listitem>
65 </varlistentry>
66
51 </variablelist> 67 </variablelist>
52 </refsection> 68 </refsection>
53 69
54 <refsection> 70 <refsection>
55 <title>Description</title> 71 <title>Description</title>
56 72
57 <para><literal>x=roots(p)</literal> returns in the complex vector 73 <para>
58 <literal>x</literal> the roots of the polynomial <literal>p</literal>. For 74 This function returns in the complex vector
59 real polynomials of degree &lt;=100 the fast RPOLY algorithm (based on 75 <literal>x</literal> the roots of the polynomial <literal>p</literal>.
60 Jenkins-Traub method) is used. In the other cases the roots are computed 76 </para>
61 as the eigenvalues of the associated companion matrix. Use 77
62 <literal>x=roots(p,'e')</literal> to force this algorithm in any 78 <para>
63 cases.</para> 79 The "e" option corresponds to method based on the eigenvalues of the
80 companion matrix.
81 </para>
82
83 <para>
84 The "f" option corresponds to the fast RPOLY algorithm, based on
85 Jenkins-Traub method.
86 </para>
87
88 <para>
89 For real polynomials of degree &lt;=100, users may consider the "f" option,
90 which might be faster in some cases.
91 On the other hand, some specific polynomials are known to be able to
92 make this option to fail.
93 </para>
94
64 </refsection> 95 </refsection>
65 96
66 <refsection> 97 <refsection>
67 <title>Examples</title> 98 <title>Examples</title>
68 99
69 <programlisting role="example"><![CDATA[ 100 <para>
101 In the following examples, we compute roots of polynomials.
102 </para>
103
104 <programlisting role="example"><![CDATA[
105// Roots given a real polynomial
106p = poly([1 2 3],"x")
107roots(p)
108// Roots, given the real coefficients
109p = [3 2 1]
110roots(p)
111// The roots of a complex polynomial
70p=poly([0,10,1+%i,1-%i],'x'); 112p=poly([0,10,1+%i,1-%i],'x');
71roots(p) 113roots(p)
72A=rand(3,3);roots(poly(A,'x')) // Evals by characteristic polynomial 114// The roots of the polynomial of a matrix
115A=rand(3,3);
116p = poly(A,'x')
117roots(p)
73spec(A) 118spec(A)
74 ]]></programlisting> 119 ]]></programlisting>
120
121 <para>
122 The polynomial representation can have a significant
123 impact on the roots.
124 In the following example, suggested by Wilkinson in the 60s and presented by Moler,
125 we consider a diagonal matrix with diagonal entries equal to 1, 2, ..., 20.
126 The eigenvalues are obviously equal to 1, 2, ..., 20.
127 If we compute the associated characteristic polynomial and compute its roots,
128 we can see that the eigenvalues are significantly different from the expected
129 ones.
130 This implies that just representing the coefficients as IEEE doubles changes the
131 roots.
132 </para>
133
134 <programlisting role="example"><![CDATA[
135A = diag(1:20);
136p = poly(A,'x')
137roots(p)
138 ]]></programlisting>
139
140 <para>
141 The "f" option produces an error if the polynomial is complex or
142 if the degree is greater than 100.
143 </para>
144
145 <programlisting role="example"><![CDATA[
146 // The following case produces an error.
147 p = %i+%s;
148 roots(p,"f")
149 // The following case produces an error.
150 p = ones(101,1);
151 roots(p,"f")
152 ]]></programlisting>
75 </refsection> 153 </refsection>
76 154
77 <refsection role="see also"> 155 <refsection role="see also">
@@ -90,19 +168,26 @@ spec(A)
90 <title>Authors</title> 168 <title>Authors</title>
91 169
92 <simplelist type="vert"> 170 <simplelist type="vert">
93 <member>Serge Steer (INRIA)</member> 171 <member>
172 Serge Steer (INRIA)
173 </member>
174 <member>
175 Copyright (C) 2011 - DIGITEO - Michael Baudin
176 </member>
94 </simplelist> 177 </simplelist>
95 </refsection> 178 </refsection>
96 179
97 <refsection> 180 <refsection>
98 <title>References</title> 181 <title>References</title>
99 182
100 <para>The RPOLY algorithm is described in "Algorithm 493: Zeros of a Real 183 <para>
184 The RPOLY algorithm is described in "Algorithm 493: Zeros of a Real
101 Polynomial", ACM TOMS Volume 1, Issue 2 (June 1975), pp. 178-189</para> 185 Polynomial", ACM TOMS Volume 1, Issue 2 (June 1975), pp. 178-189</para>
102 <para>Jenkins, M. A. and Traub, J. F. (1970), A Three-Stage Algorithm for 186 <para>Jenkins, M. A. and Traub, J. F. (1970), A Three-Stage Algorithm for
103 Real Polynomials Using Quadratic Iteration, SIAM J. Numer. Anal., 7(1970), 545-566.</para> 187 Real Polynomials Using Quadratic Iteration, SIAM J. Numer. Anal., 7(1970), 545-566.</para>
104 <para>Jenkins, M. A. and Traub, J. F. (1970), Principles for Testing Polynomial Zerofinding Programs. 188 <para>Jenkins, M. A. and Traub, J. F. (1970), Principles for Testing Polynomial Zerofinding Programs.
105 ACM TOMS 1, 1 (March 1975), pp. 26-34</para> 189 ACM TOMS 1, 1 (March 1975), pp. 26-34
190 </para>
106 191
107 </refsection> 192 </refsection>
108 193
diff --git a/scilab/modules/polynomials/help/fr_FR/roots.xml b/scilab/modules/polynomials/help/fr_FR/roots.xml
index 1e53128..06e8c5e 100644
--- a/scilab/modules/polynomials/help/fr_FR/roots.xml
+++ b/scilab/modules/polynomials/help/fr_FR/roots.xml
@@ -1,4 +1,16 @@
1<?xml version="1.0" encoding="UTF-8"?> 1<?xml version="1.0" encoding="UTF-8"?>
2<!--
3 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4 * Copyright (C) ???? - 2008 - INRIA
5 * Copyright (C) 2011 - DIGITEO - Michael Baudin
6 *
7 * This file must be used under the terms of the CeCILL.
8 * This source file is licensed as described in the file COPYING, which
9 * you should have received as part of this distribution. The terms
10 * are also available at
11 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
12 *
13 -->
2<refentry version="5.0-subset Scilab" xml:id="roots" xml:lang="fr" 14<refentry version="5.0-subset Scilab" xml:id="roots" xml:lang="fr"
3 xmlns="http://docbook.org/ns/docbook" 15 xmlns="http://docbook.org/ns/docbook"
4 xmlns:xlink="http://www.w3.org/1999/xlink" 16 xmlns:xlink="http://www.w3.org/1999/xlink"
@@ -6,10 +18,6 @@
6 xmlns:ns5="http://www.w3.org/1999/xhtml" 18 xmlns:ns5="http://www.w3.org/1999/xhtml"
7 xmlns:mml="http://www.w3.org/1998/Math/MathML" 19 xmlns:mml="http://www.w3.org/1998/Math/MathML"
8 xmlns:db="http://docbook.org/ns/docbook"> 20 xmlns:db="http://docbook.org/ns/docbook">
9 <info>
10 <pubdate>$LastChangedDate: 2008-07-11 10:34:31 +0200 (ven., 11 juil. 2008)
11 $</pubdate>
12 </info>
13 21
14 <refnamediv> 22 <refnamediv>
15 <refname>roots</refname> 23 <refname>roots</refname>
@@ -20,7 +28,10 @@
20 <refsynopsisdiv> 28 <refsynopsisdiv>
21 <title>Séquence d'appel</title> 29 <title>Séquence d'appel</title>
22 30
23 <synopsis>[x]=roots(p)</synopsis> 31 <synopsis>
32 x=roots(p)
33 x=roots(p,"algo")
34 </synopsis>
24 </refsynopsisdiv> 35 </refsynopsisdiv>
25 36
26 <refsection> 37 <refsection>
@@ -31,34 +42,108 @@
31 <term>p</term> 42 <term>p</term>
32 43
33 <listitem> 44 <listitem>
34 <para>polynôme à coefficients réels ou complexes ou vecteur des 45 <para>
35 coefficients du polynŽôme ordonnés par puissance décroissante 46 un polynôme réel ou complexe, ou
36 (compatibilité avec Matlab).</para> 47 une m-par-1 or 1-par-m matrice de doubles, les coefficients du polynômes par ordre de puissance décroissante.
48 </para>
37 </listitem> 49 </listitem>
38 </varlistentry> 50 </varlistentry>
51
52 <varlistentry>
53 <term>algo</term>
54
55 <listitem>
56 <para>
57 une chaîne de caractères, l'algorithme à utiliser (défaut algo="f").
58 Si algo="e", alors les valeurs propres de la matrice compagnion sont utilisées.
59 Si algo="f", alors l'algorithme de Jenkins-Traub est utilisé (si les coefficients
60 du polynôme sont réels et que le degré du polynôme est plus petit que 100).
61 Si algo="f" et que le polynôme est complexe, alors une erreur est générée.
62 Si algo="f" et que le polynôme est de degré inférieur à 100, alors une erreur est
63 générée.
64 </para>
65 </listitem>
66 </varlistentry>
67
39 </variablelist> 68 </variablelist>
40 </refsection> 69 </refsection>
41 70
42 <refsection> 71 <refsection>
43 <title>Description</title> 72 <title>Description</title>
44 73
45 <para><literal>x=roots(p)</literal> renvoie dans le vecteur complexe 74 <para>
46 <literal>x</literal> les racines du polynôme <literal>p</literal>. Pour 75 Cette fonction retourne dans le vecteur complexe
47 les polynômes à coefficients réels et de degré &lt;=100, l'algorithme 76 <literal>x</literal> les racines du polynôme <literal>p</literal>.
48 rapide RPOLY (fondé sur la méthode de Jenkins-Traub) est utilisé. Dans les 77 </para>
49 autres cas, les racines sont calculées comme valeurs propres de la matrice 78
50 compagnon du polynôme. Pour forcer ce dernier algorithme dans tous les 79 <para>
51 cas, utilisez <literal>x=roots(p,'e')</literal>.</para> 80 L'option "f" utilise l'algorithme rapide RPOLY, fondé sur la méthode de Jenkins-Traub.
81 </para>
82
83 <para>
84 Pour les polynôme réels de degré inférieur à 100, on peut utiliser
85 l'option "f", qui peut être plus rapide dans certains cas.
86 Toutefois, certains polynômes sont susceptibles de poser des problèmes
87 de convergences pour l'algorithme associé à l'option "f".
88 </para>
52 </refsection> 89 </refsection>
53 90
54 <refsection> 91 <refsection>
55 <title>Exemples</title> 92 <title>Exemples</title>
56 93
57 <programlisting role="example"><![CDATA[ 94 <para>
95 Dans les exemples suivants, on calcule des racines de polynômes.
96 </para>
97
98 <programlisting role="example"><![CDATA[
99// Un polynôme réel.
100p = poly([1 2 3],"x")
101roots(p)
102// Les coefficients du polynôme sont donnés.
103p = [3 2 1]
104roots(p)
105// Les racines d'un polynôme complexe.
58p=poly([0,10,1+%i,1-%i],'x'); 106p=poly([0,10,1+%i,1-%i],'x');
59roots(p) 107roots(p)
60A=rand(3,3);roots(poly(A,'x')) // comparaison via le polynôme caractéristique 108// Les racines du polynôme caractéristique d'une matrice.
61spec(A) 109A=rand(3,3);
110p = poly(A,'x')
111roots(p)
112spec(A)
113 ]]></programlisting>
114
115 <para>
116 La représentation polynômiale peut avoir un impact significatif sur les
117 racines.
118 Dans l'exemple suivant, suggéré par Wilkinson dans les années 60 et présenté
119 par Moler, on considère une matrice dont les termes diagonaux sont égaux
120 à 1, 2, ..., 20.
121 Bien entendu, les racines du polynôme caractéristique sont 1, 2, ..., 20.
122 Si on calcule le polynôme caractéristique associé et qu'on calcule ses
123 racines, on peut voir qu'elles sont significativement différentes des
124 valeurs attendues.
125 Cela montre que le seul fait de représenter les coefficients dans
126 des doubles IEEE change les racines.
127 </para>
128
129 <programlisting role="example"><![CDATA[
130A = diag(1:20);
131p = poly(A,'x')
132roots(p)
133 ]]></programlisting>
134
135 <para>
136 L'option "f" produit une erreur si le polynôme est complexe
137 ou que le degré est plus grand que 100.
138 </para>
139
140 <programlisting role="example"><![CDATA[
141 // Le cas suivant produit une erreur
142 p = %i+%s;
143 roots(p,"f")
144 // Le cas suivant produit une erreur
145 p = ones(101,1);
146 roots(p,"f")
62 ]]></programlisting> 147 ]]></programlisting>
63 </refsection> 148 </refsection>
64 149
@@ -78,23 +163,32 @@ spec(A)
78 <title>Auteurs</title> 163 <title>Auteurs</title>
79 164
80 <simplelist type="vert"> 165 <simplelist type="vert">
81 <member>Serge Steer (INRIA)</member> 166 <member>
167 Serge Steer (INRIA)
168 </member>
169 <member>
170 Copyright (C) 2011 - DIGITEO - Michael Baudin
171 </member>
82 </simplelist> 172 </simplelist>
83 </refsection> 173 </refsection>
84 174
85 <refsection> 175 <refsection>
86 <title>Bibliographie</title> 176 <title>Bibliographie</title>
87 177
88 <para>La routine RPOLY est decrite dans "Algorithm 493: Zeros of a Real 178 <para>
89 Polynomial", ACM TOMS Volume 1, Issue 2 (June 1975), pp. 178-189</para> 179 La routine RPOLY est decrite dans "Algorithm 493: Zeros of a Real
180 Polynomial", ACM TOMS Volume 1, Issue 2 (June 1975), pp. 178-189
181 </para>
90 </refsection> 182 </refsection>
91 183
92 <refsection> 184 <refsection>
93 <title>Fonctions Utilisées</title> 185 <title>Fonctions Utilisées</title>
94 186
95 <para>Le code source de rpoly.f peut être trouvé dans le repertoire 187 <para>
188 Le code source de rpoly.f peut être trouvé dans le repertoire
96 SCI/modules/polynomials/src/fortran de la distribution source de Scilab. Dans le cas où la 189 SCI/modules/polynomials/src/fortran de la distribution source de Scilab. Dans le cas où la
97 matrix compagnon est utilisée, le calcul des valeurs propres est effectué 190 matrix compagnon est utilisée, le calcul des valeurs propres est effectué
98 en utilisant les routines DGEEV et ZGEEV de LAPACK.</para> 191 en utilisant les routines DGEEV et ZGEEV de LAPACK.
192 </para>
99 </refsection> 193 </refsection>
100</refentry> 194</refentry>
diff --git a/scilab/modules/polynomials/sci_gateway/fortran/sci_f_roots.f b/scilab/modules/polynomials/sci_gateway/fortran/sci_f_roots.f
index 54cb8ac..f74384d 100644
--- a/scilab/modules/polynomials/sci_gateway/fortran/sci_f_roots.f
+++ b/scilab/modules/polynomials/sci_gateway/fortran/sci_f_roots.f
@@ -1,5 +1,6 @@
1c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 1c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2c Copyright (C) ????-2008 - INRIA 2c Copyright (C) ????-2008 - INRIA - Serge Steer
3c Copyright (C) 2011 - DIGITEO - Michael Baudin
3c 4c
4c This file must be used under the terms of the CeCILL. 5c This file must be used under the terms of the CeCILL.
5c This source file is licensed as described in the file COPYING, which 6c This source file is licensed as described in the file COPYING, which
@@ -15,8 +16,10 @@ c http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
15 integer vol 16 integer vol
16 integer fail 17 integer fail
17 logical ref,eigen 18 logical ref,eigen
18 integer eig 19 integer algo_eig
19 data eig/14/ 20 data algo_eig/14/
21 integer algo_fast
22 data algo_fast/15/
20c 23c
21 iadr(l)=l+l-1 24 iadr(l)=l+l-1
22 sadr(l)=(l/2)+1 25 sadr(l)=(l/2)+1
@@ -31,7 +34,11 @@ c
31 return 34 return
32 endif 35 endif
33c 36c
34 eigen=.false. 37c if algo=="e", then eigen = .true.
38c if algo=="f", then eigen = .false.
39c
40 eigen=.true.
41c Read the algo option
35 if (rhs.eq.2) then 42 if (rhs.eq.2) then
36 ilopt=iadr(lstk(top)) 43 ilopt=iadr(lstk(top))
37 if(istk(ilopt).lt.0) ilopt=iadr(istk(ilopt+1)) 44 if(istk(ilopt).lt.0) ilopt=iadr(istk(ilopt+1))
@@ -45,12 +52,21 @@ c
45 call error(89) 52 call error(89)
46 return 53 return
47 endif 54 endif
48 if(istk(ilopt+5).ne.2.or.istk(ilopt+6).ne.eig) then 55 if(istk(ilopt+5).ne.2) then
49 err=2 56 err=2
50 call error(116) 57 call error(116)
51 return 58 return
59 else
60 if (istk(ilopt+6)==algo_eig) then
61 eigen=.true.
62 elseif (istk(ilopt+6)==algo_fast) then
63 eigen=.false.
64 else
65 err=2
66 call error(116)
67 return
68 endif
52 endif 69 endif
53 eigen=.true.
54 rhs=rhs-1 70 rhs=rhs-1
55 top=top-1 71 top=top-1
56 endif 72 endif
@@ -68,9 +84,18 @@ c
68 m1=istk(il1+1) 84 m1=istk(il1+1)
69 n1=istk(il1+2) 85 n1=istk(il1+2)
70 mn1=m1*n1 86 mn1=m1*n1
71 if(m1*n1.eq.0) return 87 if(m1*n1.eq.0) then
88 return
89 endif
72 90
73 it1=istk(il1+3) 91 it1=istk(il1+3)
92c If "fast" algo was chosen and polynomial is complex,
93c then produce an error.
94 if ( .not.eigen .and. it1 .eq. 1 ) then
95 err=2
96 call error(116)
97 return
98 endif
74 99
75 if(istk(il1).eq.1) then 100 if(istk(il1).eq.1) then
76c for Matlab compatibility root of the vector of coefficients 101c for Matlab compatibility root of the vector of coefficients
@@ -109,12 +134,26 @@ c for Matlab compatibility root of the vector of coefficients
109 endif 134 endif
110 l1=sadr(ilr+4) 135 l1=sadr(ilr+4)
111 21 n=n-1 136 21 n=n-1
112 if(n.lt.0) goto 24 137 if(n.lt.0) then
138 goto 24
139 endif
113 t=abs(stk(lc+n)) 140 t=abs(stk(lc+n))
114 if(it1.eq.1) t=t+abs(stk(lc+n+vol)) 141 if(it1.eq.1) then
115 if(t.eq.0.0d+0) goto 21 142 t=t+abs(stk(lc+n+vol))
143 endif
144 if(t.eq.0.0d+0) then
145 goto 21
146 endif
147
148 if ( .not.eigen .and.n.gt.100) then
149c If "fast" algo was chosen and polynomial has degree greater than 100,
150c then produce an error.
151 err=2
152 call error(116)
153 return
154 endif
116 155
117 if (eigen) goto 22 156 if (.not.eigen) then
118c 157c
119c real polynomial: rpoly algorithm 158c real polynomial: rpoly algorithm
120c this alg is much more speedy, but it may happens that it gives 159c this alg is much more speedy, but it may happens that it gives
@@ -124,7 +163,6 @@ C 1.355 and 6.65 and the other ones inside a circle centered in 0
124C with radius 1 163C with radius 1
125C 164C
126 165
127 if(it1.eq.0.and.n.le.100) then
128 lp=max(lw,l1+2*n) 166 lp=max(lw,l1+2*n)
129 err=lp+n+1-lstk(bot) 167 err=lp+n+1-lstk(bot)
130 if(err.gt.0) then 168 if(err.gt.0) then
@@ -151,50 +189,52 @@ C
151 if(n.eq.0) istk(ilr+2)=0 189 if(n.eq.0) istk(ilr+2)=0
152 istk(ilr+3)=1 190 istk(ilr+3)=1
153 lstk(top+1)=l1+2*n 191 lstk(top+1)=l1+2*n
154 goto 999
155 endif
156
157 22 continue
158c
159c Companion matrix method
160 lw=lw+n*n*(it1+1)
161 err=lw+n*(it1+1)-lstk(bot)
162 if(err.gt.0) then
163 call error(17)
164 return 192 return
165 endif
166 sr=stk(lc+n)
167 call unsfdcopy(n,stk(lc),-1,stk(lw),1)
168 if(it1.eq.0) then
169 call dscal(n,-1.0d+0/sr,stk(lw),1)
170 else 193 else
171 si=stk(lc+vol+n) 194c
172 t=sr*sr+si*si 195c Companion matrix method
173 sr=-sr/t 196 lw=lw+n*n*(it1+1)
174 si=si/t 197 err=lw+n*(it1+1)-lstk(bot)
175 call unsfdcopy(n,stk(lc+vol),-1,stk(lw+n),1) 198 if(err.gt.0) then
176 call wscal(n,sr,si,stk(lw),stk(lw+n),1) 199 call error(17)
177 endif 200 return
178 call dset(n*n*(it1+1),0.0d+0,stk(l1),1) 201 endif
179 call dset(n-1,1.0d+0,stk(l1+n),n+1) 202 sr=stk(lc+n)
180 call unsfdcopy(n,stk(lw),1,stk(l1),1) 203 call unsfdcopy(n,stk(lc),-1,stk(lw),1)
181 if(it1.eq.1) call unsfdcopy(n,stk(lw+n),1,stk(l1+n*n),1) 204 if(it1.eq.0) then
182 lstk(top+1)=l1+n*n*(it1+1) 205 call dscal(n,-1.0d+0/sr,stk(lw),1)
183 istk(ilr)=1 206 else
184 istk(ilr+1)=n 207 si=stk(lc+vol+n)
185 istk(ilr+2)=n 208 t=sr*sr+si*si
186 istk(ilr+3)=it1 209 sr=-sr/t
187 fin=3 210 si=si/t
188 fun=2 211 call unsfdcopy(n,stk(lc+vol),-1,stk(lw+n),1)
212 call wscal(n,sr,si,stk(lw),stk(lw+n),1)
213 endif
214 call dset(n*n*(it1+1),0.0d+0,stk(l1),1)
215 call dset(n-1,1.0d+0,stk(l1+n),n+1)
216 call unsfdcopy(n,stk(lw),1,stk(l1),1)
217 if(it1.eq.1) then
218 call unsfdcopy(n,stk(lw+n),1,stk(l1+n*n),1)
219 endif
220 lstk(top+1)=l1+n*n*(it1+1)
221 istk(ilr)=1
222 istk(ilr+1)=n
223 istk(ilr+2)=n
224 istk(ilr+3)=it1
225 fin=3
226 fun=2
189c *call* matds(r c) 227c *call* matds(r c)
190 goto 999 228 return
229 endif
230
191c polynome de degre 0 231c polynome de degre 0
192 24 istk(ilr)=1 232 24 istk(ilr)=1
193 istk(ilr+1)=0 233 istk(ilr+1)=0
194 istk(ilr+2)=0 234 istk(ilr+2)=0
195 istk(ilr+3)=0 235 istk(ilr+3)=0
196 lstk(top+1)=sadr(ilr+4) 236 lstk(top+1)=sadr(ilr+4)
197 goto 999 237
198 999 return 238 return
199 end 239 end
200c ======================================= 240c =======================================
diff --git a/scilab/modules/polynomials/tests/nonreg_tests/bug_415.tst b/scilab/modules/polynomials/tests/nonreg_tests/bug_415.tst
index 0c716a3..6efffc4 100644
--- a/scilab/modules/polynomials/tests/nonreg_tests/bug_415.tst
+++ b/scilab/modules/polynomials/tests/nonreg_tests/bug_415.tst
@@ -187,10 +187,12 @@ function flag = assert_close ( computed, expected, epsilon )
187 if flag <> 1 then pause,end 187 if flag <> 1 then pause,end
188endfunction 188endfunction
189 189
190
191// There is no problem in this test: only the order of the
192// eigenvalues change.
190t=poly(0,"t"); 193t=poly(0,"t");
191p=t^14 - 15*t^12 - t^11 + 89*t^10 + 12*t^9 - 263*t^8 - 53*t^7 + 397*t^6 + 103*t^5 - 275*t^4 - 78*t^3 + 62*t^2 + 8*t - 7; 194p=t^14 - 15*t^12 - t^11 + 89*t^10 + 12*t^9 - 263*t^8 - 53*t^7 + 397*t^6 + 103*t^5 - 275*t^4 - 78*t^3 + 62*t^2 + 8*t - 7;
192myroots=roots(p); 195myroots=roots(p);
193//computedroots = sort(myroots);
194computed = sort_merge ( myroots , compare_complexrealimag ); 196computed = sort_merge ( myroots , compare_complexrealimag );
195expected = [ 197expected = [
196- 1.9914144710587742270747 198- 1.9914144710587742270747
diff --git a/scilab/modules/polynomials/tests/nonreg_tests/bug_7101.dia.ref b/scilab/modules/polynomials/tests/nonreg_tests/bug_7101.dia.ref
new file mode 100644
index 0000000..c12d06b
--- /dev/null
+++ b/scilab/modules/polynomials/tests/nonreg_tests/bug_7101.dia.ref
@@ -0,0 +1,272 @@
1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) 2011 - DIGITEO - Michael Baudin
4//
5// This file is distributed under the same license as the Scilab package.
6// =============================================================================
7// <-- JVM NOT MANDATORY -->
8// <-- Non-regression test for bug 7101 -->
9//
10// <-- Bugzilla URL -->
11// http://bugzilla.scilab.org/show_bug.cgi?id=7101
12//
13// <-- Short Description -->
14// The roots / Jenkins-Traub algorithm does not produce the roots, sometimes.
15// The solution is to use the eigenvalues of the companion matrix.
16//
17// sort_merge_comparison --
18// Returns -1 if x < y,
19// returns 0 if x==y,
20// returns +1 if x > y
21//
22function order = sort_merge_comparison ( x , y )
23 if x < y then
24 order = -1
25 elseif x==y then
26 order = 0
27 else
28 order = 1
29 end
30endfunction
31//
32// sort_merge --
33// Returns the sorted array x.
34// Arguments
35// x : the array to sort
36// compfun : the comparison function
37// Bruno Pincon
38// "quelques tests de rapidit´e entre diff´erents logiciels matriciels"
39// Modified by Michael Baudin to manage a comparison function
40//
41function [x] = sort_merge ( varargin )
42 [lhs,rhs]=argn();
43 if rhs<>1 & rhs<>2 then
44 errmsg = sprintf("Unexpected number of arguments : %d provided while 1 or 2 are expected.",rhs);
45 error(errmsg)
46 end
47 // Get the array x
48 x = varargin(1);
49 // Get the comparison function compfun
50 if rhs==1 then
51 compfun = sort_merge_comparison;
52 else
53 compfun = varargin(2);
54 end
55 // Proceed...
56 n = length(x)
57 if n > 1 then
58 m = floor(n/2);
59 p = n-m
60 x1 = sort_merge ( x(1:m) , compfun )
61 x2 = sort_merge ( x(m+1:n) , compfun )
62 i = 1;
63 i1 = 1;
64 i2 = 1;
65 for i = 1:n
66 order = compfun ( x1(i1) , x2(i2) );
67 if order<=0 then
68 x(i) = x1(i1)
69 i1 = i1+1
70 if (i1 > m) then
71 x(i+1:n) = x2(i2:p)
72 break
73 end
74 else
75 x(i) = x2(i2)
76 i2 = i2+1
77 if (i2 > p) then
78 x(i+1:n) = x1(i1:m)
79 break
80 end
81 end
82 end
83 end
84endfunction
85//
86// compare_complexrealimag --
87// Returns -1 if a < b,
88// returns 0 if a==b,
89// returns +1 if a > b
90// Compare first by real parts, then by imaginary parts.
91//
92function order = compare_complexrealimag ( a , b )
93 ar = real(a)
94 br = real(b)
95 if ar < br then
96 order = -1
97 elseif ar > br then
98 order = 1
99 else
100 ai = imag(a)
101 bi = imag(b)
102 if ai < bi then
103 order = -1
104 elseif ai == bi then
105 order = 0
106 else
107 order = 1
108 end
109 end
110endfunction
111//
112// assert_close --
113// Returns 1 if the two real matrices computed and expected are close,
114// i.e. if the relative distance between computed and expected is lesser than epsilon.
115// Arguments
116// computed, expected : the two matrices to compare
117// epsilon : a small number
118//
119function flag = assert_close ( computed, expected, epsilon )
120 if expected==0.0 then
121 shift = norm(computed-expected);
122 else
123 shift = norm(computed-expected)/norm(expected);
124 end
125 if shift < epsilon then
126 flag = 1;
127 else
128 flag = 0;
129 end
130 if flag <> 1 then bugmes();quit;end
131endfunction
132function y = sortmyroots(x)
133 // Sort the roots of a polynomial with a customized
134 // complex-aware sorting algorithm.
135 y = sort_merge ( x , compare_complexrealimag );
136endfunction
137// Failed on Windows 32 bits
138p = [1,1,-7,-15,1,-4,4,7,4,-53,1,-53,-8,3,3,0,9,-15];
139r = roots(p);
140e = [
141 2.9977242
142 - 2.0998215 + 1.0381514 * %i
143 - 2.0998215 - 1.0381514 * %i
144 - 1.1261224 + 0.7687233 * %i
145 - 1.1261224 - 0.7687233 * %i
146 1.1176579 + 0.5115332 * %i
147 1.1176579 - 0.5115332 * %i
148 - 0.7359417 + 0.3731641 * %i
149 - 0.7359417 - 0.3731641 * %i
150 0.2849638 + 0.9531919 * %i
151 0.2849638 - 0.9531919 * %i
152 0.0897371 + 1.0370037 * %i
153 0.0897371 - 1.0370037 * %i
154 - 0.1740455 + 0.9263179 * %i
155 - 0.1740455 - 0.9263179 * %i
156 0.6447102 + 0.2914081 * %i
157 0.6447102 - 0.2914081 * %i
158 ];
159e = sortmyroots(e);
160r = sortmyroots(r);
161assert_close ( r, e, 1.e-6 );
162// Failed on Mac, on Windows
163p=[1,1,-7,-35,1,-4,4,7,4,-88,1,-88,-8,3,3,0,9,-35];
164r = roots(p);
165e = [
166 3.6133489
167 - 2.3323533 + 2.0888127 * %i
168 - 2.3323533 - 2.0888127 * %i
169 1.0856792 + 0.5138318 * %i
170 1.0856792 - 0.5138318 * %i
171 - 1.1030013 + 0.6108696 * %i
172 - 1.1030013 - 0.6108696 * %i
173 0.3226838 + 0.9451270 * %i
174 0.3226838 - 0.9451270 * %i
175 0.0250044 + 1.0210451 * %i
176 0.0250044 - 1.0210451 * %i
177 - 0.2556563 + 0.9467085 * %i
178 - 0.2556563 - 0.9467085 * %i
179 - 0.7512303 + 0.3765797 * %i
180 - 0.7512303 - 0.3765797 * %i
181 0.7021994 + 0.3415821 * %i
182 0.7021994 - 0.3415821 * %i
183];
184e = sortmyroots(e);
185r = sortmyroots(r);
186assert_close ( r, e, 1.e-6 );
187// Failed on Linux
188p=[1,1,-7,-80,1,-4,4,7,4,-27,1,-27,-8,3,3,0,9,-80];
189r = roots(p);
190e = [
191 - 2.7595524 + 3.1924496 * %i
192 - 2.7595524 - 3.1924496 * %i
193 4.5006465
194 - 0.9689444 + 0.2683252 * %i
195 - 0.9689444 - 0.2683252 * %i
196 - 0.8111357 + 0.6166997 * %i
197 - 0.8111357 - 0.6166997 * %i
198 - 0.3893539 + 0.9194344 * %i
199 - 0.3893539 - 0.9194344 * %i
200 0.0061369 + 1.0065796 * %i
201 0.0061369 - 1.0065796 * %i
202 0.4195701 + 0.9089127 * %i
203 0.4195701 - 0.9089127 * %i
204 0.9590394 + 0.2589039 * %i
205 0.9590394 - 0.2589039 * %i
206 0.7939168 + 0.5672744 * %i
207 0.7939168 - 0.5672744 * %i
208];
209e = sortmyroots(e);
210r = sortmyroots(r);
211assert_close ( r, e, 1.e-6 );
212// Failed on Windows 32 bits
213p=[1,0,1,1,1,-1,-1,1,1,0,1,0,-1,-1,1,-2,0,0,1,-1,1];
214r = roots(p);
215e = [
216 0.5444059 + 1.3082079 * %i
217 0.5444059 - 1.3082079 * %i
218 - 1.0517348 + 0.2347104 * %i
219 - 1.0517348 - 0.2347104 * %i
220 - 0.5893898 + 0.9840032 * %i
221 - 0.5893898 - 0.9840032 * %i
222 - 0.8170407 + 0.5459189 * %i
223 - 0.8170407 - 0.5459189 * %i
224 - 0.6570402 + 0.7150468 * %i
225 - 0.6570402 - 0.7150468 * %i
226 0.0129780 + 0.9748750 * %i
227 0.0129780 - 0.9748750 * %i
228 0.9192290 + 0.4894403 * %i
229 0.9192290 - 0.4894403 * %i
230 0.8691302 + 0.0832523 * %i
231 0.8691302 - 0.0832523 * %i
232 0.4975871 + 0.6807740 * %i
233 0.4975871 - 0.6807740 * %i
234 0.2718754 + 0.7528695 * %i
235 0.2718754 - 0.7528695 * %i
236];
237e = sortmyroots(e);
238r = sortmyroots(r);
239assert_close ( r, e, 1.e-6 );
240// A loop on several polynomials
241for i=-100:100
242 if ( modulo(i,20)==0 ) then
243 mprintf("i=%d\n",i);
244 end
245 for j=-100:100
246 p=[1 1 -7 j 1 -4 4 7 4 i 1 i -8 3 3 0 9 j];
247 roots(p);
248 end;
249end;
250i=-100
251i=-80
252i=-60
253i=-40
254i=-20
255i=0
256i=20
257i=40
258i=60
259i=80
260i=100
261// A loop on random polynomials.
262// The coefficients are integers
263for i = 1:3000
264 if ( modulo(i,1000)==0 ) then
265 mprintf("i=%d\n",i);
266 end
267 p = [1 round(4*rand(1,20)-2)];
268 roots(p);
269end
270i=1000
271i=2000
272i=3000
diff --git a/scilab/modules/polynomials/tests/nonreg_tests/bug_7101.tst b/scilab/modules/polynomials/tests/nonreg_tests/bug_7101.tst
new file mode 100644
index 0000000..0e103ca
--- /dev/null
+++ b/scilab/modules/polynomials/tests/nonreg_tests/bug_7101.tst
@@ -0,0 +1,271 @@
1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) 2011 - DIGITEO - Michael Baudin
4//
5// This file is distributed under the same license as the Scilab package.
6// =============================================================================
7
8// <-- JVM NOT MANDATORY -->
9
10// <-- Non-regression test for bug 7101 -->
11//
12// <-- Bugzilla URL -->
13// http://bugzilla.scilab.org/show_bug.cgi?id=7101
14//
15// <-- Short Description -->
16// The roots / Jenkins-Traub algorithm does not produce the roots, sometimes.
17// The solution is to use the eigenvalues of the companion matrix.
18
19//
20// sort_merge_comparison --
21// Returns -1 if x < y,
22// returns 0 if x==y,
23// returns +1 if x > y
24//
25function order = sort_merge_comparison ( x , y )
26 if x < y then
27 order = -1
28 elseif x==y then
29 order = 0
30 else
31 order = 1
32 end
33endfunction
34
35//
36// sort_merge --
37// Returns the sorted array x.
38// Arguments
39// x : the array to sort
40// compfun : the comparison function
41// Bruno Pincon
42// "quelques tests de rapidite entre differents logiciels matriciels"
43// Modified by Michael Baudin to manage a comparison function
44//
45function [x] = sort_merge ( varargin )
46 [lhs,rhs]=argn();
47 if rhs<>1 & rhs<>2 then
48 errmsg = sprintf("Unexpected number of arguments : %d provided while 1 or 2 are expected.",rhs);
49 error(errmsg)
50 end
51 // Get the array x
52 x = varargin(1);
53 // Get the comparison function compfun
54 if rhs==1 then
55 compfun = sort_merge_comparison;
56 else
57 compfun = varargin(2);
58 end
59 // Proceed...
60 n = length(x)
61 if n > 1 then
62 m = floor(n/2);
63 p = n-m
64 x1 = sort_merge ( x(1:m) , compfun )
65 x2 = sort_merge ( x(m+1:n) , compfun )
66 i = 1;
67 i1 = 1;
68 i2 = 1;
69 for i = 1:n
70 order = compfun ( x1(i1) , x2(i2) );
71 if order<=0 then
72 x(i) = x1(i1)
73 i1 = i1+1
74 if (i1 > m) then
75 x(i+1:n) = x2(i2:p)
76 break
77 end
78 else
79 x(i) = x2(i2)
80 i2 = i2+1
81 if (i2 > p) then
82 x(i+1:n) = x1(i1:m)
83 break
84 end
85 end
86 end
87 end
88endfunction
89
90//
91// compare_complexrealimag --
92// Returns -1 if a < b,
93// returns 0 if a==b,
94// returns +1 if a > b
95// Compare first by real parts, then by imaginary parts.
96//
97function order = compare_complexrealimag ( a , b )
98 ar = real(a)
99 br = real(b)
100 if ar < br then
101 order = -1
102 elseif ar > br then
103 order = 1
104 else
105 ai = imag(a)
106 bi = imag(b)
107 if ai < bi then
108 order = -1
109 elseif ai == bi then
110 order = 0
111 else
112 order = 1
113 end
114 end
115endfunction
116
117//
118// assert_close --
119// Returns 1 if the two real matrices computed and expected are close,
120// i.e. if the relative distance between computed and expected is lesser than epsilon.
121// Arguments
122// computed, expected : the two matrices to compare
123// epsilon : a small number
124//
125function flag = assert_close ( computed, expected, epsilon )
126 if expected==0.0 then
127 shift = norm(computed-expected);
128 else
129 shift = norm(computed-expected)/norm(expected);
130 end
131 if shift < epsilon then
132 flag = 1;
133 else
134 flag = 0;
135 end
136 if flag <> 1 then pause,end
137endfunction
138
139function y = sortmyroots(x)
140 // Sort the roots of a polynomial with a customized
141 // complex-aware sorting algorithm.
142 y = sort_merge ( x , compare_complexrealimag );
143endfunction
144
145// Failed on Windows 32 bits
146p = [1,1,-7,-15,1,-4,4,7,4,-53,1,-53,-8,3,3,0,9,-15];
147r = roots(p);
148e = [
149 2.9977242
150 - 2.0998215 + 1.0381514 * %i
151 - 2.0998215 - 1.0381514 * %i
152 - 1.1261224 + 0.7687233 * %i
153 - 1.1261224 - 0.7687233 * %i
154 1.1176579 + 0.5115332 * %i
155 1.1176579 - 0.5115332 * %i
156 - 0.7359417 + 0.3731641 * %i
157 - 0.7359417 - 0.3731641 * %i
158 0.2849638 + 0.9531919 * %i
159 0.2849638 - 0.9531919 * %i
160 0.0897371 + 1.0370037 * %i
161 0.0897371 - 1.0370037 * %i
162 - 0.1740455 + 0.9263179 * %i
163 - 0.1740455 - 0.9263179 * %i
164 0.6447102 + 0.2914081 * %i
165 0.6447102 - 0.2914081 * %i
166 ];
167e = sortmyroots(e);
168r = sortmyroots(r);
169assert_close ( r, e, 1.e-6 );
170
171// Failed on Mac, on Windows
172p=[1,1,-7,-35,1,-4,4,7,4,-88,1,-88,-8,3,3,0,9,-35];
173r = roots(p);
174e = [
175 3.6133489
176 - 2.3323533 + 2.0888127 * %i
177 - 2.3323533 - 2.0888127 * %i
178 1.0856792 + 0.5138318 * %i
179 1.0856792 - 0.5138318 * %i
180 - 1.1030013 + 0.6108696 * %i
181 - 1.1030013 - 0.6108696 * %i
182 0.3226838 + 0.9451270 * %i
183 0.3226838 - 0.9451270 * %i
184 0.0250044 + 1.0210451 * %i
185 0.0250044 - 1.0210451 * %i
186 - 0.2556563 + 0.9467085 * %i
187 - 0.2556563 - 0.9467085 * %i
188 - 0.7512303 + 0.3765797 * %i
189 - 0.7512303 - 0.3765797 * %i
190 0.7021994 + 0.3415821 * %i
191 0.7021994 - 0.3415821 * %i
192];
193e = sortmyroots(e);
194r = sortmyroots(r);
195assert_close ( r, e, 1.e-6 );
196
197// Failed on Linux
198p=[1,1,-7,-80,1,-4,4,7,4,-27,1,-27,-8,3,3,0,9,-80];
199r = roots(p);
200e = [
201 - 2.7595524 + 3.1924496 * %i
202 - 2.7595524 - 3.1924496 * %i
203 4.5006465
204 - 0.9689444 + 0.2683252 * %i
205 - 0.9689444 - 0.2683252 * %i
206 - 0.8111357 + 0.6166997 * %i
207 - 0.8111357 - 0.6166997 * %i
208 - 0.3893539 + 0.9194344 * %i
209 - 0.3893539 - 0.9194344 * %i
210 0.0061369 + 1.0065796 * %i
211 0.0061369 - 1.0065796 * %i
212 0.4195701 + 0.9089127 * %i
213 0.4195701 - 0.9089127 * %i
214 0.9590394 + 0.2589039 * %i
215 0.9590394 - 0.2589039 * %i
216 0.7939168 + 0.5672744 * %i
217 0.7939168 - 0.5672744 * %i
218];
219e = sortmyroots(e);
220r = sortmyroots(r);
221assert_close ( r, e, 1.e-6 );
222
223// Failed on Windows 32 bits
224p=[1,0,1,1,1,-1,-1,1,1,0,1,0,-1,-1,1,-2,0,0,1,-1,1];
225r = roots(p);
226e = [
227 0.5444059 + 1.3082079 * %i
228 0.5444059 - 1.3082079 * %i
229 - 1.0517348 + 0.2347104 * %i
230 - 1.0517348 - 0.2347104 * %i
231 - 0.5893898 + 0.9840032 * %i
232 - 0.5893898 - 0.9840032 * %i
233 - 0.8170407 + 0.5459189 * %i
234 - 0.8170407 - 0.5459189 * %i
235 - 0.6570402 + 0.7150468 * %i
236 - 0.6570402 - 0.7150468 * %i
237 0.0129780 + 0.9748750 * %i
238 0.0129780 - 0.9748750 * %i
239 0.9192290 + 0.4894403 * %i
240 0.9192290 - 0.4894403 * %i
241 0.8691302 + 0.0832523 * %i
242 0.8691302 - 0.0832523 * %i
243 0.4975871 + 0.6807740 * %i
244 0.4975871 - 0.6807740 * %i
245 0.2718754 + 0.7528695 * %i
246 0.2718754 - 0.7528695 * %i
247];
248e = sortmyroots(e);
249r = sortmyroots(r);
250assert_close ( r, e, 1.e-6 );
251
252// A loop on several polynomials
253for i=-100:100
254 if ( modulo(i,20)==0 ) then
255 mprintf("i=%d\n",i);
256 end
257 for j=-100:100
258 p=[1 1 -7 j 1 -4 4 7 4 i 1 i -8 3 3 0 9 j];
259 roots(p);
260 end;
261end;
262
263// A loop on random polynomials.
264// The coefficients are integers
265for i = 1:3000
266 if ( modulo(i,1000)==0 ) then
267 mprintf("i=%d\n",i);
268 end
269 p = [1 round(4*rand(1,20)-2)];
270 roots(p);
271end
diff --git a/scilab/modules/polynomials/tests/unit_tests/poly_roots.dia.ref b/scilab/modules/polynomials/tests/unit_tests/poly_roots.dia.ref
deleted file mode 100644
index ea91886..0000000
--- a/scilab/modules/polynomials/tests/unit_tests/poly_roots.dia.ref
+++ /dev/null
@@ -1,177 +0,0 @@
1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) ????-2008 - INRIA - Michael Baudin
4//
5// This file is distributed under the same license as the Scilab package.
6// =============================================================================
7// <-- JVM NOT MANDATORY -->
8// poly_roots.tst --
9// Check the computation of the roots of a polynomial
10// with different kinds of polynomials and different
11// kinds of roots :
12// - real poly,
13// - complex poly,
14// - real roots,
15// - complex roots.
16//roots : 3 real roots -> RPOLY
17p=-6+11*%s-6*%s^2+%s^3;
18myroots=roots(p);
19computedroots = gsort(myroots);
20expectedroots = [3. ;2. ; 1.];
21if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
22//roots : 3 real roots + polynomials algebra -> RPOLY
23p=-6+11*%s-6*%s^2+%s^3;
24myroots=roots(p+0);
25computedroots = gsort(myroots);
26expectedroots = [3. ;2. ; 1.];
27if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
28//roots : 3 complex roots -> Companion matrix used
29p=-6-%i*6+(11+%i*5)*%s+(-6-%i)*%s^2+%s^3;
30myroots=roots(p);
31computedroots = gsort(myroots);
32expectedroots = [3. ;2. ; 1.+%i];
33if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
34//roots : 3 complex roots + polynomials algebra -> Companion matrix used
35p=-6-%i*6+(11+%i*5)*%s+(-6-%i)*%s^2+%s^3;
36myroots=roots(p+0);
37computedroots = gsort(myroots);
38expectedroots = [3. ;2. ; 1.+%i];
39if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
40// roots : no root at all -> RPOLY
41p=1;v=[];
42if (roots(p)<>v) then bugmes();quit;end
43if (roots(p+0)<>v) then bugmes();quit;end
44//roots : 2 complex roots -> RPOLY
45p=1+%s+%s^2;
46myroots=roots(p);
47computedroots = gsort(myroots);
48expectedroots = [-0.5 + sqrt(3.)/2.*%i ; -0.5 - sqrt(3.)/2.*%i];
49if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
50//roots : 2 roots equals 0 -> RPOLY
51p=%s^2;
52myroots=roots(p);
53computedroots = gsort(myroots);
54expectedroots = [0. ; 0. ];
55if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
56//roots : 3 real roots -> forced to companion matrix with "e" option
57p=-6+11*%s-6*%s^2+%s^3;
58myroots=roots(p,"e");
59computedroots = gsort(myroots);
60expectedroots = [3. ;2. ; 1.];
61if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
62//roots : 2 complex roots -> forced to companion matrix with "e" option
63p=1+%s+%s^2;
64myroots=roots(p,"e");
65computedroots = gsort(myroots);
66expectedroots = [-0.5 + sqrt(3.)/2.*%i ; -0.5 - sqrt(3.)/2.*%i];
67if (abs(computedroots-expectedroots)>400*%eps) then bugmes();quit;end
68// 2 real roots with a zero derivate at the root -> RPOLY
69p=(%s-%pi)^2
70 p =
71
72 2
73 9.8696044 - 6.2831853s + s
74myroots=roots(p);
75computedroots = gsort(myroots);
76expectedroots = [%pi;%pi];
77if (abs(computedroots-expectedroots)>%eps) then bugmes();quit;end
78// 2 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
79p=(%s-%pi)^2
80 p =
81
82 2
83 9.8696044 - 6.2831853s + s
84myroots=roots(p,"e");
85computedroots = gsort(myroots);
86expectedroots = [%pi;%pi];
87if (abs(computedroots-expectedroots)>10*%eps) then bugmes();quit;end
88//
89// Caution !
90// The following are difficult root-finding problems
91// with expected precision problems.
92// See "Principles for testing polynomial
93// zerofinding programs"
94// Jenkins, Traub
95// 1975
96// p.28
97// "The accuracy which one may expect to achieve in calculating
98// zeros is limited by the condition of these zeros. In particular,
99// for multiple zeros perturbations of size epsilon in the
100// coefficients cause perturbations of size epsilon^(1/m)
101// in the zeros."
102//
103//
104// 3 real roots with a zero derivate at the root -> RPOLY
105// *** PRECISION PROBLEM : only simple precision computed, instead of double precision ***
106p=(%s-%pi)^3
107 p =
108
109 2 3
110 - 31.006277 + 29.608813s - 9.424778s + s
111myroots=roots(p);
112computedroots = gsort(myroots);
113expectedroots = [%pi;%pi;%pi];
114if (abs(computedroots-expectedroots)>10^9*%eps) then bugmes();quit;end
115// 3 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
116// *** PRECISION PROBLEM : not even simple precision computed, instead of double precision expected ***
117p=(%s-%pi)^3
118 p =
119
120 2 3
121 - 31.006277 + 29.608813s - 9.424778s + s
122myroots=roots(p,"e");
123computedroots = gsort(myroots);
124expectedroots = [%pi;%pi;%pi];
125if (abs(computedroots-expectedroots)>10^11*%eps) then bugmes();quit;end
126// 4 real roots with a zero derivate at the root -> RPOLY
127// *** PRECISION PROBLEM : only simple precision computed, instead of double precision expected ***
128p=(%s-%pi)^4
129 p =
130
131 2 3 4
132 97.409091 - 124.02511s + 59.217626s - 12.566371s + s
133myroots=roots(p);
134computedroots = gsort(myroots);
135expectedroots = [%pi;%pi;%pi;%pi];
136if (abs(computedroots-expectedroots)>10^9*%eps) then bugmes();quit;end
137// 4 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
138// *** PRECISION PROBLEM : not even simple precision computed, instead of double precision expected ***
139p=(%s-%pi)^4
140 p =
141
142 2 3 4
143 97.409091 - 124.02511s + 59.217626s - 12.566371s + s
144myroots=roots(p,"e");
145computedroots = gsort(myroots);
146expectedroots = [%pi;%pi;%pi;%pi];
147if (abs(computedroots-expectedroots)>10^13*%eps) then bugmes();quit;end
148// 10 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
149// *** PRECISION PROBLEM : only one one true digit ***
150p=(%s-%pi)^10
151 p =
152
153 2 3 4
154 93648.047 - 298090.99s + 426983.9s - 362435.19s + 201891.73s
155 5 6 7 8
156 - 77116.961s + 20455.909s - 3720.7532s + 444.1322s
157 9 10
158 - 31.415927s + s
159myroots=roots(p);
160computedroots = gsort(myroots);
161expectedroots = [%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi];
162if (abs(computedroots-expectedroots)>10^10*%eps) then bugmes();quit;end
163// 10 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
164// *** PRECISION PROBLEM : only one one true digit ***
165p=(%s-%pi)^10
166 p =
167
168 2 3 4
169 93648.047 - 298090.99s + 426983.9s - 362435.19s + 201891.73s
170 5 6 7 8
171 - 77116.961s + 20455.909s - 3720.7532s + 444.1322s
172 9 10
173 - 31.415927s + s
174myroots=roots(p,"e");
175computedroots = gsort(myroots);
176expectedroots = [%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi];
177if (abs(computedroots-expectedroots)>10^15*%eps) then bugmes();quit;end
diff --git a/scilab/modules/polynomials/tests/unit_tests/poly_roots.tst b/scilab/modules/polynomials/tests/unit_tests/poly_roots.tst
deleted file mode 100644
index da35f66..0000000
--- a/scilab/modules/polynomials/tests/unit_tests/poly_roots.tst
+++ /dev/null
@@ -1,145 +0,0 @@
1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) ????-2008 - INRIA - Michael Baudin
4//
5// This file is distributed under the same license as the Scilab package.
6// =============================================================================
7
8// <-- JVM NOT MANDATORY -->
9
10// poly_roots.tst --
11// Check the computation of the roots of a polynomial
12// with different kinds of polynomials and different
13// kinds of roots :
14// - real poly,
15// - complex poly,
16// - real roots,
17// - complex roots.
18
19
20//roots : 3 real roots -> RPOLY
21p=-6+11*%s-6*%s^2+%s^3;
22myroots=roots(p);
23computedroots = gsort(myroots);
24expectedroots = [3. ;2. ; 1.];
25if (abs(computedroots-expectedroots)>400*%eps) then pause,end
26//roots : 3 real roots + polynomials algebra -> RPOLY
27p=-6+11*%s-6*%s^2+%s^3;
28myroots=roots(p+0);
29computedroots = gsort(myroots);
30expectedroots = [3. ;2. ; 1.];
31if (abs(computedroots-expectedroots)>400*%eps) then pause,end
32//roots : 3 complex roots -> Companion matrix used
33p=-6-%i*6+(11+%i*5)*%s+(-6-%i)*%s^2+%s^3;
34myroots=roots(p);
35computedroots = gsort(myroots);
36expectedroots = [3. ;2. ; 1.+%i];
37if (abs(computedroots-expectedroots)>400*%eps) then pause,end
38//roots : 3 complex roots + polynomials algebra -> Companion matrix used
39p=-6-%i*6+(11+%i*5)*%s+(-6-%i)*%s^2+%s^3;
40myroots=roots(p+0);
41computedroots = gsort(myroots);
42expectedroots = [3. ;2. ; 1.+%i];
43if (abs(computedroots-expectedroots)>400*%eps) then pause,end
44// roots : no root at all -> RPOLY
45p=1;v=[];
46if (roots(p)<>v) then pause,end
47if (roots(p+0)<>v) then pause,end
48//roots : 2 complex roots -> RPOLY
49p=1+%s+%s^2;
50myroots=roots(p);
51computedroots = gsort(myroots);
52expectedroots = [-0.5 + sqrt(3.)/2.*%i ; -0.5 - sqrt(3.)/2.*%i];
53if (abs(computedroots-expectedroots)>400*%eps) then pause,end
54//roots : 2 roots equals 0 -> RPOLY
55p=%s^2;
56myroots=roots(p);
57computedroots = gsort(myroots);
58expectedroots = [0. ; 0. ];
59if (abs(computedroots-expectedroots)>400*%eps) then pause,end
60//roots : 3 real roots -> forced to companion matrix with "e" option
61p=-6+11*%s-6*%s^2+%s^3;
62myroots=roots(p,"e");
63computedroots = gsort(myroots);
64expectedroots = [3. ;2. ; 1.];
65if (abs(computedroots-expectedroots)>400*%eps) then pause,end
66//roots : 2 complex roots -> forced to companion matrix with "e" option
67p=1+%s+%s^2;
68myroots=roots(p,"e");
69computedroots = gsort(myroots);
70expectedroots = [-0.5 + sqrt(3.)/2.*%i ; -0.5 - sqrt(3.)/2.*%i];
71if (abs(computedroots-expectedroots)>400*%eps) then pause,end
72// 2 real roots with a zero derivate at the root -> RPOLY
73p=(%s-%pi)^2
74myroots=roots(p);
75computedroots = gsort(myroots);
76expectedroots = [%pi;%pi];
77if (abs(computedroots-expectedroots)>%eps) then pause,end
78// 2 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
79p=(%s-%pi)^2
80myroots=roots(p,"e");
81computedroots = gsort(myroots);
82expectedroots = [%pi;%pi];
83if (abs(computedroots-expectedroots)>10*%eps) then pause,end
84//
85// Caution !
86// The following are difficult root-finding problems
87// with expected precision problems.
88// See "Principles for testing polynomial
89// zerofinding programs"
90// Jenkins, Traub
91// 1975
92// p.28
93// "The accuracy which one may expect to achieve in calculating
94// zeros is limited by the condition of these zeros. In particular,
95// for multiple zeros perturbations of size epsilon in the
96// coefficients cause perturbations of size epsilon^(1/m)
97// in the zeros."
98//
99//
100// 3 real roots with a zero derivate at the root -> RPOLY
101// *** PRECISION PROBLEM : only simple precision computed, instead of double precision ***
102p=(%s-%pi)^3
103myroots=roots(p);
104computedroots = gsort(myroots);
105expectedroots = [%pi;%pi;%pi];
106if (abs(computedroots-expectedroots)>10^9*%eps) then pause,end
107// 3 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
108// *** PRECISION PROBLEM : not even simple precision computed, instead of double precision expected ***
109p=(%s-%pi)^3
110myroots=roots(p,"e");
111computedroots = gsort(myroots);
112expectedroots = [%pi;%pi;%pi];
113if (abs(computedroots-expectedroots)>10^11*%eps) then pause,end
114// 4 real roots with a zero derivate at the root -> RPOLY
115// *** PRECISION PROBLEM : only simple precision computed, instead of double precision expected ***
116p=(%s-%pi)^4
117myroots=roots(p);
118computedroots = gsort(myroots);
119expectedroots = [%pi;%pi;%pi;%pi];
120if (abs(computedroots-expectedroots)>10^9*%eps) then pause,end
121// 4 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
122// *** PRECISION PROBLEM : not even simple precision computed, instead of double precision expected ***
123p=(%s-%pi)^4
124myroots=roots(p,"e");
125computedroots = gsort(myroots);
126expectedroots = [%pi;%pi;%pi;%pi];
127if (abs(computedroots-expectedroots)>10^13*%eps) then pause,end
128// 10 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
129// *** PRECISION PROBLEM : only one one true digit ***
130p=(%s-%pi)^10
131myroots=roots(p);
132computedroots = gsort(myroots);
133expectedroots = [%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi];
134if (abs(computedroots-expectedroots)>10^10*%eps) then pause,end
135// 10 real roots with a zero derivate at the root -> forced to companion matrix with "e" option
136// *** PRECISION PROBLEM : only one one true digit ***
137p=(%s-%pi)^10
138myroots=roots(p,"e");
139computedroots = gsort(myroots);
140expectedroots = [%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi];
141if (abs(computedroots-expectedroots)>10^15*%eps) then pause,end
142
143
144
145
diff --git a/scilab/modules/polynomials/tests/unit_tests/roots.dia.ref b/scilab/modules/polynomials/tests/unit_tests/roots.dia.ref
new file mode 100644
index 0000000..b42a909
--- /dev/null
+++ b/scilab/modules/polynomials/tests/unit_tests/roots.dia.ref
@@ -0,0 +1,624 @@
1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) 2008 - 2009 - INRIA - Michael Baudin
4// Copyright (C) 2011 - DIGITEO - Michael Baudin
5//
6// This file is distributed under the same license as the Scilab package.
7// =============================================================================
8// <-- JVM NOT MANDATORY -->
9//
10// sort_merge_comparison --
11// Returns -1 if x < y,
12// returns 0 if x==y,
13// returns +1 if x > y
14//
15function order = sort_merge_comparison ( x , y )
16 if x < y then
17 order = -1
18 elseif x==y then
19 order = 0
20 else
21 order = 1
22 end
23endfunction
24//
25// sort_merge --
26// Returns the sorted array x.
27// Arguments
28// x : the array to sort
29// compfun : the comparison function
30// Bruno Pincon
31// "quelques tests de rapidit´e entre diff´erents logiciels matriciels"
32// Modified by Michael Baudin to manage a comparison function
33//
34function [x] = sort_merge ( varargin )
35 [lhs,rhs]=argn();
36 if rhs<>1 & rhs<>2 then
37 errmsg = sprintf("Unexpected number of arguments : %d provided while 1 or 2 are expected.",rhs);
38 error(errmsg)
39 end
40 // Get the array x
41 x = varargin(1);
42 // Get the comparison function compfun
43 if rhs==1 then
44 compfun = sort_merge_comparison;
45 else
46 compfun = varargin(2);
47 end
48 // Proceed...
49 n = length(x)
50 if n > 1 then
51 m = floor(n/2);
52 p = n-m
53 x1 = sort_merge ( x(1:m) , compfun )
54 x2 = sort_merge ( x(m+1:n) , compfun )
55 i = 1;
56 i1 = 1;
57 i2 = 1;
58 if ( typeof(compfun)=="list" ) then
59 __compfun_fun__ = compfun(1)
60 end
61 for i = 1:n
62 if ( typeof(compfun)=="function" ) then
63 order = compfun ( x1(i1) , x2(i2) );
64 else
65 order = __compfun_fun__ ( x1(i1) , x2(i2), compfun(2:$) );
66 end
67 if order<=0 then
68 x(i) = x1(i1)
69 i1 = i1+1
70 if (i1 > m) then
71 x(i+1:n) = x2(i2:p)
72 break
73 end
74 else
75 x(i) = x2(i2)
76 i2 = i2+1
77 if (i2 > p) then
78 x(i+1:n) = x1(i1:m)
79 break
80 end
81 end
82 end
83 end
84endfunction
85//
86// compare_complexrealimag --
87// Returns -1 if a < b,
88// returns 0 if a==b,
89// returns +1 if a > b
90// Compare first by real parts, then by imaginary parts.
91//
92function order = compare_complexrealimag ( a , b )
93 ar = real(a)
94 br = real(b)
95 if ar < br then
96 order = -1
97 elseif ar > br then
98 order = 1
99 else
100 ai = imag(a)
101 bi = imag(b)
102 if ai < bi then
103 order = -1
104 elseif ai == bi then
105 order = 0
106 else
107 order = 1
108 end
109 end
110endfunction
111function order = assert_mycompcompl ( varargin )
112 // Copyright (C) 2008-2009 - INRIA - Michael Baudin
113 // Copyright (C) 2010 - 2011 - DIGITEO - Michael Baudin
114 // Compare complex numbers with a tolerance.
115 // order = assert_mycompcompl(a,b)
116 // order = assert_mycompcompl(a,b,rtol)
117 // order = assert_mycompcompl(a,b,rtol,atol)
118 [lhs,rhs]=argn()
119 if ( and(rhs <> [2 3 4] ) ) then
120 errmsg = sprintf ( gettext ( "%s: Wrong number of input arguments: %d to %d expected.") , "assert_mycompcompl" , 2 , 4 )
121 error(errmsg)
122 end
123 //
124 // Get arguments
125 a = varargin(1)
126 b = varargin(2)
127 reltol = argindefault ( rhs , varargin , 3 , sqrt(%eps) )
128 abstol = argindefault ( rhs , varargin , 4 , 0 )
129 //
130 // Check types of variables
131 if ( typeof(a) <> "constant" ) then
132 errmsg = sprintf ( gettext ( "%s: Wrong type for argument %d: Matrix expected.\n") , "assert_mycompcompl" , 1 )
133 error(errmsg)
134 end
135 if ( typeof(a) <> "constant" ) then
136 errmsg = sprintf ( gettext ( "%s: Wrong type for argument %d: Matrix expected.\n") , "assert_mycompcompl" , 2 )
137 error(errmsg)
138 end
139 if ( typeof(reltol) <> "constant" ) then
140 errmsg = sprintf ( gettext ( "%s: Wrong type for argument %d: Matrix expected.\n") , "assert_mycompcompl" , 3 )
141 error(errmsg)
142 end
143 if ( typeof(abstol) <> "constant" ) then
144 errmsg = sprintf ( gettext ( "%s: Wrong type for argument %d: Matrix expected.\n") , "assert_mycompcompl" , 4 )
145 error(errmsg)
146 end
147 //
148 // Check sizes of variables
149 if ( size(a,"*") <> 1 ) then
150 errmsg = sprintf ( gettext ( "%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n") , "assert_mycompcompl" , 1 , 1 , 1 )
151 error(errmsg)
152 end
153 if ( size(b,"*") <> 1 ) then
154 errmsg = sprintf ( gettext ( "%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n") , "assert_mycompcompl" , 2 , 1 , 1 )
155 error(errmsg)
156 end
157 if ( size(reltol,"*") <> 1 ) then
158 errmsg = sprintf ( gettext ( "%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n") , "assert_mycompcompl" , 3 , 1 , 1 )
159 error(errmsg)
160 end
161 if ( size(abstol,"*") <> 1 ) then
162 errmsg = sprintf ( gettext ( "%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n") , "assert_mycompcompl" , 4 , 1 , 1 )
163 error(errmsg)
164 end
165 //
166 // Check values of variables
167 if ( reltol < 0 ) then
168 errmsg = sprintf ( gettext ( "%s: Wrong value for input argument #%d: Must be > %d.\n") , "assert_mycompcompl" , 3 , 0 )
169 error(errmsg)
170 end
171 if ( abstol < 0 ) then
172 errmsg = sprintf ( gettext ( "%s: Wrong value for input argument #%d: Must be > %d.\n") , "assert_mycompcompl" , 4 , 0 )
173 error(errmsg)
174 end
175 //
176 ord_re = assert_compdata ( real(a) , real(b) , reltol , abstol )
177 if ( ord_re == 0 ) then
178 // Tie on the real part: compare imaginary parts
179 ord_im = assert_compdata ( imag(a) , imag(b) , reltol , abstol )
180 if ( ord_im == 0 ) then
181 // Tie on imaginary parts too: two numbers are "equal"
182 order = 0
183 elseif ( ord_im == -1 ) then
184 order = -1
185 else
186 order = 1
187 end
188 elseif ( ord_re == -1 ) then
189 order = -1
190 else
191 order = 1
192 end
193endfunction
194function order = assert_compdata ( a , b , reltol , abstol )
195 if ( a == %inf ) then
196 if ( isnan(b) ) then
197 order = -1
198 elseif ( b == %inf ) then
199 order = 0
200 else
201 order = 1
202 end
203 elseif ( a == -%inf ) then
204 if ( b == -%inf ) then
205 order = 0
206 else
207 order = -1
208 end
209 elseif ( isnan(a) ) then
210 if ( isnan(b) ) then
211 order = 0
212 else
213 order = 1
214 end
215 else
216 if ( isnan(b) ) then
217 order = -1
218 elseif ( b == -%inf ) then
219 order = 1
220 elseif ( b == %inf ) then
221 order = -1
222 else
223 areequal = abs ( a - b ) <= reltol * max ( abs(a) , abs(b) ) + abstol
224 if ( areequal ) then
225 order = 0
226 elseif ( a < b ) then
227 order = -1
228 else
229 order = 1
230 end
231 end
232 end
233endfunction
234function argin = argindefault ( rhs , vararglist , ivar , default )
235 // Returns the value of the input argument #ivar.
236 // If this argument was not provided, or was equal to the
237 // empty matrix, returns the default value.
238 if ( rhs < ivar ) then
239 argin = default
240 else
241 if ( vararglist(ivar) <> [] ) then
242 argin = vararglist(ivar)
243 else
244 argin = default
245 end
246 end
247endfunction
248function y = sortmyroots(varargin)
249 // Sort the roots of a polynomial with a customized
250 // complex-aware sorting algorithm.
251 // y = sortmyroots(x)
252 // y = sortmyroots(x,rtol)
253 // y = sortmyroots(x,rtol,atol)
254 [lhs,rhs]=argn()
255 if ( and(rhs <> [1 2 3] ) ) then
256 errmsg = sprintf ( gettext ( "%s: Wrong number of input arguments: %d to %d expected.") , "sortmyroots" , 1 , 3 )
257 error(errmsg)
258 end
259 x = varargin(1)
260 reltol = argindefault ( rhs , varargin , 2 , sqrt(%eps) )
261 abstol = argindefault ( rhs , varargin , 3 , 0 )
262 y = sort_merge ( x , list(assert_mycompcompl,reltol,abstol) );
263endfunction
264//
265// assert_close --
266// Returns 1 if the two real matrices computed and expected are close,
267// i.e. if the relative distance between computed and expected is lesser than epsilon.
268// Arguments
269// computed, expected : the two matrices to compare
270// epsilon : a small number
271//
272function flag = assert_close ( computed, expected, epsilon )
273 if expected==0.0 then
274 shift = norm(computed-expected);
275 else
276 den = norm(expected)
277 if ( den==0 ) then
278 shift = norm(computed-expected);
279 else
280 shift = norm(computed-expected)/den;
281 end
282 end
283 if shift <= epsilon then
284 flag = 1;
285 else
286 flag = 0;
287 end
288 if flag <> 1 then bugmes();quit;end
289endfunction
290//
291// assert_equal --
292// Returns 1 if the two real matrices computed and expected are equal.
293// Arguments
294// computed, expected : the two matrices to compare
295// epsilon : a small number
296//
297function flag = assert_equal ( computed , expected )
298 if computed==expected then
299 flag = 1;
300 else
301 flag = 0;
302 end
303 if flag <> 1 then bugmes();quit;end
304endfunction
305function checkroots(p,expectedroots,rtol)
306 // Checks the roots function against given roots.
307 //
308 // 1. Check default algorithm
309 myroots=roots(p);
310 computedroots = sortmyroots(myroots);
311 expectedroots = sortmyroots(expectedroots);
312 assert_close(computedroots,expectedroots,rtol);
313 //
314 // 2. Check "e" algorithm
315 myroots=roots(p,"e");
316 computedroots = sortmyroots(myroots);
317 expectedroots = sortmyroots(expectedroots);
318 assert_close(computedroots,expectedroots,rtol);
319 //
320 // 3. Check "f" algorithm
321 if ( isreal(p) ) then
322 myroots=roots(p,"f");
323 computedroots = sortmyroots(myroots);
324 expectedroots = sortmyroots(expectedroots);
325 assert_close(computedroots,expectedroots,rtol);
326 end
327endfunction
328// Check the computation of the roots of a polynomial
329// with different kinds of polynomials and different
330// kinds of roots :
331// - real poly,
332// - complex poly,
333// - real roots,
334// - complex roots.
335//roots : 3 real roots
336p=-6+11*%s-6*%s^2+%s^3;
337expectedroots = [1; 2; 3];
338checkroots(p,expectedroots,400*%eps);
339//roots : 3 real roots + polynomials algebra
340p=-6+11*%s-6*%s^2+%s^3;
341q = p+0;
342expectedroots = [1; 2; 3];
343checkroots(q,expectedroots,400*%eps);
344//roots : 3 complex roots
345p=-6-%i*6+(11+%i*5)*%s+(-6-%i)*%s^2+%s^3;
346expectedroots = [1+%i; 2 ; 3];
347checkroots(p,expectedroots,400*%eps);
348//roots : 3 complex roots + polynomials algebra
349p=-6-%i*6+(11+%i*5)*%s+(-6-%i)*%s^2+%s^3;
350q = p+0;
351expectedroots = [1+%i; 2 ; 3];
352checkroots(p,expectedroots,400*%eps);
353// roots : no root at all
354p=1;
355v=[];
356checkroots(p,[],0);
357q = p+0;
358checkroots(q,[],0);
359//roots : 2 complex roots
360p=1+%s+%s^2;
361expectedroots = [-0.5 - sqrt(3.)/2.*%i; -0.5 + sqrt(3.)/2.*%i ];
362checkroots(p,expectedroots,400*%eps);
363//roots : 2 roots equals 0
364p=%s^2;
365expectedroots = [0. ; 0. ];
366// 2 real roots with a zero derivate at the root
367p=(%s-%pi)^2;
368expectedroots = [%pi;%pi];
369checkroots(p,expectedroots,400*%eps);
370//
371// Caution !
372// The following are difficult root-finding problems
373// with expected precision problems.
374// See "Principles for testing polynomial
375// zerofinding programs"
376// Jenkins, Traub
377// 1975
378// p.28
379// "The accuracy which one may expect to achieve in calculating
380// zeros is limited by the condition of these zeros. In particular,
381// for multiple zeros perturbations of size epsilon in the
382// coefficients cause perturbations of size epsilon^(1/m)
383// in the zeros."
384//
385//
386// 3 real roots with a zero derivate at the root
387// Really difficult problem : only simple precision computed, instead of double precision ***
388p=(%s-%pi)^3;
389expectedroots = [%pi;%pi;%pi];
390checkroots(p,expectedroots,10^11*%eps);
391// 4 real roots with a zero derivate at the root
392// Really difficult problem : only simple precision
393p=(%s-%pi)^4;
394expectedroots = [%pi;%pi;%pi;%pi];
395checkroots(p,expectedroots,10^12*%eps);
396// 10 real roots with a zero derivate at the root
397// Really difficult problem : only one correct digit
398p=(%s-%pi)^10;
399expectedroots = [%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi];
400checkroots(p,expectedroots,10^15*%eps);
401// "Numerical computing with Matlab", Cleve Moler.
402A = diag(1:20);
403p = poly(A,'x');
404e = [1:20]';
405checkroots(p,e,1.e-2);
406// Tests from CPOLY
407// M. A. Jenkins and J. F. Traub. 1972.
408// Algorithm 419: zeros of a complex polynomial.
409// Commun. ACM 15, 2 (February 1972), 97-99.
410//
411// EXAMPLE 1. POLYNOMIAL WITH ZEROS 1,2,...,10.
412P=[];
413PI=[];
414P(1)=1;
415P(2)=-55;
416P(3)=1320;
417P(4)=-18150;
418P(5)=157773;
419P(6)=-902055;
420P(7) = 3416930;
421P(8)=-8409500;
422P(9)=12753576;
423P(10)=-10628640;
424P(11)=3628800;
425PI(1:11) = 0;
426P = complex(P,PI);
427E = (1:10)';
428R = roots(P);
429E = sortmyroots(E);
430R = sortmyroots(R);
431assert_close ( R, E, 1.e-10 );
432// EXAMPLE 2. ZEROS ON IMAGINARY AXIS DEGREE 3.
433// x^3-10001.0001*i*x^2-10001.0001*x+i
434P = [];
435PI=[];
436P(1)=1;
437P(2)=0;
438P(3)=-10001.0001;
439P(4)=0;
440PI(1)=0;
441PI(2)=-10001.0001;
442PI(3)=0;
443PI(4)=1;
444P = complex(P,PI);
445E = [
4460.0001*%i
447%i
44810000*%i
449];
450R = roots(P);
451E = sortmyroots(E,[],1.e-10);
452R = sortmyroots(R,[],1.e-10);
453assert_close ( R, E, 1.e-14 );
454// plot(real(R),imag(R),"bo")
455// xtitle("Roots","Real","Imaginary")
456// EXAMPLE 3. ZEROS AT 1+I,1/2*(1+I)....1/(2**-9)*(1+I)
457P = [];
458PI=[];
459P(1)=1.0;
460P(2)=-1.998046875;
461P(3)=0.0;
462P(4)=.7567065954208374D0;
463P(5)=-.2002119533717632D0;
464P(6)=1.271507365163416D-2;
465P(7)=0;
466P(8)=-1.154642632172909D-5;
467P(9)=1.584803612786345D-7;
468P(10)=-4.652065399568528D-10;
469P(11)=0;
470PI(1)=0;
471PI(2)=P(2);
472PI(3)=2.658859252929688D0;
473PI(4)=-7.567065954208374D-1;
474PI(5)=0;
475PI(6)=P(6);
476PI(7)=-7.820779428584501D-4;
477PI(8)=-P(8);
478PI(9)=0;
479PI(10)=P(10);
480PI(11)=9.094947017729282D-13;
481P = complex(P,PI);
482R = roots(P)
483 R =
484
485 1. + i
486 0.5 + 0.5i
487 0.25 + 0.25i
488 0.125 + 0.125i
489 0.0625 + 0.0625i
490 0.03125 + 0.03125i
491 0.015625 + 0.015625i
492 0.0078125 + 0.0078125i
493 0.0039062 + 0.0039062i
494 0.0019531 + 0.0019531i
495E = (1+%i)*2.^((0:-1:-9)');
496E = sortmyroots(E);
497R = sortmyroots(R);
498assert_close ( R, E, 1.e-14 );
499// EXAMPLE 4. MULTIPLE ZEROS
500// Real part:
501// 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
502// Imaginary part:
503// 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
504P = [];
505PI=[];
506P(1)=1;
507P(2)=-10;
508P(3)=3;
509P(4)=284;
510P(5)=-1293;
511P(6)=2374;
512P(7)=-1587;
513P(8)=-920;
514P(9)=2204;
515P(10)=-1344;
516P(11)=288;
517PI(1)=0;
518PI(2)=-10;
519PI(3)=100;
520PI(4)=-334;
521PI(5)=200;
522PI(6)=1394;
523PI(7) =-3836;
524PI(8)=4334;
525PI(9)=-2352;
526PI(10)=504;
527PI(11)=0;
528P = complex(P,PI);
529R = roots(P)
530 R =
531
532 - 1.776D-15 + 4.i
533 3.0000002 + 0.0000002i
534 2.9999998 - 0.0000002i
535 0.0000285 + 1.9999944i
536 - 0.0000094 + 2.0000275i
537 - 0.0000192 + 1.9999781i
538 0.9997371 + 0.0001462i
539 0.9998538 - 0.0002628i
540 1.0001463 + 0.0002629i
541 1.0002628 - 0.0001464i
542E = [
5431
5441
5451
5461
5472*%i
5482*%i
5492*%i
5503
5513
5524*%i
553];
554E = sortmyroots(E,[],1.e-3);
555R = sortmyroots(R,[],1.e-3);
556assert_close ( R, E, 1.e-4 );
557// EXAMPLE 5. 12 ZEROS EVENLY DISTRIBUTE ON A CIRCLE OF
558// RADIUS 1. CENTERED AT 0+2I
559// Real part:
560// 4095 - 67584*x^2 + 126720*x^4 - 59136*x^6 + 7920*x^8 - 264*x^10 + x^12
561// Imaginary part:
562// 24576*x - 112640x^3 + 101376x^5 - 25344x^7 + 1760x^9 - 24x^11
563P = [];
564PI=[];
565P(1)=1;
566P(2)=0;
567P(3)=-264;
568P(4)=0;
569P(5)=7920;
570P(6)=0;
571P(7)=-59136;
572P(8)=0;
573P(9)=126720;
574P(10)=0;
575P(11)=-67584;
576P(12)=0;
577P(13)=4095;
578PI(1)=0;
579PI(2)=-24;
580PI(3)=0;
581PI(4)=1760;
582PI(5)=0;
583PI(6)=-25344;
584PI(7)=0;
585PI(8)=101376;
586PI(9)=0;
587PI(10)=-112640;
588PI(11)=0;
589PI(12)=24576;
590PI(13)=0;
591P = complex(P,PI);
592R = roots(P)
593 R =
594
595 0.5 + 2.8660254i
596 0.8660254 + 2.5i
597 4.012D-11 + 3.i
598 - 0.5 + 2.8660254i
599 - 0.8660254 + 2.5i
600 1. + 2.i
601 - 1. + 2.i
602 0.8660254 + 1.5i
603 - 0.8660254 + 1.5i
604 0.5 + 1.1339746i
605 - 0.5 + 1.1339746i
606 - 4.367D-16 + i
607S3=sqrt(3);
608E = [
609-1 + 2*%i
610%i
6113*%i
6121+2*%i
613(1/2)*(-S3+3*%i)
614(1/2)*(-S3+5*%i)
615-(1/2)*%i*(S3+(-4-%i))
616(1/2)*((1+4*%i)-%i*S3)
617(1/2)*%i*(S3+(4+%i))
618(1/2)*((1+4*%i)+%i*S3)
619(1/2)*(S3+3*%i)
620(1/2)*(S3+5*%i)
621];
622E = sortmyroots(E);
623R = sortmyroots(R);
624assert_close ( R, E, 1.e-8 );
diff --git a/scilab/modules/polynomials/tests/unit_tests/roots.tst b/scilab/modules/polynomials/tests/unit_tests/roots.tst
new file mode 100644
index 0000000..6f7b574
--- /dev/null
+++ b/scilab/modules/polynomials/tests/unit_tests/roots.tst
@@ -0,0 +1,609 @@
1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) 2008 - 2009 - INRIA - Michael Baudin
4// Copyright (C) 2011 - DIGITEO - Michael Baudin
5//
6// This file is distributed under the same license as the Scilab package.
7// =============================================================================
8
9// <-- JVM NOT MANDATORY -->
10
11//
12// sort_merge_comparison --
13// Returns -1 if x < y,
14// returns 0 if x==y,
15// returns +1 if x > y
16//
17function order = sort_merge_comparison ( x , y )
18 if x < y then
19 order = -1
20 elseif x==y then
21 order = 0
22 else
23 order = 1
24 end
25endfunction
26
27//
28// sort_merge --
29// Returns the sorted array x.
30// Arguments
31// x : the array to sort
32// compfun : the comparison function
33// Bruno Pincon
34// "quelques tests de rapidit´e entre diff´erents logiciels matriciels"
35// Modified by Michael Baudin to manage a comparison function
36//
37function [x] = sort_merge ( varargin )
38 [lhs,rhs]=argn();
39 if rhs<>1 & rhs<>2 then
40 errmsg = sprintf("Unexpected number of arguments : %d provided while 1 or 2 are expected.",rhs);
41 error(errmsg)
42 end
43 // Get the array x
44 x = varargin(1);
45 // Get the comparison function compfun
46 if rhs==1 then
47 compfun = sort_merge_comparison;
48 else
49 compfun = varargin(2);
50 end
51 // Proceed...
52 n = length(x)
53 if n > 1 then
54 m = floor(n/2);
55 p = n-m
56 x1 = sort_merge ( x(1:m) , compfun )
57 x2 = sort_merge ( x(m+1:n) , compfun )
58 i = 1;
59 i1 = 1;
60 i2 = 1;
61 if ( typeof(compfun)=="list" ) then
62 __compfun_fun__ = compfun(1)
63 end
64 for i = 1:n
65 if ( typeof(compfun)=="function" ) then
66 order = compfun ( x1(i1) , x2(i2) );
67 else
68 order = __compfun_fun__ ( x1(i1) , x2(i2), compfun(2:$) );
69 end
70 if order<=0 then
71 x(i) = x1(i1)
72 i1 = i1+1
73 if (i1 > m) then
74 x(i+1:n) = x2(i2:p)
75 break
76 end
77 else
78 x(i) = x2(i2)
79 i2 = i2+1
80 if (i2 > p) then
81 x(i+1:n) = x1(i1:m)
82 break
83 end
84 end
85 end
86 end
87endfunction
88
89//
90// compare_complexrealimag --
91// Returns -1 if a < b,
92// returns 0 if a==b,
93// returns +1 if a > b
94// Compare first by real parts, then by imaginary parts.
95//
96function order = compare_complexrealimag ( a , b )
97 ar = real(a)
98 br = real(b)
99 if ar < br then
100 order = -1
101 elseif ar > br then
102 order = 1
103 else
104 ai = imag(a)
105 bi = imag(b)
106 if ai < bi then
107 order = -1
108 elseif ai == bi then
109 order = 0
110 else
111 order = 1
112 end
113 end
114endfunction
115
116function order = assert_mycompcompl ( varargin )
117 // Copyright (C) 2008-2009 - INRIA - Michael Baudin
118 // Copyright (C) 2010 - 2011 - DIGITEO - Michael Baudin
119 // Compare complex numbers with a tolerance.
120 // order = assert_mycompcompl(a,b)
121 // order = assert_mycompcompl(a,b,rtol)
122 // order = assert_mycompcompl(a,b,rtol,atol)
123
124 [lhs,rhs]=argn()
125 if ( and(rhs <> [2 3 4] ) ) then
126 errmsg = sprintf ( gettext ( "%s: Wrong number of input arguments: %d to %d expected.") , "assert_mycompcompl" , 2 , 4 )
127 error(errmsg)
128 end
129 //
130 // Get arguments
131 a = varargin(1)
132 b = varargin(2)
133 reltol = argindefault ( rhs , varargin , 3 , sqrt(%eps) )
134 abstol = argindefault ( rhs , varargin , 4 , 0 )
135 //
136 // Check types of variables
137 if ( typeof(a) <> "constant" ) then
138 errmsg = sprintf ( gettext ( "%s: Wrong type for argument %d: Matrix expected.\n") , "assert_mycompcompl" , 1 )
139 error(errmsg)
140 end
141 if ( typeof(a) <> "constant" ) then
142 errmsg = sprintf ( gettext ( "%s: Wrong type for argument %d: Matrix expected.\n") , "assert_mycompcompl" , 2 )
143 error(errmsg)
144 end
145 if ( typeof(reltol) <> "constant" ) then
146 errmsg = sprintf ( gettext ( "%s: Wrong type for argument %d: Matrix expected.\n") , "assert_mycompcompl" , 3 )
147 error(errmsg)
148 end
149 if ( typeof(abstol) <> "constant" ) then
150 errmsg = sprintf ( gettext ( "%s: Wrong type for argument %d: Matrix expected.\n") , "assert_mycompcompl" , 4 )
151 error(errmsg)
152 end
153 //
154 // Check sizes of variables
155 if ( size(a,"*") <> 1 ) then
156 errmsg = sprintf ( gettext ( "%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n") , "assert_mycompcompl" , 1 , 1 , 1 )
157 error(errmsg)
158 end
159 if ( size(b,"*") <> 1 ) then
160 errmsg = sprintf ( gettext ( "%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n") , "assert_mycompcompl" , 2 , 1 , 1 )
161 error(errmsg)
162 end
163 if ( size(reltol,"*") <> 1 ) then
164 errmsg = sprintf ( gettext ( "%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n") , "assert_mycompcompl" , 3 , 1 , 1 )
165 error(errmsg)
166 end
167 if ( size(abstol,"*") <> 1 ) then
168 errmsg = sprintf ( gettext ( "%s: Wrong size for input argument #%d: %d-by-%d matrix expected.\n") , "assert_mycompcompl" , 4 , 1 , 1 )
169 error(errmsg)
170 end
171 //
172 // Check values of variables
173 if ( reltol < 0 ) then
174 errmsg = sprintf ( gettext ( "%s: Wrong value for input argument #%d: Must be > %d.\n") , "assert_mycompcompl" , 3 , 0 )
175 error(errmsg)
176 end
177 if ( abstol < 0 ) then
178 errmsg = sprintf ( gettext ( "%s: Wrong value for input argument #%d: Must be > %d.\n") , "assert_mycompcompl" , 4 , 0 )
179 error(errmsg)
180 end
181 //
182 ord_re = assert_compdata ( real(a) , real(b) , reltol , abstol )
183 if ( ord_re == 0 ) then
184 // Tie on the real part: compare imaginary parts
185 ord_im = assert_compdata ( imag(a) , imag(b) , reltol , abstol )
186 if ( ord_im == 0 ) then
187 // Tie on imaginary parts too: two numbers are "equal"
188 order = 0
189 elseif ( ord_im == -1 ) then
190 order = -1
191 else
192 order = 1
193 end
194 elseif ( ord_re == -1 ) then
195 order = -1
196 else
197 order = 1
198 end
199endfunction
200
201function order = assert_compdata ( a , b , reltol , abstol )
202 if ( a == %inf ) then
203 if ( isnan(b) ) then
204 order = -1
205 elseif ( b == %inf ) then
206 order = 0
207 else
208 order = 1
209 end
210 elseif ( a == -%inf ) then
211 if ( b == -%inf ) then
212 order = 0
213 else
214 order = -1
215 end
216 elseif ( isnan(a) ) then
217 if ( isnan(b) ) then
218 order = 0
219 else
220 order = 1
221 end
222 else
223 if ( isnan(b) ) then
224 order = -1
225 elseif ( b == -%inf ) then
226 order = 1
227 elseif ( b == %inf ) then
228 order = -1
229 else
230 areequal = abs ( a - b ) <= reltol * max ( abs(a) , abs(b) ) + abstol
231 if ( areequal ) then
232 order = 0
233 elseif ( a < b ) then
234 order = -1
235 else
236 order = 1
237 end
238 end
239 end
240endfunction
241
242
243function argin = argindefault ( rhs , vararglist , ivar , default )
244 // Returns the value of the input argument #ivar.
245 // If this argument was not provided, or was equal to the
246 // empty matrix, returns the default value.
247 if ( rhs < ivar ) then
248 argin = default
249 else
250 if ( vararglist(ivar) <> [] ) then
251 argin = vararglist(ivar)
252 else
253 argin = default
254 end
255 end
256endfunction
257
258function y = sortmyroots(varargin)
259 // Sort the roots of a polynomial with a customized
260 // complex-aware sorting algorithm.
261 // y = sortmyroots(x)
262 // y = sortmyroots(x,rtol)
263 // y = sortmyroots(x,rtol,atol)
264 [lhs,rhs]=argn()
265 if ( and(rhs <> [1 2 3] ) ) then
266 errmsg = sprintf ( gettext ( "%s: Wrong number of input arguments: %d to %d expected.") , "sortmyroots" , 1 , 3 )
267 error(errmsg)
268 end
269 x = varargin(1)
270 reltol = argindefault ( rhs , varargin , 2 , sqrt(%eps) )
271 abstol = argindefault ( rhs , varargin , 3 , 0 )
272 y = sort_merge ( x , list(assert_mycompcompl,reltol,abstol) );
273endfunction
274
275//
276// assert_close --
277// Returns 1 if the two real matrices computed and expected are close,
278// i.e. if the relative distance between computed and expected is lesser than epsilon.
279// Arguments
280// computed, expected : the two matrices to compare
281// epsilon : a small number
282//
283function flag = assert_close ( computed, expected, epsilon )
284 if expected==0.0 then
285 shift = norm(computed-expected);
286 else
287 den = norm(expected)
288 if ( den==0 ) then
289 shift = norm(computed-expected);
290 else
291 shift = norm(computed-expected)/den;
292 end
293 end
294 if shift <= epsilon then
295 flag = 1;
296 else
297 flag = 0;
298 end
299 if flag <> 1 then pause,end
300endfunction
301
302//
303// assert_equal --
304// Returns 1 if the two real matrices computed and expected are equal.
305// Arguments
306// computed, expected : the two matrices to compare
307// epsilon : a small number
308//
309function flag = assert_equal ( computed , expected )
310 if computed==expected then
311 flag = 1;
312 else
313 flag = 0;
314 end
315 if flag <> 1 then pause,end
316endfunction
317
318function checkroots(p,expectedroots,rtol)
319 // Checks the roots function against given roots.
320 //
321 // 1. Check default algorithm
322 myroots=roots(p);
323 computedroots = sortmyroots(myroots);
324 expectedroots = sortmyroots(expectedroots);
325 assert_close(computedroots,expectedroots,rtol);
326 //
327 // 2. Check "e" algorithm
328 myroots=roots(p,"e");
329 computedroots = sortmyroots(myroots);
330 expectedroots = sortmyroots(expectedroots);
331 assert_close(computedroots,expectedroots,rtol);
332 //
333 // 3. Check "f" algorithm
334 if ( isreal(p) ) then
335 myroots=roots(p,"f");
336 computedroots = sortmyroots(myroots);
337 expectedroots = sortmyroots(expectedroots);
338 assert_close(computedroots,expectedroots,rtol);
339 end
340endfunction
341
342// Check the computation of the roots of a polynomial
343// with different kinds of polynomials and different
344// kinds of roots :
345// - real poly,
346// - complex poly,
347// - real roots,
348// - complex roots.
349
350//roots : 3 real roots
351p=-6+11*%s-6*%s^2+%s^3;
352expectedroots = [1; 2; 3];
353checkroots(p,expectedroots,400*%eps);
354//roots : 3 real roots + polynomials algebra
355p=-6+11*%s-6*%s^2+%s^3;
356q = p+0;
357expectedroots = [1; 2; 3];
358checkroots(q,expectedroots,400*%eps);
359//roots : 3 complex roots
360p=-6-%i*6+(11+%i*5)*%s+(-6-%i)*%s^2+%s^3;
361expectedroots = [1+%i; 2 ; 3];
362checkroots(p,expectedroots,400*%eps);
363//roots : 3 complex roots + polynomials algebra
364p=-6-%i*6+(11+%i*5)*%s+(-6-%i)*%s^2+%s^3;
365q = p+0;
366expectedroots = [1+%i; 2 ; 3];
367checkroots(p,expectedroots,400*%eps);
368// roots : no root at all
369p=1;
370v=[];
371checkroots(p,[],0);
372q = p+0;
373checkroots(q,[],0);
374//roots : 2 complex roots
375p=1+%s+%s^2;
376expectedroots = [-0.5 - sqrt(3.)/2.*%i; -0.5 + sqrt(3.)/2.*%i ];
377checkroots(p,expectedroots,400*%eps);
378//roots : 2 roots equals 0
379p=%s^2;
380expectedroots = [0. ; 0. ];
381// 2 real roots with a zero derivate at the root
382p=(%s-%pi)^2;
383expectedroots = [%pi;%pi];
384checkroots(p,expectedroots,400*%eps);
385//
386// Caution !
387// The following are difficult root-finding problems
388// with expected precision problems.
389// See "Principles for testing polynomial
390// zerofinding programs"
391// Jenkins, Traub
392// 1975
393// p.28
394// "The accuracy which one may expect to achieve in calculating
395// zeros is limited by the condition of these zeros. In particular,
396// for multiple zeros perturbations of size epsilon in the
397// coefficients cause perturbations of size epsilon^(1/m)
398// in the zeros."
399//
400//
401// 3 real roots with a zero derivate at the root
402// Really difficult problem : only simple precision computed, instead of double precision ***
403p=(%s-%pi)^3;
404expectedroots = [%pi;%pi;%pi];
405checkroots(p,expectedroots,10^11*%eps);
406// 4 real roots with a zero derivate at the root
407// Really difficult problem : only simple precision
408p=(%s-%pi)^4;
409expectedroots = [%pi;%pi;%pi;%pi];
410checkroots(p,expectedroots,10^12*%eps);
411// 10 real roots with a zero derivate at the root
412// Really difficult problem : only one correct digit
413p=(%s-%pi)^10;
414expectedroots = [%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi];
415checkroots(p,expectedroots,10^15*%eps);
416
417// "Numerical computing with Matlab", Cleve Moler.
418A = diag(1:20);
419p = poly(A,'x');
420e = [1:20]';
421checkroots(p,e,1.e-2);
422
423// Tests from CPOLY
424// M. A. Jenkins and J. F. Traub. 1972.
425// Algorithm 419: zeros of a complex polynomial.
426// Commun. ACM 15, 2 (February 1972), 97-99.
427//
428// EXAMPLE 1. POLYNOMIAL WITH ZEROS 1,2,...,10.
429P=[];
430PI=[];
431P(1)=1;
432P(2)=-55;
433P(3)=1320;
434P(4)=-18150;
435P(5)=157773;
436P(6)=-902055;
437P(7) = 3416930;
438P(8)=-8409500;
439P(9)=12753576;
440P(10)=-10628640;
441P(11)=3628800;
442PI(1:11) = 0;
443P = complex(P,PI);
444E = (1:10)';
445R = roots(P);
446E = sortmyroots(E);
447R = sortmyroots(R);
448assert_close ( R, E, 1.e-10 );
449
450// EXAMPLE 2. ZEROS ON IMAGINARY AXIS DEGREE 3.
451// x^3-10001.0001*i*x^2-10001.0001*x+i
452P = [];
453PI=[];
454P(1)=1;
455P(2)=0;
456P(3)=-10001.0001;
457P(4)=0;
458PI(1)=0;
459PI(2)=-10001.0001;
460PI(3)=0;
461PI(4)=1;
462P = complex(P,PI);
463E = [
4640.0001*%i
465%i
46610000*%i
467];
468R = roots(P);
469E = sortmyroots(E,[],1.e-10);
470R = sortmyroots(R,[],1.e-10);
471assert_close ( R, E, 1.e-14 );
472
473// plot(real(R),imag(R),"bo")
474// xtitle("Roots","Real","Imaginary")
475
476// EXAMPLE 3. ZEROS AT 1+I,1/2*(1+I)....1/(2**-9)*(1+I)
477P = [];
478PI=[];
479P(1)=1.0;
480P(2)=-1.998046875;
481P(3)=0.0;
482P(4)=.7567065954208374D0;
483P(5)=-.2002119533717632D0;
484P(6)=1.271507365163416D-2;
485P(7)=0;
486P(8)=-1.154642632172909D-5;
487P(9)=1.584803612786345D-7;
488P(10)=-4.652065399568528D-10;
489P(11)=0;
490PI(1)=0;
491PI(2)=P(2);
492PI(3)=2.658859252929688D0;
493PI(4)=-7.567065954208374D-1;
494PI(5)=0;
495PI(6)=P(6);
496PI(7)=-7.820779428584501D-4;
497PI(8)=-P(8);
498PI(9)=0;
499PI(10)=P(10);
500PI(11)=9.094947017729282D-13;
501P = complex(P,PI);
502R = roots(P)
503E = (1+%i)*2.^((0:-1:-9)');
504E = sortmyroots(E);
505R = sortmyroots(R);
506assert_close ( R, E, 1.e-14 );
507
508// EXAMPLE 4. MULTIPLE ZEROS
509// Real part:
510// 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
511// Imaginary part:
512// 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
513P = [];
514PI=[];
515P(1)=1;
516P(2)=-10;
517P(3)=3;
518P(4)=284;
519P(5)=-1293;
520P(6)=2374;
521P(7)=-1587;
522P(8)=-920;
523P(9)=2204;
524P(10)=-1344;
525P(11)=288;
526PI(1)=0;
527PI(2)=-10;
528PI(3)=100;
529PI(4)=-334;
530PI(5)=200;
531PI(6)=1394;
532PI(7) =-3836;
533PI(8)=4334;
534PI(9)=-2352;
535PI(10)=504;
536PI(11)=0;
537P = complex(P,PI);
538R = roots(P)
539E = [
5401
5411
5421
5431
5442*%i
5452*%i
5462*%i
5473
5483
5494*%i
550];
551E = sortmyroots(E,[],1.e-3);
552R = sortmyroots(R,[],1.e-3);
553assert_close ( R, E, 1.e-4 );
554
555// EXAMPLE 5. 12 ZEROS EVENLY DISTRIBUTE ON A CIRCLE OF
556// RADIUS 1. CENTERED AT 0+2I
557// Real part:
558// 4095 - 67584*x^2 + 126720*x^4 - 59136*x^6 + 7920*x^8 - 264*x^10 + x^12
559// Imaginary part:
560// 24576*x - 112640x^3 + 101376x^5 - 25344x^7 + 1760x^9 - 24x^11
561
562P = [];
563PI=[];
564P(1)=1;
565P(2)=0;
566P(3)=-264;
567P(4)=0;
568P(5)=7920;
569P(6)=0;
570P(7)=-59136;
571P(8)=0;
572P(9)=126720;
573P(10)=0;
574P(11)=-67584;
575P(12)=0;
576P(13)=4095;
577PI(1)=0;
578PI(2)=-24;
579PI(3)=0;
580PI(4)=1760;
581PI(5)=0;
582PI(6)=-25344;
583PI(7)=0;
584PI(8)=101376;
585PI(9)=0;
586PI(10)=-112640;
587PI(11)=0;
588PI(12)=24576;
589PI(13)=0;
590P = complex(P,PI);
591R = roots(P)
592S3=sqrt(3);
593E = [
594-1 + 2*%i
595%i
5963*%i
5971+2*%i
598(1/2)*(-S3+3*%i)
599(1/2)*(-S3+5*%i)
600-(1/2)*%i*(S3+(-4-%i))
601(1/2)*((1+4*%i)-%i*S3)
602(1/2)*%i*(S3+(4+%i))
603(1/2)*((1+4*%i)+%i*S3)
604(1/2)*(S3+3*%i)
605(1/2)*(S3+5*%i)
606];
607E = sortmyroots(E);
608R = sortmyroots(R);
609assert_close ( R, E, 1.e-8 );