elementary_functions module.
[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 - 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 "cumprod.hxx"
15
16 int cumprod(types::Double* pIn, int iOrientation, types::Double* pOut)
17 {
18     if (iOrientation == 0) // all
19     {
20         pOut->set(0, pIn->get(0));
21         for (int i = 1; i < pIn->getSize(); i++)
22         {
23             pOut->set(i, pIn->get(i) * pOut->get(i - 1));
24         }
25
26         if (pIn->isComplex())
27         {
28             pOut->setImg(0, pIn->getImg(0));
29             for (int i = 1; i < pIn->getSize(); i++)
30             {
31                 pOut->setImg(i, pIn->getImg(i) * pOut->getImg(i - 1));
32             }
33         }
34     }
35     else // cumprod on one dimension
36     {
37         int iSizeOfDimN = pIn->getDimsArray()[iOrientation - 1];
38         int iIncrement = 1;
39
40         for (int i = 0; i < iOrientation - 1; i++)
41         {
42             iIncrement *= pIn->getDimsArray()[i];
43         }
44
45         for (int j = 0; j < pIn->getSize(); j += (iIncrement * iSizeOfDimN))
46         {
47             for (int i = j; i < iIncrement + j; i++) // set the first values in out
48             {
49                 pOut->set(i, pIn->get(i));
50             }
51
52             for (int k = 1; k < iSizeOfDimN; k++) // make the cumprod for the next values
53             {
54                 for (int i = (iIncrement * k) + j; i < (iIncrement * (k + 1)) + j; i++)
55                 {
56                     pOut->set(i, pIn->get(i) * pOut->get(i - iIncrement));
57                 }
58             }
59         }
60         if (pIn->isComplex())
61         {
62             for (int j = 0; j < pIn->getSize(); j += (iIncrement * iSizeOfDimN))
63             {
64                 for (int i = j; i < iIncrement + j; i++) // set the first values in out
65                 {
66                     pOut->setImg(i, pIn->getImg(i));
67                 }
68
69                 for (int k = 1; k < iSizeOfDimN; k++) // make the cumprod for the next values
70                 {
71                     for (int i = (iIncrement * k) + j; i < (iIncrement * (k + 1)) + j; i++)
72                     {
73                         pOut->setImg(i, pIn->getImg(i) * pOut->getImg(i - iIncrement));
74                     }
75                 }
76             }
77         }
78     }
79
80     return 0;
81 }
82
83 // polynom
84 int cumprod(types::Polynom* pIn, int iOrientation, types::Polynom* pOut)
85 {
86     int iErr        = 0;
87     int iRank       = 0;
88     int iLastRank   = 0;
89     double* pdRData = NULL;
90     double* pdIData = NULL;
91     bool bComplex   = pIn->isComplex();
92
93     types::SinglePoly* pSP = NULL;
94
95     if (iOrientation == 0) // all
96     {
97         // set first element
98         iRank = pIn->get(0)->getRank();
99
100         if (bComplex)
101         {
102             pSP = new types::SinglePoly(&pdRData, &pdIData, iRank);
103         }
104         else
105         {
106             pSP = new types::SinglePoly(&pdRData, iRank);
107         }
108
109         for (int j = 0; j < iRank; j++)
110         {
111             pdRData[j] = pIn->get(0)->getCoefReal()[j];
112             if (bComplex)
113             {
114                 pdIData[j] = pIn->get(0)->getCoefImg()[j];
115             }
116         }
117
118         pOut->set(0, pSP);
119         iLastRank = iRank;
120
121         // compute next elements
122         for (int i = 1; i < pIn->getSize(); i++)
123         {
124             if (iRank > pIn->get(i)->getRank())
125             {
126                 iRank = pIn->get(i)->getRank();
127             }
128
129             if (bComplex)
130             {
131                 pSP = new types::SinglePoly(&pdRData, &pdIData, iRank);
132             }
133             else
134             {
135                 pSP = new types::SinglePoly(&pdRData, iRank);
136             }
137
138             for (int j = 0; j < iRank; j++)
139             {
140                 pdRData[j] = pIn->get(i)->getCoefReal()[j] * pOut->get(i - 1)->getCoefReal()[j];
141                 if (bComplex)
142                 {
143                     pdIData[j] = pIn->get(i)->getCoefImg()[j] * pOut->get(i - 1)->getCoefImg()[j];
144                 }
145             }
146
147             pOut->set(i, pSP);
148             iLastRank = iRank;
149         }
150     }
151     else // cumprod on one dimension
152     {
153         int iSizeOfDimN = pIn->getDimsArray()[iOrientation - 1];
154         int iIncrement = 1;
155
156         for (int i = 0; i < iOrientation - 1; i++)
157         {
158             iIncrement *= pIn->getDimsArray()[i];
159         }
160
161         for (int j = 0; j < pIn->getSize(); j += (iIncrement * iSizeOfDimN))
162         {
163             for (int i = j; i < iIncrement + j; i++) // set the first values in out
164             {
165                 iRank = pIn->get(i)->getRank();
166
167                 if (bComplex)
168                 {
169                     pSP = new types::SinglePoly(&pdRData, &pdIData, iRank);
170                 }
171                 else
172                 {
173                     pSP = new types::SinglePoly(&pdRData, iRank);
174                 }
175
176                 for (int j = 0; j < iRank; j++)
177                 {
178                     pdRData[j] = pIn->get(i)->getCoefReal()[j];
179                     if (bComplex)
180                     {
181                         pdIData[j] = pIn->get(i)->getCoefImg()[j];
182                     }
183                 }
184                 pOut->set(i, pSP);
185             }
186
187             for (int k = 1; k < iSizeOfDimN; k++) // make the cumprod for the next values
188             {
189                 for (int i = (iIncrement * k) + j; i < (iIncrement * (k + 1)) + j; i++)
190                 {
191                     iRank = pOut->get(i - iIncrement)->getRank();
192                     if (iRank > pIn->get(i)->getRank())
193                     {
194                         iRank = pIn->get(i)->getRank();
195                     }
196
197                     if (bComplex)
198                     {
199                         pSP = new types::SinglePoly(&pdRData, &pdIData, iRank);
200                     }
201                     else
202                     {
203                         pSP = new types::SinglePoly(&pdRData, iRank);
204                     }
205
206                     for (int j = 0; j < iRank; j++)
207                     {
208                         pdRData[j] = pIn->get(i)->getCoefReal()[j] * pOut->get(i - iIncrement)->getCoefReal()[j];
209                         if (bComplex)
210                         {
211                             pdIData[j] = pIn->get(i)->getCoefImg()[j] * pOut->get(i - iIncrement)->getCoefImg()[j];
212                         }
213                     }
214                     pOut->set(i, pSP);
215                 }
216             }
217         }
218     }
219
220     return iErr;
221 }
222