5180ba1138de8b7d62876a2bbb360b49ec169f38
[scilab.git] / scilab / modules / elementary_functions / src / cpp / cumsum.cpp
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2012 - DIGITEO - 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 "cumsum.hxx"
15
16 int cumsum(types::Double* pIn, int iOrientation, types::Double* pOut)
17 {
18     double* pdblInR = pIn->get();
19     double* pdblOutR = pOut->get();
20     double* pdblInI = pIn->getImg();
21     double* pdblOutI = pOut->getImg();
22
23     if (iOrientation == 0) // all
24     {
25         pdblOutR[0] = pdblInR[0];
26
27         if (pIn->isComplex())
28         {
29             pdblOutI[0] = pdblInI[0];
30             for (int i = 1; i < pIn->getSize(); i++)
31             {
32                 pdblOutR[i] = pdblOutR[i - 1] + pdblInR[i];
33                 pdblOutI[i] = pdblOutI[i - 1] + pdblInI[i];
34             }
35         }
36         else
37         {
38             for (int i = 1; i < pIn->getSize(); i++)
39             {
40                 pdblOutR[i] = pdblOutR[i - 1] + pdblInR[i];
41             }
42         }
43     }
44     else // cumsum on one dimension
45     {
46         int iSizeOfDimN = pIn->getDimsArray()[iOrientation - 1];
47         int iIncrement = 1;
48
49         for (int i = 0; i < iOrientation - 1; i++)
50         {
51             iIncrement *= pIn->getDimsArray()[i];
52         }
53
54         if (pIn->isComplex())
55         {
56             for (int j = 0; j < pIn->getSize(); j += (iIncrement * iSizeOfDimN))
57             {
58                 for (int i = j; i < iIncrement + j; i++) // set the first values in out
59                 {
60                     pdblOutR[i] = pdblInR[i];
61                     pdblOutI[i] = pdblInI[i];
62                 }
63
64                 for (int k = 1; k < iSizeOfDimN; k++) // make the cumsum for the next values
65                 {
66                     for (int i = (iIncrement * k) + j; i < (iIncrement * (k + 1)) + j; i++)
67                     {
68                         pdblOutR[i] = pdblInR[i] + pdblOutR[i - iIncrement];
69                         pdblOutI[i] = pdblInI[i] + pdblOutI[i - iIncrement];
70                     }
71                 }
72             }
73         }
74         else
75         {
76             for (int j = 0; j < pIn->getSize(); j += (iIncrement * iSizeOfDimN))
77             {
78                 for (int i = j; i < iIncrement + j; i++) // set the first values in out
79                 {
80                     pdblOutR[i] = pdblInR[i];
81                 }
82
83                 for (int k = 1; k < iSizeOfDimN; k++) // make the cumsum for the next values
84                 {
85                     for (int i = (iIncrement * k) + j; i < (iIncrement * (k + 1)) + j; i++)
86                     {
87                         pdblOutR[i] = pdblInR[i] + pdblOutR[i - iIncrement];
88                     }
89                 }
90             }
91         }
92     }
93
94     return 0;
95 }
96
97 // polynom
98 int cumsum(types::Polynom* pIn, int iOrientation, types::Polynom* pOut)
99 {
100     int iErr        = 0;
101     int iRank       = 0;
102     int iOutRank    = 0;
103     int iLastRank   = 0;
104     int iMin        = 0;
105
106     double* pdRData         = NULL;
107     double* pdIData         = NULL;
108     double* pdblReal        = NULL;
109     double* pdblImg         = NULL;
110     double* pdblLastReal    = NULL;
111     double* pdblLastImg     = NULL;
112
113     bool bComplex   = pIn->isComplex();
114
115     types::SinglePoly* pSP = NULL;
116
117     if (iOrientation == 0) // all
118     {
119         // set first element
120         iRank = pIn->get(0)->getRank();
121         pdblReal = pIn->get(0)->getCoefReal();
122         if (bComplex)
123         {
124             pSP = new types::SinglePoly(&pdRData, &pdIData, iRank);
125             pdblImg = pIn->get(0)->getCoefImg();
126             for (int j = 0; j < iRank; j++)
127             {
128                 pdRData[j] = pdblReal[j];
129                 pdIData[j] = pdblImg[j];
130             }
131         }
132         else
133         {
134             pSP = new types::SinglePoly(&pdRData, iRank);
135             for (int j = 0; j < iRank; j++)
136             {
137                 pdRData[j] = pdblReal[j];
138             }
139         }
140
141         pOut->set(0, pSP);
142         iLastRank = iRank;
143         pdblLastReal = pdblReal;
144         pdblLastImg = pdblImg;
145
146         // compute next elements
147         if (bComplex)
148         {
149             for (int i = 1; i < pIn->getSize(); i++)
150             {
151                 pdblReal = pIn->get(i)->getCoefReal();
152                 pdblImg  = pIn->get(i)->getCoefImg();
153                 iRank    = pIn->get(i)->getRank();
154
155                 iOutRank = std::max(iRank, iLastRank);
156                 iMin     = std::min(iRank, iLastRank);
157
158                 pSP = new types::SinglePoly(&pdRData, &pdIData, iOutRank);
159                 for (int j = 0; j < iMin; j++)
160                 {
161                     pdRData[j] = pdblReal[j] + pdblLastReal[j];
162                     pdIData[j] = pdblImg[j]  + pdblLastImg[j];
163                 }
164
165                 if (iOutRank == iRank)
166                 {
167                     for (int j = iMin; j < iOutRank; j++)
168                     {
169                         pdRData[j] = pdblReal[j];
170                         pdIData[j] = pdblImg[j];
171                     }
172                 }
173                 else
174                 {
175                     for (int j = iMin; j < iOutRank; j++)
176                     {
177                         pdRData[j] = pdblLastReal[j];
178                         pdIData[j] = pdblLastImg[j];
179                     }
180                 }
181
182                 pOut->set(i, pSP);
183                 pdblLastReal = pdRData;
184                 pdblLastImg = pdIData;
185                 iLastRank = iOutRank;
186                 delete pSP;
187             }
188         }
189         else
190         {
191             for (int i = 1; i < pIn->getSize(); i++)
192             {
193                 pdblReal = pIn->get(i)->getCoefReal();
194                 iRank    = pIn->get(i)->getRank();
195
196                 iOutRank = std::max(iRank, iLastRank);
197                 iMin     = std::min(iRank, iLastRank);
198
199                 pSP = new types::SinglePoly(&pdRData, iOutRank);
200                 for (int j = 0; j < iMin; j++)
201                 {
202                     pdRData[j] = pdblReal[j] + pdblLastReal[j];
203                 }
204
205                 if (iOutRank == iRank)
206                 {
207                     for (int j = iMin; j < iOutRank; j++)
208                     {
209                         pdRData[j] = pdblReal[j];
210                     }
211                 }
212                 else
213                 {
214                     for (int j = iMin; j < iOutRank; j++)
215                     {
216                         pdRData[j] = pdblLastReal[j];
217                     }
218                 }
219
220                 pOut->set(i, pSP);
221                 pdblLastReal = pdRData;
222                 iLastRank = iOutRank;
223                 delete pSP;
224             }
225         }
226     }
227     else // cumsum on one dimension
228     {
229         int iSizeOfDimN = pIn->getDimsArray()[iOrientation - 1];
230         int iIncrement = 1;
231
232         for (int i = 0; i < iOrientation - 1; i++)
233         {
234             iIncrement *= pIn->getDimsArray()[i];
235         }
236
237         if (bComplex)
238         {
239             for (int j = 0; j < pIn->getSize(); j += (iIncrement * iSizeOfDimN))
240             {
241                 for (int i = j; i < iIncrement + j; i++) // set the first values in out
242                 {
243                     pdblReal = pIn->get(i)->getCoefReal();
244                     pdblImg  = pIn->get(i)->getCoefImg();
245                     iRank    = pIn->get(i)->getRank();
246
247                     pSP = new types::SinglePoly(&pdRData, &pdIData, iRank);
248
249                     for (int j = 0; j < iRank; j++)
250                     {
251                         pdRData[j] = pdblReal[j];
252                         pdIData[j] = pdblImg[j];
253                     }
254                     pOut->set(i, pSP);
255                     delete pSP;
256                 }
257
258                 for (int k = 1; k < iSizeOfDimN; k++) // make the cumsum for the next values
259                 {
260                     for (int i = (iIncrement * k) + j; i < (iIncrement * (k + 1)) + j; i++)
261                     {
262                         iLastRank    = pOut->get(i - iIncrement)->getRank();
263                         pdblLastReal = pOut->get(i - iIncrement)->getCoefReal();
264                         pdblLastImg  = pOut->get(i - iIncrement)->getCoefImg();
265
266                         iRank    = pIn->get(i)->getRank();
267                         pdblReal = pIn->get(i)->getCoefReal();
268                         pdblImg  = pIn->get(i)->getCoefImg();
269
270                         iOutRank = std::max(iRank, iLastRank);
271                         iMin     = std::min(iRank, iLastRank);
272
273                         pSP = new types::SinglePoly(&pdRData, &pdIData, iOutRank);
274                         for (int j = 0; j < iMin; j++)
275                         {
276                             pdRData[j] = pdblReal[j] + pdblLastReal[j];
277                             pdIData[j] = pdblImg[j]  + pdblLastImg[j];
278                         }
279
280                         if (iOutRank == iRank)
281                         {
282                             for (int j = iMin; j < iOutRank; j++)
283                             {
284                                 pdRData[j] = pdblReal[j];
285                                 pdIData[j] = pdblImg[j];
286                             }
287                         }
288                         else
289                         {
290                             for (int j = iMin; j < iOutRank; j++)
291                             {
292                                 pdRData[j] = pdblLastReal[j];
293                                 pdIData[j] = pdblLastImg[j];
294                             }
295                         }
296                         pOut->set(i, pSP);
297                         delete pSP;
298                     }
299                 }
300             }
301         }
302         else
303         {
304             for (int j = 0; j < pIn->getSize(); j += (iIncrement * iSizeOfDimN))
305             {
306                 for (int i = j; i < iIncrement + j; i++) // set the first values in out
307                 {
308                     pdblReal = pIn->get(i)->getCoefReal();
309                     iRank    = pIn->get(i)->getRank();
310
311                     pSP = new types::SinglePoly(&pdRData, iRank);
312
313                     for (int j = 0; j < iRank; j++)
314                     {
315                         pdRData[j] = pdblReal[j];
316                     }
317                     pOut->set(i, pSP);
318                     delete pSP;
319                 }
320
321                 for (int k = 1; k < iSizeOfDimN; k++) // make the cumsum for the next values
322                 {
323                     for (int i = (iIncrement * k) + j; i < (iIncrement * (k + 1)) + j; i++)
324                     {
325                         iLastRank    = pOut->get(i - iIncrement)->getRank();
326                         pdblLastReal = pOut->get(i - iIncrement)->getCoefReal();
327
328                         iRank    = pIn->get(i)->getRank();
329                         pdblReal = pIn->get(i)->getCoefReal();
330
331                         iOutRank = std::max(iRank, iLastRank);
332                         iMin     = std::min(iRank, iLastRank);
333
334                         pSP = new types::SinglePoly(&pdRData, iOutRank);
335                         for (int j = 0; j < iMin; j++)
336                         {
337                             pdRData[j] = pdblReal[j] + pdblLastReal[j];
338                         }
339
340                         if (iOutRank == iRank)
341                         {
342                             for (int j = iMin; j < iOutRank; j++)
343                             {
344                                 pdRData[j] = pdblReal[j];
345                             }
346                         }
347                         else
348                         {
349                             for (int j = iMin; j < iOutRank; j++)
350                             {
351                                 pdRData[j] = pdblLastReal[j];
352                             }
353                         }
354                         pOut->set(i, pSP);
355                         delete pSP;
356                     }
357                 }
358             }
359         }
360     }
361
362     return iErr;
363 }
364