9f88f184f36d39177c5416f27b57e300553d5c24
[scilab.git] / scilab / modules / elementary_functions / src / cpp / cumprod.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 "cumprod.hxx"
15 extern "C"
16 {
17 #include "matrix_multiplication.h"
18 }
19
20 int cumprod(types::Double* pIn, int iOrientation, types::Double* pOut)
21 {
22     double* pdblInReal  = pIn->getReal();
23     double* pdblOutReal = pOut->getReal();
24     double* pdblInImg   = pIn->getImg();
25     double* pdblOutImg  = pOut->getImg();
26
27     if (iOrientation == 0) // all
28     {
29         pdblOutReal[0] = pdblInReal[0];
30         if (pIn->isComplex())
31         {
32             pdblOutImg[0] = pdblInImg[0];
33             for (int i = 1; i < pIn->getSize(); i++)
34             {
35                 pdblOutReal[i] = pdblInReal[i] * pdblOutReal[i - 1];
36                 pdblOutReal[i] -= pdblInImg[i] * pdblOutImg[i - 1];
37                 pdblOutImg[i] = (pdblInReal[i] * pdblOutImg[i - 1]) + (pdblInImg[i] * pdblOutReal[i - 1]);
38             }
39         }
40         else
41         {
42             for (int i = 1; i < pIn->getSize(); i++)
43             {
44                 pdblOutReal[i] = pdblInReal[i] * pdblOutReal[i - 1];
45             }
46         }
47     }
48     else // cumprod on one dimension
49     {
50         int iSizeOfDimN = pIn->getDimsArray()[iOrientation - 1];
51         int iIncrement = 1;
52
53         for (int i = 0; i < iOrientation - 1; i++)
54         {
55             iIncrement *= pIn->getDimsArray()[i];
56         }
57
58         if (pIn->isComplex())
59         {
60             for (int j = 0; j < pIn->getSize(); j += (iIncrement * iSizeOfDimN))
61             {
62                 for (int i = j; i < iIncrement + j; i++) // set the first values in out
63                 {
64                     pdblOutReal[i] = pdblInReal[i];
65                     pdblOutImg[i] = pdblInImg[i];
66                 }
67
68                 for (int k = 1; k < iSizeOfDimN; k++) // make the cumprod for the next values
69                 {
70                     for (int i = (iIncrement * k) + j; i < (iIncrement * (k + 1)) + j; i++)
71                     {
72                         pdblOutReal[i] = pdblInReal[i] * pdblOutReal[i - iIncrement];
73                         pdblOutReal[i] -= pdblInImg[i] * pdblOutImg[i - iIncrement];
74                         pdblOutImg[i] = (pdblInReal[i] * pdblOutImg[i - iIncrement]) + (pdblInImg[i] * pdblOutReal[i - iIncrement]);
75                     }
76                 }
77             }
78         }
79         else
80         {
81             for (int j = 0; j < pIn->getSize(); j += (iIncrement * iSizeOfDimN))
82             {
83                 for (int i = j; i < iIncrement + j; i++) // set the first values in out
84                 {
85                     pdblOutReal[i] = pdblInReal[i];
86                 }
87
88                 for (int k = 1; k < iSizeOfDimN; k++) // make the cumprod for the next values
89                 {
90                     for (int i = (iIncrement * k) + j; i < (iIncrement * (k + 1)) + j; i++)
91                     {
92                         pdblOutReal[i] = pdblInReal[i] * pdblOutReal[i - iIncrement];
93                     }
94                 }
95             }
96         }
97     }
98
99     return 0;
100 }
101
102 // polynom
103 int cumprod(types::Polynom* pIn, int iOrientation, types::Polynom* pOut)
104 {
105     int iErr        = 0;
106     int iRank       = 0;
107     int iOutRank    = 0;
108     int iLastRank   = 0;
109
110     double* pdRData         = NULL;
111     double* pdIData         = NULL;
112     double* pdblReal        = NULL;
113     double* pdblImg         = NULL;
114     double* pdblLastReal    = NULL;
115     double* pdblLastImg     = NULL;
116
117     bool bComplex   = pIn->isComplex();
118     types::SinglePoly* pSP = NULL;
119
120     if (iOrientation == 0) // all
121     {
122         // set first element
123         iRank = pIn->get(0)->getRank();
124         pdblReal = pIn->get(0)->getCoefReal();
125         if (bComplex)
126         {
127             pSP = new types::SinglePoly(&pdRData, &pdIData, iRank);
128             pdblImg = pIn->get(0)->getCoefImg();
129             for (int j = 0; j < iRank; j++)
130             {
131                 pdRData[j] = pdblReal[j];
132                 pdIData[j] = pdblImg[j];
133             }
134         }
135         else
136         {
137             pSP = new types::SinglePoly(&pdRData, iRank);
138             for (int j = 0; j < iRank; j++)
139             {
140                 pdRData[j] = pdblReal[j];
141             }
142         }
143
144         pOut->set(0, pSP);
145         iLastRank = iRank;
146         pdblLastReal = pdblReal;
147         pdblLastImg = pdblImg;
148
149         // compute next elements
150         if (bComplex)
151         {
152             for (int i = 1; i < pIn->getSize(); i++)
153             {
154                 pdblReal = pIn->get(i)->getCoefReal();
155                 pdblImg = pIn->get(i)->getCoefImg();
156                 iRank = pIn->get(i)->getRank();
157
158                 iOutRank = iLastRank + iRank - 1;
159
160                 pSP = new types::SinglePoly(&pdRData, &pdIData, iOutRank);
161                 pSP->getCoef()->setZeros();
162
163                 iMultiComplexPolyByComplexPoly( pdblReal, pdblImg, iRank,
164                                                 pdblLastReal, pdblLastImg, iLastRank,
165                                                 pdRData, pdIData, iOutRank);
166
167                 pOut->set(i, pSP);
168                 pdblLastReal = pdRData;
169                 pdblLastImg = pdIData;
170                 iLastRank = iOutRank;
171                 delete pSP;
172             }
173         }
174         else
175         {
176             for (int i = 1; i < pIn->getSize(); i++)
177             {
178                 pdblReal = pIn->get(i)->getCoefReal();
179                 iRank = pIn->get(i)->getRank();
180
181                 iOutRank = iLastRank + iRank - 1;
182
183                 pSP = new types::SinglePoly(&pdRData, iOutRank);
184                 pSP->getCoef()->setZeros();
185                 iMultiScilabPolynomByScilabPolynom(pdblReal, iRank, pdblLastReal, iLastRank, pdRData, iOutRank);
186
187                 pOut->set(i, pSP);
188                 pdblLastReal = pdRData;
189                 iLastRank = iOutRank;
190                 delete pSP;
191             }
192         }
193     }
194     else // cumprod on one dimension
195     {
196         int iSizeOfDimN = pIn->getDimsArray()[iOrientation - 1];
197         int iIncrement = 1;
198
199         for (int i = 0; i < iOrientation - 1; i++)
200         {
201             iIncrement *= pIn->getDimsArray()[i];
202         }
203
204         if (bComplex)
205         {
206             for (int j = 0; j < pIn->getSize(); j += (iIncrement * iSizeOfDimN))
207             {
208                 for (int i = j; i < iIncrement + j; i++) // set the first values in out
209                 {
210                     pdblReal = pIn->get(i)->getCoefReal();
211                     pdblImg = pIn->get(i)->getCoefImg();
212                     iRank = pIn->get(i)->getRank();
213                     pSP = new types::SinglePoly(&pdRData, &pdIData, iRank);
214
215                     for (int j = 0; j < iRank; j++)
216                     {
217                         pdRData[j] = pdblReal[j];
218                         pdIData[j] = pdblImg[j];
219                     }
220                     pOut->set(i, pSP);
221                     delete pSP;
222                 }
223
224                 for (int k = 1; k < iSizeOfDimN; k++) // make the cumprod for the next values
225                 {
226                     for (int i = (iIncrement * k) + j; i < (iIncrement * (k + 1)) + j; i++)
227                     {
228                         iLastRank = pOut->get(i - iIncrement)->getRank();
229                         pdblLastReal = pOut->get(i - iIncrement)->getCoefReal();
230                         pdblLastImg = pOut->get(i - iIncrement)->getCoefImg();
231                         iRank = pIn->get(i)->getRank();
232                         pdblReal = pIn->get(i)->getCoefReal();
233                         pdblImg = pIn->get(i)->getCoefImg();
234
235                         iOutRank = iLastRank + iRank - 1;
236                         pSP = new types::SinglePoly(&pdRData, &pdIData, iOutRank);
237
238                         iMultiComplexPolyByComplexPoly(pdblReal, pdblImg, iRank,
239                                                        pdblLastReal, pdblLastImg, iLastRank,
240                                                        pdRData, pdIData, iOutRank);
241
242                         pOut->set(i, pSP);
243                         delete pSP;
244                     }
245                 }
246             }
247         }
248         else
249         {
250             for (int j = 0; j < pIn->getSize(); j += (iIncrement * iSizeOfDimN))
251             {
252                 for (int i = j; i < iIncrement + j; i++) // set the first values in out
253                 {
254                     pdblReal = pIn->get(i)->getCoefReal();
255                     iRank = pIn->get(i)->getRank();
256                     pSP = new types::SinglePoly(&pdRData, iRank);
257
258                     for (int j = 0; j < iRank; j++)
259                     {
260                         pdRData[j] = pdblReal[j];
261                     }
262                     pOut->set(i, pSP);
263                     delete pSP;
264                 }
265
266                 for (int k = 1; k < iSizeOfDimN; k++) // make the cumprod for the next values
267                 {
268                     for (int i = (iIncrement * k) + j; i < (iIncrement * (k + 1)) + j; i++)
269                     {
270                         iLastRank = pOut->get(i - iIncrement)->getRank();
271                         pdblLastReal = pOut->get(i - iIncrement)->getCoefReal();
272                         iRank = pIn->get(i)->getRank();
273                         pdblReal = pIn->get(i)->getCoefReal();
274
275                         iOutRank = iLastRank + iRank - 1;
276                         pSP = new types::SinglePoly(&pdRData, iOutRank);
277
278                         iMultiScilabPolynomByScilabPolynom(pdblReal, iRank, pdblLastReal, iLastRank, pdRData, iOutRank);
279
280                         pOut->set(i, pSP);
281                         delete pSP;
282                     }
283                 }
284             }
285         }
286     }
287
288     pOut->updateRank();
289     return iErr;
290 }
291