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