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