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