Coverity: ast module resource leaks fixed
[scilab.git] / scilab / modules / ast / src / cpp / operations / types_power.cpp
1 /*
2 *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 *  Copyright (C) 2008-2008 - DIGITEO - Antoine ELIAS
4 *
5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
6  *
7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
9  * This file was originally licensed under the terms of the CeCILL v2.1,
10  * and continues to be available under such terms.
11  * For more information, see the COPYING file which you should have received
12  * along with this program.
13 *
14 */
15
16 #include "types_power.hxx"
17 #include "types_multiplication.hxx"
18
19 extern "C"
20 {
21 #include "operation_f.h"
22 #include "matrix_power.h"
23 #include "localization.h"
24 #include "charEncoding.h"
25 }
26
27 using namespace types;
28
29 InternalType *GenericPower(InternalType *_pLeftOperand, InternalType *_pRightOperand)
30 {
31     InternalType *pResult = NULL;
32
33     /*
34     ** DOUBLE ^ DOUBLE
35     ** DOUBLE ** DOUBLE
36     */
37     if (_pLeftOperand->isDouble() && _pRightOperand->isDouble())
38     {
39         Double *pL   = _pLeftOperand->getAs<Double>();
40         Double *pR   = _pRightOperand->getAs<Double>();
41
42         int iResult = PowerDoubleByDouble(pL, pR, (Double**)&pResult);
43         if (iResult)
44         {
45             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
46         }
47
48         return pResult;
49     }
50
51     /*
52     ** POLY ^ DOUBLE
53     ** POLY ** DOUBLE
54     */
55     if (_pLeftOperand->isPoly() && _pRightOperand->isDouble())
56     {
57         Polynom *pL   = _pLeftOperand->getAs<Polynom>();
58         Double *pR    = _pRightOperand->getAs<Double>();
59
60         int iResult = PowerPolyByDouble(pL, pR, &pResult);
61         switch (iResult)
62         {
63             case 1 :
64                 throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
65             case 2 :
66                 throw ast::InternalError(_W("Invalid exponent: expected real exponents.\n"));
67             default:
68                 //OK
69                 break;
70         }
71
72         return pResult;
73     }
74
75     /*
76     ** Default case : Return NULL will Call Overloading.
77     */
78     return NULL;
79
80 }
81
82 InternalType *GenericDotPower(InternalType *_pLeftOperand, InternalType *_pRightOperand)
83 {
84     InternalType *pResult = NULL;
85     GenericType::ScilabType TypeL = _pLeftOperand->getType();
86     GenericType::ScilabType TypeR = _pRightOperand->getType();
87
88     /*
89     ** DOUBLE .^ DOUBLE
90     ** DOUBLE .** DOUBLE
91     */
92     if (TypeL == GenericType::ScilabDouble && TypeR == GenericType::ScilabDouble)
93     {
94         Double *pL  = _pLeftOperand->getAs<Double>();
95         Double *pR  = _pRightOperand->getAs<Double>();
96
97         int iResult = DotPowerDoubleByDouble(pL, pR, (Double**)&pResult);
98         if (iResult)
99         {
100             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
101         }
102
103         return pResult;
104     }
105     /*
106     ** SPARSE  .^ DOUBLE
107     ** SPARSE .** DOUBLE
108     */
109     if (TypeL == GenericType::ScilabSparse && TypeR == GenericType::ScilabDouble)
110     {
111         types::Sparse *pL = _pLeftOperand->getAs<types::Sparse>();
112         Double *pR = _pRightOperand->getAs<Double>();
113         int iResult = DotPowerSpaseByDouble(pL, pR, &pResult);
114         if (iResult)
115         {
116             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
117         }
118         return pResult;
119
120     }
121
122     /*
123     ** POLY .^ DOUBLE
124     ** POLY .** DOUBLE
125     */
126     if (TypeL == GenericType::ScilabPolynom && TypeR == GenericType::ScilabDouble)
127     {
128         Polynom *pL   = _pLeftOperand->getAs<Polynom>();
129         Double *pR   = _pRightOperand->getAs<Double>();
130
131         int iResult = DotPowerPolyByDouble(pL, pR, &pResult);
132         switch (iResult)
133         {
134             case 1 :
135                 throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
136             case 2 :
137                 throw ast::InternalError(_W("Invalid exponent: expected real exponents.\n"));
138             default:
139                 //OK
140                 break;
141         }
142         return pResult;
143     }
144
145     /*
146     ** Default case : Return NULL will Call Overloading.
147     */
148     return NULL;
149
150 }
151
152 int PowerDoubleByDouble(Double* _pDouble1, Double* _pDouble2, Double** _pDoubleOut)
153 {
154     bool bComplex1  = _pDouble1->isComplex();
155     bool bComplex2  = _pDouble2->isComplex();
156     bool bScalar1   = _pDouble1->isScalar();
157     bool bScalar2   = _pDouble2->isScalar();
158
159     int iComplex = 1;
160
161     if (bScalar1)
162     {
163         if (bScalar2)
164         {
165             //s ^ s
166             *_pDoubleOut = new Double(1, 1, true);
167
168             if (bComplex1 == false && bComplex2 == false)
169             {
170                 iPowerRealScalarByRealScalar(_pDouble1->get(0), _pDouble2->get(0), (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), &iComplex);
171             }
172             else if (bComplex1 == false && bComplex2 == true)
173             {
174                 iPowerRealScalarByComplexScalar(_pDouble1->get(0), _pDouble2->get(0), _pDouble2->getImg(0), (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
175             }
176             else if (bComplex1 == true && bComplex2 == false)
177             {
178                 iPowerComplexScalarByRealScalar(_pDouble1->get(0), _pDouble1->getImg(0), _pDouble2->get(0), (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
179             }
180             else if (bComplex1 == true && bComplex2 == true)
181             {
182                 iPowerComplexScalarByComplexScalar(_pDouble1->get(0), _pDouble1->getImg(0), _pDouble2->get(0), _pDouble2->getImg(0), (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
183             }
184
185             if (iComplex == 0)
186             {
187                 (*_pDoubleOut)->setComplex(false);
188             }
189
190             return 0;
191         }
192         else if (_pDouble2->isVector())
193         {
194             //s ^ Vector == s .^ Vector
195             *_pDoubleOut = new Double(_pDouble2->getDims(), _pDouble2->getDimsArray(), true);
196
197             if (bComplex1 == false && bComplex2 == false)
198             {
199                 iPowerRealScalarByRealMatrix(
200                     _pDouble1->get(0),
201                     _pDouble2->get(), 1, _pDouble2->getSize(),
202                     (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), &iComplex);
203             }
204             else if (bComplex1 == false && bComplex2 == true)
205             {
206                 iPowerRealScalarByComplexMatrix(
207                     _pDouble1->get(0),
208                     _pDouble2->get(), _pDouble2->getImg(), 1, _pDouble2->getSize(),
209                     (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
210             }
211             else if (bComplex1 == true && bComplex2 == false)
212             {
213                 iPowerComplexScalarByRealMatrix(
214                     _pDouble1->get(0), _pDouble1->getImg(0),
215                     _pDouble2->get(), 1, _pDouble2->getSize(),
216                     (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
217             }
218             else if (bComplex1 == true && bComplex2 == true)
219             {
220                 iPowerComplexScalarByComplexMatrix(
221                     _pDouble1->get(0), _pDouble1->getImg(0),
222                     _pDouble2->get(), _pDouble2->getImg(), 1, _pDouble2->getSize(),
223                     (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
224             }
225
226             if (iComplex == 0)
227             {
228                 (*_pDoubleOut)->setComplex(false);
229             }
230
231             return 0;
232         }
233         else
234         {
235             // s ^ Matrix
236             // call overload %s_pow
237             return 0;
238         }
239     }
240
241     if (bScalar2)
242     {
243         if (_pDouble1->isVector())
244         {
245             //_pDouble1 is a vector and _pDouble is a scalar
246             *_pDoubleOut = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray() , true);
247
248             if (bComplex1 == false && bComplex2 == false)
249             {
250                 for (int i = 0 ; i < (*_pDoubleOut)->getSize() ; i++)
251                 {
252                     iPowerRealScalarByRealScalar(
253                         _pDouble1->get(i),
254                         _pDouble2->get(0),
255                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
256                 }
257             }
258             else if (bComplex1 == false && bComplex2 == true)
259             {
260                 for (int i = 0 ; i < (*_pDoubleOut)->getSize() ; i++)
261                 {
262                     iPowerRealScalarByComplexScalar(
263                         _pDouble1->get(i),
264                         _pDouble2->get(0), _pDouble2->getImg(0),
265                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
266                 }
267             }
268             else if (bComplex1 == true && bComplex2 == false)
269             {
270                 for (int i = 0 ; i < (*_pDoubleOut)->getSize() ; i++)
271                 {
272                     iPowerComplexScalarByRealScalar(
273                         _pDouble1->get(i), _pDouble1->getImg(i),
274                         _pDouble2->get(0),
275                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
276                 }
277             }
278             else if (bComplex1 == true && bComplex2 == true)
279             {
280                 for (int i = 0 ; i < (*_pDoubleOut)->getSize() ; i++)
281                 {
282                     iPowerComplexScalarByComplexScalar(
283                         _pDouble1->get(i), _pDouble1->getImg(i),
284                         _pDouble2->get(0), _pDouble2->getImg(0),
285                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
286                 }
287             }
288
289             if (iComplex == 0)
290             {
291                 (*_pDoubleOut)->setComplex(false);
292             }
293
294             return 0;
295         }
296         else if (_pDouble1->getRows() == _pDouble1->getCols())
297         {
298             //power of a square matrix by a scalar exponent.
299             int iRet = 0;
300             if (bComplex2)
301             {
302                 // mange by overloading
303                 // Call %s_pow
304                 return 0;
305             }
306
307             *_pDoubleOut = new Double(_pDouble1->getRows(), _pDouble1->getCols() , true);
308             if (bComplex1 == false)
309             {
310                 iRet = iPowerRealSquareMatrixByRealScalar(
311                            _pDouble1->get(), _pDouble1->getRows(), _pDouble1->getCols(),
312                            _pDouble2->get(0),
313                            (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), &iComplex);
314             }
315             else if (bComplex1 == true)
316             {
317                 iRet = iPowerComplexSquareMatrixByRealScalar(
318                            _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
319                            _pDouble2->get(0),
320                            (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
321             }
322
323             // call overload
324             if (iRet == -1)
325             {
326                 delete *_pDoubleOut;
327                 *_pDoubleOut = NULL;
328                 return 0;
329             }
330
331             if (iComplex == 0)
332             {
333                 (*_pDoubleOut)->setComplex(false);
334             }
335         }
336     }
337     return 0;
338 }
339
340 int PowerPolyByDouble(Polynom* _pPoly, Double* _pDouble, InternalType** _pOut)
341 {
342     bool bComplex1  = _pPoly->isComplex();
343     bool bScalar1   = _pPoly->isScalar();
344     double* bImg    = _pDouble->getImg();
345     bool bNumericallyComplex1 = _pDouble->isNumericallyComplex();
346
347     if (!bNumericallyComplex1)
348     {
349         return 2;
350     }
351
352     if (_pDouble->isEmpty())
353     {
354         //p ** []
355         *_pOut = Double::Empty();
356         return 0;
357     }
358
359     if (bScalar1)
360     {
361         //p ^ x or p ^ X
362         int iRank   = 0;
363         int* piRank = new int[_pDouble->getSize()];
364
365         _pPoly->getRank(&iRank);
366         for (int i = 0 ; i < _pDouble->getSize() ; i++)
367         {
368             int iInputRank = (int)_pDouble->get(i);
369             if (iInputRank < 0)
370             {
371                 //call overload
372                 _pOut = NULL;
373                 delete[] piRank;
374                 return 0;
375             }
376
377             piRank[i] = iRank * iInputRank;
378         }
379
380         Polynom* pOut = new Polynom(_pPoly->getVariableName(), _pDouble->getRows(), _pDouble->getCols(), piRank);
381         delete[] piRank;
382         pOut->setComplex(bComplex1);
383
384         for (int i = 0 ; i < _pDouble->getSize() ; i++)
385         {
386             SinglePoly* pCoeffOut = pOut->get(i);
387
388             int iCurrentRank    = 0;
389             int iLoop           = (int)_pDouble->get(i);
390
391             //initialize Out to 1
392             pCoeffOut->set(0, 1);
393             //get a copy of p
394             Polynom* pP = _pPoly->clone()->getAs<Polynom>();
395             pP->setComplex(_pPoly->isComplex());
396
397             while (iLoop)
398             {
399                 SinglePoly* ps = pP->get()[0];
400                 if (iLoop % 2)
401                 {
402                     int iRank = pP->getMaxRank();
403                     if (bComplex1)
404                     {
405                         C2F(wpmul1)(pCoeffOut->get(), pCoeffOut->getImg(), &iCurrentRank, ps->get(), ps->getImg(), &iRank, pCoeffOut->get(), pCoeffOut->getImg());
406                     }
407                     else
408                     {
409                         C2F(dpmul1)(pCoeffOut->get(), &iCurrentRank, ps->get(), &iRank, pCoeffOut->get());
410                     }
411                     iCurrentRank += iRank;
412                 }
413
414                 iLoop /= 2;
415                 if (iLoop)
416                 {
417                     //p = p * p
418                     Polynom* pTemp = NULL;
419                     MultiplyPolyByPoly(pP, pP, &pTemp);
420                     pP->killMe();
421                     pP = pTemp;
422                 }
423             }
424
425             pP->killMe();
426         }
427         *_pOut = pOut;
428     }
429     return 0;
430 }
431
432 int DotPowerSpaseByDouble(Sparse* _pSp, Double* _pDouble, InternalType** _pOut)
433 {
434     if (_pDouble->isEmpty())
435     {
436         //sp .^ []
437         *_pOut = Double::Empty();
438         return 0;
439     }
440
441     size_t iSize = _pSp->nonZeros();
442     int* Col = new int[iSize];
443     int* Row = new int[_pSp->getRows()];
444     _pSp->getColPos(Col);
445     _pSp->getNbItemByRow(Row);
446     int* iPositVal = new int[iSize];
447
448     int j = 0;
449     for (int i = 0; i < iSize;  j++)
450     {
451         for (int k = 0; k < Row[j]; k++)
452         {
453             iPositVal[i] = (Col[i] - 1) * _pSp->getRows() + j;
454             i++;
455         }
456     }
457
458     Double** pDbl = new Double*[iSize];
459     Double** pDblSp = new Double*[iSize];
460     double* pdbl = _pDouble->get();
461
462     if (_pDouble->isScalar())
463     {
464         if (_pDouble->isComplex())
465         {
466             double* pdblImg = _pDouble->getImg();
467             for (int i = 0; i < iSize; i++)
468             {
469                 pDbl[i] = new Double(pdbl[0], pdblImg[0]);
470                 pDblSp[i] = new Double(_pSp->get(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
471             }
472         }
473         else
474         {
475             for (int i = 0; i < iSize; i++)
476             {
477                 pDbl[i] = new Double(pdbl[0]);
478                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
479             }
480         }
481     }
482     else if (_pDouble->getSize() == iSize)
483     {
484         if (_pDouble->isComplex())
485         {
486             double* pdblImg = _pDouble->getImg();
487             for (int i = 0; i < iSize; i++)
488             {
489                 pDbl[i] = new Double(pdbl[i], pdblImg[i]);
490                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
491             }
492         }
493         else
494         {
495             for (int i = 0; i < iSize; i++)
496             {
497                 pDbl[i] = new Double(pdbl[i]);
498                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
499             }
500         }
501     }
502     else
503     {
504         delete[] pDblSp;
505         delete[] pDbl;
506         throw ast::InternalError(_W("Invalid exponent.\n"));
507         return 1;
508     }
509
510     Sparse* pSpTemp = new Sparse(_pSp->getRows(), _pSp->getCols(), _pSp->isComplex() || _pDouble->isComplex());
511     pSpTemp->zero_set();
512
513     Double* ppDblGet = NULL;
514     for (int i = 0; i < iSize; i++)
515     {
516         if ((pDblSp[i]->get(0) != 0) || (pDblSp[i]->getImg(0) != 0))
517         {
518             DotPowerDoubleByDouble(pDblSp[i], pDbl[i], &ppDblGet);
519             std::complex<double> cplx(ppDblGet->get(0), ppDblGet->getImg(0));
520             pSpTemp->set(iPositVal[i], cplx, false);
521         }
522     }
523
524     delete[] Col;
525     delete[] Row;
526     delete[] iPositVal;
527     for (int i = 0; i < iSize; i++)
528     {
529         delete pDbl[i];
530         delete pDblSp[i];
531     }
532     delete pDbl;
533     delete pDblSp;
534
535     pSpTemp->finalize();
536     *_pOut = pSpTemp;
537     return 0;
538
539 }
540
541
542 int DotPowerPolyByDouble(Polynom* _pPoly, Double* _pDouble, InternalType** _pOut)
543 {
544     if (_pDouble->isEmpty())
545     {
546         //p .^ []
547         *_pOut = Double::Empty();
548         return 0;
549     }
550
551     int iSize = _pPoly->getSize();
552     if (_pPoly->isScalar())
553     {
554         return PowerPolyByDouble(_pPoly, _pDouble, _pOut);
555     }
556
557     Double** pDblPower = new Double*[iSize];
558     double* pdblPower = _pDouble->get();
559     if (_pDouble->isScalar())
560     {
561         if (pdblPower[0] < 0)
562         {
563             //call overload
564             _pOut = NULL;
565             delete[] pDblPower;
566             return 0;
567         }
568
569         for (int i = 0; i < iSize; i++)
570         {
571             pDblPower[i] = new Double(pdblPower[0]);
572         }
573     }
574     else if (_pDouble->getSize() == iSize)
575     {
576         for (int i = 0; i < iSize; i++)
577         {
578             if (pdblPower[i] < 0)
579             {
580                 //call overload
581                 _pOut = NULL;
582                 delete[] pDblPower;
583                 return 0;
584             }
585
586             pDblPower[i] = new Double(pdblPower[i]);
587         }
588     }
589     else
590     {
591         delete[] pDblPower;
592         throw ast::InternalError(_W("Invalid exponent.\n"));
593     }
594
595     InternalType* pITTempOut    = NULL;
596     Polynom* pPolyTemp          = new Polynom(_pPoly->getVariableName(), 1, 1);
597     Polynom* pPolyOut           = new Polynom(_pPoly->getVariableName(), _pPoly->getDims(), _pPoly->getDimsArray());
598     SinglePoly** pSPOut         = pPolyOut->get();
599     SinglePoly** pSPTemp        = pPolyTemp->get();
600     SinglePoly** pSP            = _pPoly->get();
601
602     int iResult = 0;
603     for (int i = 0; i < iSize; i++)
604     {
605         // set singlePoly of _pPoly in pPolyTemp without copy
606         pSPTemp[0] = pSP[i];
607         iResult = PowerPolyByDouble(pPolyTemp, pDblPower[i], &pITTempOut);
608         if (iResult)
609         {
610             break;
611         }
612
613         // get singlePoly of pITTempOut and set it in pPolyOut without copy
614         SinglePoly** pSPTempOut = pITTempOut->getAs<Polynom>()->get();
615         pSPOut[i] = pSPTempOut[0];
616         // increase ref to avoid the delete of pSPTempOut[0]
617         // which are setted in pSPOut without copy.
618         pSPOut[i]->IncreaseRef();
619         delete pITTempOut;
620         pSPOut[i]->DecreaseRef();
621     }
622
623     //delete exp
624     for (int i = 0; i < iSize; i++)
625     {
626         delete pDblPower[i];
627     }
628
629     delete[] pDblPower;
630
631     // delete temporary polynom
632     // do not delete the last SinglePoly of _pPoly setted without copy in pPolyTemp
633     pSPTemp[0]->IncreaseRef();
634     delete pPolyTemp;
635     pSP[iSize - 1]->DecreaseRef();
636
637     switch (iResult)
638     {
639         case 1 :
640         {
641             delete pPolyOut;
642             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
643         }
644         case 2 :
645         {
646             delete pPolyOut;
647             throw ast::InternalError(_W("Invalid exponent.\n"));
648         }
649         default:
650             //OK
651             break;
652     }
653
654     *_pOut = pPolyOut;
655     return 0;
656 }
657
658 int DotPowerDoubleByDouble(Double* _pDouble1, Double* _pDouble2, Double** _pDoubleOut)
659 {
660     int iResultComplex = 0;
661
662     if (_pDouble1->isEmpty() || _pDouble2->isEmpty())
663     {
664         *_pDoubleOut = Double::Empty();
665     }
666     else if (_pDouble1->isIdentity())
667     {
668         return 1;
669     }
670     else if (_pDouble2->isIdentity())
671     {
672         *_pDoubleOut = dynamic_cast<Double*>(_pDouble1->clone());
673     }
674     else if (_pDouble1->isScalar())
675     {
676         //a .^ (b or B)
677         *_pDoubleOut = new Double(_pDouble2->getDims() , _pDouble2->getDimsArray(), true);
678
679         if (_pDouble1->isComplex())
680         {
681             double dblR1 = _pDouble1->get(0);
682             double dblI1 = _pDouble1->getImg(0);
683
684             if (_pDouble2->isComplex())
685             {
686                 iResultComplex = 1;
687                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
688                 {
689                     iPowerComplexScalarByComplexScalar(
690                         dblR1, dblI1,
691                         _pDouble2->get(i), _pDouble2->getImg(i),
692                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
693                 }
694             }
695             else
696             {
697                 iResultComplex = 1;
698                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
699                 {
700                     iPowerComplexScalarByRealScalar(
701                         dblR1, dblI1,
702                         _pDouble2->get(i),
703                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
704                 }
705             }
706         }
707         else
708         {
709             double dblR1 = _pDouble1->get(0);
710             if (_pDouble2->isComplex())
711             {
712                 iResultComplex = 1;
713                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
714                 {
715                     iPowerRealScalarByComplexScalar(
716                         dblR1,
717                         _pDouble2->get(i), _pDouble2->getImg(i),
718                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
719                 }
720             }
721             else
722             {
723                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
724                 {
725                     int iComplex = 1;
726                     iPowerRealScalarByRealScalar(
727                         dblR1,
728                         _pDouble2->get(i),
729                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
730                     iResultComplex |= iComplex;
731                 }
732             }
733         }
734     }
735     else if (_pDouble2->isScalar())
736     {
737         //A .^ b
738         *_pDoubleOut = new Double(_pDouble1->getDims() , _pDouble1->getDimsArray(), true);
739         if (_pDouble1->isComplex())
740         {
741             double dblR2 = _pDouble2->get(0);
742             if (_pDouble2->isComplex())
743             {
744                 double dblI2 = _pDouble2->getImg(0);
745                 iResultComplex = 1;
746                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
747                 {
748                     iPowerComplexScalarByComplexScalar(
749                         _pDouble1->get(i), _pDouble1->getImg(i),
750                         dblR2, dblI2,
751                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
752                 }
753             }
754             else
755             {
756                 double dblR2 = _pDouble2->get(0);
757                 iResultComplex = 1;
758                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
759                 {
760                     iPowerComplexScalarByRealScalar(
761                         _pDouble1->get(i), _pDouble1->getImg(i),
762                         dblR2,
763                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
764                 }
765             }
766         }
767         else
768         {
769             double dblR2 = _pDouble2->get(0);
770             if (_pDouble2->isComplex())
771             {
772                 double dblI2 = _pDouble2->getImg(0);
773                 iResultComplex = 1;
774                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
775                 {
776                     iPowerRealScalarByComplexScalar(
777                         _pDouble1->get(i),
778                         dblR2, dblI2,
779                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
780                 }
781             }
782             else
783             {
784                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
785                 {
786                     int iComplex = 1;
787                     iPowerRealScalarByRealScalar(
788                         _pDouble1->get(i),
789                         dblR2,
790                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
791                     iResultComplex |= iComplex;
792                 }
793             }
794         }
795     }
796     else
797     {
798         //A .^ B
799         //check dimension compatibilities ( same number of dimension and same size for each dimension
800         int iDims1      = _pDouble1->getDims();
801         int* piDims1    = _pDouble1->getDimsArray();
802         int iDims2      = _pDouble2->getDims();
803         int* piDims2    = _pDouble2->getDimsArray();
804
805         if (iDims1 != iDims2)
806         {
807             return 1;
808         }
809
810         for (int i = 0 ; i < iDims1 ; i++)
811         {
812             if (piDims1[i] != piDims2[i])
813             {
814                 return 1;
815             }
816         }
817
818         (*_pDoubleOut) = new Double(iDims2, piDims2, true);
819
820         if (_pDouble1->isComplex())
821         {
822             if (_pDouble2->isComplex())
823             {
824                 iResultComplex = 1;
825                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
826                 {
827                     iPowerComplexScalarByComplexScalar(
828                         _pDouble1->get(i), _pDouble1->getImg(i),
829                         _pDouble2->get(i), _pDouble2->getImg(i),
830                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
831                 }
832             }
833             else
834             {
835                 iResultComplex = 1;
836                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
837                 {
838                     iPowerComplexScalarByRealScalar(
839                         _pDouble1->get(i), _pDouble1->getImg(i),
840                         _pDouble2->get(i),
841                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
842                 }
843             }
844         }
845         else
846         {
847             if (_pDouble2->isComplex())
848             {
849                 iResultComplex = 1;
850                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
851                 {
852                     iPowerRealScalarByComplexScalar(
853                         _pDouble1->get(i),
854                         _pDouble2->get(i), _pDouble2->getImg(i),
855                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
856                 }
857             }
858             else
859             {
860                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
861                 {
862                     int iComplex = 1;
863                     iPowerRealScalarByRealScalar(
864                         _pDouble1->get(i),
865                         _pDouble2->get(i),
866                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
867                     iResultComplex |= iComplex;
868                 }
869             }
870         }
871     }
872
873     if (iResultComplex == 0)
874     {
875         (*_pDoubleOut)->setComplex(false);
876     }
877     return 0;
878 }
879