* Bug 15599 fixed: now degree of zero polynomial is -Inf
[scilab.git] / scilab / modules / ast / src / cpp / types / singlepoly.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 #include <sstream>
16 #include <cmath>
17 #include "singlepoly.hxx"
18 #include "double.hxx"
19 #include "tostring_common.hxx"
20 #include "configvariable.hxx"
21 #include "scilabWrite.hxx"
22
23 extern "C"
24 {
25 #include "log.h"
26 #include "exp.h"
27 #include "elem_common.h"
28 }
29
30 namespace types
31 {
32 SinglePoly::SinglePoly()
33 {
34     double* pdblCoefR = NULL;
35     int piDims[2] = {1, 1};
36     create(piDims, 2, &pdblCoefR, NULL);
37     pdblCoefR[0] = 0;
38 #ifndef NDEBUG
39     Inspector::addItem(this);
40 #endif
41 }
42
43 SinglePoly::SinglePoly(double** _pdblCoefR, int _iRank)
44 {
45     int piDims[2] = {_iRank + 1, 1};
46     create(piDims, 2, _pdblCoefR, NULL);
47 #ifndef NDEBUG
48     Inspector::addItem(this);
49 #endif
50 }
51
52 SinglePoly::SinglePoly(double** _pdblCoefR, double** _pdblCoefI, int _iRank)
53 {
54     int piDims[2] = {_iRank + 1, 1};
55     create(piDims, 2, _pdblCoefR, _pdblCoefI);
56 #ifndef NDEBUG
57     Inspector::addItem(this);
58 #endif
59 }
60
61 SinglePoly::~SinglePoly()
62 {
63     deleteAll();
64 #ifndef NDEBUG
65     Inspector::removeItem(this);
66 #endif
67 }
68
69 void SinglePoly::deleteAll()
70 {
71     delete[] m_pRealData;
72     m_pRealData = NULL;
73     deleteImg();
74 }
75
76 void SinglePoly::deleteImg()
77 {
78     if (m_pImgData != NULL)
79     {
80         delete[] m_pImgData;
81         m_pImgData = NULL;
82     }
83 }
84
85 double SinglePoly::getNullValue()
86 {
87     return 0;
88 }
89
90 SinglePoly* SinglePoly::createEmpty(int /*_iDims*/, int* _piDims, bool _bComplex)
91 {
92     double* pdblData = NULL;
93     SinglePoly* pSP = new SinglePoly(&pdblData, _piDims[0] - 1);
94     pSP->setComplex(_bComplex);
95
96     return pSP;
97 }
98
99 double* SinglePoly::allocData(int _iSize)
100 {
101     double* pDbl = NULL;
102     try
103     {
104         if (_iSize < 0)
105         {
106             m_pRealData = NULL;
107             m_pImgData = NULL;
108             char message[bsiz];
109             os_sprintf(message, _("Can not allocate negative size (%d).\n"),  _iSize);
110             throw ast::InternalError(message);
111         }
112         else
113         {
114             pDbl = new double[_iSize];
115         }
116     }
117     catch (std::bad_alloc &/*e*/)
118     {
119         char message[bsiz];
120         os_sprintf(message, _("Can not allocate %.2f MB memory.\n"),  (double) (_iSize * sizeof(double)) / 1.e6);
121         throw ast::InternalError(message);
122     }
123
124     return pDbl;
125 }
126
127 double SinglePoly::copyValue(double _dblData)
128 {
129     return _dblData;
130 }
131
132 int SinglePoly::getRank()
133 {
134     return m_iSize - 1;
135 }
136
137 double SinglePoly::getDegree()
138 {
139     return m_iSize==1 && m_pRealData[0]==0 && (m_pImgData == NULL || m_pImgData[0]==0) ? -INFINITY : m_iSize - 1;
140 }
141
142 bool SinglePoly::setRank(int _iRank, bool bSave)
143 {
144     double *pR = NULL;
145     double *pI = NULL;
146     if (bSave == false)
147     {
148         if (getRank() != _iRank)
149         {
150             int piDims[2] = {_iRank + 1, 1};
151             if (m_pImgData == NULL)
152             {
153                 deleteAll();
154                 create(piDims, 2, &pR, NULL);
155             }
156             else
157             {
158                 deleteAll();
159                 create(piDims, 2, &pR, &pI);
160             }
161
162             return true;
163         }
164
165         return true;
166     }
167     else
168     {
169         double* pdblOldReal = m_pRealData;
170         double* pdblOldImg  = m_pImgData;
171         int iMinSize = Min(m_iSize, _iRank + 1);
172         int piDims[2] = {_iRank + 1, 1};
173
174         if (m_pImgData == NULL)
175         {
176             create(piDims, 2, &pR, NULL);
177         }
178         else
179         {
180             create(piDims, 2, &pR, &pI);
181             memcpy(m_pImgData, pdblOldImg, iMinSize * sizeof(double));
182         }
183
184         memcpy(m_pRealData, pdblOldReal, iMinSize * sizeof(double));
185
186         if (pdblOldImg)
187         {
188             delete[] pdblOldImg;
189             pdblOldImg = NULL;
190         }
191
192         delete[] pdblOldReal;
193         pdblOldReal = NULL;
194
195         return true;
196     }
197
198     return false;
199 }
200
201 bool SinglePoly::setZeros()
202 {
203     if (m_pRealData != NULL)
204     {
205         memset(m_pRealData, 0x00, m_iSize * sizeof(double));
206     }
207     else
208     {
209         return false;
210     }
211
212     if (isComplex() == true)
213     {
214         if (m_pImgData != NULL)
215         {
216             memset(m_pImgData, 0x00, m_iSize * sizeof(double));
217         }
218         else
219         {
220             return false;
221         }
222     }
223
224     return true;
225 }
226
227 bool SinglePoly::setCoef(Double* _pdblCoefR)
228 {
229     if (m_pRealData == NULL || _pdblCoefR == NULL)
230     {
231         return false;
232     }
233
234     double *pInR    = _pdblCoefR->getReal();
235     double *pInI    = _pdblCoefR->getImg();
236
237     return setCoef(pInR, pInI);
238 }
239
240 bool SinglePoly::setCoef(const double* _pdblCoefR, const double* _pdblCoefI)
241 {
242     if (_pdblCoefI != NULL && isComplex() == false)
243     {
244         setComplex(true);
245     }
246
247     if (_pdblCoefR != NULL)
248     {
249         memcpy(m_pRealData, _pdblCoefR, m_iSize * sizeof(double));
250     }
251
252     if (_pdblCoefI != NULL)
253     {
254         memcpy(m_pImgData, _pdblCoefI, m_iSize * sizeof(double));
255     }
256
257     return true;
258 }
259
260 void SinglePoly::whoAmI()
261 {
262     std::cout << "types::SinglePoly";
263 }
264
265 bool SinglePoly::evaluate(double _dblInR, double _dblInI, double *_pdblOutR, double *_pdblOutI)
266 {
267     *_pdblOutR = 0;
268     *_pdblOutI = 0;
269     if (m_iSize == 0)
270     {
271         return true;
272     }
273
274     for (int i = 0 ; i < m_iSize ; i++)
275     {
276         //real part
277         *_pdblOutR += m_pRealData[i] * std::pow(_dblInR, i);
278         //only if variable is complex
279         if (isComplex())
280         {
281             *_pdblOutR -= m_pImgData[i] * std::pow(_dblInI, i);
282             //img part
283             *_pdblOutI += m_pRealData[i] * std::pow(_dblInR, i);
284         }
285         *_pdblOutI += m_pRealData[i] * std::pow(_dblInI, i);
286     }
287
288     return true;
289 }
290
291 void SinglePoly::updateRank(void)
292 {
293     int iNewRank = getRank();
294     if (m_pImgData)
295     {
296         for (int i = getRank(); i > 0 ; i--)
297         {
298             if (std::fabs(m_pRealData[i]) == 0.0 && std::fabs(m_pImgData[i]) == 0.0)
299             {
300                 iNewRank--;
301             }
302             else
303             {
304                 break;
305             }
306         }
307     }
308     else
309     {
310         for (int i = getRank(); i > 0 ; i--)
311         {
312             if (std::fabs(m_pRealData[i]) == 0.0)
313             {
314                 iNewRank--;
315             }
316             else
317             {
318                 break;
319             }
320         }
321     }
322
323     if (iNewRank < getRank())
324     {
325         setRank(iNewRank, true);
326     }
327 }
328
329 bool SinglePoly::toString(std::wostringstream& ostr)
330 {
331     ostr << L"FIXME : implement SinglePoly::toString" << std::endl;
332     return true;
333 }
334
335 void SinglePoly::toStringReal(const std::wstring& _szVar, std::list<std::wstring>* _pListExp, std::list<std::wstring>* _pListCoef)
336 {
337     toStringInternal(m_pRealData, _szVar, _pListExp, _pListCoef);
338 }
339
340 void SinglePoly::toStringImg(const std::wstring& _szVar, std::list<std::wstring>* _pListExp, std::list<std::wstring>* _pListCoef)
341 {
342     if (isComplex() == false)
343     {
344         _pListExp->clear();
345         _pListCoef->clear();
346         return;
347     }
348
349     toStringInternal(m_pImgData, _szVar, _pListExp, _pListCoef);
350 }
351
352 bool SinglePoly::subMatrixToString(std::wostringstream& /*ostr*/, int* /*_piDims*/, int /*_iDims*/)
353 {
354     return false;
355 }
356
357 void SinglePoly::toStringInternal(double *_pdblVal, const std::wstring& _szVar, std::list<std::wstring>* _pListExp, std::list<std::wstring>* _pListCoef)
358 {
359     int iLineLen = ConfigVariable::getConsoleWidth();
360
361     std::wostringstream ostemp;
362     std::wostringstream ostemp2;
363
364     ostemp << L" ";
365     ostemp2 << L" ";
366
367     int iLen = 0;
368     int iLastFlush = 2;
369     for (int i = 0 ; i < m_iSize ; i++)
370     {
371         if (_pdblVal[i] != 0)
372         {
373             DoubleFormat df;
374             getDoubleFormat(_pdblVal[i], &df);
375
376             if (iLen + df.iWidth + df.iSignLen >= iLineLen - 1)
377             {
378                 iLastFlush = i;
379                 _pListExp->push_back(ostemp2.str());
380                 ostemp2.str(L""); //reset stream
381                 addSpaces(&ostemp2, 11); //take from scilab ... why not ...
382
383                 _pListCoef->push_back(ostemp.str());
384                 ostemp.str(L""); //reset stream
385                 addSpaces(&ostemp, 11); //take from scilab ... why not ...
386             }
387
388             bool bFirst = ostemp.str().size() == 1;
389
390             // In scientific notation case bExp == true, so we have to print point (2.000D+10s)
391             // In other case don't print point (2s)
392             df.bPrintPoint = df.bExp;
393             df.bPrintPlusSign = bFirst == false;
394             df.bPrintOne = i == 0;
395             addDoubleValue(&ostemp, _pdblVal[i], &df);
396
397             if (i == 0)
398             {
399                 iLen = static_cast<int>(ostemp.str().size());
400             }
401             else if (i == 1)
402             {
403                 // add polynom name
404                 ostemp << _szVar;
405                 iLen = static_cast<int>(ostemp.str().size());
406             }
407             else
408             {
409                 // add polynom name and exponent
410                 ostemp << _szVar;
411                 iLen = static_cast<int>(ostemp.str().size());
412                 addSpaces(&ostemp2, iLen - static_cast<int>(ostemp2.str().size()));
413                 ostemp2 << i;
414                 int iSize = static_cast<int>(ostemp2.str().size()) - iLen;
415                 addSpaces(&ostemp, iSize);
416             }
417         }
418     }
419
420     if (iLastFlush != 0)
421     {
422         if (ostemp.str() == L" ")
423         {
424             ostemp << L"  0";
425             addSpaces(&ostemp2, static_cast<int>(ostemp.str().size()));
426         }
427
428         if (ostemp2.str() == L" ")
429         {
430             // -1 because ostemp2 have already a space
431             addSpaces(&ostemp2, static_cast<int>(ostemp.str().size()) - 1);
432         }
433
434         _pListExp->push_back(ostemp2.str());
435         _pListCoef->push_back(ostemp.str());
436     }
437
438     return;
439 }
440
441 bool SinglePoly::operator==(const InternalType& it)
442 {
443     if (const_cast<InternalType &>(it).isSinglePoly() == false)
444     {
445         return false;
446     }
447
448     SinglePoly* pP = const_cast<InternalType &>(it).getAs<SinglePoly>();
449
450     if (getRank() != pP->getRank())
451     {
452         return false;
453     }
454
455     double *pdblReal = pP->get();
456     for (int i = 0 ; i < getSize() ; i++)
457     {
458         if (m_pRealData[i] != pdblReal[i])
459         {
460             return false;
461         }
462     }
463
464     //both complex
465     if (isComplex() && pP->isComplex())
466     {
467         double *pdblImg = pP->getImg();
468         for (int i = 0 ; i < m_iSize ; i++)
469         {
470             if (m_pImgData[i] != pdblImg[i])
471             {
472                 return false;
473             }
474         }
475     }
476     //pdbl complex check all img values == 0
477     else if (pP->isComplex())
478     {
479         double *pdblImg = pP->getImg();
480         for (int i = 0 ; i < m_iSize ; i++)
481         {
482             if (pdblImg[i])
483             {
484                 return false;
485             }
486         }
487     }
488     //complex check all img values == 0
489     else if (isComplex())
490     {
491         for (int i = 0 ; i < m_iSize ; i++)
492         {
493             if (m_pImgData[i])
494             {
495                 return false;
496             }
497         }
498     }
499
500     return true;
501 }
502
503 bool SinglePoly::operator!=(const InternalType& it)
504 {
505     return !(*this == it);
506 }
507
508 SinglePoly* SinglePoly::clone()
509 {
510     SinglePoly* pPoly = NULL;
511     double *pR = NULL;
512     if (isComplex())
513     {
514         double *pI = NULL;
515         pPoly = new SinglePoly(&pR, &pI, getRank());
516         pPoly->setCoef(m_pRealData, m_pImgData);
517     }
518     else
519     {
520         pPoly = new SinglePoly(&pR, getRank());
521         pPoly->setCoef(m_pRealData, NULL);
522     }
523     return pPoly;
524 }
525
526 SinglePoly* SinglePoly::conjugate()
527 {
528     if (isComplex())
529     {
530         double *pR = NULL;
531         double *pI = NULL;
532         SinglePoly * pPoly = new SinglePoly(&pR, &pI, getRank());
533
534         Transposition::conjugate(m_iSize, m_pRealData, pR, m_pImgData, pI);
535
536         return pPoly;
537     }
538     else
539     {
540         return clone();
541     }
542 }
543
544 SinglePoly* operator*(const SinglePoly& _lhs, const SinglePoly& _rhs)
545 {
546     SinglePoly& lhs = const_cast<SinglePoly &>(_lhs);
547     SinglePoly& rhs = const_cast<SinglePoly &>(_rhs);
548     SinglePoly* pOut = NULL;
549
550     bool isComplexL = lhs.isComplex();
551     bool isComplexR = rhs.isComplex();
552     bool isComplexOut = isComplexL || isComplexR;
553
554     int iRankL = lhs.getRank();
555     int iRankR = rhs.getRank();
556     int iRankOut = lhs.getRank() + rhs.getRank();
557
558     double* pdblOutR = NULL;
559     double* pdblOutI = NULL;
560     double* pdblLR = lhs.get();
561     double* pdblLI = lhs.getImg();
562     double* pdblRR = rhs.get();
563     double* pdblRI = rhs.getImg();
564
565     if (isComplexOut)
566     {
567         pOut = new SinglePoly(&pdblOutR, &pdblOutI, iRankOut);
568         memset(pdblOutR, 0x00, sizeof(double) * (iRankOut + 1));
569         memset(pdblOutI, 0x00, sizeof(double) * (iRankOut + 1));
570     }
571     else
572     {
573         pOut = new SinglePoly(&pdblOutR, iRankOut);
574         memset(pdblOutR, 0x00, sizeof(double) * (iRankOut + 1));
575     }
576
577     if (isComplexL)
578     {
579         if (isComplexR)
580         {
581             for (int i = 0 ; i < iRankL + 1; ++i)
582             {
583                 for (int j = 0 ; j < iRankR + 1; ++j)
584                 {
585                     pdblOutR[i + j]  += pdblLR[i] * pdblRR[j] - pdblLI[i] * pdblRI[j];
586                     pdblOutI[i + j]  += pdblLI[i] * pdblRR[j] + pdblLR[i] * pdblRI[j];
587                 }
588             }
589         }
590         else
591         {
592             for (int i = 0 ; i < iRankL + 1 ; ++i)
593             {
594                 for (int j = 0 ; j < iRankR + 1 ; ++j)
595                 {
596                     pdblOutR[i + j]  += pdblLR[i] * pdblRR[j];
597                     pdblOutI[i + j]  += pdblLI[i] * pdblRR[j];
598                 }
599             }
600         }
601     }
602     else
603     {
604         if (isComplexR)
605         {
606             for (int i = 0 ; i < iRankL + 1 ; ++i)
607             {
608                 for (int j = 0 ; j < iRankR + 1 ; ++j)
609                 {
610                     pdblOutR[i + j]  += pdblLR[i] * pdblRR[j];
611                     pdblOutI[i + j]  += pdblLR[i] * pdblRI[j];
612                 }
613             }
614         }
615         else
616         {
617             for (int i = 0 ; i < iRankL + 1 ; ++i)
618             {
619                 for (int j = 0 ; j < iRankR + 1 ; ++j)
620                 {
621                     pdblOutR[i + j]  += pdblLR[i] * pdblRR[j];
622                 }
623             }
624         }
625     }
626
627     return pOut;
628 }
629 }
630
631
632