invalid read in polynomials prod fixed and bool2s memory leak fixed
[scilab.git] / scilab / modules / elementary_functions / src / cpp / prod.cpp
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2012 - Scilab Enterprises - Cedric Delamarre
4  *
5  * This file must be used under the terms of the CeCILL.
6  * This source file is licensed as described in the file COPYING, which
7  * you should have received as part of this distribution.  The terms
8  * are also available at
9  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10  *
11  */
12 /*--------------------------------------------------------------------------*/
13
14 #include "prod.hxx"
15 extern "C"
16 {
17 #include "matrix_multiplication.h"
18 }
19
20 types::Double* prod(types::Double* pIn, int iOrientation)
21 {
22     types::Double* pOut = NULL;
23     double* pdblInReal  = pIn->getReal();
24     double* pdblInImg   = pIn->getImg();
25
26     if (iOrientation == 0) // all
27     {
28         int size = pIn->getSize();
29         double dblR = pdblInReal[0];
30         if (pIn->isComplex())
31         {
32             double dblI = pdblInImg[0];
33             double dblRTmp = 0;
34             double dblITmp = 0;
35             for (int i = 1; i < size; i++)
36             {
37                 dblRTmp = dblR * pdblInReal[i] - dblI * pdblInImg[i];
38                 dblITmp = dblR * pdblInImg[i] + dblI * pdblInReal[i];
39
40                 dblR = dblRTmp;
41                 dblI = dblITmp;
42             }
43
44             pOut = new types::Double(dblR, dblI);
45         }
46         else
47         {
48             for (int i = 1; i < size; i++)
49             {
50                 dblR *= pdblInReal[i];
51             }
52
53             pOut = new types::Double(dblR);
54         }
55     }
56     else // prod on one dimension
57     {
58         // create output dimensions
59         int iDims   = pIn->getDims();
60         int* piDims = new int[iDims];
61
62         for (int i = 0 ; i < iDims ; i++)
63         {
64             piDims[i] = pIn->getDimsArray()[i];
65         }
66
67         piDims[iOrientation - 1] = 1;
68
69         // create output variable
70         pOut = new types::Double(iDims, piDims, pIn->isComplex());
71         delete[] piDims;
72         piDims = NULL;
73
74         // init output
75         double* pdblOut = pOut->get();
76         double* pdblOutImg = pOut->getImg();
77         int size = pOut->getSize();
78
79         if (pOut->isComplex())
80         {
81             for (int i = 0; i < size; i++)
82             {
83                 pdblOut[i] = 1;
84                 pdblOutImg[i] = 0;
85             }
86         }
87         else
88         {
89             for (int i = 0; i < size; i++)
90             {
91                 pdblOut[i] = 1;
92             }
93         }
94
95         // perform operations
96         int* piIndex = new int[iDims];
97         size = pIn->getSize();
98         if (pIn->isComplex())
99         {
100             for (int i = 0; i < size; i++)
101             {
102                 //get array of dim
103                 pIn->getIndexes(i, piIndex);
104
105                 //convert indexes for result
106                 piIndex[iOrientation - 1] = 0;
107                 int iIndex = pOut->getIndex(piIndex);
108
109                 double dblRTmp = pdblOut[iIndex];
110                 double dblITmp = pdblOutImg[iIndex];
111
112                 iMultiComplexMatrixByComplexMatrix( &pdblInReal[i], &pdblInImg[i], 1, 1,
113                                                     &dblRTmp, &dblITmp, 1, 1,
114                                                     &pdblOut[iIndex], &pdblOutImg[iIndex]);
115             }
116         }
117         else
118         {
119             for (int i = 0; i < size; i++)
120             {
121                 //get array of dim
122                 pIn->getIndexes(i, piIndex);
123
124                 //convert indexes for result
125                 piIndex[iOrientation - 1] = 0;
126                 int iIndex = pOut->getIndex(piIndex);
127                 pdblOut[iIndex] *= pdblInReal[i];
128             }
129         }
130
131         delete[] piIndex;
132         piIndex = NULL;
133     }
134
135     return pOut;
136 }
137
138 // polynom
139 types::Polynom* prod(types::Polynom* pIn, int iOrientation)
140 {
141     types::Polynom* pOut = NULL;
142     if (iOrientation == 0) // sum of all element
143     {
144         // get rank of all input single poly
145         int* piRanks = new int[pIn->getSize()];
146         pIn->getRank(piRanks);
147
148         // compute the output rank
149         int iRankMax = 0;
150         for (int i = 0; i < pIn->getSize(); i++)
151         {
152             iRankMax += piRanks[i];
153         }
154
155         // create output matrix of poly
156         pOut = new types::Polynom(pIn->getVariableName(), 1, 1, &iRankMax);
157         pOut->setComplex(pIn->isComplex());
158         double* pdblRealOut = pOut->get(0)->get();
159         memcpy(pdblRealOut, pIn->get(0)->get(), (piRanks[0] + 1) * sizeof(double));
160
161         // do prod
162         int iSize = pIn->getSize();
163         double* pdblTempReal = new double[iRankMax + 1];
164         int iRank = piRanks[0];
165         if (pIn->isComplex())
166         {
167             double* pdblImgOut = pOut->get(0)->getImg();
168             // alloc temporary workspace
169             double* pdblTempImg  = new double[iRankMax + 1];
170             memcpy(pdblImgOut, pIn->get(0)->getImg(), (piRanks[0] + 1) * sizeof(double));
171
172             // perform operations
173             for (int i = 1; i < iSize; i++)
174             {
175                 double* pdblRealIn = pIn->get(i)->get();
176                 double* pdblImgIn = pIn->get(i)->getImg();
177                 memcpy(pdblTempReal, pdblRealOut, (iRank + 1) * sizeof(double));
178                 memcpy(pdblTempImg,  pdblImgOut,  (iRank + 1) * sizeof(double));
179                 iMultiComplexPolyByComplexPoly(pdblTempReal, pdblTempImg, iRank + 1,
180                                                pdblRealIn, pdblImgIn, piRanks[i] + 1,
181                                                pdblRealOut, pdblImgOut, iRankMax + 1);
182                 iRank += piRanks[i];
183             }
184
185             delete[] pdblTempImg;
186         }
187         else
188         {
189             // perform operations
190             for (int i = 1; i < iSize; i++)
191             {
192                 double* pdblRealIn = pIn->get(i)->get();
193                 memcpy(pdblTempReal, pdblRealOut, (iRank + 1) * sizeof(double));
194                 iMultiScilabPolynomByScilabPolynom(pdblTempReal, iRank + 1,
195                                                    pdblRealIn, piRanks[i] + 1,
196                                                    pdblRealOut, iRankMax + 1);
197                 iRank += piRanks[i];
198             }
199         }
200
201         delete[] pdblTempReal;
202         delete[] piRanks;
203     }
204     else // sum following a dimension
205     {
206         // create output dimensions
207         int iDims = pIn->getDims();
208         int* piDims = new int[iDims];
209
210         for (int i = 0 ; i < iDims ; i++)
211         {
212             piDims[i] = pIn->getDimsArray()[i];
213         }
214
215         piDims[iOrientation - 1] = 1;
216
217         // get input ranks in int*
218         int* piRanks = new int[pIn->getSize()];
219         pIn->getRank(piRanks);
220
221         // get input ranks in types::Double
222         types::Double* pDblRanks = new types::Double(pIn->getDims(), pIn->getDimsArray());
223         for (int i = 0; i < pDblRanks->getSize(); i++)
224         {
225             pDblRanks->set(i, static_cast<double>(piRanks[i]));
226         }
227
228         // create output max ranks
229         types::Double* pDblRanksOut = new types::Double(iDims, piDims);
230         double* pdblOut = pDblRanksOut->get();
231         for (int i = 0; i < pDblRanksOut->getSize(); i++)
232         {
233             pdblOut[i] = 0;
234         }
235
236         // compute the maximum ranks for the dim n
237         int*    piIndex = new int[iDims];
238         double* pdblIn  = pDblRanks->get();
239         for (int i = 0 ; i < pDblRanks->getSize() ; i++)
240         {
241             //get array of dim
242             pDblRanks->getIndexes(i, piIndex);
243
244             //convert indexes for result
245             piIndex[iOrientation - 1] = 0;
246             int iIndex = pDblRanksOut->getIndex(piIndex);
247             pdblOut[iIndex] += pdblIn[i];
248         }
249         pDblRanks->killMe();
250
251         // move output ranks from types::Double to int*
252         int* piRankMax = new int[pDblRanksOut->getSize()];
253         int* piRankTmp = new int[pDblRanksOut->getSize()];
254         int iMaxOutputRank = 0;
255         for (int i = 0; i < pDblRanksOut->getSize(); i++)
256         {
257             piRankMax[i] = static_cast<int>(pdblOut[i]);
258             iMaxOutputRank = std::max(iMaxOutputRank, piRankMax[i]);
259         }
260
261         pDblRanksOut->killMe();
262         // create the outpout polynom
263         pOut = new types::Polynom(pIn->getVariableName(), iDims, piDims, piRankMax);
264         pOut->setComplex(pIn->isComplex());
265
266         // init output with first element of lead dimension
267         if (pIn->isComplex())
268         {
269
270             for (int i = 0; i < pOut->getSize(); i++)
271             {
272                 pOut->getIndexes(i, piIndex);
273                 int iIndex = pIn->getIndex(piIndex);
274                 memcpy(pOut->get(i)->get(), pIn->get(iIndex)->get(), (piRanks[iIndex] + 1) * sizeof(double));
275                 memcpy(pOut->get(i)->getImg(), pIn->get(iIndex)->getImg(), (piRanks[iIndex] + 1) * sizeof(double));
276                 piRankTmp[i] = piRanks[iIndex];
277             }
278         }
279         else
280         {
281             for (int i = 0; i < pOut->getSize(); i++)
282             {
283                 pOut->getIndexes(i, piIndex);
284                 int iIndex = pIn->getIndex(piIndex);
285                 memcpy(pOut->get(i)->get(), pIn->get(iIndex)->get(), (piRanks[iIndex] + 1) * sizeof(double));
286                 piRankTmp[i] = piRanks[iIndex];
287             }
288         }
289
290         // alloc temporary workspace
291         double* pdblTempReal = new double[iMaxOutputRank + 1];
292         if (pIn->isComplex())
293         {
294             // alloc temporary workspace
295             double* pdblTempImg  = new double[iMaxOutputRank + 1];
296
297             // perform operations
298             for (int i = 0 ; i < pIn->getSize() ; i++)
299             {
300                 //get array of dim
301                 pIn->getIndexes(i, piIndex);
302
303                 //convert indexes for result
304                 if (piIndex[iOrientation - 1] == 0)
305                 {
306                     // first element of lead dimension is already setted.
307                     continue;
308                 }
309
310                 piIndex[iOrientation - 1] = 0;
311                 int iIndex = pOut->getIndex(piIndex);
312
313                 // make sum on each ranks
314                 double* pdblRealIn = pIn->get(i)->get();
315                 double* pdblRealOut = pOut->get(iIndex)->get();
316                 double* pdblImgIn = pIn->get(i)->getImg();
317                 double* pdblImgOut = pOut->get(iIndex)->getImg();
318                 memcpy(pdblTempReal, pdblRealOut, (piRankTmp[iIndex] + 1) * sizeof(double));
319                 memcpy(pdblTempImg,  pdblImgOut,  (piRankTmp[iIndex] + 1) * sizeof(double));
320                 iMultiComplexPolyByComplexPoly(pdblTempReal, pdblTempImg, piRankTmp[iIndex] + 1,
321                                                pdblRealIn, pdblImgIn, piRanks[i] + 1,
322                                                pdblRealOut, pdblImgOut, piRankMax[iIndex] + 1);
323                 piRankTmp[iIndex] += piRanks[i];
324             }
325
326             delete[] pdblTempImg;
327         }
328         else
329         {
330             // perform operations
331             for (int i = 0 ; i < pIn->getSize() ; i++)
332             {
333                 //get array of dim
334                 pIn->getIndexes(i, piIndex);
335
336                 //convert indexes for result
337                 if (piIndex[iOrientation - 1] == 0)
338                 {
339                     // first element of lead dimension is already setted.
340                     continue;
341                 }
342
343                 piIndex[iOrientation - 1] = 0;
344                 int iIndex = pOut->getIndex(piIndex);
345
346                 // make sum on each ranks
347                 double* pdblRealIn = pIn->get(i)->get();
348                 double* pdblRealOut = pOut->get(iIndex)->get();
349                 memcpy(pdblTempReal, pdblRealOut, (piRankTmp[iIndex] + 1) * sizeof(double));
350                 iMultiScilabPolynomByScilabPolynom(pdblTempReal, piRankTmp[iIndex] + 1,
351                                                    pdblRealIn, piRanks[i] + 1,
352                                                    pdblRealOut, piRankMax[iIndex] + 1);
353                 piRankTmp[iIndex] += piRanks[i];
354             }
355         }
356
357         delete[] pdblTempReal;
358         delete[] piRankMax;
359         delete[] piRanks;
360         delete[] piIndex;
361         delete[] piDims;
362         delete[] piRankTmp;
363     }
364
365     pOut->updateRank();
366     return pOut;
367 }
368