fix warning in ast module.
[scilab.git] / scilab / modules / ast / src / cpp / operations / types_multiplication.cpp
1 /*
2  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  *  Copyright (C) 2008-2008 - DIGITEO - Antoine ELIAS
4  *  Copyright (C) 2010-2010 - DIGITEO - Bruno JOFRET
5  *
6  *  This file must be used under the terms of the CeCILL.
7  *  This source file is licensed as described in the file COPYING, which
8  *  you should have received as part of this distribution.  The terms
9  *  are also available at
10  *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11  *
12  */
13
14 #include "types_multiplication.hxx"
15 #include "types_addition.hxx"
16 #include "double.hxx"
17 #include "int.hxx"
18 #include "sparse.hxx"
19 #include "polynom.hxx"
20 #include "singlepoly.hxx"
21
22 #include "scilabexception.hxx"
23
24 extern "C"
25 {
26 #include "matrix_multiplication.h"
27 #include "matrix_addition.h"
28 #include "operation_f.h"
29 #include "localization.h"
30 #include "charEncoding.h"
31 #include "elem_common.h"
32 }
33
34
35 using namespace types;
36
37 InternalType *GenericTimes(InternalType *_pLeftOperand, InternalType *_pRightOperand)
38 {
39     InternalType *pResult = NULL;
40     GenericType::ScilabType TypeL = _pLeftOperand->getType();
41     GenericType::ScilabType TypeR = _pRightOperand->getType();
42
43     if (TypeL == GenericType::ScilabDouble && _pLeftOperand->getAs<Double>()->isEmpty())
44     {
45         return Double::Empty();
46     }
47
48     if (TypeR == GenericType::ScilabDouble && _pRightOperand->getAs<Double>()->isEmpty())
49     {
50         return Double::Empty();
51     }
52
53     /*
54     ** DOUBLE * DOUBLE
55     */
56     if (TypeL == GenericType::ScilabDouble && TypeR == GenericType::ScilabDouble)
57     {
58         Double *pL   = _pLeftOperand->getAs<Double>();
59         Double *pR   = _pRightOperand->getAs<Double>();
60
61         int iResult = MultiplyDoubleByDouble(pL, pR, (Double**)&pResult);
62         if (iResult)
63         {
64             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
65         }
66
67         return pResult;
68     }
69
70     /*
71     ** DOUBLE * POLY
72     */
73     else if (TypeL == InternalType::ScilabDouble && TypeR == InternalType::ScilabPolynom)
74     {
75         Double *pL   = _pLeftOperand->getAs<Double>();
76         Polynom *pR     = _pRightOperand->getAs<types::Polynom>();
77
78         int iResult = MultiplyDoubleByPoly(pL, pR, (Polynom**)&pResult);
79         if (iResult)
80         {
81             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
82         }
83
84         return pResult;
85     }
86
87     /*
88     ** POLY * DOUBLE
89     */
90     else if (TypeL == InternalType::ScilabPolynom && TypeR == InternalType::ScilabDouble)
91     {
92         Polynom *pL          = _pLeftOperand->getAs<types::Polynom>();
93         Double *pR              = _pRightOperand->getAs<Double>();
94
95         int iResult = MultiplyPolyByDouble(pL, pR, (Polynom**)&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 * POLY
106     */
107     else if (TypeL == InternalType::ScilabPolynom && TypeR == InternalType::ScilabPolynom)
108     {
109         Polynom *pL          = _pLeftOperand->getAs<types::Polynom>();
110         Polynom *pR          = _pRightOperand->getAs<types::Polynom>();
111
112         int iResult = MultiplyPolyByPoly(pL, pR, (Polynom**)&pResult);
113         if (iResult)
114         {
115             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
116         }
117
118         return pResult;
119     }
120
121     /*
122     ** SPARSE * SPARSE
123     */
124     if (TypeL == GenericType::ScilabSparse && TypeR == GenericType::ScilabSparse)
125     {
126         Sparse *pL   = _pLeftOperand->getAs<Sparse>();
127         Sparse *pR   = _pRightOperand->getAs<Sparse>();
128
129         int iResult = MultiplySparseBySparse(pL, pR, (Sparse**)&pResult);
130         if (iResult)
131         {
132             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
133         }
134
135         return pResult;
136     }
137
138     /*
139     ** DOUBLE * SPARSE
140     */
141     if (TypeL == GenericType::ScilabDouble && TypeR == GenericType::ScilabSparse)
142     {
143         Double *pL   = _pLeftOperand->getAs<Double>();
144         Sparse *pR   = _pRightOperand->getAs<Sparse>();
145
146         int iResult = MultiplyDoubleBySparse(pL, pR, (GenericType**)&pResult);
147         if (iResult)
148         {
149             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
150         }
151
152         return pResult;
153     }
154
155     /*
156     ** SPARSE * DOUBLE
157     */
158     if (TypeL == GenericType::ScilabSparse && TypeR == GenericType::ScilabDouble)
159     {
160         Sparse *pL   = _pLeftOperand->getAs<Sparse>();
161         Double *pR   = _pRightOperand->getAs<Double>();
162
163         int iResult = MultiplySparseByDouble(pL, pR, (GenericType**)&pResult);
164         if (iResult)
165         {
166             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
167         }
168
169         return pResult;
170     }
171
172     /*
173     ** Default case : Return NULL will Call Overloading.
174     */
175     return NULL;
176
177 }
178
179 int MultiplyDoubleByDouble(Double* _pDouble1, Double* _pDouble2, Double** _pDoubleOut)
180 {
181     if (_pDouble1->isScalar())
182     {
183         bool bComplex1  = _pDouble1->isComplex();
184         bool bComplex2  = _pDouble2->isComplex();
185
186         (*_pDoubleOut) = new Double(_pDouble2->getDims(), _pDouble2->getDimsArray(), bComplex1 | bComplex2);
187
188         if (bComplex1 == false && bComplex2 == false)
189         {
190             iMultiRealScalarByRealMatrix(_pDouble1->get(0), _pDouble2->get(), _pDouble2->getSize(), 1, (*_pDoubleOut)->get());
191         }
192         else if (bComplex1 == false && bComplex2 == true)
193         {
194             iMultiRealScalarByComplexMatrix(_pDouble1->get(0), _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getSize(), 1, (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
195         }
196         else if (bComplex1 == true && bComplex2 == false)
197         {
198             iMultiComplexScalarByRealMatrix(_pDouble1->get(0), _pDouble1->getImg(0), _pDouble2->get(), _pDouble2->getSize(), 1, (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
199         }
200         else //if(bComplex1 == true && bComplex2 == true)
201         {
202             iMultiComplexScalarByComplexMatrix(_pDouble1->get(0), _pDouble1->getImg(0), _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getSize(), 1, (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
203         }
204
205         return 0;
206     }
207
208     if (_pDouble2->isScalar())
209     {
210         bool bComplex1  = _pDouble1->isComplex();
211         bool bComplex2  = _pDouble2->isComplex();
212
213         (*_pDoubleOut) = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray(), bComplex1 | bComplex2);
214
215         if (bComplex1 == false && bComplex2 == false)
216         {
217             //Real Matrix by Real Scalar
218             iMultiRealScalarByRealMatrix(_pDouble2->get(0), _pDouble1->get(), _pDouble1->getSize(), 1, (*_pDoubleOut)->get());
219         }
220         else if (bComplex1 == false && bComplex2 == true)
221         {
222             //Real Matrix by Scalar Complex
223             iMultiComplexScalarByRealMatrix(_pDouble2->get(0), _pDouble2->getImg(0), _pDouble1->get(), _pDouble1->getSize(), 1, (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
224         }
225         else if (bComplex1 == true && bComplex2 == false)
226         {
227             iMultiRealScalarByComplexMatrix(_pDouble2->get(0, 0), _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getSize(), 1, (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
228         }
229         else //if(bComplex1 == true && bComplex2 == true)
230         {
231             iMultiComplexScalarByComplexMatrix(_pDouble2->get(0, 0), _pDouble2->getImg(0, 0), _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getSize(), 1, (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
232         }
233
234         return 0;
235     }
236
237     if (_pDouble1->getDims() > 2 || _pDouble2->getDims() > 2 || _pDouble1->getCols() != _pDouble2->getRows())
238     {
239         return 1;
240     }
241
242     bool bComplex1  = _pDouble1->isComplex();
243     bool bComplex2  = _pDouble2->isComplex();
244     (*_pDoubleOut) = new Double(_pDouble1->getRows(), _pDouble2->getCols(), bComplex1 | bComplex2);
245
246     if (bComplex1 == false && bComplex2 == false)
247     {
248         //Real Matrix by Real Matrix
249         iMultiRealMatrixByRealMatrix(
250             _pDouble1->get(), _pDouble1->getRows(), _pDouble1->getCols(),
251             _pDouble2->get(), _pDouble2->getRows(), _pDouble2->getCols(),
252             (*_pDoubleOut)->get());
253     }
254     else if (bComplex1 == false && bComplex2 == true)
255     {
256         //Real Matrix by Matrix Complex
257         iMultiRealMatrixByComplexMatrix(
258             _pDouble1->get(), _pDouble1->getRows(), _pDouble1->getCols(),
259             _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
260             (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
261     }
262     else if (bComplex1 == true && bComplex2 == false)
263     {
264         //Complex Matrix by Real Matrix
265         iMultiComplexMatrixByRealMatrix(
266             _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
267             _pDouble2->get(), _pDouble2->getRows(), _pDouble2->getCols(),
268             (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
269     }
270     else //if(bComplex1 == true && bComplex2 == true)
271     {
272         //Complex Matrix by Complex Matrix
273         iMultiComplexMatrixByComplexMatrix(
274             _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
275             _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
276             (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
277     }
278     return 0;
279 }
280
281 int DotMultiplyDoubleByDouble(Double* _pDouble1, Double* _pDouble2, Double**  _pDoubleOut)
282 {
283     bool bComplex1  = _pDouble1->isComplex();
284     bool bComplex2  = _pDouble2->isComplex();
285     bool bScalar1   = _pDouble1->isScalar();
286     bool bScalar2   = _pDouble2->isScalar();
287
288     if (bScalar1)
289     {
290         (*_pDoubleOut) = new Double(_pDouble2->getDims(), _pDouble2->getDimsArray(), _pDouble1->isComplex() | _pDouble2->isComplex());
291         if (bComplex1 == false && bComplex2 == false)
292         {
293             iMultiRealScalarByRealMatrix(_pDouble1->get(0), _pDouble2->get(), _pDouble2->getSize(), 1, (*_pDoubleOut)->get());
294         }
295         else if (bComplex1 == false && bComplex2 == true)
296         {
297             iMultiRealScalarByComplexMatrix(_pDouble1->get(0), _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getSize(), 1, (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
298         }
299         else if (bComplex1 == true && bComplex2 == false)
300         {
301             iMultiComplexScalarByRealMatrix(_pDouble1->get(0), _pDouble1->getImg(0), _pDouble2->get(), _pDouble2->getSize(), 1, (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
302         }
303         else //if(bComplex1 == true && bComplex2 == true)
304         {
305             iMultiComplexScalarByComplexMatrix(_pDouble1->get(0), _pDouble1->getImg(0), _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getSize(), 1, (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
306         }
307
308         return 0;
309     }
310
311     if (bScalar2)
312     {
313         (*_pDoubleOut) = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray(), _pDouble1->isComplex() | _pDouble2->isComplex());
314         if (bComplex1 == false && bComplex2 == false)
315         {
316             //Real Matrix by Real Scalar
317             iMultiRealScalarByRealMatrix(_pDouble2->get(0), _pDouble1->get(), _pDouble1->getSize(), 1, (*_pDoubleOut)->get());
318         }
319         else if (bComplex1 == false && bComplex2 == true)
320         {
321             //Real Matrix by Scalar Complex
322             iMultiComplexScalarByRealMatrix(_pDouble2->get(0), _pDouble2->getImg(0), _pDouble1->get(), _pDouble1->getSize(), 1, (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
323         }
324         else if (bComplex1 == true && bComplex2 == false)
325         {
326             iMultiRealScalarByComplexMatrix(_pDouble2->get(0), _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getSize(), 1, (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
327         }
328         else //if(bComplex1 == true && bComplex2 == true)
329         {
330             iMultiComplexScalarByComplexMatrix(_pDouble2->get(0), _pDouble2->getImg(0), _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getSize(), 1, (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
331         }
332
333         return 0;
334     }
335
336     if (_pDouble1->getDims() != _pDouble2->getDims())
337     {
338         return 1;
339     }
340
341     int* piDims1 = _pDouble1->getDimsArray();
342     int* piDims2 = _pDouble2->getDimsArray();
343
344     for (int i = 0 ; i < _pDouble1->getDims() ; i++)
345     {
346         if (piDims1[i] != piDims2[i])
347         {
348             return 0;
349         }
350     }
351
352     (*_pDoubleOut) = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray(), _pDouble1->isComplex() | _pDouble2->isComplex());
353     if (bComplex1 == false && bComplex2 == false)
354     {
355         iDotMultiplyRealMatrixByRealMatrix(_pDouble1->get(), _pDouble2->get(), (*_pDoubleOut)->get(), _pDouble1->getSize(), 1);
356     }
357     else if (bComplex1 == false && bComplex2 == true)
358     {
359         iDotMultiplyRealMatrixByComplexMatrix(_pDouble1->get(), _pDouble2->get(), _pDouble2->getImg(), (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), _pDouble1->getSize(), 1);
360     }
361     else if (bComplex1 == true && bComplex2 == false)
362     {
363         iDotMultiplyComplexMatrixByRealMatrix(_pDouble1->get(), _pDouble1->getImg(), _pDouble2->get(), (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), _pDouble1->getSize(), 1);
364     }
365     else //if(bComplex1 == true && bComplex2 == true)
366     {
367         iDotMultiplyComplexMatrixByComplexMatrix(_pDouble1->get(), _pDouble1->getImg(), _pDouble2->get(), _pDouble2->getImg(), (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), _pDouble1->getSize(), 1);
368     }
369
370     return 0;
371 }
372
373 int MultiplyDoubleByPoly(Double* _pDouble, Polynom* _pPoly, Polynom** _pPolyOut)
374 {
375     bool bComplex1  = _pDouble->isComplex();
376     bool bComplex2  = _pPoly->isComplex();
377
378     if (_pDouble->isScalar())
379     {
380         int* piRank = new int[_pPoly->getSize()];
381         for (int i = 0 ; i < _pPoly->getSize() ; i++)
382         {
383             piRank[i] = _pPoly->get(i)->getRank();
384         }
385
386         (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pPoly->getDims(), _pPoly->getDimsArray(), piRank);
387         delete[] piRank;
388         if (bComplex1 || bComplex2)
389         {
390             (*_pPolyOut)->setComplex(true);
391         }
392
393         for (int i = 0 ; i < _pPoly->getSize() ; i++)
394         {
395             SinglePoly *pPolyIn     = _pPoly->get(i);
396             double* pRealIn         = pPolyIn->get();
397             double* pImgIn          = pPolyIn->getImg();
398
399             SinglePoly *pPolyOut    = (*_pPolyOut)->get(i);
400             double* pRealOut        = pPolyOut->get();
401             double* pImgOut         = pPolyOut->getImg();
402
403             if (bComplex1 == false && bComplex2 == false)
404             {
405                 iMultiRealScalarByRealMatrix(_pDouble->get(0), pRealIn, 1, pPolyIn->getSize(), pRealOut);
406             }
407             else if (bComplex1 == false && bComplex2 == true)
408             {
409                 iMultiRealScalarByComplexMatrix(_pDouble->get(0), pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
410             }
411             else if (bComplex1 == true && bComplex2 == false)
412             {
413                 iMultiComplexScalarByRealMatrix(_pDouble->get(0), _pDouble->getImg(0), pRealIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
414             }
415             else if (bComplex1 == true && bComplex2 == true)
416             {
417                 iMultiComplexScalarByComplexMatrix(_pDouble->get(0), _pDouble->getImg(0), pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
418             }
419         }
420         (*_pPolyOut)->updateRank();
421         return 0;
422     }
423
424     if (_pPoly->isScalar())
425     {
426         int* piRank = new int[_pDouble->getSize()];
427         for (int i = 0 ; i < _pDouble->getSize() ; i++)
428         {
429             piRank[i] = _pPoly->get(0)->getRank();
430         }
431
432         (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pDouble->getDims(), _pDouble->getDimsArray(), piRank);
433         delete[] piRank;
434         if (bComplex1 || bComplex2)
435         {
436             (*_pPolyOut)->setComplex(true);
437         }
438
439         double *pDoubleR    = _pDouble->get();
440         double *pDoubleI    = _pDouble->getImg();
441
442         SinglePoly *pPolyIn = _pPoly->get(0);
443         double* pRealIn     = pPolyIn->get();
444         double* pImgIn      = pPolyIn->getImg();
445
446         for (int i = 0 ; i < _pDouble->getSize() ; i++)
447         {
448             SinglePoly *pPolyOut    = (*_pPolyOut)->get(i);
449             double* pRealOut        = pPolyOut->get();
450             double* pImgOut         = pPolyOut->getImg();
451
452             if (bComplex1 == false && bComplex2 == false)
453             {
454                 iMultiRealScalarByRealMatrix(pDoubleR[i], pRealIn, 1, pPolyIn->getSize(), pRealOut);
455             }
456             else if (bComplex1 == false && bComplex2 == true)
457             {
458                 iMultiRealScalarByComplexMatrix(pDoubleR[i], pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
459             }
460             else if (bComplex1 == true && bComplex2 == false)
461             {
462                 iMultiComplexScalarByRealMatrix(pDoubleR[i], pDoubleI[i], pRealIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
463             }
464             else if (bComplex1 == true && bComplex2 == true)
465             {
466                 iMultiComplexScalarByComplexMatrix(pDoubleR[i], pDoubleI[i], pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
467             }
468         }
469
470         (*_pPolyOut)->updateRank();
471         return 0;
472     }
473
474     if (_pPoly->getDims() > 2 || _pDouble->getDims() > 2 || _pDouble->getCols() != _pPoly->getRows())
475     {
476         return 1;
477     }
478
479     int* piRank = new int[_pDouble->getRows() * _pPoly->getCols()];
480     int iMaxRank = _pPoly->getMaxRank();
481     for (int i = 0 ; i < _pDouble->getRows() * _pPoly->getCols() ; i++)
482     {
483         piRank[i] = iMaxRank;
484     }
485
486     (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pDouble->getRows(), _pPoly->getCols(), piRank);
487     delete[] piRank;
488     if (bComplex1 || bComplex2)
489     {
490         (*_pPolyOut)->setComplex(true);
491     }
492
493     Double *pCoef = _pPoly->getCoef();
494     Double *pTemp = new Double(_pDouble->getRows(), pCoef->getCols(), bComplex1 || bComplex2);
495
496     if (bComplex1 == false && bComplex2 == false)
497     {
498         iMultiRealMatrixByRealMatrix(_pDouble->get(), _pDouble->getRows(), _pDouble->getCols(),
499                                      pCoef->get(), pCoef->getRows(), pCoef->getCols(),
500                                      pTemp->get());
501     }
502     else if (bComplex1 == false && bComplex2 == true)
503     {
504         iMultiRealMatrixByComplexMatrix(_pDouble->get(), _pDouble->getRows(), _pDouble->getCols(),
505                                         pCoef->get(), pCoef->getImg(), pCoef->getRows(), pCoef->getCols(),
506                                         pTemp->get(), pTemp->getImg());
507
508     }
509     else if (bComplex1 == true && bComplex2 == false)
510     {
511         iMultiComplexMatrixByRealMatrix(_pDouble->get(), _pDouble->getImg(), _pDouble->getRows(), _pDouble->getCols(),
512                                         pCoef->get(), pCoef->getRows(), pCoef->getCols(),
513                                         pTemp->get(), pTemp->getImg());
514     }
515     else //if(bComplex1 == true && bComplex2 == true)
516     {
517         iMultiComplexMatrixByComplexMatrix(_pDouble->get(), _pDouble->getImg(), _pDouble->getRows(), _pDouble->getCols(),
518                                            pCoef->get(), pCoef->getImg(), pCoef->getRows(), pCoef->getCols(),
519                                            pTemp->get(), pTemp->getImg());
520     }
521
522     (*_pPolyOut)->setCoef(pTemp);
523     (*_pPolyOut)->updateRank();
524     delete pTemp;
525     return 0;
526 }
527
528 int MultiplyPolyByDouble(Polynom* _pPoly, Double* _pDouble, Polynom **_pPolyOut)
529 {
530     bool bComplex1  = _pPoly->isComplex();
531     bool bComplex2  = _pDouble->isComplex();
532     bool bScalar1   = _pPoly->isScalar();
533     bool bScalar2   = _pDouble->isScalar();
534
535     if (bScalar1)
536     {
537         int* piRank = new int[_pDouble->getSize()];
538         for (int i = 0 ; i < _pDouble->getSize() ; i++)
539         {
540             piRank[i] = _pPoly->get(0)->getRank();
541         }
542
543         (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pDouble->getDims(), _pDouble->getDimsArray(), piRank);
544         delete[] piRank;
545         if (bComplex1 || bComplex2)
546         {
547             (*_pPolyOut)->setComplex(true);
548         }
549
550         double *pDoubleR    = _pDouble->get();
551         double *pDoubleI    = _pDouble->getImg();
552
553         SinglePoly *pPolyIn = _pPoly->get(0);
554         double* pRealIn     = pPolyIn->get();
555         double* pImgIn      = pPolyIn->getImg();
556
557         for (int i = 0 ; i < _pDouble->getSize() ; i++)
558         {
559             SinglePoly *pPolyOut    = (*_pPolyOut)->get(i);
560             double* pRealOut        = pPolyOut->get();
561             double* pImgOut         = pPolyOut->getImg();
562
563             if (bComplex1 == false && bComplex2 == false)
564             {
565                 iMultiRealScalarByRealMatrix(pDoubleR[i], pRealIn, 1, pPolyIn->getSize(), pRealOut);
566             }
567             else if (bComplex1 == false && bComplex2 == true)
568             {
569                 iMultiComplexScalarByRealMatrix(pDoubleR[i], pDoubleI[i], pRealIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
570             }
571             else if (bComplex1 == true && bComplex2 == false)
572             {
573                 iMultiRealScalarByComplexMatrix(pDoubleR[i], pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
574             }
575             else if (bComplex1 == true && bComplex2 == true)
576             {
577                 iMultiComplexScalarByComplexMatrix(pDoubleR[i], pDoubleI[i], pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
578             }
579         }
580
581         (*_pPolyOut)->updateRank();
582         return 0;
583     }
584     else if (bScalar2)
585     {
586         int* piRank = new int[_pPoly->getSize()];
587         for (int i = 0 ; i < _pPoly->getSize() ; i++)
588         {
589             piRank[i] = _pPoly->get(i)->getRank();
590         }
591
592         (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pPoly->getDims(), _pPoly->getDimsArray(), piRank);
593         delete[] piRank;
594         if (bComplex1 || bComplex2)
595         {
596             (*_pPolyOut)->setComplex(true);
597         }
598
599         for (int i = 0 ; i < _pPoly->getSize() ; i++)
600         {
601             SinglePoly *pPolyIn = _pPoly->get(i);
602             double* pRealIn     = pPolyIn->get();
603             double* pImgIn      = pPolyIn->getImg();
604
605             SinglePoly *pPolyOut    = (*_pPolyOut)->get(i);
606             double* pRealOut        = pPolyOut->get();
607             double* pImgOut         = pPolyOut->getImg();
608
609             if (bComplex1 == false && bComplex2 == false)
610             {
611                 iMultiRealScalarByRealMatrix(_pDouble->get(0), pRealIn, 1, pPolyIn->getSize(), pRealOut);
612             }
613             else if (bComplex1 == false && bComplex2 == true)
614             {
615                 iMultiComplexScalarByRealMatrix(_pDouble->get(0), _pDouble->getImg(0), pRealIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
616             }
617             else if (bComplex1 == true && bComplex2 == false)
618             {
619                 iMultiRealScalarByComplexMatrix(_pDouble->get(0), pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
620             }
621             else if (bComplex1 == true && bComplex2 == true)
622             {
623                 iMultiComplexScalarByComplexMatrix(_pDouble->get(0), _pDouble->getImg(0), pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
624             }
625         }
626
627         (*_pPolyOut)->updateRank();
628         return 0;
629     }
630
631     if (_pDouble->getDims() > 2 || _pPoly->getDims() > 2 || _pPoly->getCols() != _pDouble->getRows())
632     {
633         return 1;
634     }
635
636     int* piRank = new int[_pPoly->getRows() * _pDouble->getCols()];
637     int iMaxRank = _pPoly->getMaxRank();
638     for (int i = 0 ; i < _pPoly->getRows() * _pDouble->getCols() ; i++)
639     {
640         piRank[i] = iMaxRank;
641     }
642
643     (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pPoly->getRows(), _pDouble->getCols(), piRank);
644     delete[] piRank;
645     if (bComplex1 || bComplex2)
646     {
647         (*_pPolyOut)->setComplex(true);
648     }
649
650     //Distribution a la mano par appels a des sous-fonctions ( iMulti...ScalarBy...Scalar ) plus iAdd...To... )
651
652     //for each line of _pPoly
653     for (int iRow1 = 0 ; iRow1 < _pPoly->getRows() ; iRow1++)
654     {
655         //for each col of _pDouble
656         for (int iCol2 = 0 ; iCol2 < _pDouble->getCols() ; iCol2++)
657         {
658             SinglePoly* pSPOut = (*_pPolyOut)->get(iRow1, iCol2);
659             pSPOut->setZeros();
660
661             //for each rows of _pDouble / cols of _pPoly
662             for (int iRow2 = 0 ; iRow2 < _pDouble->getRows() ; iRow2++)
663             {
664                 // SinglePoly(iRow1, iRow2) * Double(iRow2, iCol2)
665                 SinglePoly* pSPIn = _pPoly->get(iRow1, iRow2);
666                 int iSize = pSPIn->getSize();
667                 double* pdblMult = new double[iSize];
668
669                 if (bComplex1 == false && bComplex2 == false)
670                 {
671                     //Real Matrix by Real Scalar
672                     iMultiRealScalarByRealMatrix(_pDouble->get(iRow2, iCol2), pSPIn->get(), iSize, 1, pdblMult);
673                     add(pSPOut->get(), (long long)iSize, pdblMult, pSPOut->get());
674                 }
675                 else if (bComplex1 == false && bComplex2 == true)
676                 {
677                     //Real Matrix by Scalar Complex
678                     double* pdblMultImg = new double[iSize];
679                     iMultiComplexScalarByRealMatrix(_pDouble->get(iRow2, iCol2), _pDouble->getImg(iRow2, iCol2), pSPIn->get(), pSPIn->getSize(), 1, pdblMult, pdblMultImg);
680                     add(pSPOut->get(), pSPOut->getImg(), (long long)iSize, pdblMult, pdblMultImg, pSPOut->get(), pSPOut->getImg());
681                     delete[] pdblMultImg;
682                 }
683                 else if (bComplex1 == true && bComplex2 == false)
684                 {
685                     double* pdblMultImg = new double[iSize];
686                     iMultiRealScalarByComplexMatrix(_pDouble->get(iRow2, iCol2), pSPIn->get(), pSPIn->getImg(), pSPIn->getSize(), 1, pdblMult, pdblMultImg);
687                     add(pSPOut->get(), pSPOut->getImg(), (long long)iSize, pdblMult, pdblMultImg, pSPOut->get(), pSPOut->getImg());
688                     delete[] pdblMultImg;
689                 }
690                 else //if(bComplex1 == true && bComplex2 == true)
691                 {
692                     double* pdblMultImg = new double[iSize];
693                     iMultiComplexScalarByComplexMatrix(_pDouble->get(iRow2, iCol2), _pDouble->getImg(iRow2, iCol2), pSPIn->get(), pSPIn->getImg(), pSPIn->getSize(), 1, pdblMult, pdblMultImg);
694                     add(pSPOut->get(), pSPOut->getImg(), (long long)iSize, pdblMult, pdblMultImg, pSPOut->get(), pSPOut->getImg());
695                     delete[] pdblMultImg;
696                 }
697
698                 delete[] pdblMult;
699             }//for(int iRow2 = 0 ; iRow2 < _pDouble->getRows() ; iRow2++)
700         }//for(int iCol2 = 0 ; iCol2 < _pDouble->getCols() ; iCol2++)
701     }//for(int iRow1 = 0 ; iRow1 < _pPoly->getRows() ; iRow1++)
702
703     (*_pPolyOut)->updateRank();
704     return 0;
705 }
706
707 int MultiplyPolyByPoly(Polynom* _pPoly1, Polynom* _pPoly2, Polynom** _pPolyOut)
708 {
709     bool bComplex1  = _pPoly1->isComplex();
710     bool bComplex2  = _pPoly2->isComplex();
711
712     if (_pPoly1->isScalar() && _pPoly2->isScalar())
713     {
714         //poly1(0) * poly2(0)
715         int iRank = _pPoly1->get(0)->getRank() + _pPoly2->get(0)->getRank();
716         (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), 1, 1, &iRank);
717         if (bComplex1 || bComplex2)
718         {
719             (*_pPolyOut)->setComplex(true);
720         }
721
722         if (bComplex1 == false && bComplex2 == false)
723         {
724             SinglePoly *pPoly1  = _pPoly1->get(0);
725             SinglePoly *pPoly2  = _pPoly2->get(0);
726             SinglePoly *pPolyOut = (*_pPolyOut)->get(0);
727
728             pPolyOut->setZeros();
729
730             iMultiScilabPolynomByScilabPolynom(
731                 pPoly1->get(), pPoly1->getSize(),
732                 pPoly2->get(), pPoly2->getSize(),
733                 pPolyOut->get(), pPolyOut->getSize());
734         }
735         else if (bComplex1 == false && bComplex2 == true)
736         {
737             SinglePoly *pPoly1  = _pPoly1->get(0);
738             SinglePoly *pPoly2  = _pPoly2->get(0);
739             SinglePoly *pPolyOut = (*_pPolyOut)->get(0);
740
741             pPolyOut->setZeros();
742
743             iMultiScilabPolynomByComplexPoly(
744                 pPoly1->get(), pPoly1->getSize(),
745                 pPoly2->get(), pPoly2->getImg(), pPoly2->getSize(),
746                 pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
747         }
748         else if (bComplex1 == true && bComplex2 == false)
749         {
750             SinglePoly *pPoly1  = _pPoly1->get(0);
751             SinglePoly *pPoly2  = _pPoly2->get(0);
752             SinglePoly *pPolyOut = (*_pPolyOut)->get(0);
753
754             pPolyOut->setZeros();
755
756             iMultiComplexPolyByScilabPolynom(
757                 pPoly1->get(), pPoly1->getImg(), pPoly1->getSize(),
758                 pPoly2->get(), pPoly2->getSize(),
759                 pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
760         }
761         else if (bComplex1 == true && bComplex2 == true)
762         {
763             SinglePoly *pPoly1   = _pPoly1->get(0);
764             SinglePoly *pPoly2   = _pPoly2->get(0);
765             SinglePoly *pPolyOut  = (*_pPolyOut)->get(0);
766
767             pPolyOut->setZeros();
768
769             iMultiComplexPolyByComplexPoly(
770                 pPoly1->get(), pPoly1->getImg(), pPoly1->getSize(),
771                 pPoly2->get(), pPoly2->getImg(), pPoly2->getSize(),
772                 pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
773         }
774
775         (*_pPolyOut)->updateRank();
776         return 0;
777     }
778
779     if (_pPoly1->isScalar())
780     {
781         //poly1(0) * poly2(n)
782         int* piRank = new int[_pPoly2->getSize()];
783         for (int i = 0 ; i < _pPoly2->getSize() ; i++)
784         {
785             piRank[i] = _pPoly1->get(0)->getRank() + _pPoly2->get(i)->getRank();
786         }
787
788         (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), _pPoly2->getDims(), _pPoly2->getDimsArray(), piRank);
789         if (bComplex1 || bComplex2)
790         {
791             (*_pPolyOut)->setComplex(true);
792         }
793         delete[] piRank;
794
795
796         SinglePoly *pPoly1  = _pPoly1->get(0);
797         if (bComplex1 == false && bComplex2 == false)
798         {
799             for (int iPoly = 0 ; iPoly < _pPoly2->getSize() ; iPoly++)
800             {
801                 SinglePoly *pPoly2  = _pPoly2->get(iPoly);
802                 SinglePoly *pPolyOut = (*_pPolyOut)->get(iPoly);
803
804                 pPolyOut->setZeros();
805
806                 iMultiScilabPolynomByScilabPolynom(
807                     pPoly1->get(), pPoly1->getSize(),
808                     pPoly2->get(), pPoly2->getSize(),
809                     pPolyOut->get(), pPolyOut->getSize());
810             }
811         }
812         else if (bComplex1 == false && bComplex2 == true)
813         {
814             for (int iPoly = 0 ; iPoly < _pPoly2->getSize() ; iPoly++)
815             {
816                 SinglePoly *pPoly2  = _pPoly2->get(iPoly);
817                 SinglePoly *pPolyOut = (*_pPolyOut)->get(iPoly);
818
819                 pPolyOut->setZeros();
820
821                 iMultiScilabPolynomByComplexPoly(
822                     pPoly1->get(), pPoly1->getSize(),
823                     pPoly2->get(), pPoly2->getImg(), pPoly2->getSize(),
824                     pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
825             }
826         }
827         else if (bComplex1 == true && bComplex2 == false)
828         {
829             for (int iPoly = 0 ; iPoly < _pPoly2->getSize() ; iPoly++)
830             {
831                 SinglePoly *pPoly2  = _pPoly2->get(iPoly);
832                 SinglePoly *pPolyOut = (*_pPolyOut)->get(iPoly);
833
834                 pPolyOut->setZeros();
835
836                 iMultiComplexPolyByScilabPolynom(
837                     pPoly1->get(), pPoly1->getImg(), pPoly1->getSize(),
838                     pPoly2->get(), pPoly2->getSize(),
839                     pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
840             }
841         }
842         else if (bComplex1 == true && bComplex2 == true)
843         {
844             for (int iPoly = 0 ; iPoly < _pPoly2->getSize() ; iPoly++)
845             {
846                 SinglePoly *pPoly2   = _pPoly2->get(iPoly);
847                 SinglePoly *pPolyOut  = (*_pPolyOut)->get(iPoly);
848
849                 pPolyOut->setZeros();
850
851                 iMultiComplexPolyByComplexPoly(
852                     pPoly1->get(), pPoly1->getImg(), pPoly1->getSize(),
853                     pPoly2->get(), pPoly2->getImg(), pPoly2->getSize(),
854                     pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
855             }
856         }
857
858         (*_pPolyOut)->updateRank();
859         return 0;
860     }
861
862     if (_pPoly2->isScalar())
863     {
864         //poly1(n) * poly2(0)
865         int* piRank = new int[_pPoly1->getSize()];
866         for (int i = 0 ; i < _pPoly1->getSize() ; i++)
867         {
868             piRank[i] = _pPoly2->get(0)->getRank() + _pPoly1->get(i)->getRank();
869         }
870
871         (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), _pPoly1->getDims(), _pPoly1->getDimsArray(), piRank);
872         if (bComplex1 || bComplex2)
873         {
874             (*_pPolyOut)->setComplex(true);
875         }
876         delete[] piRank;
877
878         SinglePoly *pPoly2  = _pPoly2->get(0);
879         if (bComplex1 == false && bComplex2 == false)
880         {
881             for (int iPoly = 0 ; iPoly < _pPoly1->getSize() ; iPoly++)
882             {
883                 SinglePoly *pPoly1  = _pPoly1->get(iPoly);
884                 SinglePoly *pPolyOut = (*_pPolyOut)->get(iPoly);
885
886                 pPolyOut->setZeros();
887
888                 iMultiScilabPolynomByScilabPolynom(
889                     pPoly1->get(), pPoly1->getSize(),
890                     pPoly2->get(), pPoly2->getSize(),
891                     pPolyOut->get(), pPolyOut->getSize());
892             }
893         }
894         else if (bComplex1 == false && bComplex2 == true)
895         {
896             for (int iPoly = 0 ; iPoly < _pPoly1->getSize() ; iPoly++)
897             {
898                 SinglePoly *pPoly1  = _pPoly1->get(iPoly);
899                 SinglePoly *pPolyOut = (*_pPolyOut)->get(iPoly);
900
901                 pPolyOut->setZeros();
902
903                 iMultiScilabPolynomByComplexPoly(
904                     pPoly1->get(), pPoly1->getSize(),
905                     pPoly2->get(), pPoly2->getImg(), pPoly2->getSize(),
906                     pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
907             }
908         }
909         else if (bComplex1 == true && bComplex2 == false)
910         {
911             for (int iPoly = 0 ; iPoly < _pPoly1->getSize() ; iPoly++)
912             {
913                 SinglePoly *pPoly1  = _pPoly1->get(iPoly);
914                 SinglePoly *pPolyOut = (*_pPolyOut)->get(iPoly);
915
916                 pPolyOut->setZeros();
917
918                 iMultiComplexPolyByScilabPolynom(
919                     pPoly1->get(), pPoly1->getImg(), pPoly1->getSize(),
920                     pPoly2->get(), pPoly2->getSize(),
921                     pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
922             }
923         }
924         else if (bComplex1 == true && bComplex2 == true)
925         {
926             for (int iPoly = 0 ; iPoly < _pPoly1->getSize() ; iPoly++)
927             {
928                 SinglePoly *pPoly1   = _pPoly1->get(iPoly);
929                 SinglePoly *pPolyOut  = (*_pPolyOut)->get(iPoly);
930
931                 pPolyOut->setZeros();
932
933                 iMultiComplexPolyByComplexPoly(
934                     pPoly1->get(), pPoly1->getImg(), pPoly1->getSize(),
935                     pPoly2->get(), pPoly2->getImg(), pPoly2->getSize(),
936                     pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
937             }
938         }
939
940         (*_pPolyOut)->updateRank();
941         return 0;
942     }
943
944     if (_pPoly1->getDims() > 2 || _pPoly2->getDims() > 2 || _pPoly1->getCols() != _pPoly2->getRows())
945     {
946         return 1;
947     }
948
949     // matrix by matrix
950     int* piRank = new int[_pPoly1->getRows() * _pPoly2->getCols()];
951     int iMaxRank = _pPoly1->getMaxRank() + _pPoly2->getMaxRank();
952     for (int i = 0 ; i < _pPoly1->getRows() * _pPoly2->getCols() ; i++)
953     {
954         piRank[i] = iMaxRank;
955     }
956
957     (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), _pPoly1->getRows(), _pPoly2->getCols(), piRank);
958     if (bComplex1 || bComplex2)
959     {
960         (*_pPolyOut)->setComplex(true);
961     }
962
963     delete[] piRank;
964
965
966     if (bComplex1 == false && bComplex2 == false)
967     {
968         double *pReal = NULL;
969         SinglePoly *pTemp  = new SinglePoly(&pReal, (*_pPolyOut)->getMaxRank());
970
971         for (int iRow = 0 ; iRow < _pPoly1->getRows() ; iRow++)
972         {
973             for (int iCol = 0 ; iCol < _pPoly2->getCols() ; iCol++)
974             {
975                 SinglePoly *pResult = (*_pPolyOut)->get(iRow, iCol);
976                 pResult->setZeros();
977
978                 for (int iCommon = 0 ; iCommon < _pPoly1->getCols() ; iCommon++)
979                 {
980                     SinglePoly *pL   = _pPoly1->get(iRow, iCommon);
981                     SinglePoly *pR   = _pPoly2->get(iCommon, iCol);
982
983                     pTemp->setZeros();
984
985                     iMultiScilabPolynomByScilabPolynom(
986                         pL->get(), pL->getSize(),
987                         pR->get(), pR->getSize(),
988                         pTemp->get(), pL->getRank() + pR->getRank() + 1);
989
990                     iAddScilabPolynomToScilabPolynom(
991                         pResult->get(), pResult->getSize(),
992                         pTemp->get(), pResult->getSize(),
993                         pResult->get(), pResult->getSize());
994                 }
995             }
996         }
997     }
998     else if (bComplex1 == false && bComplex2 == true)
999     {
1000         double *pReal = NULL;
1001         double *pImg = NULL;
1002         SinglePoly *pTemp  = new SinglePoly(&pReal, &pImg, (*_pPolyOut)->getMaxRank());
1003
1004         for (int iRow = 0 ; iRow < _pPoly1->getRows() ; iRow++)
1005         {
1006             for (int iCol = 0 ; iCol < _pPoly2->getCols() ; iCol++)
1007             {
1008                 SinglePoly *pResult = (*_pPolyOut)->get(iRow, iCol);
1009                 pResult->setZeros();
1010
1011                 for (int iCommon = 0 ; iCommon < _pPoly1->getCols() ; iCommon++)
1012                 {
1013                     SinglePoly *pL   = _pPoly1->get(iRow, iCommon);
1014                     SinglePoly *pR   = _pPoly2->get(iCommon, iCol);
1015
1016                     pTemp->setZeros();
1017
1018                     iMultiScilabPolynomByComplexPoly(
1019                         pL->get(), pL->getSize(),
1020                         pR->get(), pR->getImg(), pR->getSize(),
1021                         pTemp->get(), pTemp->getImg(), pL->getRank() + pR->getRank() + 1);
1022
1023                     iAddComplexPolyToComplexPoly(
1024                         pResult->get(), pResult->getImg(), pResult->getSize(),
1025                         pTemp->get(), pTemp->getImg(), pResult->getSize(),
1026                         pResult->get(), pResult->getImg(), pResult->getSize());
1027                 }
1028             }
1029         }
1030     }
1031     else if (bComplex1 == true && bComplex2 == false)
1032     {
1033         double *pReal = NULL;
1034         double *pImg = NULL;
1035         SinglePoly *pTemp  = new SinglePoly(&pReal, &pImg, (*_pPolyOut)->getMaxRank());
1036
1037         for (int iRow = 0 ; iRow < _pPoly1->getRows() ; iRow++)
1038         {
1039             for (int iCol = 0 ; iCol < _pPoly2->getCols() ; iCol++)
1040             {
1041                 SinglePoly *pResult = (*_pPolyOut)->get(iRow, iCol);
1042                 pResult->setZeros();
1043
1044                 for (int iCommon = 0 ; iCommon < _pPoly1->getCols() ; iCommon++)
1045                 {
1046                     SinglePoly *pL   = _pPoly1->get(iRow, iCommon);
1047                     SinglePoly *pR   = _pPoly2->get(iCommon, iCol);
1048
1049                     pTemp->setZeros();
1050
1051                     iMultiScilabPolynomByComplexPoly(
1052                         pR->get(), pR->getSize(),
1053                         pL->get(), pL->getImg(), pL->getSize(),
1054                         pTemp->get(), pTemp->getImg(), pL->getRank() + pR->getRank() + 1);
1055
1056                     iAddComplexPolyToComplexPoly(
1057                         pResult->get(), pResult->getImg(), pResult->getSize(),
1058                         pTemp->get(), pTemp->getImg(), pResult->getSize(),
1059                         pResult->get(), pResult->getImg(), pResult->getSize());
1060                 }
1061             }
1062         }
1063     }
1064     else if (bComplex1 == true && bComplex2 == true)
1065     {
1066         double *pReal = NULL;
1067         double *pImg = NULL;
1068         SinglePoly *pTemp  = new SinglePoly(&pReal, &pImg, (*_pPolyOut)->getMaxRank());
1069
1070         for (int iRow = 0 ; iRow < _pPoly1->getRows() ; iRow++)
1071         {
1072             for (int iCol = 0 ; iCol < _pPoly2->getCols() ; iCol++)
1073             {
1074                 SinglePoly *pResult = (*_pPolyOut)->get(iRow, iCol);
1075                 pResult->setZeros();
1076
1077                 for (int iCommon = 0 ; iCommon < _pPoly1->getCols() ; iCommon++)
1078                 {
1079                     SinglePoly *pL   = _pPoly1->get(iRow, iCommon);
1080                     SinglePoly *pR   = _pPoly2->get(iCommon, iCol);
1081
1082                     pTemp->setZeros();
1083
1084                     iMultiComplexPolyByComplexPoly(
1085                         pL->get(), pL->getImg(), pL->getSize(),
1086                         pR->get(), pR->getImg(), pR->getSize(),
1087                         pTemp->get(), pTemp->getImg(), pL->getRank() + pR->getRank() + 1);
1088
1089                     iAddComplexPolyToComplexPoly(
1090                         pResult->get(), pResult->getImg(), pResult->getSize(),
1091                         pTemp->get(), pTemp->getImg(), pResult->getSize(),
1092                         pResult->get(), pResult->getImg(), pResult->getSize());
1093                 }
1094             }
1095         }
1096     }
1097     (*_pPolyOut)->updateRank();
1098
1099     return 0;
1100 }
1101
1102 int MultiplySparseBySparse(Sparse* _pSparse1, Sparse* _pSparse2, Sparse** _pSparseOut)
1103 {
1104     if (_pSparse1->isScalar())
1105     {
1106         //scalar * sp
1107         Double* pDbl = NULL;
1108         if (_pSparse1->isComplex())
1109         {
1110             std::complex<double> dbl = _pSparse1->getImg(0, 0);
1111             pDbl = new Double(dbl.real(), dbl.imag());
1112         }
1113         else
1114         {
1115             pDbl = new Double(_pSparse1->get(0, 0));
1116         }
1117
1118         MultiplyDoubleBySparse(pDbl, _pSparse2, (GenericType**)_pSparseOut);
1119         delete pDbl;
1120         return 0;
1121     }
1122
1123     if (_pSparse2->isScalar())
1124     {
1125         //sp * scalar
1126         Double* pDbl = NULL;
1127         if (_pSparse2->isComplex())
1128         {
1129             std::complex<double> dbl = _pSparse2->getImg(0, 0);
1130             pDbl = new Double(dbl.real(), dbl.imag());
1131         }
1132         else
1133         {
1134             pDbl = new Double(_pSparse2->get(0, 0));
1135         }
1136
1137         MultiplySparseByDouble(_pSparse1, pDbl, (GenericType**)_pSparseOut);
1138         delete pDbl;
1139         return 0;
1140     }
1141
1142     if (_pSparse1->getCols() != _pSparse2->getRows())
1143     {
1144         return 1;
1145     }
1146
1147     *_pSparseOut = _pSparse1->multiply(*_pSparse2);
1148     return 0;
1149 }
1150
1151 int MultiplyDoubleBySparse(Double* _pDouble, Sparse *_pSparse, GenericType** _pOut)
1152 {
1153     //D * SP
1154     if (_pDouble->isScalar())
1155     {
1156         //d * SP -> SP
1157         Sparse* pOut = NULL;
1158         if (_pDouble->isComplex())
1159         {
1160             std::complex<double> dbl(_pDouble->get(0), _pDouble->getImg(0));
1161             pOut = _pSparse->multiply(dbl);
1162         }
1163         else
1164         {
1165             pOut = _pSparse->multiply(_pDouble->get(0));
1166         }
1167         *_pOut = pOut;
1168         return 0;
1169     }
1170
1171     if (_pSparse->isScalar())
1172     {
1173         //D * sp -> D .* d
1174         Double* pD = NULL;
1175
1176         if (_pSparse->isComplex())
1177         {
1178             std::complex<double> dbl(_pSparse->getImg(0, 0));
1179             pD = new Double(dbl.real(), dbl.imag());
1180         }
1181         else
1182         {
1183             pD = new Double(_pSparse->get(0, 0));
1184         }
1185
1186         InternalType* pIT = GenericDotTimes(_pDouble, pD);
1187         *_pOut = pIT->getAs<GenericType>();
1188         delete pD;
1189         return 0;
1190     }
1191
1192     if (_pDouble->getCols() != _pSparse->getRows())
1193     {
1194         return 1;
1195     }
1196
1197     //try to be smart and only compute for non zero values
1198
1199     //get some information
1200     int iNonZeros = static_cast<int>(_pSparse->nonZeros());
1201     int* pRows = new int[iNonZeros * 2];
1202     _pSparse->outputRowCol(pRows);
1203     int* pCols = pRows + iNonZeros;
1204     double* pValR = new double[iNonZeros];
1205     double* pValI = new double[iNonZeros];
1206     _pSparse->outputValues(pValR, pValI);
1207
1208     Double* pOut = new Double(_pDouble->getRows(), _pSparse->getCols(), _pDouble->isComplex() | _pSparse->isComplex());
1209     pOut->setZeros();
1210
1211     if (_pDouble->isComplex() == false && _pSparse->isComplex() == false)
1212     {
1213         for (int i = 0 ; i < iNonZeros ; i++)
1214         {
1215             int iRow = static_cast<int>(pRows[i]) - 1;
1216             int iCol = static_cast<int>(pCols[i]) - 1;
1217             double dbl = pValR[i];
1218
1219             for (int j = 0 ; j < _pDouble->getRows() ; j++)
1220             {
1221                 double dblVal = _pDouble->get(j, iRow) * dbl;
1222                 pOut->set(j, iCol, pOut->get(j, iCol) + dblVal);
1223             }
1224         }
1225     }
1226     else if (_pDouble->isComplex() == false && _pSparse->isComplex() == true)
1227     {
1228         //a * (b ci) -> ab ac
1229         for (int i = 0 ; i < iNonZeros ; i++)
1230         {
1231             int iRow = static_cast<int>(pRows[i]) - 1;
1232             int iCol = static_cast<int>(pCols[i]) - 1;
1233             double dblR = pValR[i];
1234             double dblI = pValI[i];
1235
1236             for (int j = 0 ; j < _pDouble->getRows() ; j++)
1237             {
1238                 double dblValR = _pDouble->get(j, iRow) * dblR;
1239                 double dblValI = _pDouble->get(j, iRow) * dblI;
1240                 pOut->set(j, iCol, pOut->get(j, iCol) + dblValR);
1241                 pOut->setImg(j, iCol, pOut->getImg(j, iCol) + dblValI);
1242             }
1243         }
1244     }
1245     else if (_pDouble->isComplex() == true && _pSparse->isComplex() == false)
1246     {
1247         //(a bi) * c -> ac + bc
1248         for (int i = 0 ; i < iNonZeros ; i++)
1249         {
1250             int iRow = static_cast<int>(pRows[i]) - 1;
1251             int iCol = static_cast<int>(pCols[i]) - 1;
1252             double dblR = pValR[i];
1253
1254             for (int j = 0 ; j < _pDouble->getRows() ; j++)
1255             {
1256                 double dblValR = _pDouble->get(j, iRow) * dblR;
1257                 double dblValI = _pDouble->getImg(j, iRow) * dblR;
1258                 pOut->set(j, iCol, pOut->get(j, iCol) + dblValR);
1259                 pOut->setImg(j, iCol, pOut->getImg(j, iCol) + dblValI);
1260             }
1261         }
1262     }
1263     else if (_pDouble->isComplex() == true && _pSparse->isComplex() == true)
1264     {
1265         for (int i = 0 ; i < iNonZeros ; i++)
1266         {
1267             int iRow = static_cast<int>(pRows[i]) - 1;
1268             int iCol = static_cast<int>(pCols[i]) - 1;
1269             double dblR = pValR[i];
1270             double dblI = pValI[i];
1271
1272             for (int j = 0 ; j < _pDouble->getRows() ; j++)
1273             {
1274                 double dblValR = _pDouble->get(j, iRow) * dblR - _pDouble->getImg(j, iRow) * dblI;
1275                 double dblValI = _pDouble->get(j, iRow) * dblI + _pDouble->getImg(j, iRow) * dblR;
1276                 pOut->set(j, iCol, pOut->get(j, iCol) + dblValR);
1277                 pOut->setImg(j, iCol, pOut->getImg(j, iCol) + dblValI);
1278             }
1279         }
1280     }
1281
1282     *_pOut = pOut;
1283     delete[] pRows;
1284     delete[] pValR;
1285     delete[] pValI;
1286
1287     return 0;
1288 }
1289
1290 int MultiplySparseByDouble(Sparse *_pSparse, Double*_pDouble, GenericType** _pOut)
1291 {
1292     if (_pDouble->isScalar())
1293     {
1294         //SP * d -> SP
1295         Sparse* pOut = NULL;
1296         if (_pDouble->isComplex())
1297         {
1298             std::complex<double> dbl(_pDouble->get(0), _pDouble->getImg(0));
1299             pOut = _pSparse->multiply(dbl);
1300         }
1301         else
1302         {
1303             pOut = _pSparse->multiply(_pDouble->get(0));
1304         }
1305         *_pOut = pOut;
1306         return 0;
1307     }
1308
1309     if (_pSparse->isScalar())
1310     {
1311         //D * sp -> D .* d
1312         Double* pD = NULL;
1313
1314         if (_pSparse->isComplex())
1315         {
1316             std::complex<double> dbl(_pSparse->getImg(0, 0));
1317             pD = new Double(dbl.real(), dbl.imag());
1318         }
1319         else
1320         {
1321             pD = new Double(_pSparse->get(0, 0));
1322         }
1323
1324         InternalType* pIT = GenericDotTimes(_pDouble, pD);
1325         *_pOut = pIT->getAs<GenericType>();
1326         delete pD;
1327         return 0;
1328     }
1329
1330     if (_pSparse->getCols() != _pDouble->getRows())
1331     {
1332         return 1;
1333     }
1334
1335     //try to be smart and only compute for non zero values
1336
1337     //get some information
1338     int iNonZeros = static_cast<int>(_pSparse->nonZeros());
1339     int* pRows = new int[iNonZeros * 2];
1340     _pSparse->outputRowCol(pRows);
1341     int* pCols = pRows + iNonZeros;
1342     double* pValR = new double[iNonZeros];
1343     double* pValI = new double[iNonZeros];
1344     _pSparse->outputValues(pValR, pValI);
1345
1346     Double* pOut = new Double(_pSparse->getRows(), _pDouble->getCols(), _pDouble->isComplex() | _pSparse->isComplex());
1347     pOut->setZeros();
1348
1349     if (_pDouble->isComplex() == false && _pSparse->isComplex() == false)
1350     {
1351         for (int i = 0 ; i < iNonZeros ; i++)
1352         {
1353             int iRow    = static_cast<int>(pRows[i]) - 1;
1354             int iCol    = static_cast<int>(pCols[i]) - 1;
1355             double dbl  = pValR[i];
1356
1357             for (int j = 0 ; j < _pDouble->getCols() ; j++)
1358             {
1359                 double dblVal = _pDouble->get(iCol, j) * dbl;
1360                 pOut->set(iRow, j, pOut->get(iRow, j) + dblVal);
1361             }
1362         }
1363     }
1364     else if (_pDouble->isComplex() == false && _pSparse->isComplex() == true)
1365     {
1366         //a * (b ci) -> ab ac
1367         for (int i = 0 ; i < iNonZeros ; i++)
1368         {
1369             int iRow = static_cast<int>(pRows[i]) - 1;
1370             int iCol = static_cast<int>(pCols[i]) - 1;
1371             double dblR = pValR[i];
1372             double dblI = pValI[i];
1373
1374             for (int j = 0 ; j < _pDouble->getCols() ; j++)
1375             {
1376                 double dblValR = _pDouble->get(iCol, j) * dblR;
1377                 double dblValI = _pDouble->get(iCol, j) * dblI;
1378                 pOut->set(iRow, j, pOut->get(iRow, j) + dblValR);
1379                 pOut->setImg(iRow, j, pOut->getImg(iRow, j) + dblValI);
1380             }
1381         }
1382     }
1383     else if (_pDouble->isComplex() == true && _pSparse->isComplex() == false)
1384     {
1385         //(a bi) * c -> ac + bc
1386         for (int i = 0 ; i < iNonZeros ; i++)
1387         {
1388             int iRow = static_cast<int>(pRows[i]) - 1;
1389             int iCol = static_cast<int>(pCols[i]) - 1;
1390             double dblR = pValR[i];
1391
1392             for (int j = 0 ; j < _pDouble->getCols() ; j++)
1393             {
1394                 double dblValR = _pDouble->get(iCol, j) * dblR;
1395                 double dblValI = _pDouble->getImg(iCol, j) * dblR;
1396                 pOut->set(iRow, j, pOut->get(iRow, j) + dblValR);
1397                 pOut->setImg(iRow, j, pOut->getImg(iRow, j) + dblValI);
1398             }
1399         }
1400     }
1401     else if (_pDouble->isComplex() == true && _pSparse->isComplex() == true)
1402     {
1403         for (int i = 0 ; i < iNonZeros ; i++)
1404         {
1405             int iRow = static_cast<int>(pRows[i]) - 1;
1406             int iCol = static_cast<int>(pCols[i]) - 1;
1407             double dblR = pValR[i];
1408             double dblI = pValI[i];
1409
1410             for (int j = 0 ; j < _pDouble->getCols() ; j++)
1411             {
1412                 double dblValR = _pDouble->get(iCol, j) * dblR - _pDouble->getImg(iCol, j) * dblI;
1413                 double dblValI = _pDouble->get(iCol, j) * dblI + _pDouble->getImg(iCol, j) * dblR;
1414                 pOut->set(iRow, j, pOut->get(iRow, j) + dblValR);
1415                 pOut->setImg(iRow, j, pOut->getImg(iRow, j) + dblValI);
1416             }
1417         }
1418     }
1419
1420     *_pOut = pOut;
1421     delete[] pRows;
1422     delete[] pValR;
1423     delete[] pValI;
1424
1425     return 0;
1426 }
1427
1428 int DotMultiplySparseBySparse(Sparse* _pSparse1, Sparse* _pSparse2, Sparse** _pOut)
1429 {
1430     if (_pSparse1->isScalar() || _pSparse2->isScalar())
1431     {
1432         //SP .* sp or sp .* SP
1433         return MultiplySparseBySparse(_pSparse1, _pSparse2, _pOut);
1434     }
1435
1436     if (_pSparse1->getRows() != _pSparse2->getRows() || _pSparse1->getCols() != _pSparse2->getCols())
1437     {
1438         return 1;
1439     }
1440
1441     *_pOut = _pSparse1->dotMultiply(*_pSparse2);
1442
1443     return 0;
1444 }
1445
1446 int DotMultiplyDoubleBySparse(Double* _pDouble, Sparse* _pSparse, GenericType**  _pOut)
1447 {
1448     if (_pDouble->isScalar())
1449     {
1450         return MultiplyDoubleBySparse(_pDouble, _pSparse, _pOut);
1451     }
1452
1453     if (_pSparse->isScalar())
1454     {
1455         return MultiplyDoubleBySparse(_pDouble, _pSparse, _pOut);
1456     }
1457
1458     if (_pSparse->getRows() != _pDouble->getRows() || _pSparse->getCols() != _pDouble->getCols())
1459     {
1460         return 1;
1461     }
1462
1463     Sparse* pOut = new Sparse(_pDouble->getRows(), _pDouble->getCols(), _pSparse->isComplex() || _pDouble->isComplex());
1464     //get some information
1465     int iNonZeros = static_cast<int>(_pSparse->nonZeros());
1466     int* pRows = new int[iNonZeros * 2];
1467     _pSparse->outputRowCol(pRows);
1468     int* pCols = pRows + iNonZeros;
1469
1470     if (_pDouble->isComplex() == false && _pSparse->isComplex() == false)
1471     {
1472         for (int i = 0 ; i < iNonZeros ; i++)
1473         {
1474             int iRow = static_cast<int>(pRows[i]) - 1;
1475             int iCol = static_cast<int>(pCols[i]) - 1;
1476             pOut->set(iRow, iCol, _pSparse->get(iRow, iCol) * _pDouble->get(iRow, iCol));
1477         }
1478     }
1479     else if (_pDouble->isComplex() == false && _pSparse->isComplex() == true)
1480     {
1481         for (int i = 0 ; i < iNonZeros ; i++)
1482         {
1483             int iRow = static_cast<int>(pRows[i]) - 1;
1484             int iCol = static_cast<int>(pCols[i]) - 1;
1485             std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1486             std::complex<double> newVal(dbl.real() * _pDouble->get(iRow, iCol), dbl.imag() * _pDouble->get(iRow, iCol));
1487             pOut->set(iRow, iCol, newVal);
1488         }
1489     }
1490     else if (_pDouble->isComplex() == true && _pSparse->isComplex() == false)
1491     {
1492         for (int i = 0 ; i < iNonZeros ; i++)
1493         {
1494             int iRow = static_cast<int>(pRows[i]) - 1;
1495             int iCol = static_cast<int>(pCols[i]) - 1;
1496             std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1497             std::complex<double> newVal(dbl.real() * _pDouble->get(iRow, iCol), dbl.real() * _pDouble->getImg(iRow, iCol));
1498             pOut->set(iRow, iCol, newVal);
1499         }
1500     }
1501     else if (_pDouble->isComplex() == true && _pSparse->isComplex() == true)
1502     {
1503         for (int i = 0 ; i < iNonZeros ; i++)
1504         {
1505             int iRow = static_cast<int>(pRows[i]) - 1;
1506             int iCol = static_cast<int>(pCols[i]) - 1;
1507             std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1508             double dblR = _pDouble->get(iRow, iCol) * dbl.real() - _pDouble->getImg(iRow, iCol) * dbl.imag();
1509             double dblI = _pDouble->getImg(iRow, iCol) * dbl.real() + _pDouble->get(iRow, iCol) * dbl.imag();
1510
1511             std::complex<double> newVal(dblR, dblI);
1512             pOut->set(iRow, iCol, newVal);
1513         }
1514     }
1515
1516     *_pOut = pOut;
1517     delete[] pRows;
1518
1519     return 0;
1520 }
1521
1522 int DotMultiplySparseByDouble(Sparse* _pSparse, Double* _pDouble, GenericType** _pOut)
1523 {
1524     return DotMultiplyDoubleBySparse(_pDouble, _pSparse, _pOut);
1525 }
1526
1527 int DotMultiplyPolyByDouble(Polynom* _pPoly, Double* _pDouble, Polynom** _pPolyOut)
1528 {
1529     return DotMultiplyDoubleByPoly(_pDouble, _pPoly, _pPolyOut);
1530 }
1531
1532 int DotMultiplyDoubleByPoly(Double* _pDouble, Polynom* _pPoly, Polynom** _pPolyOut)
1533 {
1534     int iSize = _pDouble->getSize();
1535     if (_pDouble->isScalar() == false &&
1536             _pPoly->isScalar() == false &&
1537             iSize != _pPoly->getSize())
1538     {
1539         return 1;
1540     }
1541
1542     int* piRanks = new int[iSize];
1543     memset(piRanks, 0x00, iSize * sizeof(int));
1544     Polynom* pPolyTemp = new Polynom(_pPoly->getVariableName(), _pDouble->getDims(), _pDouble->getDimsArray(), piRanks);
1545     delete[] piRanks;
1546     pPolyTemp->setCoef(_pDouble);
1547     int iErr = DotMultiplyPolyByPoly(pPolyTemp, _pPoly, _pPolyOut);
1548     delete pPolyTemp;
1549     return iErr;
1550 }
1551
1552 int DotMultiplyPolyByPoly(Polynom* _pPoly1, Polynom* _pPoly2, Polynom** _pPolyOut)
1553 {
1554     if (_pPoly1->isScalar() || _pPoly2->isScalar())
1555     {
1556         return MultiplyPolyByPoly(_pPoly1, _pPoly2, _pPolyOut);
1557     }
1558     else
1559     {
1560         if (_pPoly1->getSize() != _pPoly2->getSize())
1561         {
1562             return 1;
1563         }
1564
1565         int* piRank = new int[_pPoly1->getSize()];
1566         for (int i = 0 ; i < _pPoly1->getSize() ; i++)
1567         {
1568             piRank[i] = _pPoly1->get(i)->getRank() + _pPoly2->get(i)->getRank();
1569         }
1570
1571         (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), _pPoly1->getDims(), _pPoly1->getDimsArray(), piRank);
1572
1573         if (_pPoly1->isComplex() && _pPoly2->isComplex())
1574         {
1575             (*_pPolyOut)->setComplex(true);
1576             for (int i = 0; i < _pPoly1->getSize(); i++)
1577             {
1578                 SinglePoly *pSP1    = _pPoly1->get(i);
1579                 SinglePoly *pSP2    = _pPoly2->get(i);
1580                 SinglePoly *pSPOut  = (*_pPolyOut)->get(i);
1581
1582                 pSPOut->setZeros();
1583
1584                 iMultiComplexPolyByComplexPoly(
1585                     pSP1->get(), pSP1->getImg(), pSP1->getSize(),
1586                     pSP2->get(), pSP2->getImg(), pSP2->getSize(),
1587                     pSPOut->get(), pSPOut->getImg(), pSPOut->getSize());
1588
1589             }
1590         }
1591         else if (_pPoly1->isComplex())
1592         {
1593             (*_pPolyOut)->setComplex(true);
1594             for (int i = 0; i < _pPoly1->getSize(); i++)
1595             {
1596                 SinglePoly *pSP1   = _pPoly1->get(i);
1597                 SinglePoly *pSP2   = _pPoly2->get(i);
1598                 SinglePoly *pSPOut = (*_pPolyOut)->get(i);
1599
1600                 pSPOut->setZeros();
1601
1602                 iMultiComplexPolyByScilabPolynom(
1603                     pSP1->get(), pSP1->getImg(), pSP1->getSize(),
1604                     pSP2->get(), pSP2->getSize(),
1605                     pSPOut->get(), pSPOut->getImg(), pSPOut->getSize());
1606             }
1607         }
1608         else if (_pPoly2->isComplex())
1609         {
1610             (*_pPolyOut)->setComplex(true);
1611             for (int i = 0; i < _pPoly1->getSize(); i++)
1612             {
1613                 SinglePoly *pSP1   = _pPoly1->get(i);
1614                 SinglePoly *pSP2   = _pPoly2->get(i);
1615                 SinglePoly *pSPOut = (*_pPolyOut)->get(i);
1616
1617                 pSPOut->setZeros();
1618
1619                 iMultiScilabPolynomByComplexPoly(
1620                     pSP1->get(), pSP1->getSize(),
1621                     pSP2->get(), pSP2->getImg(), pSP2->getSize(),
1622                     pSPOut->get(), pSPOut->getImg(), pSPOut->getSize());
1623             }
1624         }
1625         else
1626         {
1627             for (int i = 0; i < _pPoly1->getSize(); i++)
1628             {
1629                 SinglePoly *pSP1   = _pPoly1->get(i);
1630                 SinglePoly *pSP2   = _pPoly2->get(i);
1631                 SinglePoly *pSPOut = (*_pPolyOut)->get(i);
1632
1633                 pSPOut->setZeros();
1634
1635                 iMultiScilabPolynomByScilabPolynom(
1636                     pSP1->get(), pSP1->getSize(),
1637                     pSP2->get(), pSP2->getSize(),
1638                     pSPOut->get(), pSPOut->getSize());
1639             }
1640         }
1641     }
1642
1643     return 0;
1644 }
1645