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