Typo fixes
[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         return 0;
343     }
344
345     if (bScalar1)
346     {
347         //in this case, we have to create a temporary square polynomial matrix
348         iRowResult = _pDouble->getCols();
349         iColResult = _pDouble->getRows();
350
351         piRank = new int[iRowResult * iRowResult];
352         int iMaxRank = _pPoly->getMaxRank();
353         for (int i = 0 ; i < iRowResult * iRowResult ; i++)
354         {
355             piRank[i] = iMaxRank;
356         }
357
358         pTemp = new Polynom(_pPoly->getVariableName(), iRowResult, iRowResult, piRank);
359         if (bComplex1 || bComplex2)
360         {
361             pTemp->setComplex(true);
362         }
363
364         SinglePoly *pdblData = _pPoly->get(0);
365         for (int i = 0 ; i < iRowResult ; i++)
366         {
367             pTemp->set(i, i, pdblData);
368         }
369     }
370
371     (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), iRowResult, iColResult, piRank);
372     delete[] piRank;
373     if (bComplex1 || bComplex2)
374     {
375         (*_pPolyOut)->setComplex(true);
376     }
377
378     if (bScalar1)
379     {
380         for (int i = 0 ; i < pTemp->get(0)->getSize() ; i++)
381         {
382             Double *pCoef    = pTemp->extractCoef(i);
383             Double *pResultCoef = new Double(iRowResult, iColResult, pCoef->isComplex());
384             double *pReal    = pResultCoef->getReal();
385             double *pImg    = pResultCoef->getImg();
386
387             if (bComplex1 == false && bComplex2 == false)
388             {
389                 double dblRcond = 0;
390                 iRightDivisionOfRealMatrix(
391                     pCoef->getReal(), iRowResult, iRowResult,
392                     _pDouble->getReal(), _pDouble->getRows(), _pDouble->getCols(),
393                     pReal, iRowResult, iColResult, &dblRcond);
394             }
395             else
396             {
397                 double dblRcond = 0;
398                 iRightDivisionOfComplexMatrix(
399                     pCoef->getReal(), pCoef->getImg(), iRowResult, iRowResult,
400                     _pDouble->getReal(), _pDouble->getImg(), _pDouble->getRows(), _pDouble->getCols(),
401                     pReal, pImg, iRowResult, iColResult, &dblRcond);
402             }
403
404             (*_pPolyOut)->insertCoef(i, pResultCoef);
405             delete pCoef;
406             delete pResultCoef;
407         }
408     }
409     return 0;
410 }
411
412 int RDivideDoubleByPoly(Double* /*_pDouble*/, Polynom* /*_pPoly*/, Polynom** /*_pPolyOut*/)
413 {
414     return 0;
415 }
416
417 int RDivideSparseByDouble(types::Sparse* _pSp, types::Double* _pDouble, InternalType** _pSpOut)
418 {
419     if (_pDouble->isEmpty())
420     {
421         //sp / []
422         *_pSpOut = Double::Empty();
423         return 0;
424     }
425
426     if (_pDouble->isIdentity())
427     {
428         *_pSpOut = new Sparse(*_pSp);
429         return 0;
430     }
431
432     size_t iSize = _pSp->nonZeros();
433     int* Col = new int[iSize];
434     int* Row = new int[_pSp->getRows()];
435     _pSp->getColPos(Col);
436     _pSp->getNbItemByRow(Row);
437     int* iPositVal = new int[iSize];
438
439     int idx = 0;
440     for (int i = 0; i < _pSp->getRows(); i++)
441     {
442         for (int j = 0; j < Row[i]; j++)
443         {
444             iPositVal[idx] = (Col[idx] - 1) * _pSp->getRows() + i;
445             ++idx;
446         }
447     }
448
449     delete[] Col;
450     delete[] Row;
451
452     Double** pDbl = new Double*[iSize];
453     Double** pDblSp = new Double*[iSize];
454     double* pdbl = _pDouble->get();
455
456     if (_pDouble->isScalar())
457     {
458         if (_pDouble->isComplex())
459         {
460             double* pdblImg = _pDouble->getImg();
461             for (int i = 0; i < iSize; i++)
462             {
463                 pDbl[i] = new Double(pdbl[0], pdblImg[0]);
464                 pDblSp[i] = new Double(_pSp->get(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
465             }
466         }
467         else
468         {
469             for (int i = 0; i < iSize; i++)
470             {
471                 pDbl[i] = new Double(pdbl[0]);
472                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
473             }
474         }
475     }
476     else if (_pDouble->getSize() == iSize)
477     {
478         if (_pDouble->isComplex())
479         {
480             double* pdblImg = _pDouble->getImg();
481             for (int i = 0; i < iSize; i++)
482             {
483                 pDbl[i] = new Double(pdbl[i], pdblImg[i]);
484                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
485             }
486         }
487         else
488         {
489             for (int i = 0; i < iSize; i++)
490             {
491                 pDbl[i] = new Double(pdbl[i]);
492                 pDblSp[i] = new Double(_pSp->getReal(iPositVal[i]), _pSp->getImg(iPositVal[i]).imag());
493             }
494         }
495     }
496     else
497     {
498         delete[] pDbl;
499         delete[] pDblSp;
500         delete[] iPositVal;
501         return 0;
502     }
503
504     Sparse* pSpTemp = new Sparse(_pSp->getRows(), _pSp->getCols(), _pSp->isComplex() || _pDouble->isComplex());
505     pSpTemp->zero_set();
506
507     Double* ppDblGet = NULL;
508     int iResultat;
509     for (int i = 0; i < iSize; i++)
510     {
511         if ((pDblSp[i]->get(0) != 0) || (pDblSp[i]->getImg(0) != 0))
512         {
513             iResultat = RDivideDoubleByDouble(pDblSp[i], pDbl[i], &ppDblGet);
514             if (iResultat != 0)
515             {
516                 for (int i = 0; i < iSize; ++i)
517                 {
518                     delete pDbl[i];
519                     delete pDblSp[i];
520                 }
521                 delete[] pDbl;
522                 delete[] pDblSp;
523                 delete[] iPositVal;
524                 delete pSpTemp;
525                 delete ppDblGet;
526                 return iResultat;
527             }
528             std::complex<double> cplx(ppDblGet->get(0), ppDblGet->getImg(0));
529             pSpTemp->set(iPositVal[i], cplx, false);
530             delete ppDblGet;
531         }
532     }
533
534     pSpTemp->finalize();
535
536     delete[] iPositVal;
537
538     for (int i = 0; i < iSize; ++i)
539     {
540         delete pDbl[i];
541         delete pDblSp[i];
542     }
543
544     delete[] pDbl;
545     delete[] pDblSp;
546
547     *_pSpOut = pSpTemp;
548     return 0;
549 }