fix dotpower for the Sparse
[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         pOut->setComplex(bComplex1);
366
367         for (int i = 0 ; i < _pDouble->getSize() ; i++)
368         {
369             SinglePoly* pCoeffOut = pOut->get(i);
370
371             int iCurrentRank    = 0;
372             int iLoop           = (int)_pDouble->get(i);
373
374             //initialize Out to 1
375             pCoeffOut->set(0, 1);
376             //get a copy of p
377             Polynom* pP = _pPoly->clone()->getAs<Polynom>();
378             pP->setComplex(_pPoly->isComplex());
379
380             while (iLoop)
381             {
382                 if (iLoop % 2)
383                 {
384                     int iRank = pP->getMaxRank();
385                     if (bComplex1)
386                     {
387                         C2F(wpmul1)(pCoeffOut->get(), pCoeffOut->getImg(), &iCurrentRank, pP->getCoef()->get(), pP->getCoef()->getImg(), &iRank, pCoeffOut->get(), pCoeffOut->getImg());
388                     }
389                     else
390                     {
391                         C2F(dpmul1)(pCoeffOut->get(), &iCurrentRank, pP->getCoef()->get(), &iRank, pCoeffOut->get());
392                     }
393                     iCurrentRank += iRank;
394                 }
395
396                 iLoop /= 2;
397                 if (iLoop)
398                 {
399                     //p = p * p
400                     Polynom* pTemp = NULL;
401                     MultiplyPolyByPoly(pP, pP, &pTemp);
402                     delete pP;
403                     pP = pTemp;
404                 }
405             }
406         }
407         *_pOut = pOut;
408     }
409     return 0;
410 }
411
412 int DotPowerSpaseByDouble(Sparse* _pSp, Double* _pDouble, InternalType** _pOut)
413 {
414     if (_pDouble->isEmpty())
415     {
416         //sp .^ []
417         *_pOut = Double::Empty();
418         return 0;
419     }
420
421     size_t iSize = _pSp->nonZeros();
422     int* Col = new int[iSize];
423     int* Row = new int[iSize];
424     _pSp->getColPos(Col);
425     _pSp->getNbItemByRow(Row);
426     int* iPositVal = new int[iSize];
427
428     int j = 0;
429     for (int i = 0; i < iSize;  j++)
430     {
431         for (int k = 0; k < Row[j]; k++)
432         {
433
434             iPositVal[i] = (Col[i] - 1) * _pSp->getRows() + j;
435             i++;
436         }
437     }
438
439     Double** pDbl = new Double*[iSize];
440     Double** pDblSp = new Double*[iSize];
441     double* pdbl = _pDouble->get();
442
443     if (_pDouble->isScalar())
444     {
445         if (_pDouble->isComplex())
446         {
447             double* pdblImg = _pDouble->getImg();
448             for (int i = 0; i < iSize; i++)
449             {
450                 pDbl[i] = new Double(pdbl[0], pdblImg[0]);
451                 pDblSp[i] = new Double(_pSp->get(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
452             }
453         }
454         else
455         {
456             for (int i = 0; i < iSize; i++)
457             {
458                 pDbl[i] = new Double(pdbl[0]);
459                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
460             }
461         }
462     }
463     else if (_pDouble->getSize() == iSize)
464     {
465         if (_pDouble->isComplex())
466         {
467             double* pdblImg = _pDouble->getImg();
468             for (int i = 0; i < iSize; i++)
469             {
470                 pDbl[i] = new Double(pdbl[i], pdblImg[i]);
471                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
472             }
473         }
474         else
475         {
476             for (int i = 0; i < iSize; i++)
477             {
478                 pDbl[i] = new Double(pdbl[i]);
479                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
480             }
481         }
482     }
483     else
484     {
485         delete[] pDblSp;
486         throw ast::ScilabError(_W("Invalid exponent.\n"));
487         return 1;
488     }
489
490     Sparse* pSpTemp = new Sparse(_pSp->getRows(), _pSp->getCols(), _pSp->isComplex() || _pDouble->isComplex());
491     pSpTemp->zero_set();
492
493     Double* ppDblGet = NULL;
494     for (int i = 0; i < iSize; i++)
495     {
496         if ((pDblSp[i]->get(0) != 0) || (pDblSp[i]->getImg(0) != 0))
497         {
498             DotPowerDoubleByDouble(pDblSp[i], pDbl[i], &ppDblGet);
499             std::complex<double> cplx(ppDblGet->get(0), ppDblGet->getImg(0));
500             pSpTemp->set(iPositVal[i], cplx, true);
501         }
502     }
503
504     delete Col;
505     delete Row;
506     delete iPositVal;
507
508     *_pOut = pSpTemp;
509     return 0;
510
511 }
512
513
514 int DotPowerPolyByDouble(Polynom* _pPoly, Double* _pDouble, InternalType** _pOut)
515 {
516     if (_pDouble->isEmpty())
517     {
518         //p .^ []
519         *_pOut = Double::Empty();
520         return 0;
521     }
522
523     int iSize = _pPoly->getSize();
524     Double** pDblPower  = new Double*[iSize];
525     double* pdblPower   = _pDouble->get();
526     if (_pPoly->isScalar())
527     {
528         return PowerPolyByDouble(_pPoly, _pDouble, _pOut);
529     }
530     else if (_pDouble->isScalar())
531     {
532         if (pdblPower[0] < 0)
533         {
534             //call overload
535             _pOut = NULL;
536             delete[] pDblPower;
537             return 0;
538         }
539
540         for (int i = 0; i < iSize; i++)
541         {
542             pDblPower[i] = new Double(pdblPower[0]);
543         }
544     }
545     else if (_pDouble->getSize() == iSize)
546     {
547         for (int i = 0; i < iSize; i++)
548         {
549             if (pdblPower[i] < 0)
550             {
551                 //call overload
552                 _pOut = NULL;
553                 delete[] pDblPower;
554                 return 0;
555             }
556
557             pDblPower[i] = new Double(pdblPower[i]);
558         }
559     }
560     else
561     {
562         delete[] pDblPower;
563         throw ast::ScilabError(_W("Invalid exponent.\n"));
564     }
565
566     InternalType* pITTempOut    = NULL;
567     Polynom* pPolyTemp          = new Polynom(_pPoly->getVariableName(), 1, 1);
568     Polynom* pPolyOut           = new Polynom(_pPoly->getVariableName(), _pPoly->getDims(), _pPoly->getDimsArray());
569     SinglePoly** pSPOut         = pPolyOut->get();
570     SinglePoly** pSPTemp        = pPolyTemp->get();
571     SinglePoly** pSP            = _pPoly->get();
572
573     int iResult = 0;
574     for (int i = 0; i < iSize; i++)
575     {
576         // set singlePoly of _pPoly in pPolyTemp without copy
577         pSPTemp[0] = pSP[i];
578         iResult = PowerPolyByDouble(pPolyTemp, pDblPower[i], &pITTempOut);
579         if (iResult)
580         {
581             break;
582         }
583
584         // get singlePoly of pITTempOut and set it in pPolyOut without copy
585         SinglePoly** pSPTempOut = pITTempOut->getAs<Polynom>()->get();
586         pSPOut[i] = pSPTempOut[0];
587         pSPTempOut[0] = NULL;
588         delete pITTempOut;
589     }
590
591     // delete exp
592     for (int i = 0; i < iSize; i++)
593     {
594         delete pDblPower[i];
595     }
596
597     delete pDblPower;
598
599     // delete temporary polynom
600     // do not delete the last SinglePoly of _pPoly setted without copy in pPolyTemp
601     pSPTemp[0] = NULL;
602     delete pPolyTemp;
603
604     switch (iResult)
605     {
606         case 1 :
607         {
608             delete pPolyOut;
609             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
610         }
611         case 2 :
612         {
613             delete pPolyOut;
614             throw ast::ScilabError(_W("Invalid exponent.\n"));
615         }
616         default:
617             //OK
618             break;
619     }
620
621     *_pOut = pPolyOut;
622     return 0;
623 }
624
625 int DotPowerDoubleByDouble(Double* _pDouble1, Double* _pDouble2, Double** _pDoubleOut)
626 {
627     int iResultComplex = 0;
628
629     if (_pDouble1->isEmpty() || _pDouble2->isEmpty())
630     {
631         *_pDoubleOut = Double::Empty();
632     }
633     else if (_pDouble1->isIdentity())
634     {
635         return 1;
636     }
637     else if (_pDouble2->isIdentity())
638     {
639         *_pDoubleOut = dynamic_cast<Double*>(_pDouble1->clone());
640     }
641     else if (_pDouble1->isScalar())
642     {
643         //a .^ (b or B)
644         *_pDoubleOut = new Double(_pDouble2->getDims() , _pDouble2->getDimsArray(), true);
645
646         if (_pDouble1->isComplex())
647         {
648             double dblR1 = _pDouble1->get(0);
649             double dblI1 = _pDouble1->getImg(0);
650
651             if (_pDouble2->isComplex())
652             {
653                 iResultComplex = 1;
654                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
655                 {
656                     iPowerComplexScalarByComplexScalar(
657                         dblR1, dblI1,
658                         _pDouble2->get(i), _pDouble2->getImg(i),
659                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
660                 }
661             }
662             else
663             {
664                 iResultComplex = 1;
665                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
666                 {
667                     iPowerComplexScalarByRealScalar(
668                         dblR1, dblI1,
669                         _pDouble2->get(i),
670                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
671                 }
672             }
673         }
674         else
675         {
676             double dblR1 = _pDouble1->get(0);
677             if (_pDouble2->isComplex())
678             {
679                 iResultComplex = 1;
680                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
681                 {
682                     iPowerRealScalarByComplexScalar(
683                         dblR1,
684                         _pDouble2->get(i), _pDouble2->getImg(i),
685                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
686                 }
687             }
688             else
689             {
690                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
691                 {
692                     int iComplex = 1;
693                     iPowerRealScalarByRealScalar(
694                         dblR1,
695                         _pDouble2->get(i),
696                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
697                     iResultComplex |= iComplex;
698                 }
699             }
700         }
701     }
702     else if (_pDouble2->isScalar())
703     {
704         //A .^ b
705         *_pDoubleOut = new Double(_pDouble1->getDims() , _pDouble1->getDimsArray(), true);
706         if (_pDouble1->isComplex())
707         {
708             double dblR2 = _pDouble2->get(0);
709             if (_pDouble2->isComplex())
710             {
711                 double dblI2 = _pDouble2->getImg(0);
712                 iResultComplex = 1;
713                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
714                 {
715                     iPowerComplexScalarByComplexScalar(
716                         _pDouble1->get(i), _pDouble1->getImg(i),
717                         dblR2, dblI2,
718                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
719                 }
720             }
721             else
722             {
723                 double dblR2 = _pDouble2->get(0);
724                 iResultComplex = 1;
725                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
726                 {
727                     iPowerComplexScalarByRealScalar(
728                         _pDouble1->get(i), _pDouble1->getImg(i),
729                         dblR2,
730                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
731                 }
732             }
733         }
734         else
735         {
736             double dblR2 = _pDouble2->get(0);
737             if (_pDouble2->isComplex())
738             {
739                 double dblI2 = _pDouble2->getImg(0);
740                 iResultComplex = 1;
741                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
742                 {
743                     iPowerRealScalarByComplexScalar(
744                         _pDouble1->get(i),
745                         dblR2, dblI2,
746                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
747                 }
748             }
749             else
750             {
751                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
752                 {
753                     int iComplex = 1;
754                     iPowerRealScalarByRealScalar(
755                         _pDouble1->get(i),
756                         dblR2,
757                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
758                     iResultComplex |= iComplex;
759                 }
760             }
761         }
762     }
763     else
764     {
765         //A .^ B
766         //check dimension compatibilities ( same number of dimension and same size for each dimension
767         int iDims1      = _pDouble1->getDims();
768         int* piDims1    = _pDouble1->getDimsArray();
769         int iDims2      = _pDouble2->getDims();
770         int* piDims2    = _pDouble2->getDimsArray();
771
772         if (iDims1 != iDims2)
773         {
774             return 1;
775         }
776
777         for (int i = 0 ; i < iDims1 ; i++)
778         {
779             if (piDims1[i] != piDims2[i])
780             {
781                 return 1;
782             }
783         }
784
785         (*_pDoubleOut) = new Double(iDims2, piDims2, true);
786
787         if (_pDouble1->isComplex())
788         {
789             if (_pDouble2->isComplex())
790             {
791                 iResultComplex = 1;
792                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
793                 {
794                     iPowerComplexScalarByComplexScalar(
795                         _pDouble1->get(i), _pDouble1->getImg(i),
796                         _pDouble2->get(i), _pDouble2->getImg(i),
797                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
798                 }
799             }
800             else
801             {
802                 iResultComplex = 1;
803                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
804                 {
805                     iPowerComplexScalarByRealScalar(
806                         _pDouble1->get(i), _pDouble1->getImg(i),
807                         _pDouble2->get(i),
808                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
809                 }
810             }
811         }
812         else
813         {
814             if (_pDouble2->isComplex())
815             {
816                 iResultComplex = 1;
817                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
818                 {
819                     iPowerRealScalarByComplexScalar(
820                         _pDouble1->get(i),
821                         _pDouble2->get(i), _pDouble2->getImg(i),
822                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
823                 }
824             }
825             else
826             {
827                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
828                 {
829                     int iComplex = 1;
830                     iPowerRealScalarByRealScalar(
831                         _pDouble1->get(i),
832                         _pDouble2->get(i),
833                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
834                     iResultComplex |= iComplex;
835                 }
836             }
837         }
838     }
839
840     if (iResultComplex == 0)
841     {
842         (*_pDoubleOut)->setComplex(false);
843     }
844     return 0;
845 }
846