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