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