Coverity #1321243, #1321246, #1321242, #1321236, #1321239, #1321238, #1097871, #13212...
[scilab.git] / scilab / modules / polynomials / sci_gateway / cpp / sci_pppdiv.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
26 extern "C"
27 {
28 #include "Scierror.h"
29 #include "localization.h"
30
31     extern double C2F(dasum)(int*, double*, int*);
32     extern int C2F(dpodiv)(double*, double*, int*, int*, int*);
33     extern int C2F(wpodiv)(double*, double*, double*, double*, int*, int*, int*);
34 }
35 /*--------------------------------------------------------------------------*/
36 types::Function::ReturnValue sci_pppdiv(types::typed_list &in, int _iRetCount, types::typed_list &out)
37 {
38     double* pdblInR[2]  = {NULL, NULL};// real part of denominator and numerator
39     double* pdblInI[2]  = {NULL, NULL};// rimaginary part
40     bool pbComplex[2]   = {false, false};
41     int piSize[2]       = {0, 0}; // rank+1 of denominator and numerator
42     int iErr            = 0;
43     int iOne            = 1;
44
45     std::wstring wstrName = L"";
46
47     if (in.size() != 2)
48     {
49         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "pppdiv", 2);
50         return types::Function::Error;
51     }
52
53     if (_iRetCount > 2)
54     {
55         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "pppdiv", 1, 2);
56         return types::Function::Error;
57     }
58
59     // get numerator and denominator
60     for (int i = 0; i < 2; i++)
61     {
62         if (in[i]->isDouble())
63         {
64             types::Double* pDblIn = in[i]->getAs<types::Double>();
65             if (pDblIn->isScalar() == false)
66             {
67                 Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "pppdiv", i + 1);
68                 return types::Function::Error;
69             }
70
71             piSize[i] = 1;
72             pdblInR[i] = pDblIn->get();
73             if (pDblIn->isComplex())
74             {
75                 pbComplex[i] = true;
76                 pdblInI[i] = pDblIn->getImg();
77             }
78         }
79         else if (in[i]->isPoly())
80         {
81             types::Polynom* pPolyIn = in[i]->getAs<types::Polynom>();
82             if (pPolyIn->isScalar() == false)
83             {
84                 Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "pppdiv", i + 1);
85                 return types::Function::Error;
86             }
87
88             if (wstrName != L"" && wstrName != pPolyIn->getVariableName())
89             {
90                 Scierror(999, _("%s: Wrong value for input argument #%d: A polynomial '%ls' expected.\n"), "pppdiv", i + 1, wstrName.c_str());
91                 return types::Function::Error;
92             }
93
94             wstrName = pPolyIn->getVariableName();
95             piSize[i] = pPolyIn->getMaxRank() + 1;
96             pdblInR[i] = pPolyIn->get(0)->get();
97             if (pPolyIn->isComplex())
98             {
99                 pbComplex[i] = true;
100                 pdblInI[i] = pPolyIn->get(0)->getImg();
101             }
102         }
103         else
104         {
105             std::wstring wstFuncName = L"%" + in[i]->getShortTypeStr() + L"_pppdiv";
106             return Overload::call(wstFuncName, in, _iRetCount, out);
107         }
108     }
109
110     // manage the case where the rank of the denominator
111     // is less than the rank of numerator. x / x² = 0 , rest = s.
112     if (piSize[0] < piSize[1])
113     {
114         if (_iRetCount == 2)
115         {
116             out.push_back(in[0]);
117         }
118         out.push_back(new types::Double(0));
119         return types::Function::OK;
120     }
121
122     // perform operation and set the result in pdblInR/I[0]
123     double* temp = pdblInR[0];
124     pdblInR[0] = new double[piSize[0]];
125     memcpy(pdblInR[0], temp, piSize[0] * sizeof(double));
126     if (pbComplex[0])
127     {
128         temp = pdblInI[0];
129         pdblInI[0] = new double[piSize[0]];
130         memcpy(pdblInI[0], temp, piSize[0] * sizeof(double));
131     }
132
133     int iDegreeNum = piSize[0] - 1;
134     int iDegreeDen = piSize[1] - 1;
135     if (pbComplex[0] == false && pbComplex[1] == false)
136     {
137         C2F(dpodiv)(pdblInR[0], pdblInR[1], &iDegreeNum, &iDegreeDen, &iErr);
138     }
139     else
140     {
141         if (pbComplex[0] == false)
142         {
143             pdblInI[0] = new double[piSize[0]];
144             memset(pdblInI[0], 0x00, piSize[0] * sizeof(double));
145         }
146         else if (pbComplex[1] == false)
147         {
148             pdblInI[1] = new double[piSize[1]];
149             memset(pdblInI[1], 0x00, piSize[1] * sizeof(double));
150         }
151
152         C2F(wpodiv)(pdblInR[0], pdblInI[0], pdblInR[1], pdblInI[1], &iDegreeNum, &iDegreeDen, &iErr);
153     }
154
155     // after execution : pdblInR[0](1 -> iDegreeDen) is the rest
156     //                   pdblInR[0](iDegreeDen -> end) is the coeff
157     int iSizeRest      = iDegreeDen;
158     int iSizeCoeff     = iDegreeNum - iDegreeDen + 1;
159     double* pdblRestR  = pdblInR[0];
160     double* pdblCoeffR = pdblInR[0] + iDegreeDen;
161     double* pdblRestI  = NULL;
162     double* pdblCoeffI = NULL;
163     bool bComplex      = pbComplex[0] || pbComplex[1];
164
165     if (bComplex)
166     {
167         pdblRestI  = pdblInI[0];
168         pdblCoeffI = pdblInI[0] + iDegreeDen;
169     }
170
171     // compute the real rank
172     if (bComplex) // complex case
173     {
174         for (int i = iSizeCoeff - 1; i >= 0; i--)
175         {
176             iSizeCoeff--;
177             if (std::fabs(pdblCoeffR[i]) + std::fabs(pdblCoeffI[i]))
178             {
179                 break;
180             }
181         }
182
183         for (int i = iSizeRest - 1; i >= 0; i--)
184         {
185             iSizeRest--;
186             if (std::fabs(pdblRestR[i]) + std::fabs(pdblRestI[i]))
187             {
188                 break;
189             }
190         }
191     }
192     else // real case
193     {
194         for (int i = iSizeCoeff - 1; i >= 0; i--)
195         {
196             iSizeCoeff--;
197             if (pdblCoeffR[i])
198             {
199                 break;
200             }
201         }
202
203         for (int i = iSizeRest - 1; i >= 0; i--)
204         {
205             iSizeRest--;
206             if (pdblRestR[i])
207             {
208                 break;
209             }
210         }
211     }
212
213     // set the rest of divide in output
214     if (_iRetCount == 2)
215     {
216         if (iSizeRest == 0) // return a types::Double
217         {
218             if (bComplex)
219             {
220                 out.push_back(new types::Double(pdblRestR[0], pdblRestI[0]));
221             }
222             else
223             {
224                 out.push_back(new types::Double(pdblRestR[0]));
225             }
226         }
227         else // return a types::Polynom
228         {
229             double* pdblReal = NULL;
230             types::Polynom* pPolyOut = new types::Polynom(wstrName, 1, 1);
231             types::SinglePoly* pSP = NULL;
232             int iSize = iSizeRest + 1;
233             if (bComplex && C2F(dasum)(&iSize, pdblRestI, &iOne) != 0)
234             {
235                 double* pdblImg = NULL;
236                 pSP = new types::SinglePoly(&pdblReal, &pdblImg, iSizeRest);
237                 memcpy(pdblImg, pdblRestI, iSize * sizeof(double));
238             }
239             else
240             {
241                 pSP = new types::SinglePoly(&pdblReal, iSizeRest);
242             }
243
244             memcpy(pdblReal, pdblRestR, iSize * sizeof(double));
245             pPolyOut->set(0, pSP);
246             delete pSP;
247             out.push_back(pPolyOut);
248         }
249     }
250
251     // set the result of divide in output
252     if (iSizeCoeff == 0) // return a types::Double
253     {
254         if (bComplex)
255         {
256             out.push_back(new types::Double(pdblCoeffR[0], pdblCoeffI[0]));
257         }
258         else
259         {
260             out.push_back(new types::Double(pdblCoeffR[0]));
261         }
262     }
263     else // return a types::Polynom
264     {
265         double* pdblReal = NULL;
266         types::Polynom* pPolyOut = new types::Polynom(wstrName, 1, 1);
267         types::SinglePoly* pSP = NULL;
268         int iSize = iSizeCoeff + 1;
269         if (bComplex && C2F(dasum)(&iSize, pdblCoeffI, &iOne) != 0)
270         {
271             double* pdblImg = NULL;
272             pSP = new types::SinglePoly(&pdblReal, &pdblImg, iSizeCoeff);
273             memcpy(pdblImg, pdblCoeffI, iSize * sizeof(double));
274         }
275         else
276         {
277             pSP = new types::SinglePoly(&pdblReal, iSizeCoeff);
278         }
279
280         memcpy(pdblReal, pdblCoeffR, iSize * sizeof(double));
281         pPolyOut->set(0, pSP);
282         delete pSP;
283         out.push_back(pPolyOut);
284     }
285
286     delete[] pdblInR[0];
287     if (pbComplex[0] || pbComplex[1])
288     {
289         delete[] pdblInI[0];
290     }
291
292     if (pbComplex[1] == false)
293     {
294         delete[] pdblInI[1];
295     }
296
297     return types::Function::OK;
298 }
299 /*--------------------------------------------------------------------------*/
300