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