License Header change: Removed the LICENSE_END before beta
[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 bool SinglePoly::setRank(int _iRank, bool bSave)
138 {
139     double *pR = NULL;
140     double *pI = NULL;
141     if (bSave == false)
142     {
143         if (getRank() != _iRank)
144         {
145             int piDims[2] = {_iRank + 1, 1};
146             if (m_pImgData == NULL)
147             {
148                 deleteAll();
149                 create(piDims, 2, &pR, NULL);
150             }
151             else
152             {
153                 deleteAll();
154                 create(piDims, 2, &pR, &pI);
155             }
156
157             return true;
158         }
159
160         return true;
161     }
162     else
163     {
164         double* pdblOldReal = m_pRealData;
165         double* pdblOldImg  = m_pImgData;
166         int iMinSize = Min(m_iSize, _iRank + 1);
167         int piDims[2] = {_iRank + 1, 1};
168
169         if (m_pImgData == NULL)
170         {
171             create(piDims, 2, &pR, NULL);
172         }
173         else
174         {
175             create(piDims, 2, &pR, &pI);
176             memcpy(m_pImgData, pdblOldImg, iMinSize * sizeof(double));
177         }
178
179         memcpy(m_pRealData, pdblOldReal, iMinSize * sizeof(double));
180
181         if (pdblOldImg)
182         {
183             delete[] pdblOldImg;
184             pdblOldImg = NULL;
185         }
186
187         delete[] pdblOldReal;
188         pdblOldReal = NULL;
189
190         return true;
191     }
192
193     return false;
194 }
195
196 bool SinglePoly::setZeros()
197 {
198     if (m_pRealData != NULL)
199     {
200         memset(m_pRealData, 0x00, m_iSize * sizeof(double));
201     }
202     else
203     {
204         return false;
205     }
206
207     if (isComplex() == true)
208     {
209         if (m_pImgData != NULL)
210         {
211             memset(m_pImgData, 0x00, m_iSize * sizeof(double));
212         }
213         else
214         {
215             return false;
216         }
217     }
218
219     return true;
220 }
221
222 bool SinglePoly::setCoef(Double* _pdblCoefR)
223 {
224     if (m_pRealData == NULL || _pdblCoefR == NULL)
225     {
226         return false;
227     }
228
229     double *pInR    = _pdblCoefR->getReal();
230     double *pInI    = _pdblCoefR->getImg();
231
232     return setCoef(pInR, pInI);
233 }
234
235 bool SinglePoly::setCoef(const double* _pdblCoefR, const double* _pdblCoefI)
236 {
237     if (_pdblCoefI != NULL && isComplex() == false)
238     {
239         setComplex(true);
240     }
241
242     if (_pdblCoefR != NULL)
243     {
244         memcpy(m_pRealData, _pdblCoefR, m_iSize * sizeof(double));
245     }
246
247     if (_pdblCoefI != NULL)
248     {
249         memcpy(m_pImgData, _pdblCoefI, m_iSize * sizeof(double));
250     }
251
252     return true;
253 }
254
255 void SinglePoly::whoAmI()
256 {
257     std::cout << "types::SinglePoly";
258 }
259
260 bool SinglePoly::evaluate(double _dblInR, double _dblInI, double *_pdblOutR, double *_pdblOutI)
261 {
262     *_pdblOutR = 0;
263     *_pdblOutI = 0;
264     if (m_iSize == 0)
265     {
266         return true;
267     }
268
269     for (int i = 0 ; i < m_iSize ; i++)
270     {
271         //real part
272         *_pdblOutR += m_pRealData[i] * std::pow(_dblInR, i);
273         //only if variable is complex
274         if (isComplex())
275         {
276             *_pdblOutR -= m_pImgData[i] * std::pow(_dblInI, i);
277             //img part
278             *_pdblOutI += m_pRealData[i] * std::pow(_dblInR, i);
279         }
280         *_pdblOutI += m_pRealData[i] * std::pow(_dblInI, i);
281     }
282
283     return true;
284 }
285
286 void SinglePoly::updateRank(void)
287 {
288     int iNewRank = getRank();
289     if (m_pImgData)
290     {
291         for (int i = getRank(); i > 0 ; i--)
292         {
293             if (std::fabs(m_pRealData[i]) == 0.0 && std::fabs(m_pImgData[i]) == 0.0)
294             {
295                 iNewRank--;
296             }
297             else
298             {
299                 break;
300             }
301         }
302     }
303     else
304     {
305         for (int i = getRank(); i > 0 ; i--)
306         {
307             if (std::fabs(m_pRealData[i]) == 0.0)
308             {
309                 iNewRank--;
310             }
311             else
312             {
313                 break;
314             }
315         }
316     }
317
318     if (iNewRank < getRank())
319     {
320         setRank(iNewRank, true);
321     }
322 }
323
324 bool SinglePoly::toString(std::wostringstream& ostr)
325 {
326     ostr << L"FIXME : implement SinglePoly::toString" << std::endl;
327     return true;
328 }
329
330 void SinglePoly::toStringReal(const std::wstring& _szVar, std::list<std::wstring>* _pListExp, std::list<std::wstring>* _pListCoef)
331 {
332     toStringInternal(m_pRealData, _szVar, _pListExp, _pListCoef);
333 }
334
335 void SinglePoly::toStringImg(const std::wstring& _szVar, std::list<std::wstring>* _pListExp, std::list<std::wstring>* _pListCoef)
336 {
337     if (isComplex() == false)
338     {
339         _pListExp->clear();
340         _pListCoef->clear();
341         return;
342     }
343
344     toStringInternal(m_pImgData, _szVar, _pListExp, _pListCoef);
345 }
346
347 bool SinglePoly::subMatrixToString(std::wostringstream& /*ostr*/, int* /*_piDims*/, int /*_iDims*/)
348 {
349     return false;
350 }
351
352 void SinglePoly::toStringInternal(double *_pdblVal, const std::wstring& _szVar, std::list<std::wstring>* _pListExp, std::list<std::wstring>* _pListCoef)
353 {
354     int iLineLen = ConfigVariable::getConsoleWidth();
355
356     std::wostringstream ostemp;
357     std::wostringstream ostemp2;
358
359     ostemp << L" ";
360     ostemp2 << L" ";
361
362     int iLen = 0;
363     int iLastFlush = 2;
364     for (int i = 0 ; i < m_iSize ; i++)
365     {
366         if (_pdblVal[i] != 0)
367         {
368             DoubleFormat df;
369             getDoubleFormat(_pdblVal[i], &df);
370
371             if (iLen + df.iWidth + df.iSignLen >= iLineLen - 1)
372             {
373                 iLastFlush = i;
374                 _pListExp->push_back(ostemp2.str());
375                 ostemp2.str(L""); //reset stream
376                 addSpaces(&ostemp2, 11); //take from scilab ... why not ...
377
378                 _pListCoef->push_back(ostemp.str());
379                 ostemp.str(L""); //reset stream
380                 addSpaces(&ostemp, 11); //take from scilab ... why not ...
381             }
382
383             bool bFirst = ostemp.str().size() == 1;
384
385             // In scientific notation case bExp == true, so we have to print point (2.000D+10s)
386             // In other case don't print point (2s)
387             df.bPrintPoint = df.bExp;
388             df.bPrintPlusSign = bFirst == false;
389             df.bPrintOne = i == 0;
390             addDoubleValue(&ostemp, _pdblVal[i], &df);
391
392             if (i == 0)
393             {
394                 iLen = static_cast<int>(ostemp.str().size());
395             }
396             else if (i == 1)
397             {
398                 // add polynom name
399                 ostemp << _szVar;
400                 iLen = static_cast<int>(ostemp.str().size());
401             }
402             else
403             {
404                 // add polynom name and exponent
405                 ostemp << _szVar;
406                 iLen = static_cast<int>(ostemp.str().size());
407                 addSpaces(&ostemp2, iLen - static_cast<int>(ostemp2.str().size()));
408                 ostemp2 << i;
409                 int iSize = static_cast<int>(ostemp2.str().size()) - iLen;
410                 addSpaces(&ostemp, iSize);
411             }
412         }
413     }
414
415     if (iLastFlush != 0)
416     {
417         if (ostemp.str() == L" ")
418         {
419             ostemp << L"  0";
420             addSpaces(&ostemp2, static_cast<int>(ostemp.str().size()));
421         }
422
423         if (ostemp2.str() == L" ")
424         {
425             // -1 because ostemp2 have already a space
426             addSpaces(&ostemp2, static_cast<int>(ostemp.str().size()) - 1);
427         }
428
429         _pListExp->push_back(ostemp2.str());
430         _pListCoef->push_back(ostemp.str());
431     }
432
433     return;
434 }
435
436 bool SinglePoly::operator==(const InternalType& it)
437 {
438     if (const_cast<InternalType &>(it).isSinglePoly() == false)
439     {
440         return false;
441     }
442
443     SinglePoly* pP = const_cast<InternalType &>(it).getAs<SinglePoly>();
444
445     if (getRank() != pP->getRank())
446     {
447         return false;
448     }
449
450     double *pdblReal = pP->get();
451     for (int i = 0 ; i < getSize() ; i++)
452     {
453         if (m_pRealData[i] != pdblReal[i])
454         {
455             return false;
456         }
457     }
458
459     //both complex
460     if (isComplex() && pP->isComplex())
461     {
462         double *pdblImg = pP->getImg();
463         for (int i = 0 ; i < m_iSize ; i++)
464         {
465             if (m_pImgData[i] != pdblImg[i])
466             {
467                 return false;
468             }
469         }
470     }
471     //pdbl complex check all img values == 0
472     else if (pP->isComplex())
473     {
474         double *pdblImg = pP->getImg();
475         for (int i = 0 ; i < m_iSize ; i++)
476         {
477             if (pdblImg[i])
478             {
479                 return false;
480             }
481         }
482     }
483     //complex check all img values == 0
484     else if (isComplex())
485     {
486         for (int i = 0 ; i < m_iSize ; i++)
487         {
488             if (m_pImgData[i])
489             {
490                 return false;
491             }
492         }
493     }
494
495     return true;
496 }
497
498 bool SinglePoly::operator!=(const InternalType& it)
499 {
500     return !(*this == it);
501 }
502
503 SinglePoly* SinglePoly::clone()
504 {
505     SinglePoly* pPoly = NULL;
506     double *pR = NULL;
507     if (isComplex())
508     {
509         double *pI = NULL;
510         pPoly = new SinglePoly(&pR, &pI, getRank());
511         pPoly->setCoef(m_pRealData, m_pImgData);
512     }
513     else
514     {
515         pPoly = new SinglePoly(&pR, getRank());
516         pPoly->setCoef(m_pRealData, NULL);
517     }
518     return pPoly;
519 }
520
521 SinglePoly* SinglePoly::conjugate()
522 {
523     if (isComplex())
524     {
525         double *pR = NULL;
526         double *pI = NULL;
527         SinglePoly * pPoly = new SinglePoly(&pR, &pI, getRank());
528
529         Transposition::conjugate(m_iSize, m_pRealData, pR, m_pImgData, pI);
530
531         return pPoly;
532     }
533     else
534     {
535         return clone();
536     }
537 }
538
539 SinglePoly* operator*(const SinglePoly& _lhs, const SinglePoly& _rhs)
540 {
541     SinglePoly& lhs = const_cast<SinglePoly &>(_lhs);
542     SinglePoly& rhs = const_cast<SinglePoly &>(_rhs);
543     SinglePoly* pOut = NULL;
544
545     bool isComplexL = lhs.isComplex();
546     bool isComplexR = rhs.isComplex();
547     bool isComplexOut = isComplexL || isComplexR;
548
549     int iRankL = lhs.getRank();
550     int iRankR = rhs.getRank();
551     int iRankOut = lhs.getRank() + rhs.getRank();
552
553     double* pdblOutR = NULL;
554     double* pdblOutI = NULL;
555     double* pdblLR = lhs.get();
556     double* pdblLI = lhs.getImg();
557     double* pdblRR = rhs.get();
558     double* pdblRI = rhs.getImg();
559
560     if (isComplexOut)
561     {
562         pOut = new SinglePoly(&pdblOutR, &pdblOutI, iRankOut);
563         memset(pdblOutR, 0x00, sizeof(double) * (iRankOut + 1));
564         memset(pdblOutI, 0x00, sizeof(double) * (iRankOut + 1));
565     }
566     else
567     {
568         pOut = new SinglePoly(&pdblOutR, iRankOut);
569         memset(pdblOutR, 0x00, sizeof(double) * (iRankOut + 1));
570     }
571
572     if (isComplexL)
573     {
574         if (isComplexR)
575         {
576             for (int i = 0 ; i < iRankL + 1; ++i)
577             {
578                 for (int j = 0 ; j < iRankR + 1; ++j)
579                 {
580                     pdblOutR[i + j]  += pdblLR[i] * pdblRR[j] - pdblLI[i] * pdblRI[j];
581                     pdblOutI[i + j]  += pdblLI[i] * pdblRR[j] + pdblLR[i] * pdblRI[j];
582                 }
583             }
584         }
585         else
586         {
587             for (int i = 0 ; i < iRankL + 1 ; ++i)
588             {
589                 for (int j = 0 ; j < iRankR + 1 ; ++j)
590                 {
591                     pdblOutR[i + j]  += pdblLR[i] * pdblRR[j];
592                     pdblOutI[i + j]  += pdblLI[i] * pdblRR[j];
593                 }
594             }
595         }
596     }
597     else
598     {
599         if (isComplexR)
600         {
601             for (int i = 0 ; i < iRankL + 1 ; ++i)
602             {
603                 for (int j = 0 ; j < iRankR + 1 ; ++j)
604                 {
605                     pdblOutR[i + j]  += pdblLR[i] * pdblRR[j];
606                     pdblOutI[i + j]  += pdblLR[i] * pdblRI[j];
607                 }
608             }
609         }
610         else
611         {
612             for (int i = 0 ; i < iRankL + 1 ; ++i)
613             {
614                 for (int j = 0 ; j < iRankR + 1 ; ++j)
615                 {
616                     pdblOutR[i + j]  += pdblLR[i] * pdblRR[j];
617                 }
618             }
619         }
620     }
621
622     return pOut;
623 }
624 }
625
626
627