9ecbeb395c052987e65fa6e60855ea50f9e2d19c
[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             Double* pC = (*_pPolyOut)->get(i)->getCoef();
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         for (int i = 0 ; i < iRowResult * iRowResult ; i++)
399         {
400             piRank[i] = _pPoly->getMaxRank();
401         }
402
403         pTemp = new Polynom(_pPoly->getVariableName(), iRowResult, iRowResult, piRank);
404         if (bComplex1 || bComplex2)
405         {
406             pTemp->setComplex(true);
407         }
408
409         Double *pdblData = _pPoly->get(0)->getCoef();
410         for (int i = 0 ; i < iRowResult ; i++)
411         {
412             pTemp->setCoef(i, i, pdblData);
413         }
414     }
415
416     (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), iRowResult, iColResult, piRank);
417     delete[] piRank;
418     if (bComplex1 || bComplex2)
419     {
420         (*_pPolyOut)->setComplex(true);
421     }
422
423     if (bScalar2)
424     {
425         //[p] * cst
426         for (int i = 0 ; i < _pPoly->getSize() ; i++)
427         {
428             SinglePoly *pPolyIn   = _pPoly->get(i);
429             double* pRealIn  = pPolyIn->getCoef()->getReal();
430             double* pImgIn  = pPolyIn->getCoef()->getImg();
431
432             SinglePoly *pPolyOut  = (*_pPolyOut)->get(i);
433             double* pRealOut = pPolyOut->getCoef()->getReal();
434             double* pImgOut  = pPolyOut->getCoef()->getImg();
435
436             if (bComplex1 == false && bComplex2 == false)
437             {
438                 iRightDivisionRealMatrixByRealMatrix(
439                     pRealIn, 1,
440                     _pDouble->getReal(), 0,
441                     pRealOut, 1, pPolyOut->getRank());
442             }
443             else if (bComplex1 == false && bComplex2 == true)
444             {
445                 iRightDivisionRealMatrixByComplexMatrix(
446                     pRealIn, 1,
447                     _pDouble->getReal(), _pDouble->getImg(), 0,
448                     pRealOut, pImgOut, 1, pPolyOut->getRank());
449             }
450             else if (bComplex1 == true && bComplex2 == false)
451             {
452                 iRightDivisionComplexMatrixByRealMatrix(
453                     pRealIn, pImgIn, 1,
454                     _pDouble->getReal(), 0,
455                     pRealOut, pImgOut, 1, pPolyOut->getRank());
456             }
457             else if (bComplex1 == true && bComplex2 == true)
458             {
459                 iRightDivisionComplexMatrixByComplexMatrix(
460                     pRealIn, pImgIn, 1,
461                     _pDouble->getReal(), _pDouble->getImg(), 0,
462                     pRealOut, pImgOut, 1, pPolyOut->getRank());
463             }
464         }
465     }
466     else if (bScalar1)
467     {
468         Double *pResultCoef = new Double(iRowResult, iColResult, (*_pPolyOut)->isComplex());
469         double *pReal    = pResultCoef->getReal();
470         double *pImg    = pResultCoef->getImg();
471
472         for (int i = 0 ; i < pTemp->get(0)->getRank() ; i++)
473         {
474             Double *pCoef    = pTemp->extractCoef(i);
475             Double *pResultCoef = new Double(iRowResult, iColResult, pCoef->isComplex());
476             double *pReal    = pResultCoef->getReal();
477             double *pImg    = pResultCoef->getImg();
478
479             if (bComplex1 == false && bComplex2 == false)
480             {
481                 double dblRcond = 0;
482                 iRightDivisionOfRealMatrix(
483                     pCoef->getReal(), iRowResult, iRowResult,
484                     _pDouble->getReal(), _pDouble->getRows(), _pDouble->getCols(),
485                     pReal, iRowResult, iColResult, &dblRcond);
486             }
487             else
488             {
489                 double dblRcond = 0;
490                 iRightDivisionOfComplexMatrix(
491                     pCoef->getReal(), pCoef->getImg(), iRowResult, iRowResult,
492                     _pDouble->getReal(), _pDouble->getImg(), _pDouble->getRows(), _pDouble->getCols(),
493                     pReal, pImg, iRowResult, iColResult, &dblRcond);
494             }
495
496             (*_pPolyOut)->insertCoef(i, pResultCoef);
497             delete pCoef;
498         }
499         delete pResultCoef;
500     }
501     return 0;
502 }
503
504 int RDivideDoubleByPoly(Double* _pDouble, Polynom* _pPoly, Polynom** _pPolyOut)
505 {
506     return 0;
507 }
508
509 int DotRDivideDoubleByDouble(Double* _pDouble1, Double* _pDouble2, Double** _pDoubleOut)
510 {
511     int iErr        = 0;
512     bool bComplex1  = _pDouble1->isComplex();
513     bool bComplex2  = _pDouble2->isComplex();
514     bool bScalar1   = _pDouble1->isScalar();
515     bool bScalar2   = _pDouble2->isScalar();
516
517     if (bScalar1)
518     {
519         //x ./ Y
520         int iInc1       = 0;
521         int iInc2       = 1;
522         int iIncOut     = 1;
523         int iSizeResult = _pDouble2->getSize();
524
525         *_pDoubleOut    = new Double(_pDouble2->getDims(), _pDouble2->getDimsArray(), bComplex1 || bComplex2);
526
527         if (bComplex1 == false && bComplex2 == false)
528         {
529             // r ./ R
530             iErr = iRightDivisionRealMatrixByRealMatrix(
531                        _pDouble1->getReal(), iInc1,
532                        _pDouble2->getReal(), iInc2,
533                        (*_pDoubleOut)->getReal(), iIncOut, iSizeResult);
534         }
535         else if (bComplex1 == false && bComplex2 == true)
536         {
537             // r ./ C
538             iErr = iRightDivisionRealMatrixByComplexMatrix(
539                        _pDouble1->getReal(), iInc1,
540                        _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
541                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
542         }
543         else if (bComplex1 == true && bComplex2 == false)
544         {
545             // c ./ R
546             iErr = iRightDivisionComplexMatrixByRealMatrix(
547                        _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
548                        _pDouble2->getReal(), iInc2,
549                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
550         }
551         else if (bComplex1 == true && bComplex2 == true)
552         {
553             // c ./ C
554             iErr = iRightDivisionComplexMatrixByComplexMatrix(
555                        _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
556                        _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
557                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
558         }
559     }
560     else if (bScalar2)
561     {
562         //X ./ y
563         int iInc1       = 1;
564         int iInc2       = 0;
565         int iIncOut     = 1;
566         int iSizeResult = _pDouble1->getSize();
567
568         *_pDoubleOut    = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray(), bComplex1 || bComplex2);
569
570         if (bComplex1 == false && bComplex2 == false)
571         {
572             // r ./ R
573             iErr = iRightDivisionRealMatrixByRealMatrix(
574                        _pDouble1->getReal(), iInc1,
575                        _pDouble2->getReal(), iInc2,
576                        (*_pDoubleOut)->getReal(), iIncOut, iSizeResult);
577         }
578         else if (bComplex1 == false && bComplex2 == true)
579         {
580             // r ./ C
581             iErr = iRightDivisionRealMatrixByComplexMatrix(
582                        _pDouble1->getReal(), iInc1,
583                        _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
584                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
585         }
586         else if (bComplex1 == true && bComplex2 == false)
587         {
588             // c ./ R
589             iErr = iRightDivisionComplexMatrixByRealMatrix(
590                        _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
591                        _pDouble2->getReal(), iInc2,
592                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
593         }
594         else if (bComplex1 == true && bComplex2 == true)
595         {
596             // c ./ C
597             iErr = iRightDivisionComplexMatrixByComplexMatrix(
598                        _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
599                        _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
600                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
601         }
602     }
603     else
604     {
605         //X ./ Y
606         //check dimension compatibilities ( same number of dimension and same size for each dimension
607         int iDims1      = _pDouble1->getDims();
608         int* piDims1    = _pDouble1->getDimsArray();
609         int iDims2      = _pDouble2->getDims();
610         int* piDims2    = _pDouble2->getDimsArray();
611
612         if (iDims1 != iDims2)
613         {
614             return 1;
615         }
616
617         for (int i = 0 ; i < iDims1 ; i++)
618         {
619             if (piDims1[i] != piDims2[i])
620             {
621                 return 1;
622             }
623         }
624
625         (*_pDoubleOut) = new Double(iDims2, piDims2, bComplex1 || bComplex2);
626
627         int iErr        = 0;
628         int iInc1       = 1;
629         int iInc2       = 1;
630         int iIncOut     = 1;
631         int iSizeResult = _pDouble1->getSize();
632
633         if (bComplex1 == false && bComplex2 == false)
634         {
635             // r ./ R
636             iErr = iRightDivisionRealMatrixByRealMatrix(
637                        _pDouble1->getReal(), iInc1,
638                        _pDouble2->getReal(), iInc2,
639                        (*_pDoubleOut)->getReal(), iIncOut, iSizeResult);
640         }
641         else if (bComplex1 == false && bComplex2 == true)
642         {
643             // r ./ C
644             iErr = iRightDivisionRealMatrixByComplexMatrix(
645                        _pDouble1->getReal(), iInc1,
646                        _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
647                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
648         }
649         else if (bComplex1 == true && bComplex2 == false)
650         {
651             // c ./ R
652             iErr = iRightDivisionComplexMatrixByRealMatrix(
653                        _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
654                        _pDouble2->getReal(), iInc2,
655                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
656         }
657         else if (bComplex1 == true && bComplex2 == true)
658         {
659             // c ./ C
660             iErr = iRightDivisionComplexMatrixByComplexMatrix(
661                        _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
662                        _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
663                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
664         }
665     }
666     return iErr;
667 }
668
669 int DotRDividePolyByDouble(Polynom* _pPoly1, Double* _pDouble2, Polynom** _pPolyOut)
670 {
671     int iErr        = 0;
672     bool bComplex1  = _pPoly1->isComplex();
673     bool bComplex2  = _pDouble2->isComplex();
674
675     //X ./ Y
676     //check dimension compatibilities ( same number of dimension and same size for each dimension
677     int iDims1      = _pPoly1->getDims();
678     int* piDims1    = _pPoly1->getDimsArray();
679     int iMaxRank    = _pPoly1->getMaxRank();
680     int iSizePoly   = _pPoly1->getSize();
681     int iDims2      = _pDouble2->getDims();
682     int* piDims2    = _pDouble2->getDimsArray();
683
684     if (iDims1 != iDims2)
685     {
686         return 1;
687     }
688
689     for (int i = 0 ; i < iDims1 ; i++)
690     {
691         if (piDims1[i] != piDims2[i])
692         {
693             return 1;
694         }
695     }
696
697     // compute output ranks
698     int* piRanks = new int[iSizePoly];
699     for (int i = 0; i < iSizePoly; i++)
700     {
701         piRanks[i] = iMaxRank;
702     }
703
704     // create output and working table
705     (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), iDims2, piDims2, piRanks);
706     delete[] piRanks;
707     Double* pDblCoefOut = new Double(_pPoly1->getRows(), _pPoly1->getCols() * iMaxRank, bComplex1 || bComplex2);
708     double* pdblCoef2   = new double[_pPoly1->getRows() * _pPoly1->getCols() * iMaxRank];
709     Double* pDblCoef1   = _pPoly1->getCoef();
710
711     int iZero = 0;
712     double* pdbl = _pDouble2->get();
713     for (int i = 0; i < iSizePoly; i++)
714     {
715         C2F(dcopy)(&iMaxRank, pdbl + i, &iZero, pdblCoef2 + i, &iSizePoly);
716     }
717
718     int iInc1       = 1;
719     int iInc2       = 1;
720     int iIncOut     = 1;
721     int iSizeResult = pDblCoefOut->getSize();
722
723     if (bComplex1 == false && bComplex2 == false)
724     {
725         // r ./ R
726         iErr = iRightDivisionRealMatrixByRealMatrix(
727                    pDblCoef1->getReal(), iInc1,
728                    pdblCoef2, iInc2,
729                    pDblCoefOut->getReal(), iIncOut, iSizeResult);
730     }
731     else if (bComplex1 == false && bComplex2 == true)
732     {
733         // r ./ C
734         //        iErr = iRightDivisionRealMatrixByComplexMatrix(
735         //                   _pDouble1->getReal(), iInc1,
736         //                   _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
737         //                   (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
738
739         // waiting for polynom rewrite about complex
740         iErr = 10;
741     }
742     else if (bComplex1 == true && bComplex2 == false)
743     {
744         // c ./ R
745         //        iErr = iRightDivisionComplexMatrixByRealMatrix(
746         //                   _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
747         //                   _pDouble2->getReal(), iInc2,
748         //                   (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
749
750         // waiting for polynom rewrite about complex
751         iErr = 10;
752     }
753     else if (bComplex1 == true && bComplex2 == true)
754     {
755         // c ./ C
756         //        iErr = iRightDivisionComplexMatrixByComplexMatrix(
757         //                   _pDouble1->getReal(), _pDouble1->getImg(), iInc1,
758         //                   _pDouble2->getReal(), _pDouble2->getImg(), iInc2,
759         //                   (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iIncOut, iSizeResult);
760
761         // waiting for polynom rewrite about complex
762         iErr = 10;
763     }
764
765     (*_pPolyOut)->setCoef(pDblCoefOut);
766     (*_pPolyOut)->updateRank();
767
768     delete pDblCoefOut;
769     delete pDblCoef1;
770     delete[] pdblCoef2;
771
772     return iErr;
773 }