fix divide Sparse by Double
[scilab.git] / scilab / modules / ast / src / cpp / operations / types_divide.cpp
1 /*
2  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  *  Copyright (C) 2008-2008 - DIGITEO - Antoine ELIAS
4  *
5  *  This file must be used under the terms of the CeCILL.
6  *  This source file is licensed as described in the file COPYING, which
7  *  you should have received as part of this distribution.  The terms
8  *  are also available at
9  *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10  *
11  */
12
13 #include "types_divide.hxx"
14 #include "types_finite.hxx"
15
16 #include "scilabexception.hxx"
17
18 extern "C"
19 {
20 #include "matrix_right_division.h"
21 #include "sciprint.h"
22 #include "localization.h"
23 #include "charEncoding.h"
24 #include "configvariable_interface.h"
25 #include "elem_common.h"
26 }
27
28 using namespace types;
29
30 InternalType *GenericRDivide(InternalType *_pLeftOperand, InternalType *_pRightOperand)
31 {
32     InternalType *pResult       = NULL;
33     GenericType::ScilabType TypeL = _pLeftOperand->getType();
34     GenericType::ScilabType TypeR = _pRightOperand->getType();
35
36     int iResult = 0;
37
38     if (_pLeftOperand->isDouble() && _pLeftOperand->getAs<Double>()->isEmpty())
39     {
40         return Double::Empty();
41     }
42
43     if (_pRightOperand->isDouble() && _pRightOperand->getAs<Double>()->isEmpty())
44     {
45         return Double::Empty();
46     }
47
48     /*
49     ** DOUBLE / DOUBLE
50     */
51     if (TypeL == GenericType::ScilabDouble && TypeR == GenericType::ScilabDouble)
52     {
53         Double *pL  = _pLeftOperand->getAs<Double>();
54         Double *pR  = _pRightOperand->getAs<Double>();
55
56         iResult = RDivideDoubleByDouble(pL, pR, (Double**)&pResult);
57     }
58
59     /*
60     ** POLY / DOUBLE
61     */
62     else if (TypeL == GenericType::ScilabPolynom && TypeR == GenericType::ScilabDouble)
63     {
64         Polynom *pL = _pLeftOperand->getAs<types::Polynom>();
65         Double *pR  = _pRightOperand->getAs<Double>();
66
67         iResult = RDividePolyByDouble(pL, pR, (Polynom**)&pResult);
68     }
69
70     /*
71     ** DOUBLE / POLY
72     */
73     else if (TypeL == GenericType::ScilabDouble && TypeR == GenericType::ScilabPolynom)
74     {
75         Double *pL  = _pLeftOperand->getAs<Double>();
76         Polynom *pR = _pRightOperand->getAs<types::Polynom>();
77
78         iResult = RDivideDoubleByPoly(pL, pR, (Polynom**)&pResult);
79     }
80
81     /*
82     ** SPARSE / DOUBLE
83     */
84     else if (TypeL == GenericType::ScilabSparse && TypeR == GenericType::ScilabDouble)
85     {
86         Sparse *pL = _pLeftOperand->getAs<types::Sparse>();
87         Double *pR = _pRightOperand->getAs<Double>();
88
89         iResult = RDivideSparseByDouble(pL, pR, &pResult);
90     }
91
92     //manage errors
93     if (iResult)
94     {
95         switch (iResult)
96         {
97             case 1 :
98                 throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
99             case 2 :
100                 throw ast::ScilabError(_W("With NaN or Inf a division by scalar expected.\n"));
101             case 3 :
102                 throw ast::ScilabError(_W("Division by zero...\n"));
103             case 4 :
104                 if (getWarningMode())
105                 {
106                     sciprint(_("Warning : Division by zero...\n"));
107                 }
108                 break;
109                 //            default : throw ast::ScilabError(_W("Operator / : Error %d not yet managed.\n"), iResult);
110             default :
111                 sciprint(_("Operator / : Error %d not yet managed.\n"), iResult);
112         }
113     }
114
115     /*
116     ** Default case : Return NULL will Call Overloading.
117     */
118     return pResult;
119 }
120
121
122 int RDivideDoubleByDouble(Double *_pDouble1, Double *_pDouble2, Double **_pDoubleOut)
123 {
124     int iErr = 0;
125
126     //check finite values of _pDouble1 and _pDouble2
127     if (isDoubleFinite(_pDouble1) == false || isDoubleFinite(_pDouble2) == false)
128     {
129         if (_pDouble2->isScalar() == false)
130         {
131             return 2;
132         }
133     }
134
135     if (_pDouble2->isScalar())
136     {
137         //Y / x
138         int iInc1       = 1;
139         int iInc2       = 0;
140         bool bComplex1  = _pDouble1->isComplex();
141         bool bComplex2  = _pDouble2->isComplex();
142
143         *_pDoubleOut    = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray(), bComplex1 || bComplex2);
144
145         if (bComplex1 == false && bComplex2 == false)
146         {
147             // Real1 \ Real2 -> Real2 / Real1
148             iErr = iRightDivisionRealMatrixByRealMatrix(
149                        _pDouble1->get(), iInc1,
150                        _pDouble2->get(), iInc2,
151                        (*_pDoubleOut)->get(), 1, _pDouble1->getSize());
152         }
153         else if (bComplex1 == false && bComplex2 == true)
154         {
155             // Real \ Complex -> Complex / Real
156             iErr = iRightDivisionRealMatrixByComplexMatrix(
157                        _pDouble1->get(), iInc1,
158                        _pDouble2->get(), _pDouble2->getImg(), iInc2,
159                        (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), 1, _pDouble1->getSize());
160         }
161         else if (bComplex1 == true && bComplex2 == false)
162         {
163             // Complex \ Real -> Real / Complex
164             iErr = iRightDivisionComplexMatrixByRealMatrix(
165                        _pDouble1->get(), _pDouble1->getImg(), iInc1,
166                        _pDouble2->get(), iInc2,
167                        (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), 1, _pDouble1->getSize());
168         }
169         else if (bComplex1 == true && bComplex2 == true)
170         {
171             // Complex \ Complex
172             iErr = iRightDivisionComplexMatrixByComplexMatrix(
173                        _pDouble1->get(), _pDouble1->getImg(), iInc1,
174                        _pDouble2->get(), _pDouble2->getImg(), iInc2,
175                        (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), 1, _pDouble1->getSize());
176         }
177
178         return iErr;
179     }
180
181     if (_pDouble1->isScalar())
182     {
183         if (_pDouble2->getDims() > 2)
184         {
185             //not managed, call overload
186             return 0;
187         }
188
189         double dblSavedR = 0;
190         double dblSavedI = 0;
191         Double *pdblTemp = NULL;
192
193         int iRowResult = _pDouble2->getCols();
194         int iColResult = _pDouble2->getRows();
195
196         //in this case, we have to create a temporary square matrix
197         pdblTemp = new Double(iRowResult, iRowResult, _pDouble1->isComplex());
198         pdblTemp->setZeros();
199
200         if (_pDouble1->isComplex())
201         {
202             dblSavedR = _pDouble1->getReal()[0];
203             dblSavedI = _pDouble1->getImg()[0];
204             for (int i = 0 ; i < iRowResult ; i++)
205             {
206                 pdblTemp->set(i, i, dblSavedR);
207                 pdblTemp->setImg(i, i, dblSavedI);
208             }
209         }
210         else
211         {
212             dblSavedR = _pDouble1->getReal()[0];
213             for (int i = 0 ; i < iRowResult ; i++)
214             {
215                 pdblTemp->set(i, i, dblSavedR);
216             }
217         }
218
219         *_pDoubleOut = new Double(iRowResult, iColResult, _pDouble1->isComplex() || _pDouble2->isComplex());
220
221         if ((*_pDoubleOut)->isComplex())
222         {
223             double dblRcond = 0;
224             iErr = iRightDivisionOfComplexMatrix(
225                        pdblTemp->getReal(), pdblTemp->getImg(), pdblTemp->getRows(), pdblTemp->getCols(),
226                        _pDouble2->getReal(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
227                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iRowResult, iColResult, &dblRcond);
228         }
229         else
230         {
231             double dblRcond = 0;
232             iErr = iRightDivisionOfRealMatrix(
233                        pdblTemp->getReal(), pdblTemp->getRows(), pdblTemp->getCols(),
234                        _pDouble2->getReal(), _pDouble2->getRows(), _pDouble2->getCols(),
235                        (*_pDoubleOut)->getReal(), iRowResult, iColResult, &dblRcond);
236         }
237         delete pdblTemp;
238         return iErr;
239     }
240
241     if (_pDouble1->getDims() > 2 || _pDouble2->getDims() > 2 || _pDouble1->getCols() != _pDouble2->getCols())
242     {
243         //not managed
244         return 1;
245     }
246
247     *_pDoubleOut = new Double(_pDouble1->getRows(), _pDouble2->getRows(), _pDouble1->isComplex() || _pDouble2->isComplex());
248     if ((*_pDoubleOut)->isComplex())
249     {
250         double dblRcond = 0;
251         iErr = iRightDivisionOfComplexMatrix(
252                    _pDouble1->getReal(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
253                    _pDouble2->getReal(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
254                    (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), _pDouble1->getRows(), _pDouble2->getRows(), &dblRcond);
255     }
256     else
257     {
258         double dblRcond = 0;
259         iErr = iRightDivisionOfRealMatrix(
260                    _pDouble1->getReal(), _pDouble1->getRows(), _pDouble1->getCols(),
261                    _pDouble2->getReal(), _pDouble2->getRows(), _pDouble2->getCols(),
262                    (*_pDoubleOut)->getReal(), _pDouble1->getRows(), _pDouble2->getRows(), &dblRcond);
263     }
264
265     return iErr;
266 }
267
268 int RDividePolyByDouble(Polynom* _pPoly, Double* _pDouble, Polynom** _pPolyOut)
269 {
270     bool bComplex1  = _pPoly->isComplex();
271     bool bComplex2  = _pDouble->isComplex();
272     bool bScalar1   = _pPoly->getRows() == 1  && _pPoly->getCols() == 1;
273     bool bScalar2   = _pDouble->getRows() == 1 && _pDouble->getCols() == 1;
274
275     Polynom *pTemp = NULL; //use only if _pPoly is scalar and _pDouble not.
276
277     int iRowResult  = 0;
278     int iColResult = 0;
279     int *piRank   = NULL;
280
281     /* if(bScalar1 && bScalar2)
282     {
283     iRowResult = 1;
284     iColResult = 1;
285
286     piRank = new int[1];
287     piRank[0] = _pPoly->get(0)->getRank();
288     }
289     else */
290
291     if (bScalar2)
292     {
293         double dblDivR = _pDouble->get(0);
294         double dblDivI = _pDouble->getImg(0);
295
296         (*_pPolyOut) = _pPoly->clone()->getAs<Polynom>();
297         if (_pDouble->isComplex())
298         {
299             (*_pPolyOut)->setComplex(true);
300         }
301
302         for (int i = 0 ; i < _pPoly->getSize() ; i++)
303         {
304             bool bComplex1 = _pPoly->isComplex();
305             bool bComplex2 = _pDouble->isComplex();
306
307             SinglePoly* pC = (*_pPolyOut)->get(i);
308
309             if (bComplex1 == false && bComplex2 == false)
310             {
311                 iRightDivisionRealMatrixByRealMatrix(pC->get(), 1, &dblDivR, 0, pC->get(), 1, pC->getSize());
312             }
313             else if (bComplex1 == true && bComplex2 == false)
314             {
315                 iRightDivisionComplexMatrixByRealMatrix(pC->get(), pC->getImg(), 1, &dblDivR, 0, pC->get(), pC->getImg(), 1, pC->getSize());
316             }
317             else if (bComplex1 == false && bComplex2 == true)
318             {
319                 iRightDivisionRealMatrixByComplexMatrix(pC->get(), 1, &dblDivR, &dblDivI, 0, pC->get(), pC->getImg(), 1, pC->getSize());
320             }
321             else if (bComplex1 == true && bComplex2 == true)
322             {
323                 iRightDivisionComplexMatrixByComplexMatrix(pC->get(), pC->getImg(), 1, &dblDivR, &dblDivI, 0, pC->get(), pC->getImg(), 1, pC->getSize());
324             }
325         }
326
327         return 0;
328     }
329
330     if (bScalar1)
331     {
332         //in this case, we have to create a temporary square polinomial matrix
333         iRowResult = _pDouble->getCols();
334         iColResult = _pDouble->getRows();
335
336         piRank = new int[iRowResult * iRowResult];
337         int iMaxRank = _pPoly->getMaxRank();
338         for (int i = 0 ; i < iRowResult * iRowResult ; i++)
339         {
340             piRank[i] = iMaxRank;
341         }
342
343         pTemp = new Polynom(_pPoly->getVariableName(), iRowResult, iRowResult, piRank);
344         if (bComplex1 || bComplex2)
345         {
346             pTemp->setComplex(true);
347         }
348
349         SinglePoly *pdblData = _pPoly->get(0);
350         for (int i = 0 ; i < iRowResult ; i++)
351         {
352             pTemp->set(i, i, pdblData);
353         }
354     }
355
356     (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), iRowResult, iColResult, piRank);
357     delete[] piRank;
358     if (bComplex1 || bComplex2)
359     {
360         (*_pPolyOut)->setComplex(true);
361     }
362
363     if (bScalar2)
364     {
365         //[p] * cst
366         for (int i = 0 ; i < _pPoly->getSize() ; i++)
367         {
368             SinglePoly *pPolyIn   = _pPoly->get(i);
369             double* pRealIn  = pPolyIn->get();
370             double* pImgIn  = pPolyIn->getImg();
371
372             SinglePoly *pPolyOut  = (*_pPolyOut)->get(i);
373             double* pRealOut = pPolyOut->get();
374             double* pImgOut  = pPolyOut->getImg();
375
376             if (bComplex1 == false && bComplex2 == false)
377             {
378                 iRightDivisionRealMatrixByRealMatrix(
379                     pRealIn, 1,
380                     _pDouble->getReal(), 0,
381                     pRealOut, 1, pPolyOut->getSize());
382             }
383             else if (bComplex1 == false && bComplex2 == true)
384             {
385                 iRightDivisionRealMatrixByComplexMatrix(
386                     pRealIn, 1,
387                     _pDouble->getReal(), _pDouble->getImg(), 0,
388                     pRealOut, pImgOut, 1, pPolyOut->getSize());
389             }
390             else if (bComplex1 == true && bComplex2 == false)
391             {
392                 iRightDivisionComplexMatrixByRealMatrix(
393                     pRealIn, pImgIn, 1,
394                     _pDouble->getReal(), 0,
395                     pRealOut, pImgOut, 1, pPolyOut->getSize());
396             }
397             else if (bComplex1 == true && bComplex2 == true)
398             {
399                 iRightDivisionComplexMatrixByComplexMatrix(
400                     pRealIn, pImgIn, 1,
401                     _pDouble->getReal(), _pDouble->getImg(), 0,
402                     pRealOut, pImgOut, 1, pPolyOut->getSize());
403             }
404         }
405     }
406     else if (bScalar1)
407     {
408         for (int i = 0 ; i < pTemp->get(0)->getSize() ; i++)
409         {
410             Double *pCoef    = pTemp->extractCoef(i);
411             Double *pResultCoef = new Double(iRowResult, iColResult, pCoef->isComplex());
412             double *pReal    = pResultCoef->getReal();
413             double *pImg    = pResultCoef->getImg();
414
415             if (bComplex1 == false && bComplex2 == false)
416             {
417                 double dblRcond = 0;
418                 iRightDivisionOfRealMatrix(
419                     pCoef->getReal(), iRowResult, iRowResult,
420                     _pDouble->getReal(), _pDouble->getRows(), _pDouble->getCols(),
421                     pReal, iRowResult, iColResult, &dblRcond);
422             }
423             else
424             {
425                 double dblRcond = 0;
426                 iRightDivisionOfComplexMatrix(
427                     pCoef->getReal(), pCoef->getImg(), iRowResult, iRowResult,
428                     _pDouble->getReal(), _pDouble->getImg(), _pDouble->getRows(), _pDouble->getCols(),
429                     pReal, pImg, iRowResult, iColResult, &dblRcond);
430             }
431
432             (*_pPolyOut)->insertCoef(i, pResultCoef);
433             delete pCoef;
434             delete pResultCoef;
435         }
436     }
437     return 0;
438 }
439
440 int RDivideDoubleByPoly(Double* /*_pDouble*/, Polynom* /*_pPoly*/, Polynom** /*_pPolyOut*/)
441 {
442     return 0;
443 }
444
445 int RDivideSparseByDouble(types::Sparse* _pSp, types::Double* _pDouble, InternalType** _pSpOut)
446 {
447     if (_pDouble->isEmpty())
448     {
449         //sp / []
450         *_pSpOut = Double::Empty();
451         return 0;
452     }
453
454     size_t iSize = _pSp->nonZeros();
455     int* Col = new int[iSize];
456     int* Row = new int[iSize];
457     _pSp->getColPos(Col);
458     _pSp->getNbItemByRow(Row);
459     int* iPositVal = new int[iSize];
460
461     int j = 0;
462     for (int i = 0; i < iSize; j++)
463     {
464         for (int k = 0; k < Row[j]; k++)
465         {
466             iPositVal[i] = (Col[i] - 1) * _pSp->getRows() + j;
467             i++;
468         }
469     }
470
471     Double** pDbl = new Double*[iSize];
472     Double** pDblSp = new Double*[iSize];
473     double* pdbl = _pDouble->get();
474
475     if (_pDouble->isScalar())
476     {
477         if (_pDouble->isComplex())
478         {
479             double* pdblImg = _pDouble->getImg();
480             for (int i = 0; i < iSize; i++)
481             {
482                 pDbl[i] = new Double(pdbl[0], pdblImg[0]);
483                 pDblSp[i] = new Double(_pSp->get(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
484             }
485         }
486         else
487         {
488             for (int i = 0; i < iSize; i++)
489             {
490                 pDbl[i] = new Double(pdbl[0]);
491                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
492             }
493         }
494     }
495     else if (_pDouble->getSize() == iSize)
496     {
497         if (_pDouble->isComplex())
498         {
499             double* pdblImg = _pDouble->getImg();
500             for (int i = 0; i < iSize; i++)
501             {
502                 pDbl[i] = new Double(pdbl[i], pdblImg[i]);
503                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
504             }
505         }
506         else
507         {
508             for (int i = 0; i < iSize; i++)
509             {
510                 pDbl[i] = new Double(pdbl[i]);
511                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
512             }
513         }
514     }
515     else
516     {
517         delete[] pDblSp;
518         throw ast::ScilabError(_W("Invalid exponent.\n"));
519         return 1;
520     }
521
522     Sparse* pSpTemp = new Sparse(_pSp->getRows(), _pSp->getCols(), _pSp->isComplex() || _pDouble->isComplex());
523     pSpTemp->zero_set();
524
525     Double* ppDblGet = NULL;
526     int iResultat;
527     for (int i = 0; i < iSize; i++)
528     {
529         if ((pDblSp[i]->get(0) != 0) || (pDblSp[i]->getImg(0) != 0))
530         {
531             iResultat = RDivideDoubleByDouble(pDblSp[i], pDbl[i], &ppDblGet);
532             if (iResultat != 0)
533             {
534                 return iResultat;
535             }
536             std::complex<double> cplx(ppDblGet->get(0), ppDblGet->getImg(0));
537             pSpTemp->set(iPositVal[i], cplx, true);
538         }
539     }
540     delete Col;
541     delete iPositVal;
542
543     *_pSpOut = pSpTemp;
544     return 0;
545 }
546
547
548
549
550 int DotRDivideDoubleByDouble(Double* _pDouble1, Double* _pDouble2, Double** _pDoubleOut)
551 {
552     int iErr        = 0;
553     bool bComplex1  = _pDouble1->isComplex();
554     bool bComplex2  = _pDouble2->isComplex();
555     bool bScalar1   = _pDouble1->isScalar();
556     bool bScalar2   = _pDouble2->isScalar();
557
558     if (bScalar1)
559     {
560         //x ./ Y
561         int iInc1       = 0;
562         int iInc2       = 1;
563         int iIncOut     = 1;
564         int iSizeResult = _pDouble2->getSize();
565
566         *_pDoubleOut    = new Double(_pDouble2->getDims(), _pDouble2->getDimsArray(), bComplex1 || bComplex2);
567
568         if (bComplex1 == false && bComplex2 == false)
569         {
570             // r ./ R
571             iErr = iRightDivisionRealMatrixByRealMatrix(
572                        _pDouble1->getReal(), iInc1,
573                        _pDouble2->getReal(), iInc2,
574                        (*_pDoubleOut)->getReal(), iIncOut, iSizeResult);
575         }
576         else if (bComplex1 == false && bComplex2 == true)
577         {
578             // r ./ C
579             iErr = iRightDivisionRealMatrixByComplexMatrix(
580                        _pDouble1->getReal(), iInc1,
581                        _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
582                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
583         }
584         else if (bComplex1 == true && bComplex2 == false)
585         {
586             // c ./ R
587             iErr = iRightDivisionComplexMatrixByRealMatrix(
588                        _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
589                        _pDouble2->getReal(), iInc2,
590                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
591         }
592         else if (bComplex1 == true && bComplex2 == true)
593         {
594             // c ./ C
595             iErr = iRightDivisionComplexMatrixByComplexMatrix(
596                        _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
597                        _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
598                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
599         }
600     }
601     else if (bScalar2)
602     {
603         //X ./ y
604         int iInc1       = 1;
605         int iInc2       = 0;
606         int iIncOut     = 1;
607         int iSizeResult = _pDouble1->getSize();
608
609         *_pDoubleOut    = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray(), bComplex1 || bComplex2);
610
611         if (bComplex1 == false && bComplex2 == false)
612         {
613             // r ./ R
614             iErr = iRightDivisionRealMatrixByRealMatrix(
615                        _pDouble1->getReal(), iInc1,
616                        _pDouble2->getReal(), iInc2,
617                        (*_pDoubleOut)->getReal(), iIncOut, iSizeResult);
618         }
619         else if (bComplex1 == false && bComplex2 == true)
620         {
621             // r ./ C
622             iErr = iRightDivisionRealMatrixByComplexMatrix(
623                        _pDouble1->getReal(), iInc1,
624                        _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
625                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
626         }
627         else if (bComplex1 == true && bComplex2 == false)
628         {
629             // c ./ R
630             iErr = iRightDivisionComplexMatrixByRealMatrix(
631                        _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
632                        _pDouble2->getReal(), iInc2,
633                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
634         }
635         else if (bComplex1 == true && bComplex2 == true)
636         {
637             // c ./ C
638             iErr = iRightDivisionComplexMatrixByComplexMatrix(
639                        _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
640                        _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
641                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
642         }
643     }
644     else
645     {
646         //X ./ Y
647         //check dimension compatibilities ( same number of dimension and same size for each dimension
648         int iDims1      = _pDouble1->getDims();
649         int* piDims1    = _pDouble1->getDimsArray();
650         int iDims2      = _pDouble2->getDims();
651         int* piDims2    = _pDouble2->getDimsArray();
652
653         if (iDims1 != iDims2)
654         {
655             return 1;
656         }
657
658         for (int i = 0 ; i < iDims1 ; i++)
659         {
660             if (piDims1[i] != piDims2[i])
661             {
662                 return 1;
663             }
664         }
665
666         (*_pDoubleOut) = new Double(iDims2, piDims2, bComplex1 || bComplex2);
667
668         int iInc1       = 1;
669         int iInc2       = 1;
670         int iIncOut     = 1;
671         int iSizeResult = _pDouble1->getSize();
672
673         if (bComplex1 == false && bComplex2 == false)
674         {
675             // r ./ R
676             iErr = iRightDivisionRealMatrixByRealMatrix(
677                        _pDouble1->getReal(), iInc1,
678                        _pDouble2->getReal(), iInc2,
679                        (*_pDoubleOut)->getReal(), iIncOut, iSizeResult);
680         }
681         else if (bComplex1 == false && bComplex2 == true)
682         {
683             // r ./ C
684             iErr = iRightDivisionRealMatrixByComplexMatrix(
685                        _pDouble1->getReal(), iInc1,
686                        _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
687                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
688         }
689         else if (bComplex1 == true && bComplex2 == false)
690         {
691             // c ./ R
692             iErr = iRightDivisionComplexMatrixByRealMatrix(
693                        _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
694                        _pDouble2->getReal(), iInc2,
695                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
696         }
697         else if (bComplex1 == true && bComplex2 == true)
698         {
699             // c ./ C
700             iErr = iRightDivisionComplexMatrixByComplexMatrix(
701                        _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
702                        _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
703                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
704         }
705     }
706     return iErr;
707 }
708
709 int DotRDividePolyByDouble(Polynom* _pPoly1, Double* _pDouble2, Polynom** _pPolyOut)
710 {
711     int iErr        = 0;
712     bool bComplex1  = _pPoly1->isComplex();
713     bool bComplex2  = _pDouble2->isComplex();
714
715     //X ./ Y
716     //check dimension compatibilities ( same number of dimension and same size for each dimension
717     int iDims1      = _pPoly1->getDims();
718     int* piDims1    = _pPoly1->getDimsArray();
719     int iMaxSize    = _pPoly1->getMaxRank() + 1;
720     int iSizePoly   = _pPoly1->getSize();
721     int iDims2      = _pDouble2->getDims();
722     int* piDims2    = _pDouble2->getDimsArray();
723
724     if (iDims1 != iDims2)
725     {
726         return 1;
727     }
728
729     for (int i = 0 ; i < iDims1 ; i++)
730     {
731         if (piDims1[i] != piDims2[i])
732         {
733             return 1;
734         }
735     }
736
737     // compute output ranks
738     int* piRanks = new int[iSizePoly];
739     for (int i = 0; i < iSizePoly; i++)
740     {
741         piRanks[i] = iMaxSize - 1;
742     }
743
744     // create output and working table
745     (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), iDims2, piDims2, piRanks);
746     delete[] piRanks;
747     Double* pDblCoefOut = new Double(_pPoly1->getRows(), _pPoly1->getCols() * iMaxSize, bComplex1 || bComplex2);
748     double* pdblCoef2   = new double[_pPoly1->getRows() * _pPoly1->getCols() * iMaxSize];
749     Double* pDblCoef1   = _pPoly1->getCoef();
750
751     int iZero = 0;
752     double* pdbl = _pDouble2->get();
753     for (int i = 0; i < iSizePoly; i++)
754     {
755         C2F(dcopy)(&iMaxSize, pdbl + i, &iZero, pdblCoef2 + i, &iSizePoly);
756     }
757
758     int iInc1       = 1;
759     int iInc2       = 1;
760     int iIncOut     = 1;
761     int iSizeResult = pDblCoefOut->getSize();
762
763     if (bComplex1 == false && bComplex2 == false)
764     {
765         // r ./ R
766         iErr = iRightDivisionRealMatrixByRealMatrix(
767                    pDblCoef1->getReal(), iInc1,
768                    pdblCoef2, iInc2,
769                    pDblCoefOut->getReal(), iIncOut, iSizeResult);
770     }
771     else if (bComplex1 == false && bComplex2 == true)
772     {
773         // r ./ C
774         //        iErr = iRightDivisionRealMatrixByComplexMatrix(
775         //                   _pDouble1->getReal(), iInc1,
776         //                   _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
777         //                   (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
778
779         // waiting for polynom rewrite about complex
780         iErr = 10;
781     }
782     else if (bComplex1 == true && bComplex2 == false)
783     {
784         // c ./ R
785         //        iErr = iRightDivisionComplexMatrixByRealMatrix(
786         //                   _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
787         //                   _pDouble2->getReal(), iInc2,
788         //                   (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
789
790         // waiting for polynom rewrite about complex
791         iErr = 10;
792     }
793     else if (bComplex1 == true && bComplex2 == true)
794     {
795         // c ./ C
796         //        iErr = iRightDivisionComplexMatrixByComplexMatrix(
797         //                   _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
798         //                   _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
799         //                   (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
800
801         // waiting for polynom rewrite about complex
802         iErr = 10;
803     }
804
805     (*_pPolyOut)->setCoef(pDblCoefOut);
806     (*_pPolyOut)->updateRank();
807
808     delete pDblCoefOut;
809     delete pDblCoef1;
810     delete[] pdblCoef2;
811
812     return iErr;
813 }