reforge dot multiplication
[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     bool bScalar1   = _pDouble->isScalar();
378     bool bScalar2   = _pPoly->isScalar();
379
380     if (_pDouble->isScalar())
381     {
382         int* piRank = new int[_pPoly->getSize()];
383         for (int i = 0 ; i < _pPoly->getSize() ; i++)
384         {
385             piRank[i] = _pPoly->get(i)->getRank();
386         }
387
388         (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pPoly->getDims(), _pPoly->getDimsArray(), piRank);
389         delete[] piRank;
390         if (bComplex1 || bComplex2)
391         {
392             (*_pPolyOut)->setComplex(true);
393         }
394
395         for (int i = 0 ; i < _pPoly->getSize() ; i++)
396         {
397             SinglePoly *pPolyIn     = _pPoly->get(i);
398             double* pRealIn         = pPolyIn->get();
399             double* pImgIn          = pPolyIn->getImg();
400
401             SinglePoly *pPolyOut    = (*_pPolyOut)->get(i);
402             double* pRealOut        = pPolyOut->get();
403             double* pImgOut         = pPolyOut->getImg();
404
405             if (bComplex1 == false && bComplex2 == false)
406             {
407                 iMultiRealScalarByRealMatrix(_pDouble->get(0), pRealIn, 1, pPolyIn->getSize(), pRealOut);
408             }
409             else if (bComplex1 == false && bComplex2 == true)
410             {
411                 iMultiRealScalarByComplexMatrix(_pDouble->get(0), pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
412             }
413             else if (bComplex1 == true && bComplex2 == false)
414             {
415                 iMultiComplexScalarByRealMatrix(_pDouble->get(0), _pDouble->getImg(0), pRealIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
416             }
417             else if (bComplex1 == true && bComplex2 == true)
418             {
419                 iMultiComplexScalarByComplexMatrix(_pDouble->get(0), _pDouble->getImg(0), pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
420             }
421         }
422         (*_pPolyOut)->updateRank();
423         return 0;
424     }
425
426     if (_pPoly->isScalar())
427     {
428         int* piRank = new int[_pDouble->getSize()];
429         for (int i = 0 ; i < _pDouble->getSize() ; i++)
430         {
431             piRank[i] = _pPoly->get(0)->getRank();
432         }
433
434         (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pDouble->getDims(), _pDouble->getDimsArray(), piRank);
435         delete[] piRank;
436         if (bComplex1 || bComplex2)
437         {
438             (*_pPolyOut)->setComplex(true);
439         }
440
441         double *pDoubleR    = _pDouble->get();
442         double *pDoubleI    = _pDouble->getImg();
443
444         SinglePoly *pPolyIn = _pPoly->get(0);
445         double* pRealIn     = pPolyIn->get();
446         double* pImgIn      = pPolyIn->getImg();
447
448         for (int i = 0 ; i < _pDouble->getSize() ; i++)
449         {
450             SinglePoly *pPolyOut    = (*_pPolyOut)->get(i);
451             double* pRealOut        = pPolyOut->get();
452             double* pImgOut         = pPolyOut->getImg();
453
454             if (bComplex1 == false && bComplex2 == false)
455             {
456                 iMultiRealScalarByRealMatrix(pDoubleR[i], pRealIn, 1, pPolyIn->getSize(), pRealOut);
457             }
458             else if (bComplex1 == false && bComplex2 == true)
459             {
460                 iMultiRealScalarByComplexMatrix(pDoubleR[i], pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
461             }
462             else if (bComplex1 == true && bComplex2 == false)
463             {
464                 iMultiComplexScalarByRealMatrix(pDoubleR[i], pDoubleI[i], pRealIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
465             }
466             else if (bComplex1 == true && bComplex2 == true)
467             {
468                 iMultiComplexScalarByComplexMatrix(pDoubleR[i], pDoubleI[i], pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
469             }
470         }
471
472         (*_pPolyOut)->updateRank();
473         return 0;
474     }
475
476     if (_pPoly->getDims() > 2 || _pDouble->getDims() > 2 || _pDouble->getCols() != _pPoly->getRows())
477     {
478         return 1;
479     }
480
481     int* piRank = new int[_pDouble->getRows() * _pPoly->getCols()];
482     int iMaxRank = _pPoly->getMaxRank();
483     for (int i = 0 ; i < _pDouble->getRows() * _pPoly->getCols() ; i++)
484     {
485         piRank[i] = iMaxRank;
486     }
487
488     (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pDouble->getRows(), _pPoly->getCols(), piRank);
489     delete[] piRank;
490     if (bComplex1 || bComplex2)
491     {
492         (*_pPolyOut)->setComplex(true);
493     }
494
495     Double *pCoef = _pPoly->getCoef();
496     Double *pTemp = new Double(_pDouble->getRows(), pCoef->getCols(), bComplex1 || bComplex2);
497
498     if (bComplex1 == false && bComplex2 == false)
499     {
500         iMultiRealMatrixByRealMatrix(_pDouble->get(), _pDouble->getRows(), _pDouble->getCols(),
501                                      pCoef->get(), pCoef->getRows(), pCoef->getCols(),
502                                      pTemp->get());
503     }
504     else if (bComplex1 == false && bComplex2 == true)
505     {
506         iMultiRealMatrixByComplexMatrix(_pDouble->get(), _pDouble->getRows(), _pDouble->getCols(),
507                                         pCoef->get(), pCoef->getImg(), pCoef->getRows(), pCoef->getCols(),
508                                         pTemp->get(), pTemp->getImg());
509
510     }
511     else if (bComplex1 == true && bComplex2 == false)
512     {
513         iMultiComplexMatrixByRealMatrix(_pDouble->get(), _pDouble->getImg(), _pDouble->getRows(), _pDouble->getCols(),
514                                         pCoef->get(), pCoef->getRows(), pCoef->getCols(),
515                                         pTemp->get(), pTemp->getImg());
516     }
517     else //if(bComplex1 == true && bComplex2 == true)
518     {
519         iMultiComplexMatrixByComplexMatrix(_pDouble->get(), _pDouble->getImg(), _pDouble->getRows(), _pDouble->getCols(),
520                                            pCoef->get(), pCoef->getImg(), pCoef->getRows(), pCoef->getCols(),
521                                            pTemp->get(), pTemp->getImg());
522     }
523
524     (*_pPolyOut)->setCoef(pTemp);
525     (*_pPolyOut)->updateRank();
526     delete pTemp;
527     return 0;
528 }
529
530 int MultiplyPolyByDouble(Polynom* _pPoly, Double* _pDouble, Polynom **_pPolyOut)
531 {
532     bool bComplex1  = _pPoly->isComplex();
533     bool bComplex2  = _pDouble->isComplex();
534     bool bScalar1   = _pPoly->isScalar();
535     bool bScalar2   = _pDouble->isScalar();
536
537     if (bScalar1)
538     {
539         int* piRank = new int[_pDouble->getSize()];
540         for (int i = 0 ; i < _pDouble->getSize() ; i++)
541         {
542             piRank[i] = _pPoly->get(0)->getRank();
543         }
544
545         (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pDouble->getDims(), _pDouble->getDimsArray(), piRank);
546         delete[] piRank;
547         if (bComplex1 || bComplex2)
548         {
549             (*_pPolyOut)->setComplex(true);
550         }
551
552         double *pDoubleR    = _pDouble->get();
553         double *pDoubleI    = _pDouble->getImg();
554
555         SinglePoly *pPolyIn = _pPoly->get(0);
556         double* pRealIn     = pPolyIn->get();
557         double* pImgIn      = pPolyIn->getImg();
558
559         for (int i = 0 ; i < _pDouble->getSize() ; i++)
560         {
561             SinglePoly *pPolyOut    = (*_pPolyOut)->get(i);
562             double* pRealOut        = pPolyOut->get();
563             double* pImgOut         = pPolyOut->getImg();
564
565             if (bComplex1 == false && bComplex2 == false)
566             {
567                 iMultiRealScalarByRealMatrix(pDoubleR[i], pRealIn, 1, pPolyIn->getSize(), pRealOut);
568             }
569             else if (bComplex1 == false && bComplex2 == true)
570             {
571                 iMultiComplexScalarByRealMatrix(pDoubleR[i], pDoubleI[i], pRealIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
572             }
573             else if (bComplex1 == true && bComplex2 == false)
574             {
575                 iMultiRealScalarByComplexMatrix(pDoubleR[i], pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
576             }
577             else if (bComplex1 == true && bComplex2 == true)
578             {
579                 iMultiComplexScalarByComplexMatrix(pDoubleR[i], pDoubleI[i], pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
580             }
581         }
582
583         (*_pPolyOut)->updateRank();
584         return 0;
585     }
586     else if (bScalar2)
587     {
588         int* piRank = new int[_pPoly->getSize()];
589         for (int i = 0 ; i < _pPoly->getSize() ; i++)
590         {
591             piRank[i] = _pPoly->get(i)->getRank();
592         }
593
594         (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pPoly->getDims(), _pPoly->getDimsArray(), piRank);
595         delete[] piRank;
596         if (bComplex1 || bComplex2)
597         {
598             (*_pPolyOut)->setComplex(true);
599         }
600
601         for (int i = 0 ; i < _pPoly->getSize() ; i++)
602         {
603             SinglePoly *pPolyIn = _pPoly->get(i);
604             double* pRealIn     = pPolyIn->get();
605             double* pImgIn      = pPolyIn->getImg();
606
607             SinglePoly *pPolyOut    = (*_pPolyOut)->get(i);
608             double* pRealOut        = pPolyOut->get();
609             double* pImgOut         = pPolyOut->getImg();
610
611             if (bComplex1 == false && bComplex2 == false)
612             {
613                 iMultiRealScalarByRealMatrix(_pDouble->get(0), pRealIn, 1, pPolyIn->getSize(), pRealOut);
614             }
615             else if (bComplex1 == false && bComplex2 == true)
616             {
617                 iMultiComplexScalarByRealMatrix(_pDouble->get(0), _pDouble->getImg(0), pRealIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
618             }
619             else if (bComplex1 == true && bComplex2 == false)
620             {
621                 iMultiRealScalarByComplexMatrix(_pDouble->get(0), pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
622             }
623             else if (bComplex1 == true && bComplex2 == true)
624             {
625                 iMultiComplexScalarByComplexMatrix(_pDouble->get(0), _pDouble->getImg(0), pRealIn, pImgIn, 1, pPolyIn->getSize(), pRealOut, pImgOut);
626             }
627         }
628
629         (*_pPolyOut)->updateRank();
630         return 0;
631     }
632
633     if (_pDouble->getDims() > 2 || _pPoly->getDims() > 2 || _pPoly->getCols() != _pDouble->getRows())
634     {
635         return 1;
636     }
637
638     int* piRank = new int[_pPoly->getRows() * _pDouble->getCols()];
639     int iMaxRank = _pPoly->getMaxRank();
640     for (int i = 0 ; i < _pPoly->getRows() * _pDouble->getCols() ; i++)
641     {
642         piRank[i] = iMaxRank;
643     }
644
645     (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pPoly->getRows(), _pDouble->getCols(), piRank);
646     delete[] piRank;
647     if (bComplex1 || bComplex2)
648     {
649         (*_pPolyOut)->setComplex(true);
650     }
651
652     //Distribution a la mano par appels a des sous-fonctions ( iMulti...ScalarBy...Scalar ) plus iAdd...To... )
653
654     //for each line of _pPoly
655     for (int iRow1 = 0 ; iRow1 < _pPoly->getRows() ; iRow1++)
656     {
657         //for each col of _pDouble
658         for (int iCol2 = 0 ; iCol2 < _pDouble->getCols() ; iCol2++)
659         {
660             SinglePoly* pSPOut = (*_pPolyOut)->get(iRow1, iCol2);
661             pSPOut->setZeros();
662
663             //for each rows of _pDouble / cols of _pPoly
664             for (int iRow2 = 0 ; iRow2 < _pDouble->getRows() ; iRow2++)
665             {
666                 // SinglePoly(iRow1, iRow2) * Double(iRow2, iCol2)
667                 SinglePoly* pSPIn = _pPoly->get(iRow1, iRow2);
668                 int iSize = pSPIn->getSize();
669                 double* pdblMult = new double[iSize];
670
671                 if (bComplex1 == false && bComplex2 == false)
672                 {
673                     //Real Matrix by Real Scalar
674                     iMultiRealScalarByRealMatrix(_pDouble->get(iRow2, iCol2), pSPIn->get(), iSize, 1, pdblMult);
675                     add(pSPOut->get(), (long long)iSize, pdblMult, pSPOut->get());
676                 }
677                 else if (bComplex1 == false && bComplex2 == true)
678                 {
679                     //Real Matrix by Scalar Complex
680                     double* pdblMultImg = new double[iSize];
681                     iMultiComplexScalarByRealMatrix(_pDouble->get(iRow2, iCol2), _pDouble->getImg(iRow2, iCol2), pSPIn->get(), pSPIn->getSize(), 1, pdblMult, pdblMultImg);
682                     add(pSPOut->get(), pSPOut->getImg(), (long long)iSize, pdblMult, pdblMultImg, pSPOut->get(), pSPOut->getImg());
683                     delete[] pdblMultImg;
684                 }
685                 else if (bComplex1 == true && bComplex2 == false)
686                 {
687                     double* pdblMultImg = new double[iSize];
688                     iMultiRealScalarByComplexMatrix(_pDouble->get(iRow2, iCol2), pSPIn->get(), pSPIn->getImg(), pSPIn->getSize(), 1, pdblMult, pdblMultImg);
689                     add(pSPOut->get(), pSPOut->getImg(), (long long)iSize, pdblMult, pdblMultImg, pSPOut->get(), pSPOut->getImg());
690                     delete[] pdblMultImg;
691                 }
692                 else //if(bComplex1 == true && bComplex2 == true)
693                 {
694                     double* pdblMultImg = new double[iSize];
695                     iMultiComplexScalarByComplexMatrix(_pDouble->get(iRow2, iCol2), _pDouble->getImg(iRow2, iCol2), pSPIn->get(), pSPIn->getImg(), pSPIn->getSize(), 1, pdblMult, pdblMultImg);
696                     add(pSPOut->get(), pSPOut->getImg(), (long long)iSize, pdblMult, pdblMultImg, pSPOut->get(), pSPOut->getImg());
697                     delete[] pdblMultImg;
698                 }
699
700                 delete[] pdblMult;
701             }//for(int iRow2 = 0 ; iRow2 < _pDouble->getRows() ; iRow2++)
702         }//for(int iCol2 = 0 ; iCol2 < _pDouble->getCols() ; iCol2++)
703     }//for(int iRow1 = 0 ; iRow1 < _pPoly->getRows() ; iRow1++)
704
705     (*_pPolyOut)->updateRank();
706     return 0;
707 }
708
709 int MultiplyPolyByPoly(Polynom* _pPoly1, Polynom* _pPoly2, Polynom** _pPolyOut)
710 {
711     bool bComplex1  = _pPoly1->isComplex();
712     bool bComplex2  = _pPoly2->isComplex();
713
714     if (_pPoly1->isScalar() && _pPoly2->isScalar())
715     {
716         //poly1(0) * poly2(0)
717         int iRank = _pPoly1->get(0)->getRank() + _pPoly2->get(0)->getRank();
718         (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), 1, 1, &iRank);
719         if (bComplex1 || bComplex2)
720         {
721             (*_pPolyOut)->setComplex(true);
722         }
723
724         if (bComplex1 == false && bComplex2 == false)
725         {
726             SinglePoly *pPoly1  = _pPoly1->get(0);
727             SinglePoly *pPoly2  = _pPoly2->get(0);
728             SinglePoly *pPolyOut = (*_pPolyOut)->get(0);
729
730             pPolyOut->setZeros();
731
732             iMultiScilabPolynomByScilabPolynom(
733                 pPoly1->get(), pPoly1->getSize(),
734                 pPoly2->get(), pPoly2->getSize(),
735                 pPolyOut->get(), pPolyOut->getSize());
736         }
737         else if (bComplex1 == false && bComplex2 == true)
738         {
739             SinglePoly *pPoly1  = _pPoly1->get(0);
740             SinglePoly *pPoly2  = _pPoly2->get(0);
741             SinglePoly *pPolyOut = (*_pPolyOut)->get(0);
742
743             pPolyOut->setZeros();
744
745             iMultiScilabPolynomByComplexPoly(
746                 pPoly1->get(), pPoly1->getSize(),
747                 pPoly2->get(), pPoly2->getImg(), pPoly2->getSize(),
748                 pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
749         }
750         else if (bComplex1 == true && bComplex2 == false)
751         {
752             SinglePoly *pPoly1  = _pPoly1->get(0);
753             SinglePoly *pPoly2  = _pPoly2->get(0);
754             SinglePoly *pPolyOut = (*_pPolyOut)->get(0);
755
756             pPolyOut->setZeros();
757
758             iMultiComplexPolyByScilabPolynom(
759                 pPoly1->get(), pPoly1->getImg(), pPoly1->getSize(),
760                 pPoly2->get(), pPoly2->getSize(),
761                 pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
762         }
763         else if (bComplex1 == true && bComplex2 == true)
764         {
765             SinglePoly *pPoly1   = _pPoly1->get(0);
766             SinglePoly *pPoly2   = _pPoly2->get(0);
767             SinglePoly *pPolyOut  = (*_pPolyOut)->get(0);
768
769             pPolyOut->setZeros();
770
771             iMultiComplexPolyByComplexPoly(
772                 pPoly1->get(), pPoly1->getImg(), pPoly1->getSize(),
773                 pPoly2->get(), pPoly2->getImg(), pPoly2->getSize(),
774                 pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
775         }
776
777         (*_pPolyOut)->updateRank();
778         return 0;
779     }
780
781     if (_pPoly1->isScalar())
782     {
783         //poly1(0) * poly2(n)
784         int* piRank = new int[_pPoly2->getSize()];
785         for (int i = 0 ; i < _pPoly2->getSize() ; i++)
786         {
787             piRank[i] = _pPoly1->get(0)->getRank() + _pPoly2->get(i)->getRank();
788         }
789
790         (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), _pPoly2->getDims(), _pPoly2->getDimsArray(), piRank);
791         if (bComplex1 || bComplex2)
792         {
793             (*_pPolyOut)->setComplex(true);
794         }
795         delete[] piRank;
796
797
798         SinglePoly *pPoly1  = _pPoly1->get(0);
799         if (bComplex1 == false && bComplex2 == false)
800         {
801             for (int iPoly = 0 ; iPoly < _pPoly2->getSize() ; iPoly++)
802             {
803                 SinglePoly *pPoly2  = _pPoly2->get(iPoly);
804                 SinglePoly *pPolyOut = (*_pPolyOut)->get(iPoly);
805
806                 pPolyOut->setZeros();
807
808                 iMultiScilabPolynomByScilabPolynom(
809                     pPoly1->get(), pPoly1->getSize(),
810                     pPoly2->get(), pPoly2->getSize(),
811                     pPolyOut->get(), pPolyOut->getSize());
812             }
813         }
814         else if (bComplex1 == false && bComplex2 == true)
815         {
816             for (int iPoly = 0 ; iPoly < _pPoly2->getSize() ; iPoly++)
817             {
818                 SinglePoly *pPoly2  = _pPoly2->get(iPoly);
819                 SinglePoly *pPolyOut = (*_pPolyOut)->get(iPoly);
820
821                 pPolyOut->setZeros();
822
823                 iMultiScilabPolynomByComplexPoly(
824                     pPoly1->get(), pPoly1->getSize(),
825                     pPoly2->get(), pPoly2->getImg(), pPoly2->getSize(),
826                     pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
827             }
828         }
829         else if (bComplex1 == true && bComplex2 == false)
830         {
831             for (int iPoly = 0 ; iPoly < _pPoly2->getSize() ; iPoly++)
832             {
833                 SinglePoly *pPoly2  = _pPoly2->get(iPoly);
834                 SinglePoly *pPolyOut = (*_pPolyOut)->get(iPoly);
835
836                 pPolyOut->setZeros();
837
838                 iMultiComplexPolyByScilabPolynom(
839                     pPoly1->get(), pPoly1->getImg(), pPoly1->getSize(),
840                     pPoly2->get(), pPoly2->getSize(),
841                     pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
842             }
843         }
844         else if (bComplex1 == true && bComplex2 == true)
845         {
846             for (int iPoly = 0 ; iPoly < _pPoly2->getSize() ; iPoly++)
847             {
848                 SinglePoly *pPoly2   = _pPoly2->get(iPoly);
849                 SinglePoly *pPolyOut  = (*_pPolyOut)->get(iPoly);
850
851                 pPolyOut->setZeros();
852
853                 iMultiComplexPolyByComplexPoly(
854                     pPoly1->get(), pPoly1->getImg(), pPoly1->getSize(),
855                     pPoly2->get(), pPoly2->getImg(), pPoly2->getSize(),
856                     pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
857             }
858         }
859
860         (*_pPolyOut)->updateRank();
861         return 0;
862     }
863
864     if (_pPoly2->isScalar())
865     {
866         //poly1(n) * poly2(0)
867         int* piRank = new int[_pPoly1->getSize()];
868         for (int i = 0 ; i < _pPoly1->getSize() ; i++)
869         {
870             piRank[i] = _pPoly2->get(0)->getRank() + _pPoly1->get(i)->getRank();
871         }
872
873         (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), _pPoly1->getDims(), _pPoly1->getDimsArray(), piRank);
874         if (bComplex1 || bComplex2)
875         {
876             (*_pPolyOut)->setComplex(true);
877         }
878         delete[] piRank;
879
880         SinglePoly *pPoly2  = _pPoly2->get(0);
881         if (bComplex1 == false && bComplex2 == false)
882         {
883             for (int iPoly = 0 ; iPoly < _pPoly1->getSize() ; iPoly++)
884             {
885                 SinglePoly *pPoly1  = _pPoly1->get(iPoly);
886                 SinglePoly *pPolyOut = (*_pPolyOut)->get(iPoly);
887
888                 pPolyOut->setZeros();
889
890                 iMultiScilabPolynomByScilabPolynom(
891                     pPoly1->get(), pPoly1->getSize(),
892                     pPoly2->get(), pPoly2->getSize(),
893                     pPolyOut->get(), pPolyOut->getSize());
894             }
895         }
896         else if (bComplex1 == false && bComplex2 == true)
897         {
898             for (int iPoly = 0 ; iPoly < _pPoly1->getSize() ; iPoly++)
899             {
900                 SinglePoly *pPoly1  = _pPoly1->get(iPoly);
901                 SinglePoly *pPolyOut = (*_pPolyOut)->get(iPoly);
902
903                 pPolyOut->setZeros();
904
905                 iMultiScilabPolynomByComplexPoly(
906                     pPoly1->get(), pPoly1->getSize(),
907                     pPoly2->get(), pPoly2->getImg(), pPoly2->getSize(),
908                     pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
909             }
910         }
911         else if (bComplex1 == true && bComplex2 == false)
912         {
913             for (int iPoly = 0 ; iPoly < _pPoly1->getSize() ; iPoly++)
914             {
915                 SinglePoly *pPoly1  = _pPoly1->get(iPoly);
916                 SinglePoly *pPolyOut = (*_pPolyOut)->get(iPoly);
917
918                 pPolyOut->setZeros();
919
920                 iMultiComplexPolyByScilabPolynom(
921                     pPoly1->get(), pPoly1->getImg(), pPoly1->getSize(),
922                     pPoly2->get(), pPoly2->getSize(),
923                     pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
924             }
925         }
926         else if (bComplex1 == true && bComplex2 == true)
927         {
928             for (int iPoly = 0 ; iPoly < _pPoly1->getSize() ; iPoly++)
929             {
930                 SinglePoly *pPoly1   = _pPoly1->get(iPoly);
931                 SinglePoly *pPolyOut  = (*_pPolyOut)->get(iPoly);
932
933                 pPolyOut->setZeros();
934
935                 iMultiComplexPolyByComplexPoly(
936                     pPoly1->get(), pPoly1->getImg(), pPoly1->getSize(),
937                     pPoly2->get(), pPoly2->getImg(), pPoly2->getSize(),
938                     pPolyOut->get(), pPolyOut->getImg(), pPolyOut->getSize());
939             }
940         }
941
942         (*_pPolyOut)->updateRank();
943         return 0;
944     }
945
946     if (_pPoly1->getDims() > 2 || _pPoly2->getDims() > 2 || _pPoly1->getCols() != _pPoly2->getRows())
947     {
948         return 1;
949     }
950
951     // matrix by matrix
952     int* piRank = new int[_pPoly1->getRows() * _pPoly2->getCols()];
953     int iMaxRank = _pPoly1->getMaxRank() + _pPoly2->getMaxRank();
954     for (int i = 0 ; i < _pPoly1->getRows() * _pPoly2->getCols() ; i++)
955     {
956         piRank[i] = iMaxRank;
957     }
958
959     (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), _pPoly1->getRows(), _pPoly2->getCols(), piRank);
960     if (bComplex1 || bComplex2)
961     {
962         (*_pPolyOut)->setComplex(true);
963     }
964
965     delete[] piRank;
966
967
968     if (bComplex1 == false && bComplex2 == false)
969     {
970         double *pReal = NULL;
971         SinglePoly *pTemp  = new SinglePoly(&pReal, (*_pPolyOut)->getMaxRank());
972
973         for (int iRow = 0 ; iRow < _pPoly1->getRows() ; iRow++)
974         {
975             for (int iCol = 0 ; iCol < _pPoly2->getCols() ; iCol++)
976             {
977                 SinglePoly *pResult = (*_pPolyOut)->get(iRow, iCol);
978                 pResult->setZeros();
979
980                 for (int iCommon = 0 ; iCommon < _pPoly1->getCols() ; iCommon++)
981                 {
982                     SinglePoly *pL   = _pPoly1->get(iRow, iCommon);
983                     SinglePoly *pR   = _pPoly2->get(iCommon, iCol);
984
985                     pTemp->setZeros();
986
987                     iMultiScilabPolynomByScilabPolynom(
988                         pL->get(), pL->getSize(),
989                         pR->get(), pR->getSize(),
990                         pTemp->get(), pL->getRank() + pR->getRank() + 1);
991
992                     iAddScilabPolynomToScilabPolynom(
993                         pResult->get(), pResult->getSize(),
994                         pTemp->get(), pResult->getSize(),
995                         pResult->get(), pResult->getSize());
996                 }
997             }
998         }
999     }
1000     else if (bComplex1 == false && bComplex2 == true)
1001     {
1002         double *pReal = NULL;
1003         double *pImg = NULL;
1004         SinglePoly *pTemp  = new SinglePoly(&pReal, &pImg, (*_pPolyOut)->getMaxRank());
1005
1006         for (int iRow = 0 ; iRow < _pPoly1->getRows() ; iRow++)
1007         {
1008             for (int iCol = 0 ; iCol < _pPoly2->getCols() ; iCol++)
1009             {
1010                 SinglePoly *pResult = (*_pPolyOut)->get(iRow, iCol);
1011                 pResult->setZeros();
1012
1013                 for (int iCommon = 0 ; iCommon < _pPoly1->getCols() ; iCommon++)
1014                 {
1015                     SinglePoly *pL   = _pPoly1->get(iRow, iCommon);
1016                     SinglePoly *pR   = _pPoly2->get(iCommon, iCol);
1017
1018                     pTemp->setZeros();
1019
1020                     iMultiScilabPolynomByComplexPoly(
1021                         pL->get(), pL->getSize(),
1022                         pR->get(), pR->getImg(), pR->getSize(),
1023                         pTemp->get(), pTemp->getImg(), pL->getRank() + pR->getRank() + 1);
1024
1025                     iAddComplexPolyToComplexPoly(
1026                         pResult->get(), pResult->getImg(), pResult->getSize(),
1027                         pTemp->get(), pTemp->getImg(), pResult->getSize(),
1028                         pResult->get(), pResult->getImg(), pResult->getSize());
1029                 }
1030             }
1031         }
1032     }
1033     else if (bComplex1 == true && bComplex2 == false)
1034     {
1035         double *pReal = NULL;
1036         double *pImg = NULL;
1037         SinglePoly *pTemp  = new SinglePoly(&pReal, &pImg, (*_pPolyOut)->getMaxRank());
1038
1039         for (int iRow = 0 ; iRow < _pPoly1->getRows() ; iRow++)
1040         {
1041             for (int iCol = 0 ; iCol < _pPoly2->getCols() ; iCol++)
1042             {
1043                 SinglePoly *pResult = (*_pPolyOut)->get(iRow, iCol);
1044                 pResult->setZeros();
1045
1046                 for (int iCommon = 0 ; iCommon < _pPoly1->getCols() ; iCommon++)
1047                 {
1048                     SinglePoly *pL   = _pPoly1->get(iRow, iCommon);
1049                     SinglePoly *pR   = _pPoly2->get(iCommon, iCol);
1050
1051                     pTemp->setZeros();
1052
1053                     iMultiScilabPolynomByComplexPoly(
1054                         pR->get(), pR->getSize(),
1055                         pL->get(), pL->getImg(), pL->getSize(),
1056                         pTemp->get(), pTemp->getImg(), pL->getRank() + pR->getRank() + 1);
1057
1058                     iAddComplexPolyToComplexPoly(
1059                         pResult->get(), pResult->getImg(), pResult->getSize(),
1060                         pTemp->get(), pTemp->getImg(), pResult->getSize(),
1061                         pResult->get(), pResult->getImg(), pResult->getSize());
1062                 }
1063             }
1064         }
1065     }
1066     else if (bComplex1 == true && bComplex2 == true)
1067     {
1068         double *pReal = NULL;
1069         double *pImg = NULL;
1070         SinglePoly *pTemp  = new SinglePoly(&pReal, &pImg, (*_pPolyOut)->getMaxRank());
1071
1072         for (int iRow = 0 ; iRow < _pPoly1->getRows() ; iRow++)
1073         {
1074             for (int iCol = 0 ; iCol < _pPoly2->getCols() ; iCol++)
1075             {
1076                 SinglePoly *pResult = (*_pPolyOut)->get(iRow, iCol);
1077                 pResult->setZeros();
1078
1079                 for (int iCommon = 0 ; iCommon < _pPoly1->getCols() ; iCommon++)
1080                 {
1081                     SinglePoly *pL   = _pPoly1->get(iRow, iCommon);
1082                     SinglePoly *pR   = _pPoly2->get(iCommon, iCol);
1083
1084                     pTemp->setZeros();
1085
1086                     iMultiComplexPolyByComplexPoly(
1087                         pL->get(), pL->getImg(), pL->getSize(),
1088                         pR->get(), pR->getImg(), pR->getSize(),
1089                         pTemp->get(), pTemp->getImg(), pL->getRank() + pR->getRank() + 1);
1090
1091                     iAddComplexPolyToComplexPoly(
1092                         pResult->get(), pResult->getImg(), pResult->getSize(),
1093                         pTemp->get(), pTemp->getImg(), pResult->getSize(),
1094                         pResult->get(), pResult->getImg(), pResult->getSize());
1095                 }
1096             }
1097         }
1098     }
1099     (*_pPolyOut)->updateRank();
1100
1101     return 0;
1102 }
1103
1104 int MultiplySparseBySparse(Sparse* _pSparse1, Sparse* _pSparse2, Sparse** _pSparseOut)
1105 {
1106     if (_pSparse1->isScalar())
1107     {
1108         //scalar * sp
1109         Double* pDbl = NULL;
1110         if (_pSparse1->isComplex())
1111         {
1112             std::complex<double> dbl = _pSparse1->getImg(0, 0);
1113             pDbl = new Double(dbl.real(), dbl.imag());
1114         }
1115         else
1116         {
1117             pDbl = new Double(_pSparse1->get(0, 0));
1118         }
1119
1120         MultiplyDoubleBySparse(pDbl, _pSparse2, (GenericType**)_pSparseOut);
1121         delete pDbl;
1122         return 0;
1123     }
1124
1125     if (_pSparse2->isScalar())
1126     {
1127         //sp * scalar
1128         Double* pDbl = NULL;
1129         if (_pSparse2->isComplex())
1130         {
1131             std::complex<double> dbl = _pSparse2->getImg(0, 0);
1132             pDbl = new Double(dbl.real(), dbl.imag());
1133         }
1134         else
1135         {
1136             pDbl = new Double(_pSparse2->get(0, 0));
1137         }
1138
1139         MultiplySparseByDouble(_pSparse1, pDbl, (GenericType**)_pSparseOut);
1140         delete pDbl;
1141         return 0;
1142     }
1143
1144     if (_pSparse1->getCols() != _pSparse2->getRows())
1145     {
1146         return 1;
1147     }
1148
1149     *_pSparseOut = _pSparse1->multiply(*_pSparse2);
1150     return 0;
1151 }
1152
1153 int MultiplyDoubleBySparse(Double* _pDouble, Sparse *_pSparse, GenericType** _pOut)
1154 {
1155     //D * SP
1156     if (_pDouble->isScalar())
1157     {
1158         //d * SP -> SP
1159         Sparse* pOut = NULL;
1160         if (_pDouble->isComplex())
1161         {
1162             std::complex<double> dbl(_pDouble->get(0), _pDouble->getImg(0));
1163             pOut = _pSparse->multiply(dbl);
1164         }
1165         else
1166         {
1167             pOut = _pSparse->multiply(_pDouble->get(0));
1168         }
1169         *_pOut = pOut;
1170         return 0;
1171     }
1172
1173     if (_pSparse->isScalar())
1174     {
1175         //D * sp -> D .* d
1176         Double* pD = NULL;
1177
1178         if (_pSparse->isComplex())
1179         {
1180             std::complex<double> dbl(_pSparse->getImg(0, 0));
1181             pD = new Double(dbl.real(), dbl.imag());
1182         }
1183         else
1184         {
1185             pD = new Double(_pSparse->get(0, 0));
1186         }
1187
1188         InternalType* pIT = GenericDotTimes(_pDouble, pD);
1189         *_pOut = pIT->getAs<GenericType>();
1190         delete pD;
1191         return 0;
1192     }
1193
1194     if (_pDouble->getCols() != _pSparse->getRows())
1195     {
1196         return 1;
1197     }
1198
1199     //try to be smart and only compute for non zero values
1200
1201     //get some information
1202     int iNonZeros = static_cast<int>(_pSparse->nonZeros());
1203     int* pRows = new int[iNonZeros * 2];
1204     _pSparse->outputRowCol(pRows);
1205     int* pCols = pRows + iNonZeros;
1206     double* pValR = new double[iNonZeros];
1207     double* pValI = new double[iNonZeros];
1208     _pSparse->outputValues(pValR, pValI);
1209
1210     Double* pOut = new Double(_pDouble->getRows(), _pSparse->getCols(), _pDouble->isComplex() | _pSparse->isComplex());
1211     pOut->setZeros();
1212
1213     if (_pDouble->isComplex() == false && _pSparse->isComplex() == false)
1214     {
1215         for (int i = 0 ; i < iNonZeros ; i++)
1216         {
1217             int iRow = static_cast<int>(pRows[i]) - 1;
1218             int iCol = static_cast<int>(pCols[i]) - 1;
1219             double dbl = pValR[i];
1220
1221             for (int j = 0 ; j < _pDouble->getRows() ; j++)
1222             {
1223                 double dblVal = _pDouble->get(j, iRow) * dbl;
1224                 pOut->set(j, iCol, pOut->get(j, iCol) + dblVal);
1225             }
1226         }
1227     }
1228     else if (_pDouble->isComplex() == false && _pSparse->isComplex() == true)
1229     {
1230         //a * (b ci) -> ab ac
1231         for (int i = 0 ; i < iNonZeros ; i++)
1232         {
1233             int iRow = static_cast<int>(pRows[i]) - 1;
1234             int iCol = static_cast<int>(pCols[i]) - 1;
1235             double dblR = pValR[i];
1236             double dblI = pValI[i];
1237
1238             for (int j = 0 ; j < _pDouble->getRows() ; j++)
1239             {
1240                 double dblValR = _pDouble->get(j, iRow) * dblR;
1241                 double dblValI = _pDouble->get(j, iRow) * dblI;
1242                 pOut->set(j, iCol, pOut->get(j, iCol) + dblValR);
1243                 pOut->setImg(j, iCol, pOut->getImg(j, iCol) + dblValI);
1244             }
1245         }
1246     }
1247     else if (_pDouble->isComplex() == true && _pSparse->isComplex() == false)
1248     {
1249         //(a bi) * c -> ac + bc
1250         for (int i = 0 ; i < iNonZeros ; i++)
1251         {
1252             int iRow = static_cast<int>(pRows[i]) - 1;
1253             int iCol = static_cast<int>(pCols[i]) - 1;
1254             double dblR = pValR[i];
1255
1256             for (int j = 0 ; j < _pDouble->getRows() ; j++)
1257             {
1258                 double dblValR = _pDouble->get(j, iRow) * dblR;
1259                 double dblValI = _pDouble->getImg(j, iRow) * dblR;
1260                 pOut->set(j, iCol, pOut->get(j, iCol) + dblValR);
1261                 pOut->setImg(j, iCol, pOut->getImg(j, iCol) + dblValI);
1262             }
1263         }
1264     }
1265     else if (_pDouble->isComplex() == true && _pSparse->isComplex() == true)
1266     {
1267         for (int i = 0 ; i < iNonZeros ; i++)
1268         {
1269             int iRow = static_cast<int>(pRows[i]) - 1;
1270             int iCol = static_cast<int>(pCols[i]) - 1;
1271             double dblR = pValR[i];
1272             double dblI = pValI[i];
1273
1274             for (int j = 0 ; j < _pDouble->getRows() ; j++)
1275             {
1276                 double dblValR = _pDouble->get(j, iRow) * dblR - _pDouble->getImg(j, iRow) * dblI;
1277                 double dblValI = _pDouble->get(j, iRow) * dblI + _pDouble->getImg(j, iRow) * dblR;
1278                 pOut->set(j, iCol, pOut->get(j, iCol) + dblValR);
1279                 pOut->setImg(j, iCol, pOut->getImg(j, iCol) + dblValI);
1280             }
1281         }
1282     }
1283
1284     *_pOut = pOut;
1285     delete[] pRows;
1286     delete[] pValR;
1287     delete[] pValI;
1288
1289     return 0;
1290 }
1291
1292 int MultiplySparseByDouble(Sparse *_pSparse, Double*_pDouble, GenericType** _pOut)
1293 {
1294     if (_pDouble->isScalar())
1295     {
1296         //SP * d -> SP
1297         Sparse* pOut = NULL;
1298         if (_pDouble->isComplex())
1299         {
1300             std::complex<double> dbl(_pDouble->get(0), _pDouble->getImg(0));
1301             pOut = _pSparse->multiply(dbl);
1302         }
1303         else
1304         {
1305             pOut = _pSparse->multiply(_pDouble->get(0));
1306         }
1307         *_pOut = pOut;
1308         return 0;
1309     }
1310
1311     if (_pSparse->isScalar())
1312     {
1313         //D * sp -> D .* d
1314         Double* pD = NULL;
1315
1316         if (_pSparse->isComplex())
1317         {
1318             std::complex<double> dbl(_pSparse->getImg(0, 0));
1319             pD = new Double(dbl.real(), dbl.imag());
1320         }
1321         else
1322         {
1323             pD = new Double(_pSparse->get(0, 0));
1324         }
1325
1326         InternalType* pIT = GenericDotTimes(_pDouble, pD);
1327         *_pOut = pIT->getAs<GenericType>();
1328         delete pD;
1329         return 0;
1330     }
1331
1332     if (_pSparse->getCols() != _pDouble->getRows())
1333     {
1334         return 1;
1335     }
1336
1337     //try to be smart and only compute for non zero values
1338
1339     //get some information
1340     int iNonZeros = static_cast<int>(_pSparse->nonZeros());
1341     int* pRows = new int[iNonZeros * 2];
1342     _pSparse->outputRowCol(pRows);
1343     int* pCols = pRows + iNonZeros;
1344     double* pValR = new double[iNonZeros];
1345     double* pValI = new double[iNonZeros];
1346     _pSparse->outputValues(pValR, pValI);
1347
1348     Double* pOut = new Double(_pSparse->getRows(), _pDouble->getCols(), _pDouble->isComplex() | _pSparse->isComplex());
1349     pOut->setZeros();
1350
1351     if (_pDouble->isComplex() == false && _pSparse->isComplex() == false)
1352     {
1353         for (int i = 0 ; i < iNonZeros ; i++)
1354         {
1355             int iRow    = static_cast<int>(pRows[i]) - 1;
1356             int iCol    = static_cast<int>(pCols[i]) - 1;
1357             double dbl  = pValR[i];
1358
1359             for (int j = 0 ; j < _pDouble->getCols() ; j++)
1360             {
1361                 double dblVal = _pDouble->get(iCol, j) * dbl;
1362                 pOut->set(iRow, j, pOut->get(iRow, j) + dblVal);
1363             }
1364         }
1365     }
1366     else if (_pDouble->isComplex() == false && _pSparse->isComplex() == true)
1367     {
1368         //a * (b ci) -> ab ac
1369         for (int i = 0 ; i < iNonZeros ; i++)
1370         {
1371             int iRow = static_cast<int>(pRows[i]) - 1;
1372             int iCol = static_cast<int>(pCols[i]) - 1;
1373             double dblR = pValR[i];
1374             double dblI = pValI[i];
1375
1376             for (int j = 0 ; j < _pDouble->getCols() ; j++)
1377             {
1378                 double dblValR = _pDouble->get(iCol, j) * dblR;
1379                 double dblValI = _pDouble->get(iCol, j) * dblI;
1380                 pOut->set(iRow, j, pOut->get(iRow, j) + dblValR);
1381                 pOut->setImg(iRow, j, pOut->getImg(iRow, j) + dblValI);
1382             }
1383         }
1384     }
1385     else if (_pDouble->isComplex() == true && _pSparse->isComplex() == false)
1386     {
1387         //(a bi) * c -> ac + bc
1388         for (int i = 0 ; i < iNonZeros ; i++)
1389         {
1390             int iRow = static_cast<int>(pRows[i]) - 1;
1391             int iCol = static_cast<int>(pCols[i]) - 1;
1392             double dblR = pValR[i];
1393
1394             for (int j = 0 ; j < _pDouble->getCols() ; j++)
1395             {
1396                 double dblValR = _pDouble->get(iCol, j) * dblR;
1397                 double dblValI = _pDouble->getImg(iCol, j) * dblR;
1398                 pOut->set(iRow, j, pOut->get(iRow, j) + dblValR);
1399                 pOut->setImg(iRow, j, pOut->getImg(iRow, j) + dblValI);
1400             }
1401         }
1402     }
1403     else if (_pDouble->isComplex() == true && _pSparse->isComplex() == true)
1404     {
1405         for (int i = 0 ; i < iNonZeros ; i++)
1406         {
1407             int iRow = static_cast<int>(pRows[i]) - 1;
1408             int iCol = static_cast<int>(pCols[i]) - 1;
1409             double dblR = pValR[i];
1410             double dblI = pValI[i];
1411
1412             for (int j = 0 ; j < _pDouble->getCols() ; j++)
1413             {
1414                 double dblValR = _pDouble->get(iCol, j) * dblR - _pDouble->getImg(iCol, j) * dblI;
1415                 double dblValI = _pDouble->get(iCol, j) * dblI + _pDouble->getImg(iCol, j) * dblR;
1416                 pOut->set(iRow, j, pOut->get(iRow, j) + dblValR);
1417                 pOut->setImg(iRow, j, pOut->getImg(iRow, j) + dblValI);
1418             }
1419         }
1420     }
1421
1422     *_pOut = pOut;
1423     delete[] pRows;
1424     delete[] pValR;
1425     delete[] pValI;
1426
1427     return 0;
1428 }
1429
1430 int DotMultiplySparseBySparse(Sparse* _pSparse1, Sparse* _pSparse2, Sparse** _pOut)
1431 {
1432     if (_pSparse1->isScalar() || _pSparse2->isScalar())
1433     {
1434         //SP .* sp or sp .* SP
1435         return MultiplySparseBySparse(_pSparse1, _pSparse2, _pOut);
1436     }
1437
1438     if (_pSparse1->getRows() != _pSparse2->getRows() || _pSparse1->getCols() != _pSparse2->getCols())
1439     {
1440         return 1;
1441     }
1442
1443     *_pOut = _pSparse1->dotMultiply(*_pSparse2);
1444
1445     return 0;
1446 }
1447
1448 int DotMultiplyDoubleBySparse(Double* _pDouble, Sparse* _pSparse, GenericType**  _pOut)
1449 {
1450     if (_pDouble->isScalar())
1451     {
1452         return MultiplyDoubleBySparse(_pDouble, _pSparse, _pOut);
1453     }
1454
1455     if (_pSparse->isScalar())
1456     {
1457         return MultiplyDoubleBySparse(_pDouble, _pSparse, _pOut);
1458     }
1459
1460     if (_pSparse->getRows() != _pDouble->getRows() || _pSparse->getCols() != _pDouble->getCols())
1461     {
1462         return 1;
1463     }
1464
1465     Sparse* pOut = new Sparse(_pDouble->getRows(), _pDouble->getCols(), _pSparse->isComplex() || _pDouble->isComplex());
1466     //get some information
1467     int iNonZeros = static_cast<int>(_pSparse->nonZeros());
1468     int* pRows = new int[iNonZeros * 2];
1469     _pSparse->outputRowCol(pRows);
1470     int* pCols = pRows + iNonZeros;
1471
1472     if (_pDouble->isComplex() == false && _pSparse->isComplex() == false)
1473     {
1474         for (int i = 0 ; i < iNonZeros ; i++)
1475         {
1476             int iRow = static_cast<int>(pRows[i]) - 1;
1477             int iCol = static_cast<int>(pCols[i]) - 1;
1478             pOut->set(iRow, iCol, _pSparse->get(iRow, iCol) * _pDouble->get(iRow, iCol));
1479         }
1480     }
1481     else if (_pDouble->isComplex() == false && _pSparse->isComplex() == true)
1482     {
1483         for (int i = 0 ; i < iNonZeros ; i++)
1484         {
1485             int iRow = static_cast<int>(pRows[i]) - 1;
1486             int iCol = static_cast<int>(pCols[i]) - 1;
1487             std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1488             std::complex<double> newVal(dbl.real() * _pDouble->get(iRow, iCol), dbl.imag() * _pDouble->get(iRow, iCol));
1489             pOut->set(iRow, iCol, newVal);
1490         }
1491     }
1492     else if (_pDouble->isComplex() == true && _pSparse->isComplex() == false)
1493     {
1494         for (int i = 0 ; i < iNonZeros ; i++)
1495         {
1496             int iRow = static_cast<int>(pRows[i]) - 1;
1497             int iCol = static_cast<int>(pCols[i]) - 1;
1498             std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1499             std::complex<double> newVal(dbl.real() * _pDouble->get(iRow, iCol), dbl.real() * _pDouble->getImg(iRow, iCol));
1500             pOut->set(iRow, iCol, newVal);
1501         }
1502     }
1503     else if (_pDouble->isComplex() == true && _pSparse->isComplex() == true)
1504     {
1505         for (int i = 0 ; i < iNonZeros ; i++)
1506         {
1507             int iRow = static_cast<int>(pRows[i]) - 1;
1508             int iCol = static_cast<int>(pCols[i]) - 1;
1509             std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1510             double dblR = _pDouble->get(iRow, iCol) * dbl.real() - _pDouble->getImg(iRow, iCol) * dbl.imag();
1511             double dblI = _pDouble->getImg(iRow, iCol) * dbl.real() + _pDouble->get(iRow, iCol) * dbl.imag();
1512
1513             std::complex<double> newVal(dblR, dblI);
1514             pOut->set(iRow, iCol, newVal);
1515         }
1516     }
1517
1518     *_pOut = pOut;
1519     delete[] pRows;
1520
1521     return 0;
1522 }
1523
1524 int DotMultiplySparseByDouble(Sparse* _pSparse, Double* _pDouble, GenericType** _pOut)
1525 {
1526     return DotMultiplyDoubleBySparse(_pDouble, _pSparse, _pOut);
1527 }
1528
1529 int DotMultiplyPolyByDouble(Polynom* _pPoly, Double* _pDouble, Polynom** _pPolyOut)
1530 {
1531     return DotMultiplyDoubleByPoly(_pDouble, _pPoly, _pPolyOut);
1532 }
1533
1534 int DotMultiplyDoubleByPoly(Double* _pDouble, Polynom* _pPoly, Polynom** _pPolyOut)
1535 {
1536     int iSize = _pDouble->getSize();
1537     if (_pDouble->isScalar() == false &&
1538             _pPoly->isScalar() == false &&
1539             iSize != _pPoly->getSize())
1540     {
1541         return 1;
1542     }
1543
1544     int* piRanks = new int[iSize];
1545     memset(piRanks, 0x00, iSize * sizeof(int));
1546     Polynom* pPolyTemp = new Polynom(_pPoly->getVariableName(), _pDouble->getDims(), _pDouble->getDimsArray(), piRanks);
1547     delete[] piRanks;
1548     pPolyTemp->setCoef(_pDouble);
1549     int iErr = DotMultiplyPolyByPoly(pPolyTemp, _pPoly, _pPolyOut);
1550     delete pPolyTemp;
1551     return iErr;
1552 }
1553
1554 int DotMultiplyPolyByPoly(Polynom* _pPoly1, Polynom* _pPoly2, Polynom** _pPolyOut)
1555 {
1556     if (_pPoly1->isScalar() || _pPoly2->isScalar())
1557     {
1558         return MultiplyPolyByPoly(_pPoly1, _pPoly2, _pPolyOut);
1559     }
1560     else
1561     {
1562         if (_pPoly1->getSize() != _pPoly2->getSize())
1563         {
1564             return 1;
1565         }
1566
1567         int* piRank = new int[_pPoly1->getSize()];
1568         for (int i = 0 ; i < _pPoly1->getSize() ; i++)
1569         {
1570             piRank[i] = _pPoly1->get(i)->getRank() + _pPoly2->get(i)->getRank();
1571         }
1572
1573         (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), _pPoly1->getDims(), _pPoly1->getDimsArray(), piRank);
1574
1575         if (_pPoly1->isComplex() && _pPoly2->isComplex())
1576         {
1577             (*_pPolyOut)->setComplex(true);
1578             for (int i = 0; i < _pPoly1->getSize(); i++)
1579             {
1580                 SinglePoly *pSP1    = _pPoly1->get(i);
1581                 SinglePoly *pSP2    = _pPoly2->get(i);
1582                 SinglePoly *pSPOut  = (*_pPolyOut)->get(i);
1583
1584                 pSPOut->setZeros();
1585
1586                 iMultiComplexPolyByComplexPoly(
1587                     pSP1->get(), pSP1->getImg(), pSP1->getSize(),
1588                     pSP2->get(), pSP2->getImg(), pSP2->getSize(),
1589                     pSPOut->get(), pSPOut->getImg(), pSPOut->getSize());
1590
1591             }
1592         }
1593         else if (_pPoly1->isComplex())
1594         {
1595             (*_pPolyOut)->setComplex(true);
1596             for (int i = 0; i < _pPoly1->getSize(); i++)
1597             {
1598                 SinglePoly *pSP1   = _pPoly1->get(i);
1599                 SinglePoly *pSP2   = _pPoly2->get(i);
1600                 SinglePoly *pSPOut = (*_pPolyOut)->get(i);
1601
1602                 pSPOut->setZeros();
1603
1604                 iMultiComplexPolyByScilabPolynom(
1605                     pSP1->get(), pSP1->getImg(), pSP1->getSize(),
1606                     pSP2->get(), pSP2->getSize(),
1607                     pSPOut->get(), pSPOut->getImg(), pSPOut->getSize());
1608             }
1609         }
1610         else if (_pPoly2->isComplex())
1611         {
1612             (*_pPolyOut)->setComplex(true);
1613             for (int i = 0; i < _pPoly1->getSize(); i++)
1614             {
1615                 SinglePoly *pSP1   = _pPoly1->get(i);
1616                 SinglePoly *pSP2   = _pPoly2->get(i);
1617                 SinglePoly *pSPOut = (*_pPolyOut)->get(i);
1618
1619                 pSPOut->setZeros();
1620
1621                 iMultiScilabPolynomByComplexPoly(
1622                     pSP1->get(), pSP1->getSize(),
1623                     pSP2->get(), pSP2->getImg(), pSP2->getSize(),
1624                     pSPOut->get(), pSPOut->getImg(), pSPOut->getSize());
1625             }
1626         }
1627         else
1628         {
1629             for (int i = 0; i < _pPoly1->getSize(); i++)
1630             {
1631                 SinglePoly *pSP1   = _pPoly1->get(i);
1632                 SinglePoly *pSP2   = _pPoly2->get(i);
1633                 SinglePoly *pSPOut = (*_pPolyOut)->get(i);
1634
1635                 pSPOut->setZeros();
1636
1637                 iMultiScilabPolynomByScilabPolynom(
1638                     pSP1->get(), pSP1->getSize(),
1639                     pSP2->get(), pSP2->getSize(),
1640                     pSPOut->get(), pSPOut->getSize());
1641             }
1642         }
1643     }
1644
1645     return 0;
1646 }
1647