refactoring polynom.
[scilab.git] / scilab / modules / ast / src / cpp / operations / types_substraction.cpp
1 /*
2 *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 *  Copyright (C) 2008-2008 - DIGITEO - Antoine ELIAS
4 *  Copyright (C) 2010-2010 - DIGITEO - Bruno JOFRET
5 *
6 *  This file must be used under the terms of the CeCILL.
7 *  This source file is licensed as described in the file COPYING, which
8 *  you should have received as part of this distribution.  The terms
9 *  are also available at
10 *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11 *
12 */
13
14 #include <stdio.h>
15
16 #include "types_substraction.hxx"
17 #include "double.hxx"
18 #include "polynom.hxx"
19 #include "int.hxx"
20 #include "sparse.hxx"
21 #include "generic_operations.hxx"
22
23 extern "C"
24 {
25 #include "matrix_substraction.h"
26 #include "elem_common.h" //dset
27 }
28
29 using namespace types;
30
31 InternalType* GenericUnaryMinus(InternalType* _pRightOperand)
32 {
33     // - Double
34     if (_pRightOperand->isDouble())
35     {
36         Double *pR = dynamic_cast<Double*>(_pRightOperand->clone());
37         bool bComplex = pR->isComplex();
38         double* pReal = pR->get();
39
40         if (bComplex)
41         {
42             double* pImg = pR->getImg();
43             for (int i = 0 ; i < pR->getSize() ; i++)
44             {
45                 pReal[i] *= -1;
46                 pImg[i]  *= -1;
47             }
48         }
49         else
50         {
51             for (int i = 0 ; i < pR->getSize() ; i++)
52             {
53                 pReal[i] *= -1;
54             }
55         }
56
57         return pR;
58     }
59
60     // - Polynom
61     if (_pRightOperand->isPoly())
62     {
63         double* pReal = NULL;
64         double* pImg  = NULL;
65
66         SinglePoly* pSPCoef = NULL;
67
68         Polynom *pR = dynamic_cast<Polynom*>(_pRightOperand->clone());
69         bool bComplex = pR->isComplex();
70
71         if (bComplex)
72         {
73             for (int i = 0; i < pR->getSize(); i++)
74             {
75                 pSPCoef = pR->get(i);
76                 pReal = pSPCoef->get();
77                 pImg = pSPCoef->getImg();
78
79                 for (int j = 0; j < pSPCoef->getSize(); j++)
80                 {
81                     pReal[j] *= -1;
82                     pImg[j] *= -1;
83                 }
84             }
85         }
86         else
87         {
88             for (int i = 0; i < pR->getSize(); i++)
89             {
90                 pSPCoef = pR->get(i);
91                 pReal = pSPCoef->get();
92
93                 for (int j = 0; j < pSPCoef->getSize(); j++)
94                 {
95                     pReal[j] *= -1;
96                 }
97             }
98         }
99
100         return pR;
101     }
102
103     // - Int
104     if (_pRightOperand->isInt())
105     {
106         switch (_pRightOperand->getType())
107         {
108             case InternalType::ScilabInt8 :
109             {
110                 Int8 *pR = _pRightOperand->clone()->getAs<Int8>();
111                 char* pReal = pR->get();
112
113                 for (int i = 0 ; i < pR->getSize() ; i++)
114                 {
115                     pReal[i] *= -1;
116                 }
117                 return pR;
118             }
119             case InternalType::ScilabInt16 :
120             {
121                 Int16 *pR = _pRightOperand->clone()->getAs<Int16>();
122                 short* pReal = pR->get();
123
124                 for (int i = 0 ; i < pR->getSize() ; i++)
125                 {
126                     pReal[i] *= -1;
127                 }
128                 return pR;
129             }
130             case InternalType::ScilabInt32 :
131             {
132                 Int32 *pR = _pRightOperand->clone()->getAs<Int32>();
133                 int* pReal = pR->get();
134
135                 for (int i = 0 ; i < pR->getSize() ; i++)
136                 {
137                     pReal[i] *= -1;
138                 }
139                 return pR;
140             }
141             case InternalType::ScilabInt64 :
142             {
143                 Int64 *pR = _pRightOperand->clone()->getAs<Int64>();
144                 long long* pReal = pR->get();
145
146                 for (int i = 0 ; i < pR->getSize() ; i++)
147                 {
148                     pReal[i] *= -1;
149                 }
150                 return pR;
151             }
152             case InternalType::ScilabUInt8 :
153             {
154                 UInt8 *pR = _pRightOperand->clone()->getAs<UInt8>();
155                 unsigned char* pReal = pR->get();
156
157                 for (int i = 0 ; i < pR->getSize() ; i++)
158                 {
159                     pReal[i] *= -1;
160                 }
161                 return pR;
162             }
163             case InternalType::ScilabUInt16 :
164             {
165                 UInt16 *pR = _pRightOperand->clone()->getAs<UInt16>();
166                 unsigned short* pReal = pR->get();
167
168                 for (int i = 0 ; i < pR->getSize() ; i++)
169                 {
170                     pReal[i] *= -1;
171                 }
172                 return pR;
173             }
174             case InternalType::ScilabUInt32 :
175             {
176                 UInt32 *pR = _pRightOperand->clone()->getAs<UInt32>();
177                 unsigned int* pReal = pR->get();
178
179                 for (int i = 0 ; i < pR->getSize() ; i++)
180                 {
181                     pReal[i] *= -1;
182                 }
183                 return pR;
184             }
185             case InternalType::ScilabUInt64 :
186             {
187                 UInt64 *pR = _pRightOperand->clone()->getAs<UInt64>();
188                 unsigned long long* pReal = pR->get();
189
190                 for (int i = 0 ; i < pR->getSize() ; i++)
191                 {
192                     pReal[i] *= -1;
193                 }
194                 return pR;
195             }
196             default :
197                 return NULL;
198         }
199     }
200
201     // - Sparse
202     if (_pRightOperand->isSparse())
203     {
204         Sparse* pSp = dynamic_cast<Sparse*>(_pRightOperand)->multiply(-1);
205         return pSp;
206     }
207
208     /*
209     ** Default case : Return NULL will Call Overloading.
210     */
211     return NULL;
212 }
213
214 InternalType* GenericMinus(InternalType* _pLeftOperand, InternalType* _pRightOperand)
215 {
216     InternalType *pResult = NULL;
217
218     // InternalType - [] and [] - []
219     if (_pRightOperand->isDouble() && _pRightOperand->getAs<Double>()->isEmpty())
220     {
221         return _pLeftOperand->clone();
222     }
223
224     // [] - InternalType
225     if (_pLeftOperand->isDouble() && _pLeftOperand->getAs<Double>()->isEmpty())
226     {
227         if (_pRightOperand->isSparse())
228         {
229             InternalType* pSpOut = NULL;
230             Sparse* pSp = new Sparse(1, 1);
231             pSp->set(0, 0, -1.0);
232             pSpOut = GenericTimes((InternalType*)pSp, _pRightOperand);
233             delete pSp;
234             return pSpOut;
235         }
236
237         Double pDblZero = Double(0);
238         return GenericMinus((InternalType*)&pDblZero, _pRightOperand);
239     }
240
241     /*
242     ** DOUBLE - DOUBLE
243     */
244     if (_pLeftOperand->isDouble() && _pRightOperand->isDouble())
245     {
246         Double *pL = _pLeftOperand->getAs<Double>();
247         Double *pR = _pRightOperand->getAs<Double>();
248
249         int iResult = SubstractDoubleToDouble(pL, pR, (Double**)&pResult);
250         if (iResult != 0)
251         {
252             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
253         }
254         return pResult;
255     }
256
257     /*
258     ** DOUBLE - POLY
259     */
260     else if (_pLeftOperand->isDouble() && _pRightOperand->isPoly())
261     {
262         Double *pL              = _pLeftOperand->getAs<Double>();
263         Polynom *pR          = _pRightOperand->getAs<types::Polynom>();
264
265         int iResult = SubstractPolyToDouble(pL, pR, (Polynom**)&pResult);
266         if (iResult != 0)
267         {
268             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
269         }
270         return pResult;
271     }
272
273     /*
274     ** POLY - DOUBLE
275     */
276     else if (_pLeftOperand->isPoly() && _pRightOperand->isDouble())
277     {
278         Polynom *pL   = _pLeftOperand->getAs<types::Polynom>();
279         Double *pR    = _pRightOperand->getAs<Double>();
280
281         int iResult = SubstractDoubleToPoly(pL, pR, (Polynom**)&pResult);
282         if (iResult != 0)
283         {
284             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
285         }
286         return pResult;
287     }
288
289     /*
290     ** POLY - POLY
291     */
292     else if (_pLeftOperand->isPoly() && _pRightOperand->isPoly())
293     {
294         Polynom *pL   = _pLeftOperand->getAs<types::Polynom>();
295         Polynom *pR   = _pRightOperand->getAs<types::Polynom>();
296
297         int iResult = SubstractPolyToPoly(pL, pR, (Polynom**)&pResult);
298         if (iResult != 0)
299         {
300             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
301         }
302         return pResult;
303     }
304
305     /*
306     ** SPARSE - SPARSE
307     */
308     else if (_pLeftOperand->isSparse() && _pRightOperand->isSparse())
309     {
310         Sparse *pL = _pLeftOperand->getAs<Sparse>();
311         Sparse *pR = _pRightOperand->getAs<Sparse>();
312
313         int iResult = SubstractSparseToSparse(pR, pL, (GenericType**)&pResult);
314         if (iResult != 0)
315         {
316             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
317         }
318
319         return pResult;
320     }
321
322     /*
323     ** SPARSE - DOUBLE
324     */
325     else if (_pLeftOperand->isSparse() && _pRightOperand->isDouble())
326     {
327         Sparse *pL = _pLeftOperand->getAs<Sparse>();
328         Double *pR = _pRightOperand->getAs<Double>();
329
330         int iResult = SubstractDoubleToSparse(pR, pL, (GenericType**)&pResult);
331         if (iResult != 0)
332         {
333             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
334         }
335
336         return pResult;
337     }
338
339     /*
340     ** DOUBLE - SPARSE
341     */
342     else if (_pLeftOperand->isDouble() && _pRightOperand->isSparse())
343     {
344         Double *pL = _pLeftOperand->getAs<Double>();
345         Sparse *pR = _pRightOperand->getAs<Sparse>();
346
347         int iResult = SubstractSparseToDouble(pR, pL, (GenericType**)&pResult);
348         if (iResult != 0)
349         {
350             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
351         }
352
353         return pResult;
354     }
355     /*
356     ** Default case : Return NULL will Call Overloading.
357     */
358     return NULL;
359 }
360
361 int SubstractDoubleToDouble(Double* _pDouble1, Double* _pDouble2, Double** _pDoubleOut)
362 {
363     bool bComplex1  = _pDouble1->isComplex();
364     bool bComplex2  = _pDouble2->isComplex();
365
366     if (_pDouble1->isIdentity())
367     {
368         if (_pDouble2->getDims() > 2)
369         {
370             //unable to substract identity matrix and greater than 2 dimensions matrix
371             return 1;
372         }
373
374         (*_pDoubleOut) = new Double(_pDouble2->getRows(), _pDouble2->getCols(), bComplex1 || bComplex2);
375
376         if (bComplex1 == false && bComplex2 == false)
377         {
378             iSubstractRealIdentityToRealMatrix(
379                 _pDouble1->get(0),
380                 _pDouble2->get(), _pDouble2->getRows(), _pDouble2->getCols(),
381                 (*_pDoubleOut)->get());
382         }
383         else if (bComplex1 == false && bComplex2 == true)
384         {
385             iSubstractRealIdentityToComplexMatrix(
386                 _pDouble1->get(0),
387                 _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
388                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
389         }
390         else if (bComplex1 == true && bComplex2 == false)
391         {
392             iSubstractComplexIdentityToRealMatrix(
393                 _pDouble1->get(0), _pDouble1->getImg(0),
394                 _pDouble2->get(), _pDouble2->getRows(), _pDouble2->getCols(),
395                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
396         }
397         else if (bComplex1 == true && bComplex2 == true)
398         {
399             iSubstractComplexIdentityToComplexMatrix(
400                 _pDouble1->get(0), _pDouble1->getImg(0),
401                 _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
402                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
403         }
404
405         //return opposite
406         double* pDbl = (*_pDoubleOut)->get();
407         if ((*_pDoubleOut)->isComplex())
408         {
409             double* pDblImg = (*_pDoubleOut)->getImg();
410             for (int i = 0 ; i < _pDouble2->getSize() ; i++)
411             {
412                 pDbl[i] = -pDbl[i];
413                 pDblImg[i] = -pDblImg[i];
414             }
415         }
416         else
417         {
418             for (int i = 0 ; i < _pDouble2->getSize() ; i++)
419             {
420                 pDbl[i] = -pDbl[i];
421             }
422         }
423         return 0;
424     }
425
426     if (_pDouble2->isIdentity())
427     {
428         if (_pDouble1->getDims() > 2)
429         {
430             //unable to substract identity matrix and greater than 2 dimensions matrix
431             return 1;
432         }
433
434         (*_pDoubleOut) = new Double(_pDouble1->getRows(), _pDouble1->getCols(), bComplex1 || bComplex2);
435         if (bComplex1 == false && bComplex2 == false)
436         {
437             iSubstractRealIdentityToRealMatrix(
438                 _pDouble2->get(0),
439                 _pDouble1->get(), _pDouble1->getRows(), _pDouble1->getCols(),
440                 (*_pDoubleOut)->get());
441         }
442         else if (bComplex1 == false && bComplex2 == true)
443         {
444             iSubstractComplexIdentityToRealMatrix(
445                 _pDouble2->get(0), _pDouble2->getImg(0),
446                 _pDouble1->get(), _pDouble1->getRows(), _pDouble1->getCols(),
447                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
448         }
449         else if (bComplex1 == true && bComplex2 == false)
450         {
451             iSubstractRealIdentityToComplexMatrix(
452                 _pDouble2->get(0),
453                 _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
454                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
455         }
456         else if (bComplex1 == true && bComplex2 == true)
457         {
458             iSubstractComplexIdentityToComplexMatrix(
459                 _pDouble2->get(0), _pDouble2->getImg(0),
460                 _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
461                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
462         }
463
464         return 0;
465     }
466
467     if (_pDouble1->isScalar())
468     {
469         //add pL with each element of pR
470         (*_pDoubleOut) = new Double(_pDouble2->getDims(), _pDouble2->getDimsArray(), bComplex1 || bComplex2);
471         if (bComplex1 == false && bComplex2 == false)
472         {
473             iSubstractRealMatrixToRealScalar(
474                 _pDouble2->get(), _pDouble2->getDimsArray(), _pDouble2->getDims(),
475                 _pDouble1->get(0),
476                 (*_pDoubleOut)->get());
477         }
478         else if (bComplex1 == false && bComplex2 == true)
479         {
480             iSubstractComplexMatrixToRealScalar(
481                 _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getDimsArray(), _pDouble2->getDims(),
482                 _pDouble1->get(0),
483                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
484         }
485         else if (bComplex1 == true && bComplex2 == false)
486         {
487             iSubstractRealMatrixToComplexScalar(
488                 _pDouble2->get(), _pDouble2->getDimsArray(), _pDouble2->getDims(),
489                 _pDouble1->get(0), _pDouble1->getImg(0),
490                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
491         }
492         else if (bComplex1 == true && bComplex2 == true)
493         {
494             iSubstractComplexMatrixToComplexScalar(
495                 _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getDimsArray(), _pDouble2->getDims(),
496                 _pDouble1->get(0), _pDouble1->getImg(0),
497                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
498         }
499
500         return 0;
501     }
502
503     if (_pDouble2->isScalar())
504     {
505         //add pL with each element of pR
506         (*_pDoubleOut) = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray(), bComplex1 || bComplex2);
507         if (bComplex1 == false && bComplex2 == false)
508         {
509             iSubstractRealScalarToRealMatrix(
510                 _pDouble2->get(0),
511                 _pDouble1->get(), _pDouble1->getDimsArray(), _pDouble1->getDims(),
512                 (*_pDoubleOut)->get());
513         }
514         else if (bComplex1 == false && bComplex2 == true)
515         {
516             iSubstractComplexScalarToRealMatrix(
517                 _pDouble2->get(0), _pDouble2->getImg(0),
518                 _pDouble1->get(), _pDouble1->getDimsArray(), _pDouble1->getDims(),
519                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
520         }
521         else if (bComplex1 == true && bComplex2 == false)
522         {
523             iSubstractRealScalarToComplexMatrix(
524                 _pDouble2->get(0),
525                 _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getDimsArray(), _pDouble1->getDims(),
526                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
527         }
528         else if (bComplex1 == true && bComplex2 == true)
529         {
530             iSubstractComplexScalarToComplexMatrix(
531                 _pDouble2->get(0), _pDouble2->getImg(0),
532                 _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getDimsArray(), _pDouble1->getDims(),
533                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
534         }
535
536         return 0;
537     }
538
539
540     //add pL and pR element wise
541
542     int iDims1 = _pDouble1->getDims();
543     int iDims2 = _pDouble2->getDims();
544
545     if (iDims1 != iDims2)
546     {
547         return 1;
548     }
549
550     int* piDims1    = _pDouble1->getDimsArray();
551     int* piDims2    = _pDouble2->getDimsArray();
552     for (int i = 0 ; i < _pDouble1->getDims() ; i++)
553     {
554         if (piDims1[i] != piDims2[i])
555         {
556             return 1;
557         }
558     }
559
560     (*_pDoubleOut) = new Double(iDims2, piDims2, bComplex1 || bComplex2);
561     if (bComplex1 == false && bComplex2 == false)
562     {
563         iSubstractRealMatrixToRealMatrix(
564             _pDouble2->get(),
565             _pDouble1->get(), piDims1, iDims1,
566             (*_pDoubleOut)->get());
567     }
568     else if (bComplex1 == false && bComplex2 == true)
569     {
570         iSubstractComplexMatrixToRealMatrix(
571             _pDouble2->get(), _pDouble2->getImg(),
572             _pDouble1->get(), piDims1, iDims1,
573             (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
574     }
575     else if (bComplex1 == true && bComplex2 == false)
576     {
577         iSubstractRealMatrixToComplexMatrix(
578             _pDouble2->get(),
579             _pDouble1->get(), _pDouble1->getImg(), piDims1, iDims1,
580             (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
581     }
582     else if (bComplex1 == true && bComplex2 == true)
583     {
584         iSubstractComplexMatrixToComplexMatrix(
585             _pDouble2->get(), _pDouble2->getImg(),
586             _pDouble1->get(), _pDouble1->getImg(), piDims1, iDims1,
587             (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
588     }
589
590     return 0;
591 }
592
593 //1 - %s
594 int SubstractPolyToDouble(Double *_pDouble, Polynom *_pPoly, Polynom** _pPolyOut)
595 {
596     double *pInDblR = _pDouble->get();
597     double *pInDblI = _pDouble->getImg();
598
599     if (_pDouble->isScalar())
600     {
601         //Clone original poly
602         (*_pPolyOut) = _pPoly->clone()->getAs<Polynom>();
603
604         if (_pDouble->isComplex())
605         {
606             (*_pPolyOut)->setComplex(true);
607         }
608
609         for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
610         {
611             SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
612             double *pOutPolyR       = pOutPoly->get();
613
614             pOutPolyR[0] = pInDblR[0] - pOutPolyR[0];
615
616             for (int j = 1 ; j < pOutPoly->getSize() ; j++)
617             {
618                 pOutPolyR[j] = - pOutPolyR[j];
619             }
620         }
621
622         if ((*_pPolyOut)->isComplex())
623         {
624             for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
625             {
626                 SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
627                 double *pOutPolyI       = pOutPoly->getImg();
628
629                 pOutPolyI[0]            = (pInDblI == NULL ? 0 : pInDblI[0]) - (pOutPolyI == NULL ? 0 : pOutPolyI[0]);
630
631                 for (int j = 1 ; j < pOutPoly->getSize() ; j++)
632                 {
633                     pOutPolyI[j] = - pOutPolyI[j];
634                 }
635             }
636         }
637
638         return 0;
639     }
640
641     if (_pPoly->isScalar())
642     {
643         //[] - %s
644         int *piRank = new int[_pDouble->getSize()];
645         for (int i = 0 ; i < _pDouble->getSize() ; i++)
646         {
647             piRank[i] = _pPoly->get(0)->getRank();
648         }
649
650         (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pDouble->getDims(), _pDouble->getDimsArray(), piRank);
651
652         for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
653         {
654             SinglePoly *pInPoly   = _pPoly->get(0);
655             SinglePoly *pOutPoly  = (*_pPolyOut)->get(i);
656
657             double *pInPolyR = pInPoly->get();
658             double *pOutPolyR = pOutPoly->get();
659
660             pOutPolyR[0] = pInDblR[i] - pInPolyR[0];
661
662             for (int j = 1 ; j < pOutPoly->getSize() ; j++)
663             {
664                 pOutPolyR[j] = - pInPolyR[j];
665             }
666         }
667
668         if (_pPoly->isComplex() || _pDouble->isComplex())
669         {
670             (*_pPolyOut)->setComplex(true);
671             for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
672             {
673                 SinglePoly *pInPoly   = _pPoly->get(0);
674                 SinglePoly *pOutPoly  = (*_pPolyOut)->get(i);
675
676                 double *pInPolyI = pInPoly->getImg();
677                 double *pOutPolyI = pOutPoly->getImg();
678
679                 pOutPolyI[0] = (pInDblI != NULL ? pInDblI[i] : 0) - (pInPolyI != NULL ? pInPolyI[0] : 0);
680
681                 for (int j = 1 ; j < pOutPoly->getSize() ; j++)
682                 {
683                     pOutPolyI[j] = -pInPolyI[j];
684                 }
685             }
686         }
687
688         return 0;
689     }
690
691     //check dimension compatibilities
692     int iLeftDims = _pDouble->getDims();
693     int iRightDims = _pPoly->getDims();
694
695     if (iLeftDims != iRightDims)
696     {
697         return 1;
698     }
699
700     int* piLeftDims = _pDouble->getDimsArray();
701     int* piRightDims = _pPoly->getDimsArray();
702
703     for (int i = 0 ; i < iLeftDims ; i++)
704     {
705         if (piLeftDims[i] != piRightDims[i])
706         {
707             return 1;
708         }
709     }
710
711
712     //Clone original poly
713     (*_pPolyOut) = _pPoly->clone()->getAs<Polynom>();
714
715     for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
716     {
717         SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
718         double *pOutPolyR       = pOutPoly->get();
719
720         pOutPolyR[0] = pInDblR[i] - pOutPolyR[0];
721
722         for (int j = 1 ; j < pOutPoly->getSize() ; j++)
723         {
724             pOutPolyR[j] = - pOutPolyR[j];
725         }
726     }
727
728     if (_pPoly->isComplex() || _pDouble->isComplex())
729     {
730         (*_pPolyOut)->setComplex(true);
731         for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
732         {
733             SinglePoly *pOutPoly  = (*_pPolyOut)->get(i);
734             double *pOutPolyI = pOutPoly->getImg();
735
736             pOutPolyI[0] = (pInDblI != NULL ? pInDblI[i] : 0) - (pOutPolyI != NULL ? pOutPolyI[0] : 0);
737
738             for (int j = 1 ; j < pOutPoly->getSize() ; j++)
739             {
740                 pOutPolyI[j] = - pOutPolyI[j];
741             }
742         }
743     }
744
745     return 0;
746 }
747
748 //%s - 1
749 int SubstractDoubleToPoly(Polynom *_pPoly, Double *_pDouble, Polynom **_pPolyOut)
750 {
751     double *pInDblR   = _pDouble->get();
752     double *pInDblI   = _pDouble->getImg();
753
754     if (_pDouble->isScalar())
755     {
756         //Clone original poly
757         (*_pPolyOut) = _pPoly->clone()->getAs<Polynom>();
758         for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
759         {
760             SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
761             double *pOutPolyR       = pOutPoly->get();
762             pOutPolyR[0]   -= pInDblR[0];
763         }
764
765         if ((*_pPolyOut)->isComplex() || _pDouble->isComplex())
766         {
767             (*_pPolyOut)->setComplex(true);
768             for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
769             {
770                 SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
771                 double *pOutPolyI       = pOutPoly->getImg();
772
773                 pOutPolyI[0]            -=  (pInDblI == NULL ? 0 : pInDblI[0]);
774             }
775         }
776     }
777     else if (_pPoly->isScalar())
778     {
779         int *piRank = new int[_pDouble->getSize()];
780         for (int i = 0 ; i < _pDouble->getSize() ; i++)
781         {
782             piRank[i] = _pPoly->get(0)->getRank();
783         }
784
785         (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pDouble->getDims(), _pDouble->getDimsArray(), piRank);
786         (*_pPolyOut)->setComplex(_pPoly->isComplex() || _pDouble->isComplex());
787
788         if (_pPoly->isComplex())
789         {
790             SinglePoly *pInPoly = _pPoly->get(0);
791             double *pInPolyR    = pInPoly->get();
792             double *pInPolyI    = pInPoly->getImg();
793             for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
794             {
795                 SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
796                 double *pOutPolyR       = pOutPoly->get();
797                 double *pOutPolyI       = pOutPoly->getImg();
798
799                 pOutPolyR[0]            = pInPolyR[0] - pInDblR[i];
800                 pOutPolyI[0]            = pInPolyI[0];
801                 for (int j = 1 ; j < pOutPoly->getSize() ; j++)
802                 {
803                     pOutPolyR[j]        = pInPolyR[j];
804                     pOutPolyI[j]        = pInPolyI[j];
805                 }
806             }
807         }
808         else if (_pDouble->isComplex())
809         {
810             SinglePoly *pInPoly = _pPoly->get(0);
811             double *pInPolyR    = pInPoly->get();
812             for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
813             {
814                 SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
815                 double *pOutPolyR       = pOutPoly->get();
816                 double *pOutPolyI       = pOutPoly->getImg();
817
818                 pOutPolyR[0]            = pInPolyR[0] - pInDblR[i];
819                 pOutPolyI[0]            = - pInDblI[i];
820                 for (int j = 1 ; j < pOutPoly->getSize() ; j++)
821                 {
822                     pOutPolyR[j]        = pInPolyR[j];
823                     pOutPolyI[j]        = 0;
824                 }
825             }
826         }
827         else
828         {
829             SinglePoly *pInPoly = _pPoly->get(0);
830             double *pInPolyR    = pInPoly->get();
831             for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
832             {
833                 SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
834                 double *pOutPolyR       = pOutPoly->get();
835
836                 pOutPolyR[0]            = pInPolyR[0] - pInDblR[i];
837                 for (int j = 1 ; j < pOutPoly->getSize() ; j++)
838                 {
839                     pOutPolyR[j]        = pInPolyR[j];
840                 }
841             }
842         }
843
844         delete[] piRank;
845     }
846     else
847     {
848         //check dimension compatibilities
849         int iLeftDims = _pDouble->getDims();
850         int iRightDims = _pPoly->getDims();
851
852         if (iLeftDims != iRightDims)
853         {
854             return 1;
855         }
856
857         int* piLeftDims = _pDouble->getDimsArray();
858         int* piRightDims = _pPoly->getDimsArray();
859
860         for (int i = 0 ; i < iLeftDims ; i++)
861         {
862             if (piLeftDims[i] != piRightDims[i])
863             {
864                 return 1;
865             }
866         }
867
868         //Clone original poly
869         (*_pPolyOut) = _pPoly->clone()->getAs<Polynom>();
870         for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
871         {
872             SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
873             double *pOutPolyR       = pOutPoly->get();
874             pOutPolyR[0]            -= pInDblR[i];
875         }
876
877         if (_pPoly->isComplex() || _pDouble->isComplex())
878         {
879             (*_pPolyOut)->setComplex(true);
880             for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
881             {
882                 SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
883                 double *pOutPolyI       = pOutPoly->getImg();
884
885                 pOutPolyI[0]            -= (pInDblI != NULL ? pInDblI[i] : 0);
886             }
887         }
888     }
889     return 0;
890 }
891
892 int SubstractPolyToPoly(Polynom *_pPoly1, Polynom *_pPoly2, Polynom **_pPolyOut)
893 {
894     if (_pPoly1->isScalar())
895     {
896         SinglePoly* p1Coef  = _pPoly1->get(0);
897         int iRank1          = p1Coef->getRank();
898         int* pRank2         = new int[_pPoly2->getSize()];
899         int* pRankOut       = new int[_pPoly2->getSize()];
900
901         _pPoly2->getRank(pRank2);
902         for (int i = 0 ; i < _pPoly2->getSize() ; i++)
903         {
904             pRankOut[i] = std::max(iRank1, pRank2[i]);
905         }
906
907         (*_pPolyOut) = new Polynom(_pPoly2->getVariableName(), _pPoly2->getDims(), _pPoly2->getDimsArray(), pRankOut);
908         (*_pPolyOut)->setComplex(_pPoly1->isComplex() || _pPoly2->isComplex());
909
910         //Result P1(0) + P2(i)
911         double* p1R = p1Coef->get();
912         if ((*_pPolyOut)->isComplex())
913         {
914             double* p1I = p1Coef->getImg();
915             for (int i = 0 ; i < _pPoly2->getSize() ; i++)
916             {
917                 SinglePoly* p2Coef   = _pPoly2->get(i);
918                 double* p2R          = p2Coef->get();
919                 double* p2I          = p2Coef->getImg();
920
921                 SinglePoly* pOutCoef = (*_pPolyOut)->get(i);
922                 double* pOutR        = pOutCoef->get();
923                 double* pOutI        = pOutCoef->getImg();
924
925                 for (int j = 0 ; j < pRankOut[i] + 1 ; j++)
926                 {
927                     if (j > iRank1)
928                     {
929                         pOutR[j] = - p2R[j];
930                         pOutI[j] = - (p2I ? p2I[j] : 0);
931                     }
932                     else if (j > pRank2[i])
933                     {
934                         pOutR[j] = p1R[j];
935                         pOutI[j] = (p1I ? p1I[j] : 0);
936                     }
937                     else
938                     {
939                         pOutR[j] = p1R[j] - p2R[j];
940                         pOutI[j] = (p1I ? p1I[j] : 0) - (p2I ? p2I[j] : 0);
941                     }
942                 }
943             }
944         }
945         else
946         {
947             for (int i = 0 ; i < _pPoly2->getSize() ; i++)
948             {
949                 SinglePoly* p2Coef   = _pPoly2->get(i);
950                 double* p2R          = p2Coef->get();
951
952                 SinglePoly* pOutCoef = (*_pPolyOut)->get(i);
953                 double* pOutR        = pOutCoef->get();
954
955                 for (int j = 0 ; j < pRankOut[i] + 1 ; j++)
956                 {
957                     if (j > iRank1)
958                     {
959                         pOutR[j] = - p2R[j];
960                     }
961                     else if (j > pRank2[i])
962                     {
963                         pOutR[j] = p1R[j];
964                     }
965                     else
966                     {
967                         pOutR[j] = p1R[j] - p2R[j];
968                     }
969                 }
970             }
971         }
972
973         delete[] pRankOut;
974         delete[] pRank2;
975         (*_pPolyOut)->updateRank();
976         return 0;
977     }
978
979     if (_pPoly2->isScalar())
980     {
981         //size(p2) == 1
982         int *pRank1     = new int[_pPoly1->getSize()];
983         int iRank2      = _pPoly2->get(0)->getRank();
984         int *pRankOut   = new int[_pPoly1->getSize()];
985
986         _pPoly1->getRank(pRank1);
987         for (int i = 0 ; i < _pPoly1->getSize() ; i++)
988         {
989             pRankOut[i] = std::max(pRank1[i], iRank2);
990         }
991
992         (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), _pPoly1->getDims(), _pPoly1->getDimsArray(), pRankOut);
993         (*_pPolyOut)->setComplex(_pPoly1->isComplex() || _pPoly2->isComplex());
994
995         //Result P1(i) + P2(0)
996         SinglePoly *p2Coef          = _pPoly2->get(0);
997         double *p2R                 = p2Coef->get();
998
999         if ((*_pPolyOut)->isComplex())
1000         {
1001             double *p2I                 = p2Coef->getImg();
1002             for (int i = 0 ; i < _pPoly1->getSize() ; i++)
1003             {
1004                 SinglePoly *p1Coef      = _pPoly1->get(i);
1005                 double *p1R             = p1Coef->get();
1006                 double *p1I             = p1Coef->getImg();
1007
1008                 SinglePoly *pOutCoef    = (*_pPolyOut)->get(i);
1009                 double *pOutR           = pOutCoef->get();
1010                 double *pOutI           = pOutCoef->getImg();
1011
1012                 for (int j = 0 ; j < pRankOut[i] + 1 ; j++)
1013                 {
1014                     if (j > pRank1[j])
1015                     {
1016                         pOutR[j] = - p2R[j];
1017                         pOutI[j] = - (p2I ? p2I[j] : 0);
1018                     }
1019                     else if (j > iRank2)
1020                     {
1021                         pOutR[j] = p1R[j];
1022                         pOutI[j] = (p1I ? p1I[j] : 0);
1023                     }
1024                     else
1025                     {
1026                         pOutR[j] = p1R[j] - p2R[j];
1027                         pOutI[j] = (p1I ? p1I[j] : 0) - (p2I ? p2I[j] : 0);
1028                     }
1029                 }
1030             }
1031         }
1032         else
1033         {
1034             for (int i = 0 ; i < _pPoly1->getSize() ; i++)
1035             {
1036                 SinglePoly *p1Coef      = _pPoly1->get(i);
1037                 double *p1R             = p1Coef->get();
1038
1039                 SinglePoly *pOutCoef    = (*_pPolyOut)->get(i);
1040                 double *pOutR           = pOutCoef->get();
1041
1042                 for (int j = 0 ; j < pRankOut[i] + 1 ; j++)
1043                 {
1044                     if (j > pRank1[j])
1045                     {
1046                         pOutR[j] = - p2R[j];
1047                     }
1048                     else if (j > iRank2)
1049                     {
1050                         pOutR[j] = p1R[j];
1051                     }
1052                     else
1053                     {
1054                         pOutR[j] = p1R[j] - p2R[j];
1055                     }
1056                 }
1057             }
1058         }
1059
1060         delete[] pRankOut;
1061         delete[] pRank1;
1062         (*_pPolyOut)->updateRank();
1063         return 0;
1064     }
1065
1066     //check dimension compatibilities
1067     int iDims1  = _pPoly1->getDims();
1068     int iDims2  = _pPoly2->getDims();
1069
1070     if (iDims1 != iDims2)
1071     {
1072         return 1;
1073     }
1074
1075     int* piDims1    = _pPoly1->getDimsArray();
1076     int* piDims2    = _pPoly2->getDimsArray();
1077
1078     for (int i = 0 ; i < iDims1 ; i++)
1079     {
1080         if (piDims1[i] != piDims2[i])
1081         {
1082             return 1;
1083         }
1084     }
1085
1086     int *pRankOut   = new int[_pPoly1->getSize()];
1087     int *pRank1     = new int[_pPoly1->getSize()];
1088     int *pRank2     = new int[_pPoly2->getSize()];
1089
1090     _pPoly1->getRank(pRank1);
1091     _pPoly2->getRank(pRank2);
1092     for (int i = 0 ; i < _pPoly1->getSize() ; i++)
1093     {
1094         pRankOut[i] = std::max(pRank1[i], pRank2[i]);
1095     }
1096
1097     (*_pPolyOut) = new Polynom(_pPoly2->getVariableName(), iDims1, piDims1, pRankOut);
1098     (*_pPolyOut)->setComplex(_pPoly1->isComplex() || _pPoly2->isComplex());
1099
1100     //Result P1(i) + P2(i)
1101     for (int i = 0 ; i < _pPoly1->getSize() ; i++)
1102     {
1103         double *p1R     = _pPoly1->get(i)->get();
1104         double *p2R     = _pPoly2->get(i)->get();
1105         double *pOutR   = (*_pPolyOut)->get(i)->get();
1106         int iMin        = std::min(pRank1[i], pRank2[i]);
1107         int iMax        = std::max(pRank1[i], pRank2[i]);
1108
1109         for (int j = 0 ; j < iMin + 1 ; j++)
1110         {
1111             pOutR[j]    = p1R[j] - p2R[j];
1112         }
1113
1114         double *pTemp   = NULL;
1115         int iCoef       = 1;
1116         if (pRank1[i] > pRank2[i])
1117         {
1118             pTemp       = p1R;
1119             iCoef       = 1;
1120         }
1121         else
1122         {
1123             pTemp       = p2R;
1124             iCoef       = -1;
1125         }
1126
1127         for (int j = iMin + 1 ; j < iMax + 1 ; j++)
1128         {
1129             pOutR[j]    = pTemp[j] * iCoef;
1130         }
1131
1132         if ((*_pPolyOut)->isComplex())
1133         {
1134             double *p1I     = _pPoly1->get(i)->getImg();
1135             double *p2I     = _pPoly2->get(i)->getImg();
1136             double *pOutI   = (*_pPolyOut)->get(i)->getImg();
1137
1138             for (int j = 0 ; j < iMin + 1 ; j++)
1139             {
1140                 pOutI[j]    = (p1I == NULL ? 0 : p1I[j]) - (p2I == NULL ? 0 : p2I[j]);
1141             }
1142
1143             for (int j = iMin + 1 ; j < iMax + 1 ; j++)
1144             {
1145                 pOutI[j]  = pTemp[j] * iCoef;
1146             }
1147         }
1148     }
1149
1150     delete[] pRankOut;
1151     delete[] pRank1;
1152     delete[] pRank2;
1153
1154     (*_pPolyOut)->updateRank();
1155     return 0;
1156 }
1157
1158 //_pSparse2 - _pSparse1
1159 int SubstractSparseToSparse(Sparse* _pSparse1, Sparse* _pSparse2, GenericType **_pSparseOut)
1160 {
1161     //check scalar hidden in a sparse ;)
1162     if (_pSparse1->getRows() == 1 && _pSparse1->getCols() == 1)
1163     {
1164         //do scalar - sp
1165         Double* pDbl = NULL;
1166         if (_pSparse1->isComplex())
1167         {
1168             std::complex<double> dbl = _pSparse1->getImg(0, 0);
1169             pDbl = new Double(dbl.real(), dbl.imag());
1170         }
1171         else
1172         {
1173             pDbl = new Double(_pSparse1->get(0, 0));
1174         }
1175
1176         SubstractSparseToDouble(_pSparse2, pDbl, (GenericType**)_pSparseOut);
1177         delete pDbl;
1178         return 0;
1179     }
1180
1181     if (_pSparse2->getRows() == 1 && _pSparse2->getCols() == 1)
1182     {
1183         //do sp - scalar
1184         Double* pDbl = NULL;
1185         if (_pSparse2->isComplex())
1186         {
1187             std::complex<double> dbl = _pSparse2->getImg(0, 0);
1188             pDbl = new Double(dbl.real(), dbl.imag());
1189         }
1190         else
1191         {
1192             pDbl = new Double(_pSparse2->get(0, 0));
1193         }
1194
1195         SubstractDoubleToSparse(pDbl, _pSparse1, (GenericType**)_pSparseOut);
1196         delete pDbl;
1197         return 0;
1198     }
1199
1200     if (_pSparse1->getRows() != _pSparse2->getRows() || _pSparse1->getCols() != _pSparse2->getCols())
1201     {
1202         //dimensions not match
1203         return 1;
1204     }
1205
1206     if (_pSparse1->nonZeros() == 0)
1207     {
1208         //sp - sp([])
1209         *_pSparseOut = new Sparse(*_pSparse2);
1210         return 0;
1211     }
1212
1213     if (_pSparse2->nonZeros() == 0)
1214     {
1215         //sp([]) - sp
1216         Sparse* pOut = new Sparse(*_pSparse1);
1217         pOut->opposite();
1218         *_pSparseOut = pOut;
1219         return 0;
1220     }
1221
1222     //copy _pSparse1 in _pSparseOut
1223     *_pSparseOut = _pSparse2->substract(*_pSparse1);
1224     return 0;
1225 }
1226
1227 int SubstractSparseToDouble(Sparse* _pSparse, Double* _pDouble, GenericType **_pOut)
1228 {
1229     //D - SP
1230     int iOne = 1; //fortran
1231     bool bComplex1 = _pSparse->isComplex();
1232     bool bComplex2 = _pDouble->isComplex();
1233
1234     if (_pDouble->isIdentity())
1235     {
1236         //convert to sp
1237         Sparse* pS = new Sparse(_pSparse->getRows(), _pSparse->getCols(), _pDouble->isComplex());
1238         if (pS->isComplex())
1239         {
1240             for (int i = 0 ; i < std::min(_pSparse->getRows() , _pSparse->getCols()) ; i++)
1241             {
1242                 pS->set(i, i, std::complex<double>(_pDouble->get(0), _pDouble->getImg(0)));
1243             }
1244         }
1245         else
1246         {
1247             for (int i = 0 ; i < std::min(_pSparse->getRows() , _pSparse->getCols()) ; i++)
1248             {
1249                 pS->set(i, i, _pDouble->get(0));
1250             }
1251         }
1252
1253         SubstractSparseToSparse(_pSparse, pS, _pOut);
1254         delete pS;
1255         return 0;
1256     }
1257
1258     if (_pSparse->isScalar() && _pDouble->isScalar())
1259     {
1260         //d - sp
1261         Double* pRes = _pDouble->clone()->getAs<Double>();
1262         pRes->setComplex(bComplex1 | bComplex2);
1263         if (bComplex1)
1264         {
1265             std::complex<double> dbl = _pSparse->getImg(0, 0);
1266             pRes->set(0, pRes->get(0) - dbl.real());
1267             pRes->setImg(0, pRes->getImg(0) - dbl.imag());
1268         }
1269         else
1270         {
1271             pRes->set(0, pRes->get(0) - _pSparse->get(0, 0));
1272         }
1273
1274         *_pOut = pRes;
1275         return 0;
1276     }
1277
1278     if (_pDouble->isScalar())
1279     {
1280         //d - SP
1281         Double* pRes = new Double(_pSparse->getRows(), _pSparse->getCols(), bComplex1 | bComplex2);
1282         int iSize = _pSparse->getSize();
1283         double dblVal = _pDouble->get(0);
1284         C2F(dset)(&iSize, &dblVal, pRes->get(), &iOne);
1285         if (bComplex2)
1286         {
1287             double dblValI = _pDouble->getImg(0);
1288             C2F(dset)(&iSize, &dblValI, pRes->getImg(), &iOne);
1289         }
1290         else if (bComplex1)
1291         {
1292             //initialize imag part at 0
1293             double dblValI = 0;
1294             C2F(dset)(&iSize, &dblValI, pRes->getImg(), &iOne);
1295         }
1296
1297         int nonZeros = static_cast<int>(_pSparse->nonZeros());
1298         int* pRows = new int[nonZeros * 2];
1299         _pSparse->outputRowCol(pRows);
1300         int* pCols = pRows + nonZeros;
1301
1302         if (bComplex1)
1303         {
1304             for (int i = 0 ; i < nonZeros ; i++)
1305             {
1306                 int iRow = static_cast<int>(pRows[i]) - 1;
1307                 int iCol = static_cast<int>(pCols[i]) - 1;
1308                 std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1309                 pRes->set(iRow, iCol, pRes->get(iRow, iCol) - dbl.real());
1310                 pRes->setImg(iRow, iCol, pRes->getImg(iRow, iCol) - dbl.imag());
1311             }
1312         }
1313         else
1314         {
1315             for (int i = 0 ; i < nonZeros ; i++)
1316             {
1317                 int iRow = static_cast<int>(pRows[i]) - 1;
1318                 int iCol = static_cast<int>(pCols[i]) - 1;
1319                 pRes->set(iRow, iCol, pRes->get(iRow, iCol) - _pSparse->get(iRow, iCol));
1320             }
1321         }
1322         *_pOut = pRes;
1323
1324         //clear
1325         delete[] pRows;
1326         return 0;
1327     }
1328
1329     if (_pSparse->isScalar())
1330     {
1331         //D - sp
1332         Double* pRes = _pDouble->clone()->getAs<Double>();
1333         pRes->setComplex(bComplex1 | bComplex2);
1334
1335         if (bComplex1)
1336         {
1337             double* pReal = pRes->get();
1338             double* pImg = pRes->getImg();
1339             for (int i = 0 ; i < pRes->getSize() ; i++)
1340             {
1341                 std::complex<double> dbl = _pSparse->getImg(0, 0);
1342                 pReal[i] -= dbl.real();
1343                 pImg[i] -= dbl.imag();
1344             }
1345         }
1346         else
1347         {
1348             double* pReal = pRes->get();
1349             for (int i = 0 ; i < pRes->getSize() ; i++)
1350             {
1351                 pReal[i] -= _pSparse->get(0, 0);
1352             }
1353         }
1354         *_pOut = pRes;
1355         return 0;
1356     }
1357
1358
1359     if (_pSparse->getRows() == _pDouble->getRows() && _pSparse->getCols() == _pDouble->getCols())
1360     {
1361         //D - SP
1362         Double* pRes = _pDouble->clone()->getAs<Double>();
1363         pRes->setComplex(bComplex1 | bComplex2);
1364
1365         int nonZeros = static_cast<int>(_pSparse->nonZeros());
1366         int* pRows = new int[nonZeros * 2];
1367         _pSparse->outputRowCol(pRows);
1368         int* pCols = pRows + nonZeros;
1369
1370         if (bComplex1)
1371         {
1372             for (int i = 0 ; i < nonZeros ; i++)
1373             {
1374                 int iRow = static_cast<int>(pRows[i]) - 1;
1375                 int iCol = static_cast<int>(pCols[i]) - 1;
1376                 std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1377                 pRes->set(iRow, iCol, pRes->get(iRow, iCol) - dbl.real());
1378                 pRes->setImg(iRow, iCol, pRes->getImg(iRow, iCol) - dbl.imag());
1379             }
1380         }
1381         else
1382         {
1383             for (int i = 0 ; i < nonZeros ; i++)
1384             {
1385                 int iRow = static_cast<int>(pRows[i]) - 1;
1386                 int iCol = static_cast<int>(pCols[i]) - 1;
1387                 pRes->set(iRow, iCol, pRes->get(iRow, iCol) - _pSparse->get(iRow, iCol));
1388             }
1389         }
1390
1391         //clear
1392         delete[] pRows;
1393         *_pOut = pRes;
1394         return 0;
1395     }
1396     return 1;
1397 }
1398
1399 int SubstractDoubleToSparse(Double* _pDouble, Sparse* _pSparse, GenericType **_pOut)
1400 {
1401     //SP - D => -(D - SP)
1402     int iRet = SubstractSparseToDouble(_pSparse, _pDouble, _pOut);
1403     if (iRet)
1404     {
1405         return iRet;
1406     }
1407
1408     //compute opposite of result
1409     if ((*_pOut)->isDouble())
1410     {
1411         Double* pD          = (*_pOut)->getAs<Double>();
1412         int iSize           = pD->getSize();
1413         double dblMinusOne  = -1;
1414         int iOne            = 1;
1415
1416         C2F(dscal)(&iSize, &dblMinusOne, pD->get(), &iOne);
1417         if (pD->isComplex())
1418         {
1419             C2F(dscal)(&iSize, &dblMinusOne, pD->getImg(), &iOne);
1420         }
1421     }
1422     else if ((*_pOut)->isSparse())
1423     {
1424         Sparse* pS = (*_pOut)->getAs<Sparse>();
1425         pS->opposite();
1426     }
1427     else
1428     {
1429         return 1;
1430     }
1431
1432     return 0;
1433 }
1434