1bdc70bba5ffa6c1b585da58e25f5ef06d9d2f3f
[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         throw ast::InternalError(_W("Invalid exponent.\n"));
506         return 1;
507     }
508
509     Sparse* pSpTemp = new Sparse(_pSp->getRows(), _pSp->getCols(), _pSp->isComplex() || _pDouble->isComplex());
510     pSpTemp->zero_set();
511
512     Double* ppDblGet = NULL;
513     for (int i = 0; i < iSize; i++)
514     {
515         if ((pDblSp[i]->get(0) != 0) || (pDblSp[i]->getImg(0) != 0))
516         {
517             DotPowerDoubleByDouble(pDblSp[i], pDbl[i], &ppDblGet);
518             std::complex<double> cplx(ppDblGet->get(0), ppDblGet->getImg(0));
519             pSpTemp->set(iPositVal[i], cplx, false);
520         }
521     }
522
523     delete[] Col;
524     delete[] Row;
525     delete[] iPositVal;
526
527     pSpTemp->finalize();
528     *_pOut = pSpTemp;
529     return 0;
530
531 }
532
533
534 int DotPowerPolyByDouble(Polynom* _pPoly, Double* _pDouble, InternalType** _pOut)
535 {
536     if (_pDouble->isEmpty())
537     {
538         //p .^ []
539         *_pOut = Double::Empty();
540         return 0;
541     }
542
543     int iSize = _pPoly->getSize();
544     if (_pPoly->isScalar())
545     {
546         return PowerPolyByDouble(_pPoly, _pDouble, _pOut);
547     }
548
549     Double** pDblPower = new Double*[iSize];
550     double* pdblPower = _pDouble->get();
551     if (_pDouble->isScalar())
552     {
553         if (pdblPower[0] < 0)
554         {
555             //call overload
556             _pOut = NULL;
557             delete[] pDblPower;
558             return 0;
559         }
560
561         for (int i = 0; i < iSize; i++)
562         {
563             pDblPower[i] = new Double(pdblPower[0]);
564         }
565     }
566     else if (_pDouble->getSize() == iSize)
567     {
568         for (int i = 0; i < iSize; i++)
569         {
570             if (pdblPower[i] < 0)
571             {
572                 //call overload
573                 _pOut = NULL;
574                 delete[] pDblPower;
575                 return 0;
576             }
577
578             pDblPower[i] = new Double(pdblPower[i]);
579         }
580     }
581     else
582     {
583         delete[] pDblPower;
584         throw ast::InternalError(_W("Invalid exponent.\n"));
585     }
586
587     InternalType* pITTempOut    = NULL;
588     Polynom* pPolyTemp          = new Polynom(_pPoly->getVariableName(), 1, 1);
589     Polynom* pPolyOut           = new Polynom(_pPoly->getVariableName(), _pPoly->getDims(), _pPoly->getDimsArray());
590     SinglePoly** pSPOut         = pPolyOut->get();
591     SinglePoly** pSPTemp        = pPolyTemp->get();
592     SinglePoly** pSP            = _pPoly->get();
593
594     int iResult = 0;
595     for (int i = 0; i < iSize; i++)
596     {
597         // set singlePoly of _pPoly in pPolyTemp without copy
598         pSPTemp[0] = pSP[i];
599         iResult = PowerPolyByDouble(pPolyTemp, pDblPower[i], &pITTempOut);
600         if (iResult)
601         {
602             break;
603         }
604
605         // get singlePoly of pITTempOut and set it in pPolyOut without copy
606         SinglePoly** pSPTempOut = pITTempOut->getAs<Polynom>()->get();
607         pSPOut[i] = pSPTempOut[0];
608         // increase ref to avoid the delete of pSPTempOut[0]
609         // which are setted in pSPOut without copy.
610         pSPOut[i]->IncreaseRef();
611         delete pITTempOut;
612         pSPOut[i]->DecreaseRef();
613     }
614
615     //delete exp
616     for (int i = 0; i < iSize; i++)
617     {
618         delete pDblPower[i];
619     }
620
621     delete[] pDblPower;
622
623     // delete temporary polynom
624     // do not delete the last SinglePoly of _pPoly setted without copy in pPolyTemp
625     pSPTemp[0]->IncreaseRef();
626     delete pPolyTemp;
627     pSP[iSize - 1]->DecreaseRef();
628
629     switch (iResult)
630     {
631         case 1 :
632         {
633             delete pPolyOut;
634             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
635         }
636         case 2 :
637         {
638             delete pPolyOut;
639             throw ast::InternalError(_W("Invalid exponent.\n"));
640         }
641         default:
642             //OK
643             break;
644     }
645
646     *_pOut = pPolyOut;
647     return 0;
648 }
649
650 int DotPowerDoubleByDouble(Double* _pDouble1, Double* _pDouble2, Double** _pDoubleOut)
651 {
652     int iResultComplex = 0;
653
654     if (_pDouble1->isEmpty() || _pDouble2->isEmpty())
655     {
656         *_pDoubleOut = Double::Empty();
657     }
658     else if (_pDouble1->isIdentity())
659     {
660         return 1;
661     }
662     else if (_pDouble2->isIdentity())
663     {
664         *_pDoubleOut = dynamic_cast<Double*>(_pDouble1->clone());
665     }
666     else if (_pDouble1->isScalar())
667     {
668         //a .^ (b or B)
669         *_pDoubleOut = new Double(_pDouble2->getDims() , _pDouble2->getDimsArray(), true);
670
671         if (_pDouble1->isComplex())
672         {
673             double dblR1 = _pDouble1->get(0);
674             double dblI1 = _pDouble1->getImg(0);
675
676             if (_pDouble2->isComplex())
677             {
678                 iResultComplex = 1;
679                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
680                 {
681                     iPowerComplexScalarByComplexScalar(
682                         dblR1, dblI1,
683                         _pDouble2->get(i), _pDouble2->getImg(i),
684                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
685                 }
686             }
687             else
688             {
689                 iResultComplex = 1;
690                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
691                 {
692                     iPowerComplexScalarByRealScalar(
693                         dblR1, dblI1,
694                         _pDouble2->get(i),
695                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
696                 }
697             }
698         }
699         else
700         {
701             double dblR1 = _pDouble1->get(0);
702             if (_pDouble2->isComplex())
703             {
704                 iResultComplex = 1;
705                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
706                 {
707                     iPowerRealScalarByComplexScalar(
708                         dblR1,
709                         _pDouble2->get(i), _pDouble2->getImg(i),
710                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
711                 }
712             }
713             else
714             {
715                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
716                 {
717                     int iComplex = 1;
718                     iPowerRealScalarByRealScalar(
719                         dblR1,
720                         _pDouble2->get(i),
721                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
722                     iResultComplex |= iComplex;
723                 }
724             }
725         }
726     }
727     else if (_pDouble2->isScalar())
728     {
729         //A .^ b
730         *_pDoubleOut = new Double(_pDouble1->getDims() , _pDouble1->getDimsArray(), true);
731         if (_pDouble1->isComplex())
732         {
733             double dblR2 = _pDouble2->get(0);
734             if (_pDouble2->isComplex())
735             {
736                 double dblI2 = _pDouble2->getImg(0);
737                 iResultComplex = 1;
738                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
739                 {
740                     iPowerComplexScalarByComplexScalar(
741                         _pDouble1->get(i), _pDouble1->getImg(i),
742                         dblR2, dblI2,
743                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
744                 }
745             }
746             else
747             {
748                 double dblR2 = _pDouble2->get(0);
749                 iResultComplex = 1;
750                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
751                 {
752                     iPowerComplexScalarByRealScalar(
753                         _pDouble1->get(i), _pDouble1->getImg(i),
754                         dblR2,
755                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
756                 }
757             }
758         }
759         else
760         {
761             double dblR2 = _pDouble2->get(0);
762             if (_pDouble2->isComplex())
763             {
764                 double dblI2 = _pDouble2->getImg(0);
765                 iResultComplex = 1;
766                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
767                 {
768                     iPowerRealScalarByComplexScalar(
769                         _pDouble1->get(i),
770                         dblR2, dblI2,
771                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
772                 }
773             }
774             else
775             {
776                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
777                 {
778                     int iComplex = 1;
779                     iPowerRealScalarByRealScalar(
780                         _pDouble1->get(i),
781                         dblR2,
782                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
783                     iResultComplex |= iComplex;
784                 }
785             }
786         }
787     }
788     else
789     {
790         //A .^ B
791         //check dimension compatibilities ( same number of dimension and same size for each dimension
792         int iDims1      = _pDouble1->getDims();
793         int* piDims1    = _pDouble1->getDimsArray();
794         int iDims2      = _pDouble2->getDims();
795         int* piDims2    = _pDouble2->getDimsArray();
796
797         if (iDims1 != iDims2)
798         {
799             return 1;
800         }
801
802         for (int i = 0 ; i < iDims1 ; i++)
803         {
804             if (piDims1[i] != piDims2[i])
805             {
806                 return 1;
807             }
808         }
809
810         (*_pDoubleOut) = new Double(iDims2, piDims2, true);
811
812         if (_pDouble1->isComplex())
813         {
814             if (_pDouble2->isComplex())
815             {
816                 iResultComplex = 1;
817                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
818                 {
819                     iPowerComplexScalarByComplexScalar(
820                         _pDouble1->get(i), _pDouble1->getImg(i),
821                         _pDouble2->get(i), _pDouble2->getImg(i),
822                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
823                 }
824             }
825             else
826             {
827                 iResultComplex = 1;
828                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
829                 {
830                     iPowerComplexScalarByRealScalar(
831                         _pDouble1->get(i), _pDouble1->getImg(i),
832                         _pDouble2->get(i),
833                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
834                 }
835             }
836         }
837         else
838         {
839             if (_pDouble2->isComplex())
840             {
841                 iResultComplex = 1;
842                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
843                 {
844                     iPowerRealScalarByComplexScalar(
845                         _pDouble1->get(i),
846                         _pDouble2->get(i), _pDouble2->getImg(i),
847                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
848                 }
849             }
850             else
851             {
852                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
853                 {
854                     int iComplex = 1;
855                     iPowerRealScalarByRealScalar(
856                         _pDouble1->get(i),
857                         _pDouble2->get(i),
858                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
859                     iResultComplex |= iComplex;
860                 }
861             }
862         }
863     }
864
865     if (iResultComplex == 0)
866     {
867         (*_pDoubleOut)->setComplex(false);
868     }
869     return 0;
870 }
871