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