932fc517b8c0d8abd2fe4fd5b8c738a2106ce325
[scilab.git] / scilab / modules / elementary_functions / src / cpp / sum.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 "sum.hxx"
15
16 types::Double* sum(types::Double* pIn, int iOrientation)
17 {
18     types::Double* pOut = NULL;
19     double* pdblInReal  = pIn->getReal();
20     double* pdblInImg   = pIn->getImg();
21
22     if (iOrientation == 0) // all
23     {
24         double dblR = 0;
25         double dblI = 0;
26
27         if (pIn->isComplex())
28         {
29             for (int i = 0 ; i < pIn->getSize() ; i++)
30             {
31                 dblR += pdblInReal[i];
32                 dblI += pdblInImg[i];
33             }
34
35             pOut = new types::Double(dblR, dblI);
36         }
37         else
38         {
39             for (int i = 0 ; i < pIn->getSize() ; i++)
40             {
41                 dblR += pdblInReal[i];
42             }
43
44             pOut = new types::Double(dblR);
45         }
46     }
47     else // sum on one dimension
48     {
49         // create output dimensions
50         int iDims   = pIn->getDims();
51         int* piDims = new int[iDims];
52
53         for (int i = 0 ; i < iDims ; i++)
54         {
55             piDims[i] = pIn->getDimsArray()[i];
56         }
57
58         piDims[iOrientation - 1] = 1;
59
60         // create output variable
61         pOut = new types::Double(iDims, piDims, pIn->isComplex());
62         // init output
63         pOut->setZeros();
64
65         delete piDims;
66         piDims = NULL;
67
68         // perform operations
69         double* pdblOut = pOut->get();
70         double* pdblOutImg = NULL;
71         int* piIndex = new int[iDims];
72
73         if (pIn->isComplex())
74         {
75             pdblOutImg = pOut->getImg();
76             for (int i = 0 ; i < pIn->getSize() ; i++)
77             {
78                 //get array of dim
79                 pIn->getIndexes(i, piIndex);
80
81                 //convert indexes for result
82                 piIndex[iOrientation - 1] = 0;
83                 int iIndex = pOut->getIndex(piIndex);
84                 pdblOutImg[iIndex]  += pdblInImg[i];
85                 pdblOut[iIndex]     += pdblInReal[i];
86             }
87         }
88         else
89         {
90             for (int i = 0 ; i < pIn->getSize() ; i++)
91             {
92                 //get array of dim
93                 pIn->getIndexes(i, piIndex);
94
95                 //convert indexes for result
96                 piIndex[iOrientation - 1] = 0;
97                 int iIndex = pOut->getIndex(piIndex);
98                 pdblOut[iIndex] += pdblInReal[i];
99             }
100         }
101
102         delete piIndex;
103         piIndex = NULL;
104     }
105
106     return pOut;
107 }
108
109 // polynom
110 types::Polynom* sum(types::Polynom* pIn, int iOrientation)
111 {
112     types::Polynom* pOut = NULL;
113     if (iOrientation == 0) // sum of all element
114     {
115         // get rank of all input single poly
116         int* piRanks = new int[pIn->getSize()];
117         pIn->getRank(piRanks);
118
119         // create output matrix of poly
120         int iRankMax = pIn->getMaxRank();
121         pOut = new types::Polynom(pIn->getVariableName(), 1, 1, &iRankMax);
122         pOut->setComplex(pIn->isComplex());
123         pOut->setZeros();
124
125         // do sum
126         double* dblRealOut = pOut->get(0)->get();
127         if (pIn->isComplex())
128         {
129             double* dblImgOut = pOut->get(0)->get();
130             // perform operations
131             for (int i = 0; i < pIn->getSize(); i++)
132             {
133                 double* dblRealIn = pIn->get(i)->get();
134                 double* dblImgIn = pIn->get(i)->getImg();
135                 for (int iRankN = 0; iRankN < piRanks[i] + 1; iRankN++)
136                 {
137                     dblRealOut[iRankN] += dblRealIn[iRankN];
138                     dblImgOut[iRankN]  += dblImgIn[iRankN];
139                 }
140             }
141         }
142         else
143         {
144             // perform operations
145             for (int i = 0; i < pIn->getSize(); i++)
146             {
147                 double* dblRealIn = pIn->get(i)->get();
148                 for (int iRankN = 0; iRankN < piRanks[i] + 1; iRankN++)
149                 {
150                     dblRealOut[iRankN] += dblRealIn[iRankN];
151                 }
152             }
153         }
154     }
155     else // sum on one dimension
156     {
157         // create output dimensions
158         int iDims = pIn->getDims();
159         int* piDims = new int[iDims];
160
161         for (int i = 0 ; i < iDims ; i++)
162         {
163             piDims[i] = pIn->getDimsArray()[i];
164         }
165
166         piDims[iOrientation - 1] = 1;
167
168         // get ranks of all polynom
169         int* piRanks = new int[pIn->getSize()];
170         pIn->getRank(piRanks);
171
172         // get input ranks in types::Double
173         types::Double* pDblRanks = new types::Double(pIn->getDims(), pIn->getDimsArray());
174         for (int i = 0; i < pDblRanks->getSize(); i++)
175         {
176             pDblRanks->set(i, static_cast<double>(piRanks[i]));
177         }
178
179         // create output max ranks
180         types::Double* pDblRanksOut = new types::Double(iDims, piDims);
181         pDblRanksOut->setZeros();
182
183         // compute the maximum ranks for the dim n
184         int*    piIndex = new int[iDims];
185         double* pdblIn  = pDblRanks->get();
186         double* pdblOut = pDblRanksOut->get();
187         for (int i = 0 ; i < pDblRanks->getSize() ; i++)
188         {
189             //get array of dim
190             pDblRanks->getIndexes(i, piIndex);
191
192             //convert indexes for result
193             piIndex[iOrientation - 1] = 0;
194             int iIndex = pDblRanksOut->getIndex(piIndex);
195             pdblOut[iIndex] = std::max(pdblOut[iIndex], pdblIn[i]);
196         }
197
198         // move output ranks from types::Double to int*
199         int* piRankMax = new int[pDblRanksOut->getSize()];
200         for (int i = 0; i < pDblRanksOut->getSize(); i++)
201         {
202             piRankMax[i] = static_cast<int>(pdblOut[i]);
203         }
204
205         // create the outpout polynom
206         pOut = new types::Polynom(pIn->getVariableName(), iDims, piDims, piRankMax);
207         pOut->setComplex(pIn->isComplex());
208         pOut->setZeros();
209
210         if (pIn->isComplex())
211         {
212             for (int i = 0 ; i < pIn->getSize() ; i++)
213             {
214                 //get array of dim
215                 pIn->getIndexes(i, piIndex);
216
217                 //convert indexes for result
218                 piIndex[iOrientation - 1] = 0;
219                 int iIndex = pOut->getIndex(piIndex);
220
221                 // make the sum for each ranks
222                 double* dblRealIn = pIn->get(i)->get();
223                 double* dblRealOut = pOut->get(iIndex)->get();
224                 double* dblImgIn = pIn->get(i)->getImg();
225                 double* dblImgOut = pOut->get(iIndex)->getImg();
226                 for (int iRankN = 0; iRankN < piRanks[i] + 1; iRankN++)
227                 {
228                     dblRealOut[iRankN] += dblRealIn[iRankN];
229                     dblImgOut[iRankN]  += dblImgIn[iRankN];
230                 }
231             }
232         }
233         else
234         {
235             for (int i = 0 ; i < pIn->getSize() ; i++)
236             {
237                 //get array of dim
238                 pIn->getIndexes(i, piIndex);
239
240                 //convert indexes for result
241                 piIndex[iOrientation - 1] = 0;
242                 int iIndex = pOut->getIndex(piIndex);
243
244                 // make sum on each ranks
245                 double* dblRealIn = pIn->get(i)->get();
246                 double* dblRealOut = pOut->get(iIndex)->get();
247                 for (int iRankN = 0; iRankN < piRanks[i] + 1; iRankN++)
248                 {
249                     dblRealOut[iRankN] += dblRealIn[iRankN];
250                 }
251             }
252         }
253     }
254
255     return pOut;
256 }