4056d728bb78fb57aafdaca754d454f3dd205da8
[scilab.git] / scilab / modules / polynomials / sci_gateway / cpp / sci_sfact.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 #include "polynomials_gw.hxx"
17 #include "function.hxx"
18 #include "double.hxx"
19 #include "polynom.hxx"
20 #include "overload.hxx"
21
22 extern "C"
23 {
24 #include "Scierror.h"
25 #include "localization.h"
26 #include "sciprint.h"
27 #include "elem_common.h"
28     extern int C2F(sfact1)(double*, int*, double*, int*, int*);
29     extern int C2F(sfact2)(double*, int*, int*, double*, int*, int*);
30 }
31 /*--------------------------------------------------------------------------*/
32 types::Function::ReturnValue sci_sfact(types::typed_list &in, int _iRetCount, types::typed_list &out)
33 {
34     types::Polynom* pPolyIn  = NULL;
35     types::Polynom* pPolyOut = NULL;
36
37     int iMaxIt = 100;
38     int iErr   = 0;
39     int iOne   = 1;
40
41     if (in.size() != 1)
42     {
43         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "sfact", 1);
44         return types::Function::Error;
45     }
46
47     if (_iRetCount > 1)
48     {
49         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "sfact", 1);
50         return types::Function::Error;
51     }
52
53     if (in[0]->isPoly() == false)
54     {
55         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_sfact";
56         return Overload::call(wstFuncName, in, _iRetCount, out);
57     }
58
59     pPolyIn = in[0]->getAs<types::Polynom>();
60
61     if (pPolyIn->isComplex())
62     {
63         Scierror(999, _("%s: Wrong value for input argument #%d: A real matrix expected.\n"), "sfact", 1);
64         return types::Function::Error;
65     }
66
67     if (pPolyIn->isScalar())
68     {
69         double* pdblCoef = pPolyIn->get(0)->get();
70
71         // check symmetry
72         int iDegree = pPolyIn->get(0)->getRank();
73         int iDegD2  = (int)(iDegree / 2);
74         int n       = 1 + iDegD2;
75
76         if (2 * iDegD2 != iDegree)
77         {
78             Scierror(999, _("%s: Wrong value for input argument #%d: A symmetric polynomial expected.\n"), "sfact", 1);
79             return types::Function::Error;
80         }
81
82         for (int i = 0; i < n; i++)
83         {
84             if (pdblCoef[i] != pdblCoef[iDegree - i])
85             {
86                 Scierror(999, _("%s: Wrong value for input argument #%d: A symmetric polynomial expected.\n"), "sfact", 1);
87                 return types::Function::Error;
88             }
89         }
90
91         // create result
92         pPolyOut = new types::Polynom(pPolyIn->getVariableName(), 1, 1);
93         double* pdblCoefOut = NULL;
94         types::SinglePoly* pSP = new types::SinglePoly(&pdblCoefOut, iDegD2);
95         C2F(dcopy)(&n, pdblCoef, &iOne, pdblCoefOut, &iOne);
96
97         // perform operation
98         double* pdblWork = new double[7 * n];
99         C2F(sfact1)(pdblCoefOut, &iDegD2, pdblWork, &iMaxIt, &iErr);
100         delete pdblWork;
101         if (iErr == 2)
102         {
103             delete pSP;
104             delete pPolyOut;
105             Scierror(999, _("%s: Wrong value for input argument #%d: Non negative or null value expected at degree %d.\n"), "sfact", 1, n);
106             return types::Function::Error;
107         }
108         else if (iErr == 1)
109         {
110             delete pSP;
111             delete pPolyOut;
112             Scierror(999, _("%s: Wrong value for input argument #%d: Convergence problem.\n"), "sfact", 1);
113             return types::Function::Error;
114         }
115         else if (iErr < 0)
116         {
117             sciprint("%s: warning: Convergence at 10^%d near.\n", "sfact", iErr);
118         }
119
120         // return result
121         pPolyOut->set(0, pSP);
122         delete pSP;
123     }
124     else
125     {
126         if (pPolyIn->getRows() != pPolyIn->getCols())
127         {
128             Scierror(999, _("%s: Wrong size for input argument #%d: Square matrix expected.\n"), "sfact", 1);
129             return types::Function::Error;
130         }
131
132         int iSize   = pPolyIn->getSize();
133         int iRows   = pPolyIn->getRows();
134         int iDegree = pPolyIn->getMaxRank();
135         int iDegD2  = (int)(iDegree / 2);
136         int n       = 1 + iDegD2;
137
138         double* pdblOut  = new double[iSize * n];
139         double* pdblWork = new double[(n + 1) * iRows * ((n + 1)*iRows) + 1];
140
141         memset(pdblOut, 0x00, iSize * n * sizeof(double));
142
143         for (int i = 0; i < iSize; i++)
144         {
145             double* pdblIn = pPolyIn->get(i)->get();
146             int iSizeToCpy = 2 + pPolyIn->get(i)->getSize() - 1 - n;
147             if (iSizeToCpy > 0)
148             {
149                 C2F(dcopy)(&iSizeToCpy, pdblIn + n - 1, &iOne, pdblOut + i, &iSize);
150             }
151         }
152
153         int nm1 = n - 1;
154         iMaxIt += n;
155         C2F(sfact2)(pdblOut, &iRows, &nm1, pdblWork, &iMaxIt, &iErr);
156         delete pdblWork;
157
158         if (iErr < 0)
159         {
160             delete[] pdblOut;
161             Scierror(999, _("%s: Wrong value for input argument #%d: Convergence problem.\n"), "sfact", 1);
162             return types::Function::Error;
163         }
164         else if (iErr > 0)
165         {
166             delete[] pdblOut;
167             Scierror(999, _("%s: Wrong value for input argument #%d: singular or asymmetric problem.\n"), "sfact", 1);
168             return types::Function::Error;
169         }
170
171         pPolyOut = new types::Polynom(pPolyIn->getVariableName(), pPolyIn->getDims(), pPolyIn->getDimsArray());
172         for (int i = 0; i < iSize; i++)
173         {
174             double* pdblCoefOut = NULL;
175             types::SinglePoly* pSP = new types::SinglePoly(&pdblCoefOut, nm1);
176             C2F(dcopy)(&n, pdblOut + i, &iSize, pdblCoefOut, &iOne);
177             pPolyOut->set(i, pSP);
178             delete pSP;
179         }
180
181         delete[] pdblOut;
182     }
183
184     out.push_back(pPolyOut);
185     return types::Function::OK;
186 }
187 /*--------------------------------------------------------------------------*/
188