5f689618066436ad1879a8622d3cca0c5790a189
[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     GenericType::ScilabType TypeL = _pLeftOperand->getType();
29     GenericType::ScilabType TypeR = _pRightOperand->getType();
30
31     /*
32     ** DOUBLE ^ DOUBLE
33     ** DOUBLE ** DOUBLE
34     */
35     if (_pLeftOperand->isDouble() && _pRightOperand->isDouble())
36     {
37         Double *pL   = _pLeftOperand->getAs<Double>();
38         Double *pR   = _pRightOperand->getAs<Double>();
39
40         int iResult = PowerDoubleByDouble(pL, pR, (Double**)&pResult);
41         if (iResult)
42         {
43             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
44         }
45
46         return pResult;
47     }
48
49     /*
50     ** POLY ^ DOUBLE
51     ** POLY ** DOUBLE
52     */
53     if (_pLeftOperand->isPoly() && _pRightOperand->isDouble())
54     {
55         Polynom *pL   = _pLeftOperand->getAs<Polynom>();
56         Double *pR    = _pRightOperand->getAs<Double>();
57
58         int iResult = PowerPolyByDouble(pL, pR, &pResult);
59         switch (iResult)
60         {
61             case 1 :
62                 throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
63             case 2 :
64                 throw ast::ScilabError(_W("Invalid exponent.\n"));
65             default:
66                 //OK
67                 break;
68         }
69
70         return pResult;
71     }
72
73     /*
74     ** Default case : Return NULL will Call Overloading.
75     */
76     return NULL;
77
78 }
79
80 InternalType *GenericDotPower(InternalType *_pLeftOperand, InternalType *_pRightOperand)
81 {
82     InternalType *pResult = NULL;
83     GenericType::ScilabType TypeL = _pLeftOperand->getType();
84     GenericType::ScilabType TypeR = _pRightOperand->getType();
85
86     /*
87     ** DOUBLE .^ DOUBLE
88     ** DOUBLE .** DOUBLE
89     */
90     if (TypeL == GenericType::ScilabDouble && TypeR == GenericType::ScilabDouble)
91     {
92         Double *pL  = _pLeftOperand->getAs<Double>();
93         Double *pR  = _pRightOperand->getAs<Double>();
94
95         int iResult = DotPowerDoubleByDouble(pL, pR, (Double**)&pResult);
96         if (iResult)
97         {
98             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
99         }
100
101         return pResult;
102     }
103
104     /*
105     ** POLY .^ DOUBLE
106     ** POLY .** DOUBLE
107     */
108     if (TypeL == GenericType::ScilabPolynom && TypeR == GenericType::ScilabDouble)
109     {
110         Polynom *pL   = _pLeftOperand->getAs<Polynom>();
111         Double *pR   = _pRightOperand->getAs<Double>();
112
113         int iResult = PowerPolyByDouble(pL, pR, &pResult);
114         switch (iResult)
115         {
116             case 1 :
117                 throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
118             case 2 :
119                 throw ast::ScilabError(_W("Invalid exponent.\n"));
120             default:
121                 //OK
122                 break;
123         }
124         return pResult;
125     }
126
127     /*
128     ** Default case : Return NULL will Call Overloading.
129     */
130     return NULL;
131
132 }
133
134 int PowerDoubleByDouble(Double* _pDouble1, Double* _pDouble2, Double** _pDoubleOut)
135 {
136     bool bComplex1  = _pDouble1->isComplex();
137     bool bComplex2  = _pDouble2->isComplex();
138     bool bScalar1   = _pDouble1->isScalar();
139     bool bScalar2   = _pDouble2->isScalar();
140
141     int iComplex = 1;
142
143     if (bScalar1 && bScalar2)
144     {
145         //s ^ s
146         *_pDoubleOut = new Double(1, 1, true);
147
148         if (bComplex1 == false && bComplex2 == false)
149         {
150             iPowerRealScalarByRealScalar(_pDouble1->get(0), _pDouble2->get(0), (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), &iComplex);
151         }
152         else if (bComplex1 == false && bComplex2 == true)
153         {
154             iPowerRealScalarByComplexScalar(_pDouble1->get(0), _pDouble2->get(0), _pDouble2->getImg(0), (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
155         }
156         else if (bComplex1 == true && bComplex2 == false)
157         {
158             iPowerComplexScalarByRealScalar(_pDouble1->get(0), _pDouble1->getImg(0), _pDouble2->get(0), (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
159         }
160         else if (bComplex1 == true && bComplex2 == true)
161         {
162             iPowerComplexScalarByComplexScalar(_pDouble1->get(0), _pDouble1->getImg(0), _pDouble2->get(0), _pDouble2->getImg(0), (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
163         }
164
165         if (iComplex == 0)
166         {
167             (*_pDoubleOut)->setComplex(false);
168         }
169
170         return 0;
171     }
172     else if (bScalar1 && _pDouble2->getDims() == 2)
173     {
174         //s ^ []
175         *_pDoubleOut = new Double(_pDouble2->getRows(), _pDouble2->getCols(), true);
176
177         if (bComplex1 == false && bComplex2 == false)
178         {
179             iPowerRealScalarByRealMatrix(
180                 _pDouble1->get(0),
181                 _pDouble2->get(), _pDouble2->getRows(), _pDouble2->getCols(),
182                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), &iComplex);
183         }
184         else if (bComplex1 == false && bComplex2 == true)
185         {
186             iPowerRealScalarByComplexMatrix(
187                 _pDouble1->get(0),
188                 _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
189                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
190         }
191         else if (bComplex1 == true && bComplex2 == false)
192         {
193             iPowerComplexScalarByRealMatrix(
194                 _pDouble1->get(0), _pDouble1->getImg(0),
195                 _pDouble2->get(), _pDouble2->getRows(), _pDouble2->getCols(),
196                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
197         }
198         else if (bComplex1 == true && bComplex2 == true)
199         {
200             iPowerComplexScalarByComplexMatrix(
201                 _pDouble1->get(0), _pDouble1->getImg(0),
202                 _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
203                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
204         }
205
206         if (iComplex == 0)
207         {
208             (*_pDoubleOut)->setComplex(false);
209         }
210
211         return 0;
212     }
213
214     if (bScalar2 && _pDouble1->getDims() == 2 && _pDouble1->isVector() )
215     {
216         //_pDouble1 is a vector and _pDouble is a scalar
217         *_pDoubleOut = new Double(_pDouble1->getRows(), _pDouble1->getCols() , true);
218
219         if (bComplex1 == false && bComplex2 == false)
220         {
221             for (int i = 0 ; i < (*_pDoubleOut)->getSize() ; i++)
222             {
223                 iPowerRealScalarByRealScalar(
224                     _pDouble1->get(i),
225                     _pDouble2->get(0),
226                     &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
227             }
228         }
229         else if (bComplex1 == false && bComplex2 == true)
230         {
231             for (int i = 0 ; i < (*_pDoubleOut)->getSize() ; i++)
232             {
233                 iPowerRealScalarByComplexScalar(
234                     _pDouble1->get(i),
235                     _pDouble2->get(0), _pDouble2->getImg(0),
236                     &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
237             }
238         }
239         else if (bComplex1 == true && bComplex2 == false)
240         {
241             for (int i = 0 ; i < (*_pDoubleOut)->getSize() ; i++)
242             {
243                 iPowerComplexScalarByRealScalar(
244                     _pDouble1->get(i), _pDouble1->getImg(i),
245                     _pDouble2->get(0),
246                     &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
247             }
248         }
249         else if (bComplex1 == true && bComplex2 == true)
250         {
251             for (int i = 0 ; i < (*_pDoubleOut)->getSize() ; i++)
252             {
253                 iPowerComplexScalarByComplexScalar(
254                     _pDouble1->get(i), _pDouble1->getImg(i),
255                     _pDouble2->get(0), _pDouble2->getImg(0),
256                     &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
257             }
258         }
259
260         if (iComplex == 0)
261         {
262             (*_pDoubleOut)->setComplex(false);
263         }
264
265         return 0;
266     }
267
268     if (bScalar2 && ( _pDouble1->getRows() == _pDouble1->getCols()))
269     {
270         //power of a square matrix by a scalar exponent.
271         int iRet = 0;
272         if (bComplex2)
273         {
274             //mange by overloading
275             return 0;
276         }
277
278         *_pDoubleOut = new Double(_pDouble1->getRows(), _pDouble1->getCols() , true);
279         if (bComplex1 == false)
280         {
281             iRet = iPowerRealSquareMatrixByRealScalar(
282                        _pDouble1->get(), _pDouble1->getRows(), _pDouble1->getCols(),
283                        _pDouble2->get(0),
284                        (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), &iComplex);
285         }
286         else if (bComplex1 == true)
287         {
288             iRet = iPowerComplexSquareMatrixByRealScalar(
289                        _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
290                        _pDouble2->get(0),
291                        (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
292         }
293
294         // call overload
295         if (iRet == -1)
296         {
297             delete *_pDoubleOut;
298             *_pDoubleOut = NULL;
299             return 0;
300         }
301
302         if (iComplex == 0)
303         {
304             (*_pDoubleOut)->setComplex(false);
305         }
306     }
307     return 0;
308 }
309
310 int PowerPolyByDouble(Polynom* _pPoly, Double* _pDouble, InternalType** _pOut)
311 {
312     bool bComplex1  = _pPoly->isComplex();
313     bool bComplex2  = _pDouble->isComplex();
314     bool bScalar1   = _pPoly->isScalar();
315     bool bScalar2   = _pDouble->isScalar();
316
317     int iComplex = 1;
318
319     if (bComplex2)
320     {
321         //invalid exponent.
322         return 2;
323     }
324
325     if (_pDouble->isEmpty())
326     {
327         //p ** []
328         *_pOut = Double::Empty();
329         return 0;
330     }
331
332     if (bScalar1)
333     {
334         //p ^ x or p ^ X
335         int iRank   = 0;
336         int* piRank = new int[_pDouble->getSize()];
337
338         _pPoly->getRank(&iRank);
339         for (int i = 0 ; i < _pDouble->getSize() ; i++)
340         {
341             int iInputRank = (int)_pDouble->get(i);
342             if (iInputRank < 0)
343             {
344                 //call overload
345                 _pOut = NULL;
346                 delete[] piRank;
347                 return 0;
348             }
349
350             piRank[i] = ((iRank - 1) * iInputRank) + 1;
351         }
352
353         Polynom* pOut = new Polynom(_pPoly->getVariableName(), _pDouble->getRows(), _pDouble->getCols(), piRank);
354         pOut->setComplex(bComplex1);
355
356         for (int i = 0 ; i < _pDouble->getSize() ; i++)
357         {
358             Double* pCoeffOut = pOut->get(i)->getCoef();
359
360             int iCurrentRank    = 0;
361             int iLoop           = (int)_pDouble->get(i);
362
363             //initialize Out to 1
364             pCoeffOut->set(0, 1);
365             //get a copy of p
366             Polynom* pP = _pPoly->clone()->getAs<Polynom>();
367             pP->setComplex(_pPoly->isComplex());
368
369             while (iLoop)
370             {
371                 if (iLoop % 2)
372                 {
373                     int iRank = pP->getMaxRank() - 1;
374                     if (bComplex1)
375                     {
376                         C2F(wpmul1)(pCoeffOut->get(), pCoeffOut->getImg(), &iCurrentRank, pP->getCoef()->get(), pP->getCoef()->getImg(), &iRank, pCoeffOut->get(), pCoeffOut->getImg());
377                     }
378                     else
379                     {
380                         C2F(dpmul1)(pCoeffOut->get(), &iCurrentRank, pP->getCoef()->get(), &iRank, pCoeffOut->get());
381                     }
382                     iCurrentRank += iRank;
383                 }
384
385                 iLoop /= 2;
386                 if (iLoop)
387                 {
388                     //p = p * p
389                     Polynom* pTemp = NULL;
390                     MultiplyPolyByPoly(pP, pP, &pTemp);
391                     delete pP;
392                     pP = pTemp;
393                 }
394             }
395         }
396         *_pOut = pOut;
397     }
398     return 0;
399 }
400
401 int DotPowerDoubleByDouble(Double* _pDouble1, Double* _pDouble2, Double** _pDoubleOut)
402 {
403     int iResultComplex = 0;
404
405     if (_pDouble1->isEmpty() || _pDouble2->isEmpty())
406     {
407         *_pDoubleOut = Double::Empty();
408     }
409     else if (_pDouble1->isIdentity())
410     {
411         return 1;
412     }
413     else if (_pDouble2->isIdentity())
414     {
415         *_pDoubleOut = dynamic_cast<Double*>(_pDouble1->clone());
416     }
417     else if (_pDouble1->isScalar())
418     {
419         //a .^ (b or B)
420         *_pDoubleOut = new Double(_pDouble2->getDims() , _pDouble2->getDimsArray(), true);
421
422         if (_pDouble1->isComplex())
423         {
424             double dblR1 = _pDouble1->get(0);
425             double dblI1 = _pDouble1->getImg(0);
426
427             if (_pDouble2->isComplex())
428             {
429                 iResultComplex = 1;
430                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
431                 {
432                     iPowerComplexScalarByComplexScalar(
433                         dblR1, dblI1,
434                         _pDouble2->get(i), _pDouble2->getImg(i),
435                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
436                 }
437             }
438             else
439             {
440                 iResultComplex = 1;
441                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
442                 {
443                     iPowerComplexScalarByRealScalar(
444                         dblR1, dblI1,
445                         _pDouble2->get(i),
446                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
447                 }
448             }
449         }
450         else
451         {
452             double dblR1 = _pDouble1->get(0);
453             if (_pDouble2->isComplex())
454             {
455                 iResultComplex = 1;
456                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
457                 {
458                     iPowerRealScalarByComplexScalar(
459                         dblR1,
460                         _pDouble2->get(i), _pDouble2->getImg(i),
461                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
462                 }
463             }
464             else
465             {
466                 for (int i = 0 ; i < _pDouble2->getSize() ; i++)
467                 {
468                     int iComplex = 1;
469                     iPowerRealScalarByRealScalar(
470                         dblR1,
471                         _pDouble2->get(i),
472                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
473                     iResultComplex |= iComplex;
474                 }
475             }
476         }
477     }
478     else if (_pDouble2->isScalar())
479     {
480         //A .^ b
481         *_pDoubleOut = new Double(_pDouble1->getDims() , _pDouble1->getDimsArray(), true);
482         if (_pDouble1->isComplex())
483         {
484             double dblR2 = _pDouble2->get(0);
485             if (_pDouble2->isComplex())
486             {
487                 double dblI2 = _pDouble2->getImg(0);
488                 iResultComplex = 1;
489                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
490                 {
491                     iPowerComplexScalarByComplexScalar(
492                         _pDouble1->get(i), _pDouble1->getImg(i),
493                         dblR2, dblI2,
494                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
495                 }
496             }
497             else
498             {
499                 double dblR2 = _pDouble2->get(0);
500                 iResultComplex = 1;
501                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
502                 {
503                     iPowerComplexScalarByRealScalar(
504                         _pDouble1->get(i), _pDouble1->getImg(i),
505                         dblR2,
506                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
507                 }
508             }
509         }
510         else
511         {
512             double dblR2 = _pDouble2->get(0);
513             if (_pDouble2->isComplex())
514             {
515                 double dblI2 = _pDouble2->getImg(0);
516                 iResultComplex = 1;
517                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
518                 {
519                     iPowerRealScalarByComplexScalar(
520                         _pDouble1->get(i),
521                         dblR2, dblI2,
522                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
523                 }
524             }
525             else
526             {
527                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
528                 {
529                     int iComplex = 1;
530                     iPowerRealScalarByRealScalar(
531                         _pDouble1->get(i),
532                         dblR2,
533                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
534                     iResultComplex |= iComplex;
535                 }
536             }
537         }
538     }
539     else
540     {
541         //A .^ B
542         //check dimension compatibilities ( same number of dimension and same size for each dimension
543         int iDims1      = _pDouble1->getDims();
544         int* piDims1    = _pDouble1->getDimsArray();
545         int iDims2      = _pDouble2->getDims();
546         int* piDims2    = _pDouble2->getDimsArray();
547
548         if (iDims1 != iDims2)
549         {
550             return 1;
551         }
552
553         for (int i = 0 ; i < iDims1 ; i++)
554         {
555             if (piDims1[i] != piDims2[i])
556             {
557                 return 1;
558             }
559         }
560
561         (*_pDoubleOut) = new Double(iDims2, piDims2, true);
562
563         if (_pDouble1->isComplex())
564         {
565             if (_pDouble2->isComplex())
566             {
567                 iResultComplex = 1;
568                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
569                 {
570                     iPowerComplexScalarByComplexScalar(
571                         _pDouble1->get(i), _pDouble1->getImg(i),
572                         _pDouble2->get(i), _pDouble2->getImg(i),
573                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
574                 }
575             }
576             else
577             {
578                 iResultComplex = 1;
579                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
580                 {
581                     iPowerComplexScalarByRealScalar(
582                         _pDouble1->get(i), _pDouble1->getImg(i),
583                         _pDouble2->get(i),
584                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
585                 }
586             }
587         }
588         else
589         {
590             if (_pDouble2->isComplex())
591             {
592                 iResultComplex = 1;
593                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
594                 {
595                     iPowerRealScalarByComplexScalar(
596                         _pDouble1->get(i),
597                         _pDouble2->get(i), _pDouble2->getImg(i),
598                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i]);
599                 }
600             }
601             else
602             {
603                 for (int i = 0 ; i < _pDouble1->getSize() ; i++)
604                 {
605                     int iComplex = 1;
606                     iPowerRealScalarByRealScalar(
607                         _pDouble1->get(i),
608                         _pDouble2->get(i),
609                         &(*_pDoubleOut)->get()[i], &(*_pDoubleOut)->getImg()[i], &iComplex);
610                     iResultComplex |= iComplex;
611                 }
612             }
613         }
614     }
615
616     if (iResultComplex == 0)
617     {
618         (*_pDoubleOut)->setComplex(false);
619     }
620     return 0;
621 }
622