Coverity: ast module resource leaks fixed
[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  * Copyright (C) 2012 - 2016 - Scilab Enterprises
6  *
7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
9  * This file was originally licensed under the terms of the CeCILL v2.1,
10  * and continues to be available under such terms.
11  * For more information, see the COPYING file which you should have received
12  * along with this program.
13  *
14  */
15
16 #include "types_divide.hxx"
17 #include "types_finite.hxx"
18
19 extern "C"
20 {
21 #include "matrix_right_division.h"
22 #include "sciprint.h"
23 #include "localization.h"
24 #include "charEncoding.h"
25 #include "configvariable_interface.h"
26 #include "elem_common.h"
27 }
28
29 using namespace types;
30
31 InternalType *GenericRDivide(InternalType *_pLeftOperand, InternalType *_pRightOperand)
32 {
33     InternalType *pResult       = NULL;
34     GenericType::ScilabType TypeL = _pLeftOperand->getType();
35     GenericType::ScilabType TypeR = _pRightOperand->getType();
36
37     int iResult = 0;
38
39     if (_pLeftOperand->isDouble() && _pLeftOperand->getAs<Double>()->isEmpty())
40     {
41         return Double::Empty();
42     }
43
44     if (_pRightOperand->isDouble() && _pRightOperand->getAs<Double>()->isEmpty())
45     {
46         return Double::Empty();
47     }
48
49     /*
50     ** DOUBLE / DOUBLE
51     */
52     if (TypeL == GenericType::ScilabDouble && TypeR == GenericType::ScilabDouble)
53     {
54         Double *pL  = _pLeftOperand->getAs<Double>();
55         Double *pR  = _pRightOperand->getAs<Double>();
56
57         iResult = RDivideDoubleByDouble(pL, pR, (Double**)&pResult);
58     }
59
60     /*
61     ** POLY / DOUBLE
62     */
63     else if (TypeL == GenericType::ScilabPolynom && TypeR == GenericType::ScilabDouble)
64     {
65         Polynom *pL = _pLeftOperand->getAs<types::Polynom>();
66         Double *pR  = _pRightOperand->getAs<Double>();
67
68         iResult = RDividePolyByDouble(pL, pR, (Polynom**)&pResult);
69     }
70
71     /*
72     ** DOUBLE / POLY
73     */
74     else if (TypeL == GenericType::ScilabDouble && TypeR == GenericType::ScilabPolynom)
75     {
76         Double *pL  = _pLeftOperand->getAs<Double>();
77         Polynom *pR = _pRightOperand->getAs<types::Polynom>();
78
79         iResult = RDivideDoubleByPoly(pL, pR, (Polynom**)&pResult);
80     }
81
82     /*
83     ** SPARSE / DOUBLE
84     */
85     else if (TypeL == GenericType::ScilabSparse && TypeR == GenericType::ScilabDouble)
86     {
87         Sparse *pL = _pLeftOperand->getAs<types::Sparse>();
88         Double *pR = _pRightOperand->getAs<Double>();
89
90         iResult = RDivideSparseByDouble(pL, pR, &pResult);
91     }
92
93     //manage errors
94     if (iResult)
95     {
96         switch (iResult)
97         {
98             case 1 :
99                 throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
100             case 2 :
101                 throw ast::InternalError(_W("With NaN or Inf a division by scalar expected.\n"));
102             case 3 :
103                 throw ast::InternalError(_W("Division by zero...\n"));
104             case 4 :
105                 if (getWarningMode())
106                 {
107                     sciprint(_("Warning : Division by zero...\n"));
108                 }
109                 break;
110             //            default : throw ast::InternalError(_W("Operator / : Error %d not yet managed.\n"), iResult);
111             default :
112                 sciprint(_("Operator / : Error %d not yet managed.\n"), iResult);
113         }
114     }
115
116     /*
117     ** Default case : Return NULL will Call Overloading.
118     */
119     return pResult;
120 }
121
122
123 int RDivideDoubleByDouble(Double *_pDouble1, Double *_pDouble2, Double **_pDoubleOut)
124 {
125     int iErr = 0;
126
127     //check finite values of _pDouble1 and _pDouble2
128     if (isDoubleFinite(_pDouble1) == false || isDoubleFinite(_pDouble2) == false)
129     {
130         if (_pDouble2->isScalar() == false)
131         {
132             return 2;
133         }
134     }
135
136     if (_pDouble2->isScalar())
137     {
138         //Y / x
139         int iInc1       = 1;
140         int iInc2       = 0;
141         bool bComplex1  = _pDouble1->isComplex();
142         bool bComplex2  = _pDouble2->isComplex();
143
144         *_pDoubleOut    = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray(), bComplex1 || bComplex2);
145
146         if (bComplex1 == false && bComplex2 == false)
147         {
148             // Real1 \ Real2 -> Real2 / Real1
149             iErr = iRightDivisionRealMatrixByRealMatrix(
150                        _pDouble1->get(), iInc1,
151                        _pDouble2->get(), iInc2,
152                        (*_pDoubleOut)->get(), 1, _pDouble1->getSize());
153         }
154         else if (bComplex1 == false && bComplex2 == true)
155         {
156             // Real \ Complex -> Complex / Real
157             iErr = iRightDivisionRealMatrixByComplexMatrix(
158                        _pDouble1->get(), iInc1,
159                        _pDouble2->get(), _pDouble2->getImg(), iInc2,
160                        (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), 1, _pDouble1->getSize());
161         }
162         else if (bComplex1 == true && bComplex2 == false)
163         {
164             // Complex \ Real -> Real / Complex
165             iErr = iRightDivisionComplexMatrixByRealMatrix(
166                        _pDouble1->get(), _pDouble1->getImg(), iInc1,
167                        _pDouble2->get(), iInc2,
168                        (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), 1, _pDouble1->getSize());
169         }
170         else if (bComplex1 == true && bComplex2 == true)
171         {
172             // Complex \ Complex
173             iErr = iRightDivisionComplexMatrixByComplexMatrix(
174                        _pDouble1->get(), _pDouble1->getImg(), iInc1,
175                        _pDouble2->get(), _pDouble2->getImg(), iInc2,
176                        (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), 1, _pDouble1->getSize());
177         }
178
179         return iErr;
180     }
181
182     if (_pDouble1->isScalar())
183     {
184         if (_pDouble2->getDims() > 2)
185         {
186             //not managed, call overload
187             return 0;
188         }
189
190         // x / eye() = x
191         if (_pDouble2->isIdentity() )
192         {
193             *_pDoubleOut    = new Double(*_pDouble1);
194             return 0;
195         }
196         double dblSavedR = 0;
197         double dblSavedI = 0;
198         Double *pdblTemp = NULL;
199
200         int iRowResult = _pDouble2->getCols();
201         int iColResult = _pDouble2->getRows();
202
203
204
205
206         //in this case, we have to create a temporary square matrix
207         pdblTemp = new Double(iRowResult, iRowResult, _pDouble1->isComplex());
208         pdblTemp->setZeros();
209
210         if (_pDouble1->isComplex())
211         {
212             dblSavedR = _pDouble1->getReal()[0];
213             dblSavedI = _pDouble1->getImg()[0];
214             for (int i = 0 ; i < iRowResult ; i++)
215             {
216                 pdblTemp->set(i, i, dblSavedR);
217                 pdblTemp->setImg(i, i, dblSavedI);
218             }
219         }
220         else
221         {
222             dblSavedR = _pDouble1->getReal()[0];
223             for (int i = 0 ; i < iRowResult ; i++)
224             {
225                 pdblTemp->set(i, i, dblSavedR);
226             }
227         }
228
229         *_pDoubleOut = new Double(iRowResult, iColResult, _pDouble1->isComplex() || _pDouble2->isComplex());
230
231         if ((*_pDoubleOut)->isComplex())
232         {
233             double dblRcond = 0;
234             iErr = iRightDivisionOfComplexMatrix(
235                        pdblTemp->getReal(), pdblTemp->getImg(), pdblTemp->getRows(), pdblTemp->getCols(),
236                        _pDouble2->getReal(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
237                        (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), iRowResult, iColResult, &dblRcond);
238         }
239         else
240         {
241             double dblRcond = 0;
242             iErr = iRightDivisionOfRealMatrix(
243                        pdblTemp->getReal(), pdblTemp->getRows(), pdblTemp->getCols(),
244                        _pDouble2->getReal(), _pDouble2->getRows(), _pDouble2->getCols(),
245                        (*_pDoubleOut)->getReal(), iRowResult, iColResult, &dblRcond);
246         }
247         delete pdblTemp;
248         return iErr;
249     }
250
251     if (_pDouble1->getDims() > 2 || _pDouble2->getDims() > 2 || _pDouble1->getCols() != _pDouble2->getCols())
252     {
253         //not managed
254         return 0;
255     }
256
257     *_pDoubleOut = new Double(_pDouble1->getRows(), _pDouble2->getRows(), _pDouble1->isComplex() || _pDouble2->isComplex());
258     if ((*_pDoubleOut)->isComplex())
259     {
260         double dblRcond = 0;
261         iErr = iRightDivisionOfComplexMatrix(
262                    _pDouble1->getReal(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
263                    _pDouble2->getReal(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
264                    (*_pDoubleOut)->getReal(), (*_pDoubleOut)->getImg(), _pDouble1->getRows(), _pDouble2->getRows(), &dblRcond);
265     }
266     else
267     {
268         double dblRcond = 0;
269         iErr = iRightDivisionOfRealMatrix(
270                    _pDouble1->getReal(), _pDouble1->getRows(), _pDouble1->getCols(),
271                    _pDouble2->getReal(), _pDouble2->getRows(), _pDouble2->getCols(),
272                    (*_pDoubleOut)->getReal(), _pDouble1->getRows(), _pDouble2->getRows(), &dblRcond);
273     }
274
275     return iErr;
276 }
277
278 int RDividePolyByDouble(Polynom* _pPoly, Double* _pDouble, Polynom** _pPolyOut)
279 {
280     bool bComplex1  = _pPoly->isComplex();
281     bool bComplex2  = _pDouble->isComplex();
282     bool bScalar1   = _pPoly->getRows() == 1  && _pPoly->getCols() == 1;
283     bool bScalar2   = _pDouble->getRows() == 1 && _pDouble->getCols() == 1;
284
285     Polynom *pTemp = NULL; //use only if _pPoly is scalar and _pDouble not.
286
287     int iRowResult  = 0;
288     int iColResult = 0;
289     int *piRank   = NULL;
290
291     /* if(bScalar1 && bScalar2)
292     {
293     iRowResult = 1;
294     iColResult = 1;
295
296     piRank = new int[1];
297     piRank[0] = _pPoly->get(0)->getRank();
298     }
299     else */
300
301     if (bScalar1 == false && bScalar2 == false)
302     {
303         // call overload
304         return 0;
305     }
306
307     if (bScalar2)
308     {
309         double dblDivR = _pDouble->get(0);
310         double dblDivI = _pDouble->getImg(0);
311
312         (*_pPolyOut) = _pPoly->clone()->getAs<Polynom>();
313         if (_pDouble->isComplex())
314         {
315             (*_pPolyOut)->setComplex(true);
316         }
317
318         for (int i = 0 ; i < _pPoly->getSize() ; i++)
319         {
320             bool bComplex1 = _pPoly->isComplex();
321             bool bComplex2 = _pDouble->isComplex();
322
323             SinglePoly* pC = (*_pPolyOut)->get(i);
324
325             if (bComplex1 == false && bComplex2 == false)
326             {
327                 iRightDivisionRealMatrixByRealMatrix(pC->get(), 1, &dblDivR, 0, pC->get(), 1, pC->getSize());
328             }
329             else if (bComplex1 == true && bComplex2 == false)
330             {
331                 iRightDivisionComplexMatrixByRealMatrix(pC->get(), pC->getImg(), 1, &dblDivR, 0, pC->get(), pC->getImg(), 1, pC->getSize());
332             }
333             else if (bComplex1 == false && bComplex2 == true)
334             {
335                 iRightDivisionRealMatrixByComplexMatrix(pC->get(), 1, &dblDivR, &dblDivI, 0, pC->get(), pC->getImg(), 1, pC->getSize());
336             }
337             else if (bComplex1 == true && bComplex2 == true)
338             {
339                 iRightDivisionComplexMatrixByComplexMatrix(pC->get(), pC->getImg(), 1, &dblDivR, &dblDivI, 0, pC->get(), pC->getImg(), 1, pC->getSize());
340             }
341         }
342
343         return 0;
344     }
345
346     if (bScalar1)
347     {
348         //in this case, we have to create a temporary square polinomial matrix
349         iRowResult = _pDouble->getCols();
350         iColResult = _pDouble->getRows();
351
352         piRank = new int[iRowResult * iRowResult];
353         int iMaxRank = _pPoly->getMaxRank();
354         for (int i = 0 ; i < iRowResult * iRowResult ; i++)
355         {
356             piRank[i] = iMaxRank;
357         }
358
359         pTemp = new Polynom(_pPoly->getVariableName(), iRowResult, iRowResult, piRank);
360         if (bComplex1 || bComplex2)
361         {
362             pTemp->setComplex(true);
363         }
364
365         SinglePoly *pdblData = _pPoly->get(0);
366         for (int i = 0 ; i < iRowResult ; i++)
367         {
368             pTemp->set(i, i, pdblData);
369         }
370     }
371
372     (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), iRowResult, iColResult, piRank);
373     delete[] piRank;
374     if (bComplex1 || bComplex2)
375     {
376         (*_pPolyOut)->setComplex(true);
377     }
378
379     if (bScalar2)
380     {
381         //[p] * cst
382         for (int i = 0 ; i < _pPoly->getSize() ; i++)
383         {
384             SinglePoly *pPolyIn   = _pPoly->get(i);
385             double* pRealIn  = pPolyIn->get();
386             double* pImgIn  = pPolyIn->getImg();
387
388             SinglePoly *pPolyOut  = (*_pPolyOut)->get(i);
389             double* pRealOut = pPolyOut->get();
390             double* pImgOut  = pPolyOut->getImg();
391
392             if (bComplex1 == false && bComplex2 == false)
393             {
394                 iRightDivisionRealMatrixByRealMatrix(
395                     pRealIn, 1,
396                     _pDouble->getReal(), 0,
397                     pRealOut, 1, pPolyOut->getSize());
398             }
399             else if (bComplex1 == false && bComplex2 == true)
400             {
401                 iRightDivisionRealMatrixByComplexMatrix(
402                     pRealIn, 1,
403                     _pDouble->getReal(), _pDouble->getImg(), 0,
404                     pRealOut, pImgOut, 1, pPolyOut->getSize());
405             }
406             else if (bComplex1 == true && bComplex2 == false)
407             {
408                 iRightDivisionComplexMatrixByRealMatrix(
409                     pRealIn, pImgIn, 1,
410                     _pDouble->getReal(), 0,
411                     pRealOut, pImgOut, 1, pPolyOut->getSize());
412             }
413             else if (bComplex1 == true && bComplex2 == true)
414             {
415                 iRightDivisionComplexMatrixByComplexMatrix(
416                     pRealIn, pImgIn, 1,
417                     _pDouble->getReal(), _pDouble->getImg(), 0,
418                     pRealOut, pImgOut, 1, pPolyOut->getSize());
419             }
420         }
421     }
422     else if (bScalar1)
423     {
424         for (int i = 0 ; i < pTemp->get(0)->getSize() ; i++)
425         {
426             Double *pCoef    = pTemp->extractCoef(i);
427             Double *pResultCoef = new Double(iRowResult, iColResult, pCoef->isComplex());
428             double *pReal    = pResultCoef->getReal();
429             double *pImg    = pResultCoef->getImg();
430
431             if (bComplex1 == false && bComplex2 == false)
432             {
433                 double dblRcond = 0;
434                 iRightDivisionOfRealMatrix(
435                     pCoef->getReal(), iRowResult, iRowResult,
436                     _pDouble->getReal(), _pDouble->getRows(), _pDouble->getCols(),
437                     pReal, iRowResult, iColResult, &dblRcond);
438             }
439             else
440             {
441                 double dblRcond = 0;
442                 iRightDivisionOfComplexMatrix(
443                     pCoef->getReal(), pCoef->getImg(), iRowResult, iRowResult,
444                     _pDouble->getReal(), _pDouble->getImg(), _pDouble->getRows(), _pDouble->getCols(),
445                     pReal, pImg, iRowResult, iColResult, &dblRcond);
446             }
447
448             (*_pPolyOut)->insertCoef(i, pResultCoef);
449             delete pCoef;
450             delete pResultCoef;
451         }
452     }
453     return 0;
454 }
455
456 int RDivideDoubleByPoly(Double* /*_pDouble*/, Polynom* /*_pPoly*/, Polynom** /*_pPolyOut*/)
457 {
458     return 0;
459 }
460
461 int RDivideSparseByDouble(types::Sparse* _pSp, types::Double* _pDouble, InternalType** _pSpOut)
462 {
463     if (_pDouble->isEmpty())
464     {
465         //sp / []
466         *_pSpOut = Double::Empty();
467         return 0;
468     }
469
470     if (_pDouble->isIdentity())
471     {
472         *_pSpOut = new Sparse(*_pSp);
473         return 0;
474     }
475
476     size_t iSize = _pSp->nonZeros();
477     int* Col = new int[iSize];
478     int* Row = new int[_pSp->getRows()];
479     _pSp->getColPos(Col);
480     _pSp->getNbItemByRow(Row);
481     int* iPositVal = new int[iSize];
482
483     int idx = 0;
484     for (int i = 0; i < _pSp->getRows(); i++)
485     {
486         for (int j = 0; j < Row[i]; j++)
487         {
488             iPositVal[idx] = (Col[idx] - 1) * _pSp->getRows() + i;
489             ++idx;
490         }
491     }
492
493     delete[] Col;
494     delete[] Row;
495
496     Double** pDbl = new Double*[iSize];
497     Double** pDblSp = new Double*[iSize];
498     double* pdbl = _pDouble->get();
499
500     if (_pDouble->isScalar())
501     {
502         if (_pDouble->isComplex())
503         {
504             double* pdblImg = _pDouble->getImg();
505             for (int i = 0; i < iSize; i++)
506             {
507                 pDbl[i] = new Double(pdbl[0], pdblImg[0]);
508                 pDblSp[i] = new Double(_pSp->get(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
509             }
510         }
511         else
512         {
513             for (int i = 0; i < iSize; i++)
514             {
515                 pDbl[i] = new Double(pdbl[0]);
516                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
517             }
518         }
519     }
520     else if (_pDouble->getSize() == iSize)
521     {
522         if (_pDouble->isComplex())
523         {
524             double* pdblImg = _pDouble->getImg();
525             for (int i = 0; i < iSize; i++)
526             {
527                 pDbl[i] = new Double(pdbl[i], pdblImg[i]);
528                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
529             }
530         }
531         else
532         {
533             for (int i = 0; i < iSize; i++)
534             {
535                 pDbl[i] = new Double(pdbl[i]);
536                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
537             }
538         }
539     }
540     else
541     {
542         delete[] pDbl;
543         delete[] pDblSp;
544         delete[] iPositVal;
545         return 0;
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                 for (int i = 0; i < iSize; ++i)
561                 {
562                     delete pDbl[i];
563                     delete pDblSp[i];
564                 }
565                 delete[] pDbl;
566                 delete[] pDblSp;
567                 delete[] iPositVal;
568                 delete pSpTemp;
569                 delete ppDblGet;
570                 return iResultat;
571             }
572             std::complex<double> cplx(ppDblGet->get(0), ppDblGet->getImg(0));
573             pSpTemp->set(iPositVal[i], cplx, false);
574             delete ppDblGet;
575         }
576     }
577
578     pSpTemp->finalize();
579
580     delete[] iPositVal;
581
582     for (int i = 0; i < iSize; ++i)
583     {
584         delete pDbl[i];
585         delete pDblSp[i];
586     }
587
588     delete[] pDbl;
589     delete[] pDblSp;
590
591     *_pSpOut = pSpTemp;
592     return 0;
593 }