2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2012 - Scilab Enterprises - Cedric DELAMARRE
5 * Copyright (C) 2012 - 2016 - Scilab Enterprises
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.
15 /*--------------------------------------------------------------------------*/
19 #include "polynomials_gw.hxx"
20 #include "function.hxx"
23 #include "polynom.hxx"
24 #include "overload.hxx"
25 #include "context.hxx"
26 #include "execvisitor.hxx"
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*);
36 /*--------------------------------------------------------------------------*/
37 types::Function::ReturnValue sci_roots(types::typed_list &in, int _iRetCount, types::typed_list &out)
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;
44 double* pdblInReal = NULL;
45 double* pdblInImg = NULL;
46 double* pdblTempReal = NULL;
47 double* pdblTempImg = NULL;
52 bool bComplex = false;
53 types::Function::ReturnValue ret = types::Function::OK;
55 if (in.size() < 1 || in.size() > 2)
57 Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "roots", 1, 2);
58 return types::Function::Error;
63 Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "roots", 1);
64 return types::Function::Error;
70 if (in[1]->isString() == false)
72 Scierror(999, _("%s: Wrong type for input argument #%d : string expected.\n"), "roots", 2);
73 return types::Function::Error;
76 types::String* pStrAlgo = in[1]->getAs<types::String>();
77 if (pStrAlgo->isScalar() == false)
79 Scierror(999, _("%s: Wrong size for input argument #%d : A scalar expected.\n"), "roots", 2);
80 return types::Function::Error;
83 wstrAlgo = pStrAlgo->get(0);
84 if (wstrAlgo != L"e" && wstrAlgo != L"f")
86 Scierror(999, _("%s: Wrong value for input argument #%d : ""%s"" or ""%s"" expected.\n"), "roots", 2, "e", "f");
87 return types::Function::Error;
91 if (in[0]->isDouble())
93 // for Matlab compatibility root of the vector of coefficients
94 pDblIn = in[0]->getAs<types::Double>();
95 if (pDblIn->isEmpty())
97 out.push_back(types::Double::Empty());
98 return types::Function::OK;
101 iSize = pDblIn->getSize();
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())
110 pdblInImg = new double[iSize];
111 C2F(dcopy)(&iSize, pDblIn->getImg(), &iOne, pdblInImg, &imOne);
114 else if (in[0]->isPoly())
116 pPolyIn = in[0]->getAs<types::Polynom>();
118 if (pPolyIn->isScalar() == false)
120 Scierror(999, _("%s: Wrong size for input argument #%d : A scalar expected.\n"), "roots", 1);
121 return types::Function::Error;
124 iSize = pPolyIn->getMaxRank() + 1;
125 pdblInReal = pPolyIn->get(0)->get();
127 if (pPolyIn->isComplex())
130 pdblInImg = pPolyIn->get(0)->getImg();
135 std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_roots";
136 return Overload::call(wstFuncName, in, _iRetCount, out);
139 // If "fast" algo was chosen and polynomial is complex,
140 // then produce an error.
141 if (wstrAlgo == L"f" && bComplex)
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;
153 out.push_back(types::Double::Empty());
154 return types::Function::OK;
157 t = std::fabs(pdblInReal[iSize]);
160 t += std::fabs(pdblInImg[iSize]);
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)
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;
172 if (wstrAlgo == L"f")
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
181 pDblOut = new types::Double(iSize, 1, true);
183 int iSizeP1 = iSize + 1;
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;
194 Scierror(999, _("%s: Convergence problem...\n"), "roots");
198 Scierror(999, _("%s: Leading coefficient is zero.\n"), "roots");
202 Scierror(999, _("%s: Too high degree (max 100).\n"), "roots");
205 return types::Function::Error;
208 out.push_back(pDblOut);
210 else // wstrAlgo == L"e"
212 // Companion matrix method
213 int iSizeM1 = iSize - 1;
214 int iSizeP1 = iSize + 1;
216 double* pdblTempReal = new double[iSize];
217 double* pdblTempImg = NULL;
218 double sr = pdblInReal[iSize];
220 C2F(dcopy)(&iSize, pdblInReal, &iOne, pdblTempReal, &imOne);
224 pdblTempImg = new double[iSize];
225 C2F(dcopy)(&iSize, pdblInImg, &iOne, pdblTempImg, &imOne);
227 double si = pdblInImg[iSize];
228 double t = sr * sr + si * si;
232 C2F(wscal)(&iSize, &sr, &si, pdblTempReal, pdblTempImg, &iOne);
236 double dbl = -1 / sr;
237 C2F(dscal)(&iSize, &dbl, pdblTempReal, &iOne);
240 pDblOut = new types::Double(iSize, iSize, bComplex);
241 double* pdblOutReal = pDblOut->get();
242 double* pdblOutImg = NULL;
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;
251 pdblOutImg = pDblOut->getImg();
252 memset(pdblOutImg, 0x00, pDblOut->getSize() * sizeof(double));
253 C2F(dcopy)(&iSize, pdblTempImg, &iOne, pdblOutImg, &iOne);
254 delete[] pdblTempImg;
258 types::InternalType* pSpec = symbol::Context::getInstance()->get(symbol::Symbol(L"spec"));
259 if (pSpec && pSpec->isFunction())
261 types::Function *funcSpec = pSpec->getAs<types::Function>();
262 types::typed_list tlInput;
263 types::optional_list tlOpt;
264 tlInput.push_back(pDblOut);
266 ret = funcSpec->call(tlInput, tlOpt, 1, out);
271 Scierror(999, _("%s: unable to find spec function\n"), "roots");
272 return types::Function::Error;
287 /*--------------------------------------------------------------------------*/