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