a3181f60fb38739a15ada60c1c0edb0affbcc077
[scilab.git] / scilab / modules / polynomials / sci_gateway / cpp / sci_roots.cpp
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2012 - Scilab Enterprises - Cedric DELAMARRE
4  *
5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
6  *
7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
9  * This file was originally licensed under the terms of the CeCILL v2.1,
10  * and continues to be available under such terms.
11  * For more information, see the COPYING file which you should have received
12  * along with this program.
13  *
14  */
15 /*--------------------------------------------------------------------------*/
16
17 #include <cmath>
18
19 #include "polynomials_gw.hxx"
20 #include "function.hxx"
21 #include "double.hxx"
22 #include "string.hxx"
23 #include "polynom.hxx"
24 #include "overload.hxx"
25 #include "context.hxx"
26 #include "execvisitor.hxx"
27
28 extern "C"
29 {
30 #include "Scierror.h"
31 #include "localization.h"
32 #include "elem_common.h"
33     extern int C2F(rpoly)(double*, int*, double*, double*, int*);
34     extern int C2F(wscal)(int*, double*, double*, double*, double*, int*);
35 }
36 /*--------------------------------------------------------------------------*/
37 types::Function::ReturnValue sci_roots(types::typed_list &in, int _iRetCount, types::typed_list &out)
38 {
39     std::wstring wstrAlgo   = L"e"; // e = eigen (default), f = fast
40     types::Double* pDblIn   = NULL;
41     types::Double* pDblOut  = NULL;
42     types::Polynom* pPolyIn = NULL;
43
44     double* pdblInReal   = NULL;
45     double* pdblInImg    = NULL;
46     double* pdblTempReal = NULL;
47     double* pdblTempImg  = NULL;
48
49     int iOne = 1;
50     int imOne = -1;
51     int iSize = 0;
52     bool bComplex = false;
53     types::Function::ReturnValue ret = types::Function::OK;
54
55     if (in.size() < 1 || in.size() > 2)
56     {
57         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "roots", 1, 2);
58         return types::Function::Error;
59     }
60
61     if (_iRetCount > 1)
62     {
63         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "roots", 1);
64         return types::Function::Error;
65     }
66
67     // get algo type
68     if (in.size() == 2)
69     {
70         if (in[1]->isString() == false)
71         {
72             Scierror(999, _("%s: Wrong type for input argument #%d : string expected.\n"), "roots", 2);
73             return types::Function::Error;
74         }
75
76         types::String* pStrAlgo = in[1]->getAs<types::String>();
77         if (pStrAlgo->isScalar() == false)
78         {
79             Scierror(999, _("%s: Wrong size for input argument #%d : A scalar expected.\n"), "roots", 2);
80             return types::Function::Error;
81         }
82
83         wstrAlgo = pStrAlgo->get(0);
84         if (wstrAlgo != L"e" && wstrAlgo != L"f")
85         {
86             Scierror(999, _("%s: Wrong value for input argument #%d : ""%s"" or ""%s"" expected.\n"), "roots", 2, "e", "f");
87             return types::Function::Error;
88         }
89     }
90
91     if (in[0]->isDouble())
92     {
93         // for Matlab compatibility root of the vector of coefficients
94         pDblIn = in[0]->getAs<types::Double>();
95         if (pDblIn->isEmpty())
96         {
97             out.push_back(types::Double::Empty());
98             return types::Function::OK;
99         }
100
101         iSize = pDblIn->getSize();
102
103         // old fortran function dtild
104         // switch elements of a vector. [1 2 3] => [3 2 1]
105         pdblInReal = new double[iSize];
106         C2F(dcopy)(&iSize, pDblIn->get(), &iOne, pdblInReal, &imOne);
107         if (pDblIn->isComplex())
108         {
109             bComplex = true;
110             pdblInImg = new double[iSize];
111             C2F(dcopy)(&iSize, pDblIn->getImg(), &iOne, pdblInImg, &imOne);
112         }
113     }
114     else if (in[0]->isPoly())
115     {
116         pPolyIn = in[0]->getAs<types::Polynom>();
117
118         if (pPolyIn->isScalar() == false)
119         {
120             Scierror(999, _("%s: Wrong size for input argument #%d : A scalar expected.\n"), "roots", 1);
121             return types::Function::Error;
122         }
123
124         iSize      = pPolyIn->getMaxRank() + 1;
125         pdblInReal = pPolyIn->get(0)->get();
126
127         if (pPolyIn->isComplex())
128         {
129             bComplex = true;
130             pdblInImg = pPolyIn->get(0)->getImg();
131         }
132     }
133     else
134     {
135         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_roots";
136         return Overload::call(wstFuncName, in, _iRetCount, out);
137     }
138
139     // If "fast" algo was chosen and polynomial is complex,
140     // then produce an error.
141     if (wstrAlgo == L"f" && bComplex)
142     {
143         Scierror(999, _("%s: Wrong value for input argument #%d : If algo is ""%s"", a real is expected.\n"), "roots", 2, "f");
144         return types::Function::Error;
145     }
146
147     double t = 0;
148     while (t == 0)
149     {
150         iSize--;
151         if (iSize < 0)
152         {
153             out.push_back(types::Double::Empty());
154             return types::Function::OK;
155         }
156
157         t = std::fabs(pdblInReal[iSize]);
158         if (bComplex)
159         {
160             t += std::fabs(pdblInImg[iSize]);
161         }
162     }
163
164     // If "fast" algo was chosen and polynomial has degree greater than 100,
165     // then produce an error.
166     if (wstrAlgo == L"f" && iSize > 100)
167     {
168         Scierror(999, _("%s: Wrong value for input argument #%d : If algo is ""%s"", a degree less than %d expected.\n"), "roots", 2, "f", 100);
169         return types::Function::Error;
170     }
171
172     if (wstrAlgo == L"f")
173     {
174         // real polynomial: rpoly algorithm
175         // this alg is much more speedy, but it may happens that it gives
176         // erroneous results without messages : example
177         // roots(%s^31-8*%s^30+9*%s^29+0.995) should have two real roots near
178         // 1.355 and 6.65 and the other ones inside a circle centered in 0
179         // with radius 1
180
181         pDblOut = new types::Double(iSize, 1, true);
182         int iFail = 0;
183         int iSizeP1 = iSize + 1;
184
185         double* pdblTempReal = new double[iSize + 1];
186         C2F(dcopy)(&iSizeP1, pdblInReal, &iOne, pdblTempReal, &imOne);
187         C2F(rpoly)(pdblTempReal, &iSize, pDblOut->get(), pDblOut->getImg(), &iFail);
188         delete pdblTempReal;
189
190         if (iFail)
191         {
192             if (iFail == 1)
193             {
194                 Scierror(999, _("%s: Convergence problem...\n"), "roots");
195             }
196             else if (iFail == 2)
197             {
198                 Scierror(999, _("%s: Leading coefficient is zero.\n"), "roots");
199             }
200             else if (iFail == 3)
201             {
202                 Scierror(999, _("%s: Too high degree (max 100).\n"), "roots");
203             }
204
205             return types::Function::Error;
206         }
207
208         out.push_back(pDblOut);
209     }
210     else // wstrAlgo == L"e"
211     {
212         // Companion matrix method
213         int iSizeM1 = iSize - 1;
214         int iSizeP1 = iSize + 1;
215         double dblOne = 1;
216         double* pdblTempReal = new double[iSize];
217         double* pdblTempImg = NULL;
218         double sr = pdblInReal[iSize];
219
220         C2F(dcopy)(&iSize, pdblInReal, &iOne, pdblTempReal, &imOne);
221
222         if (bComplex)
223         {
224             pdblTempImg = new double[iSize];
225             C2F(dcopy)(&iSize, pdblInImg, &iOne, pdblTempImg, &imOne);
226
227             double si = pdblInImg[iSize];
228             double t = sr * sr + si * si;
229             sr = -sr / t;
230             si = si / t;
231
232             C2F(wscal)(&iSize, &sr, &si, pdblTempReal, pdblTempImg, &iOne);
233         }
234         else
235         {
236             double dbl = -1 / sr;
237             C2F(dscal)(&iSize, &dbl, pdblTempReal, &iOne);
238         }
239
240         pDblOut = new types::Double(iSize, iSize, bComplex);
241         double* pdblOutReal = pDblOut->get();
242         double* pdblOutImg = NULL;
243
244         memset(pdblOutReal, 0x00, pDblOut->getSize() * sizeof(double));
245         C2F(dset)(&iSizeM1, &dblOne, &pdblOutReal[iSize], &iSizeP1);
246         C2F(dcopy)(&iSize, pdblTempReal, &iOne, pdblOutReal, &iOne);
247         delete[] pdblTempReal;
248
249         if (bComplex)
250         {
251             pdblOutImg = pDblOut->getImg();
252             memset(pdblOutImg, 0x00, pDblOut->getSize() * sizeof(double));
253             C2F(dcopy)(&iSize, pdblTempImg, &iOne, pdblOutImg, &iOne);
254             delete[] pdblTempImg;
255         }
256
257         //call spec
258         types::InternalType* pSpec = symbol::Context::getInstance()->get(symbol::Symbol(L"spec"));
259         if (pSpec && pSpec->isFunction())
260         {
261             types::Function *funcSpec = pSpec->getAs<types::Function>();
262             types::typed_list tlInput;
263             types::optional_list tlOpt;
264             tlInput.push_back(pDblOut);
265
266             ret = funcSpec->call(tlInput, tlOpt, 1, out);
267             pDblOut->killMe();
268         }
269         else
270         {
271             Scierror(999, _("%s: unable to find spec function\n"), "roots");
272             return types::Function::Error;
273         }
274     }
275
276     if (pDblIn)
277     {
278         delete pdblInReal;
279         if (bComplex)
280         {
281             delete pdblInImg;
282         }
283     }
284
285     return ret;
286 }
287 /*--------------------------------------------------------------------------*/
288