* Bug 16483 fixed: substraction of complex polynomial matrix was broken
[scilab.git] / scilab / modules / ast / src / cpp / operations / types_subtraction.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  * Copyright (C) 2012 - 2016 - Scilab Enterprises
7  *
8  * This file is hereby licensed under the terms of the GNU GPL v2.0,
9  * pursuant to article 5.3.4 of the CeCILL v.2.1.
10  * This file was originally licensed under the terms of the CeCILL v2.1,
11  * and continues to be available under such terms.
12  * For more information, see the COPYING file which you should have received
13  * along with this program.
14 *
15 */
16
17 #include <algorithm>
18
19 #include "types_subtraction.hxx"
20 #include "types_opposite.hxx"
21 #include "double.hxx"
22 #include "polynom.hxx"
23 #include "int.hxx"
24 #include "sparse.hxx"
25 #include "generic_operations.hxx"
26
27 extern "C"
28 {
29 #include "elem_common.h" //dset
30 #include "Sciwarning.h" //Sciwarning
31 }
32
33 using namespace types;
34
35 //define arrays on operation functions
36 static sub_function pSubfunction[types::InternalType::IdLast][types::InternalType::IdLast] = {NULL};
37
38 void fillSubtractFunction()
39 {
40 #define scilab_fill_sub(id1, id2, func, typeIn1, typeIn2, typeOut) \
41     pSubfunction[types::InternalType::Id ## id1][types::InternalType::Id ## id2] = (sub_function)&sub_##func<typeIn1, typeIn2, typeOut>
42
43     //Double
44     //Matrix - Matrix
45     scilab_fill_sub(Double, Double, M_M, Double, Double, Double);
46     scilab_fill_sub(Double, Int8, M_M, Double, Int8, Int8);
47     scilab_fill_sub(Double, UInt8, M_M, Double, UInt8, UInt8);
48     scilab_fill_sub(Double, Int16, M_M, Double, Int16, Int16);
49     scilab_fill_sub(Double, UInt16, M_M, Double, UInt16, UInt16);
50     scilab_fill_sub(Double, Int32, M_M, Double, Int32, Int32);
51     scilab_fill_sub(Double, UInt32, M_M, Double, UInt32, UInt32);
52     scilab_fill_sub(Double, Int64, M_M, Double, Int64, Int64);
53     scilab_fill_sub(Double, UInt64, M_M, Double, UInt64, UInt64);
54     scilab_fill_sub(Double, Bool, M_M, Double, Bool, Double);
55     scilab_fill_sub(Double, Polynom, M_M, Double, Polynom, Polynom);
56     scilab_fill_sub(Double, Sparse, M_M, Double, Sparse, Double);
57
58     //Matrix - Matrix Complex
59     scilab_fill_sub(Double, DoubleComplex, M_MC, Double, Double, Double);
60     scilab_fill_sub(Double, PolynomComplex, M_M, Double, Polynom, Polynom);
61     scilab_fill_sub(Double, SparseComplex, M_M, Double, Sparse, Double);
62
63     //Matrix - Scalar
64     scilab_fill_sub(Double, ScalarDouble, M_S, Double, Double, Double);
65     scilab_fill_sub(Double, ScalarInt8, M_S, Double, Int8, Int8);
66     scilab_fill_sub(Double, ScalarUInt8, M_S, Double, UInt8, UInt8);
67     scilab_fill_sub(Double, ScalarInt16, M_S, Double, Int16, Int16);
68     scilab_fill_sub(Double, ScalarUInt16, M_S, Double, UInt16, UInt16);
69     scilab_fill_sub(Double, ScalarInt32, M_S, Double, Int32, Int32);
70     scilab_fill_sub(Double, ScalarUInt32, M_S, Double, UInt32, UInt32);
71     scilab_fill_sub(Double, ScalarInt64, M_S, Double, Int64, Int64);
72     scilab_fill_sub(Double, ScalarUInt64, M_S, Double, UInt64, UInt64);
73     scilab_fill_sub(Double, ScalarBool, M_S, Double, Bool, Double);
74     scilab_fill_sub(Double, ScalarPolynom, M_M, Double, Polynom, Polynom);
75
76     //Matrix - Scalar Complex
77     scilab_fill_sub(Double, ScalarDoubleComplex, M_SC, Double, Double, Double);
78     scilab_fill_sub(Double, ScalarPolynomComplex, M_M, Double, Polynom, Polynom);
79     //Matrix - Empty
80     scilab_fill_sub(Double, Empty, M_E, Double, Double, Double);
81
82
83     //Matrix Complex - Matrix
84     scilab_fill_sub(DoubleComplex, Double, MC_M, Double, Double, Double);
85     scilab_fill_sub(DoubleComplex, DoubleComplex, MC_MC, Double, Double, Double);
86     scilab_fill_sub(DoubleComplex, ScalarDouble, MC_S, Double, Double, Double);
87     scilab_fill_sub(DoubleComplex, ScalarDoubleComplex, MC_SC, Double, Double, Double);
88     scilab_fill_sub(DoubleComplex, Empty, MC_E, Double, Double, Double);
89     scilab_fill_sub(DoubleComplex, Polynom, M_M, Double, Polynom, Polynom);
90     scilab_fill_sub(DoubleComplex, PolynomComplex, M_M, Double, Polynom, Polynom);
91     scilab_fill_sub(DoubleComplex, ScalarPolynom, M_M, Double, Polynom, Polynom);
92     scilab_fill_sub(DoubleComplex, ScalarPolynomComplex, M_M, Double, Polynom, Polynom);
93     scilab_fill_sub(DoubleComplex, Sparse, M_M, Double, Sparse, Double);
94     scilab_fill_sub(DoubleComplex, SparseComplex, M_M, Double, Sparse, Double);
95
96     //Scalar - Matrix
97     scilab_fill_sub(ScalarDouble, Double, S_M, Double, Double, Double);
98     scilab_fill_sub(ScalarDouble, Int8, S_M, Double, Int8, Int8);
99     scilab_fill_sub(ScalarDouble, UInt8, S_M, Double, UInt8, UInt8);
100     scilab_fill_sub(ScalarDouble, Int16, S_M, Double, Int16, Int16);
101     scilab_fill_sub(ScalarDouble, UInt16, S_M, Double, UInt16, UInt16);
102     scilab_fill_sub(ScalarDouble, Int32, S_M, Double, Int32, Int32);
103     scilab_fill_sub(ScalarDouble, UInt32, S_M, Double, UInt32, UInt32);
104     scilab_fill_sub(ScalarDouble, Int64, S_M, Double, Int64, Int64);
105     scilab_fill_sub(ScalarDouble, UInt64, S_M, Double, UInt64, UInt64);
106     scilab_fill_sub(ScalarDouble, Bool, S_M, Double, Bool, Double);
107     scilab_fill_sub(ScalarDouble, Polynom, M_M, Double, Polynom, Polynom);
108     scilab_fill_sub(ScalarDouble, Sparse, M_M, Double, Sparse, Double);
109
110     //Scalar - Matrix Complex
111     scilab_fill_sub(ScalarDouble, DoubleComplex, S_MC, Double, Double, Double);
112     scilab_fill_sub(ScalarDouble, PolynomComplex, M_M, Double, Polynom, Polynom);
113     scilab_fill_sub(ScalarDouble, SparseComplex, M_M, Double, Sparse, Double);
114
115     //Scalar - Scalar
116     scilab_fill_sub(ScalarDouble, ScalarDouble, S_S, Double, Double, Double);
117     scilab_fill_sub(ScalarDouble, ScalarInt8, S_S, Double, Int8, Int8);
118     scilab_fill_sub(ScalarDouble, ScalarUInt8, S_S, Double, UInt8, UInt8);
119     scilab_fill_sub(ScalarDouble, ScalarInt16, S_S, Double, Int16, Int16);
120     scilab_fill_sub(ScalarDouble, ScalarUInt16, S_S, Double, UInt16, UInt16);
121     scilab_fill_sub(ScalarDouble, ScalarInt32, S_S, Double, Int32, Int32);
122     scilab_fill_sub(ScalarDouble, ScalarUInt32, S_S, Double, UInt32, UInt32);
123     scilab_fill_sub(ScalarDouble, ScalarInt64, S_S, Double, Int64, Int64);
124     scilab_fill_sub(ScalarDouble, ScalarUInt64, S_S, Double, UInt64, UInt64);
125     scilab_fill_sub(ScalarDouble, ScalarBool, S_S, Double, Bool, Double);
126     scilab_fill_sub(ScalarDouble, ScalarPolynom, M_M, Double, Polynom, Polynom);
127
128     //Scalar - Scalar Complex
129     scilab_fill_sub(ScalarDouble, ScalarDoubleComplex, S_SC, Double, Double, Double);
130     scilab_fill_sub(ScalarDouble, PolynomComplex, M_M, Double, Polynom, Polynom);
131     scilab_fill_sub(ScalarDouble, ScalarPolynomComplex, M_M, Double, Polynom, Polynom);
132
133     //Scalar - Empty
134     scilab_fill_sub(ScalarDouble, Empty, S_E, Double, Double, Double);
135
136     //Scalar Complex - Matrix
137     scilab_fill_sub(ScalarDoubleComplex, Double, SC_M, Double, Double, Double);
138     scilab_fill_sub(ScalarDoubleComplex, Polynom, M_M, Double, Polynom, Polynom);
139     scilab_fill_sub(ScalarDoubleComplex, Sparse, M_M, Double, Sparse, Double);
140     //Scalar Complex - Matrix Complex
141     scilab_fill_sub(ScalarDoubleComplex, DoubleComplex, SC_MC, Double, Double, Double);
142     scilab_fill_sub(ScalarDoubleComplex, PolynomComplex, M_M, Double, Polynom, Polynom);
143     scilab_fill_sub(ScalarDoubleComplex, SparseComplex, M_M, Double, Sparse, Double);
144     //Scalar Complex - Scalar
145     scilab_fill_sub(ScalarDoubleComplex, ScalarDouble, SC_S, Double, Double, Double);
146     scilab_fill_sub(ScalarDoubleComplex, ScalarPolynom, M_M, Double, Polynom, Polynom);
147     //Scalar Complex - Scalar Complex
148     scilab_fill_sub(ScalarDoubleComplex, ScalarDoubleComplex, SC_SC, Double, Double, Double);
149     scilab_fill_sub(ScalarDoubleComplex, ScalarPolynomComplex, M_M, Double, Polynom, Polynom);
150     //Scalar Complex - Empty
151     scilab_fill_sub(ScalarDoubleComplex, Empty, SC_E, Double, Double, Double);
152
153     //Empty - Matrix
154     scilab_fill_sub(Empty, Double, E_M, Double, Double, Double);
155     scilab_fill_sub(Empty, Int8, E_M, Double, Int8, Int8);
156     scilab_fill_sub(Empty, UInt8, E_M, Double, UInt8, UInt8);
157     scilab_fill_sub(Empty, Int16, E_M, Double, Int16, Int16);
158     scilab_fill_sub(Empty, UInt16, E_M, Double, UInt16, UInt16);
159     scilab_fill_sub(Empty, Int32, E_M, Double, Int32, Int32);
160     scilab_fill_sub(Empty, UInt32, E_M, Double, UInt32, UInt32);
161     scilab_fill_sub(Empty, Int64, E_M, Double, Int64, Int64);
162     scilab_fill_sub(Empty, UInt64, E_M, Double, UInt64, UInt64);
163     //scilab_fill_sub(Empty, Bool, E_M, Double, Bool, Double);
164     scilab_fill_sub(Empty, Polynom, E_M, Double, Polynom, Polynom);
165     scilab_fill_sub(Empty, PolynomComplex, E_MC, Double, Polynom, Polynom);
166     scilab_fill_sub(Empty, Sparse, E_M, Double, Sparse, Sparse);
167     scilab_fill_sub(Empty, SparseComplex, E_MC, Double, Sparse, Sparse);
168
169     //Empty - Matrix Complex
170     scilab_fill_sub(Empty, DoubleComplex, E_MC, Double, Double, Double);
171     //Empty - Scalar
172     scilab_fill_sub(Empty, ScalarDouble, E_M, Double, Double, Double);
173     scilab_fill_sub(Empty, ScalarInt8, E_M, Double, Int8, Int8);
174     scilab_fill_sub(Empty, ScalarUInt8, E_M, Double, UInt8, UInt8);
175     scilab_fill_sub(Empty, ScalarInt16, E_M, Double, Int16, Int16);
176     scilab_fill_sub(Empty, ScalarUInt16, E_M, Double, UInt16, UInt16);
177     scilab_fill_sub(Empty, ScalarInt32, E_M, Double, Int32, Int32);
178     scilab_fill_sub(Empty, ScalarUInt32, E_M, Double, UInt32, UInt32);
179     scilab_fill_sub(Empty, ScalarInt64, E_M, Double, Int64, Int64);
180     scilab_fill_sub(Empty, ScalarUInt64, E_M, Double, UInt64, UInt64);
181     scilab_fill_sub(Empty, ScalarBool, E_M, Double, Bool, Double);
182     scilab_fill_sub(Empty, ScalarPolynom, E_M, Double, Polynom, Polynom);
183
184     //Empty - Scalar Complex
185     scilab_fill_sub(Empty, ScalarDoubleComplex, E_MC, Double, Double, Double);
186     scilab_fill_sub(Empty, ScalarPolynomComplex, E_MC, Double, Polynom, Polynom);
187     //Empty - Empty
188     scilab_fill_sub(Empty, Empty, E_E, Double, Double, Double);
189     //Empty - eye
190     scilab_fill_sub(Empty, Identity, E_I, Double, Double, Double);
191     scilab_fill_sub(Empty, IdentityComplex, E_IC, Double, Double, Double);
192
193     //Matrix - Identity
194     scilab_fill_sub(Double, Identity, M_I, Double, Double, Double);
195     scilab_fill_sub(Double, IdentityComplex, M_IC, Double, Double, Double);
196     scilab_fill_sub(DoubleComplex, Identity, MC_I, Double, Double, Double);
197     scilab_fill_sub(DoubleComplex, IdentityComplex, MC_IC, Double, Double, Double);
198     scilab_fill_sub(ScalarDouble, Identity, S_I, Double, Double, Double);
199     scilab_fill_sub(ScalarDouble, IdentityComplex, S_IC, Double, Double, Double);
200     scilab_fill_sub(ScalarDoubleComplex, Identity, SC_I, Double, Double, Double);
201     scilab_fill_sub(ScalarDoubleComplex, IdentityComplex, SC_IC, Double, Double, Double);
202
203     //Int8
204     //Matrix - Matrix
205     scilab_fill_sub(Int8, Double, M_M, Int8, Double, Int8);
206     scilab_fill_sub(Int8, Int8, M_M, Int8, Int8, Int8);
207     scilab_fill_sub(Int8, UInt8, M_M, Int8, UInt8, UInt8);
208     scilab_fill_sub(Int8, Int16, M_M, Int8, Int16, Int16);
209     scilab_fill_sub(Int8, UInt16, M_M, Int8, UInt16, UInt16);
210     scilab_fill_sub(Int8, Int32, M_M, Int8, Int32, Int32);
211     scilab_fill_sub(Int8, UInt32, M_M, Int8, UInt32, UInt32);
212     scilab_fill_sub(Int8, Int64, M_M, Int8, Int64, Int64);
213     scilab_fill_sub(Int8, UInt64, M_M, Int8, UInt64, UInt64);
214     scilab_fill_sub(Int8, Bool, M_M, Int8, Bool, Int8);
215     scilab_fill_sub(Int8, Empty, M_E, Int8, Double, Int8);
216
217     //Matrix - Scalar
218     scilab_fill_sub(Int8, ScalarDouble, M_S, Int8, Double, Int8);
219     scilab_fill_sub(Int8, ScalarInt8, M_S, Int8, Int8, Int8);
220     scilab_fill_sub(Int8, ScalarUInt8, M_S, Int8, UInt8, UInt8);
221     scilab_fill_sub(Int8, ScalarInt16, M_S, Int8, Int16, Int16);
222     scilab_fill_sub(Int8, ScalarUInt16, M_S, Int8, UInt16, UInt16);
223     scilab_fill_sub(Int8, ScalarInt32, M_S, Int8, Int32, Int32);
224     scilab_fill_sub(Int8, ScalarUInt32, M_S, Int8, UInt32, UInt32);
225     scilab_fill_sub(Int8, ScalarInt64, M_S, Int8, Int64, Int64);
226     scilab_fill_sub(Int8, ScalarUInt64, M_S, Int8, UInt64, UInt64);
227     scilab_fill_sub(Int8, ScalarBool, M_S, Int8, Bool, Int8);
228
229     //Scalar - Matrix
230     scilab_fill_sub(ScalarInt8, Double, S_M, Int8, Double, Int8);
231     scilab_fill_sub(ScalarInt8, Int8, S_M, Int8, Int8, Int8);
232     scilab_fill_sub(ScalarInt8, UInt8, S_M, Int8, UInt8, UInt8);
233     scilab_fill_sub(ScalarInt8, Int16, S_M, Int8, Int16, Int16);
234     scilab_fill_sub(ScalarInt8, UInt16, S_M, Int8, UInt16, UInt16);
235     scilab_fill_sub(ScalarInt8, Int32, S_M, Int8, Int32, Int32);
236     scilab_fill_sub(ScalarInt8, UInt32, S_M, Int8, UInt32, UInt32);
237     scilab_fill_sub(ScalarInt8, Int64, S_M, Int8, Int64, Int64);
238     scilab_fill_sub(ScalarInt8, UInt64, S_M, Int8, UInt64, UInt64);
239     scilab_fill_sub(ScalarInt8, Bool, S_M, Int8, Bool, Int8);
240     scilab_fill_sub(ScalarInt8, Empty, S_E, Int8, Double, Int8);
241
242     //Scalar - Scalar
243     scilab_fill_sub(ScalarInt8, ScalarDouble, S_S, Int8, Double, Int8);
244     scilab_fill_sub(ScalarInt8, ScalarInt8, S_S, Int8, Int8, Int8);
245     scilab_fill_sub(ScalarInt8, ScalarUInt8, S_S, Int8, UInt8, UInt8);
246     scilab_fill_sub(ScalarInt8, ScalarInt16, S_S, Int8, Int16, Int16);
247     scilab_fill_sub(ScalarInt8, ScalarUInt16, S_S, Int8, UInt16, UInt16);
248     scilab_fill_sub(ScalarInt8, ScalarInt32, S_S, Int8, Int32, Int32);
249     scilab_fill_sub(ScalarInt8, ScalarUInt32, S_S, Int8, UInt32, UInt32);
250     scilab_fill_sub(ScalarInt8, ScalarInt64, S_S, Int8, Int64, Int64);
251     scilab_fill_sub(ScalarInt8, ScalarUInt64, S_S, Int8, UInt64, UInt64);
252     scilab_fill_sub(ScalarInt8, ScalarBool, S_S, Int8, Bool, Int8);
253
254     //UInt8
255     //Matrix - Matrix
256     scilab_fill_sub(UInt8, Double, M_M, UInt8, Double, UInt8);
257     scilab_fill_sub(UInt8, Int8, M_M, UInt8, Int8, UInt8);
258     scilab_fill_sub(UInt8, UInt8, M_M, UInt8, UInt8, UInt8);
259     scilab_fill_sub(UInt8, Int16, M_M, UInt8, Int16, UInt16);
260     scilab_fill_sub(UInt8, UInt16, M_M, UInt8, UInt16, UInt16);
261     scilab_fill_sub(UInt8, Int32, M_M, UInt8, Int32, UInt32);
262     scilab_fill_sub(UInt8, UInt32, M_M, UInt8, UInt32, UInt32);
263     scilab_fill_sub(UInt8, Int64, M_M, UInt8, Int64, UInt64);
264     scilab_fill_sub(UInt8, UInt64, M_M, UInt8, UInt64, UInt64);
265     scilab_fill_sub(UInt8, Bool, M_M, UInt8, Bool, UInt8);
266     scilab_fill_sub(UInt8, Empty, M_E, UInt8, Double, UInt8);
267
268     //Matrix - Scalar
269     scilab_fill_sub(UInt8, ScalarDouble, M_S, UInt8, Double, UInt8);
270     scilab_fill_sub(UInt8, ScalarInt8, M_S, UInt8, Int8, UInt8);
271     scilab_fill_sub(UInt8, ScalarUInt8, M_S, UInt8, UInt8, UInt8);
272     scilab_fill_sub(UInt8, ScalarInt16, M_S, UInt8, Int16, UInt16);
273     scilab_fill_sub(UInt8, ScalarUInt16, M_S, UInt8, UInt16, UInt16);
274     scilab_fill_sub(UInt8, ScalarInt32, M_S, UInt8, Int32, UInt32);
275     scilab_fill_sub(UInt8, ScalarUInt32, M_S, UInt8, UInt32, UInt32);
276     scilab_fill_sub(UInt8, ScalarInt64, M_S, UInt8, Int64, UInt64);
277     scilab_fill_sub(UInt8, ScalarUInt64, M_S, UInt8, UInt64, UInt64);
278     scilab_fill_sub(UInt8, ScalarBool, M_S, UInt8, Bool, UInt8);
279
280     //Scalar - Matrix
281     scilab_fill_sub(ScalarUInt8, Double, S_M, UInt8, Double, UInt8);
282     scilab_fill_sub(ScalarUInt8, Int8, S_M, UInt8, Int8, UInt8);
283     scilab_fill_sub(ScalarUInt8, UInt8, S_M, UInt8, UInt8, UInt8);
284     scilab_fill_sub(ScalarUInt8, Int16, S_M, UInt8, Int16, UInt16);
285     scilab_fill_sub(ScalarUInt8, UInt16, S_M, UInt8, UInt16, UInt16);
286     scilab_fill_sub(ScalarUInt8, Int32, S_M, UInt8, Int32, UInt32);
287     scilab_fill_sub(ScalarUInt8, UInt32, S_M, UInt8, UInt32, UInt32);
288     scilab_fill_sub(ScalarUInt8, Int64, S_M, UInt8, Int64, UInt64);
289     scilab_fill_sub(ScalarUInt8, UInt64, S_M, UInt8, UInt64, UInt64);
290     scilab_fill_sub(ScalarUInt8, Bool, S_M, UInt8, Bool, UInt8);
291     scilab_fill_sub(ScalarUInt8, Empty, S_E, UInt8, Double, UInt8);
292
293     //Scalar - Scalar
294     scilab_fill_sub(ScalarUInt8, ScalarDouble, S_S, UInt8, Double, UInt8);
295     scilab_fill_sub(ScalarUInt8, ScalarInt8, S_S, UInt8, Int8, UInt8);
296     scilab_fill_sub(ScalarUInt8, ScalarUInt8, S_S, UInt8, UInt8, UInt8);
297     scilab_fill_sub(ScalarUInt8, ScalarInt16, S_S, UInt8, Int16, UInt16);
298     scilab_fill_sub(ScalarUInt8, ScalarUInt16, S_S, UInt8, UInt16, UInt16);
299     scilab_fill_sub(ScalarUInt8, ScalarInt32, S_S, UInt8, Int32, UInt32);
300     scilab_fill_sub(ScalarUInt8, ScalarUInt32, S_S, UInt8, UInt32, UInt32);
301     scilab_fill_sub(ScalarUInt8, ScalarInt64, S_S, UInt8, Int64, UInt64);
302     scilab_fill_sub(ScalarUInt8, ScalarUInt64, S_S, UInt8, UInt64, UInt64);
303     scilab_fill_sub(ScalarUInt8, ScalarBool, S_S, UInt8, Bool, UInt8);
304
305     //Int16
306     //Matrix - Matrix
307     scilab_fill_sub(Int16, Double, M_M, Int16, Double, Int16);
308     scilab_fill_sub(Int16, Int8, M_M, Int16, Int8, Int16);
309     scilab_fill_sub(Int16, UInt8, M_M, Int16, UInt8, UInt16);
310     scilab_fill_sub(Int16, Int16, M_M, Int16, Int16, Int16);
311     scilab_fill_sub(Int16, UInt16, M_M, Int16, UInt16, UInt16);
312     scilab_fill_sub(Int16, Int32, M_M, Int16, Int32, Int32);
313     scilab_fill_sub(Int16, UInt32, M_M, Int16, UInt32, UInt32);
314     scilab_fill_sub(Int16, Int64, M_M, Int16, Int64, Int64);
315     scilab_fill_sub(Int16, UInt64, M_M, Int16, UInt64, UInt64);
316     scilab_fill_sub(Int16, Bool, M_M, Int16, Bool, Int16);
317     scilab_fill_sub(Int16, Empty, M_E, Int16, Double, Int16);
318
319     //Matrix - Scalar
320     scilab_fill_sub(Int16, ScalarDouble, M_S, Int16, Double, Int16);
321     scilab_fill_sub(Int16, ScalarInt8, M_S, Int16, Int8, Int16);
322     scilab_fill_sub(Int16, ScalarUInt8, M_S, Int16, UInt8, UInt16);
323     scilab_fill_sub(Int16, ScalarInt16, M_S, Int16, Int16, Int16);
324     scilab_fill_sub(Int16, ScalarUInt16, M_S, Int16, UInt16, UInt16);
325     scilab_fill_sub(Int16, ScalarInt32, M_S, Int16, Int32, Int32);
326     scilab_fill_sub(Int16, ScalarUInt32, M_S, Int16, UInt32, UInt32);
327     scilab_fill_sub(Int16, ScalarInt64, M_S, Int16, Int64, Int64);
328     scilab_fill_sub(Int16, ScalarUInt64, M_S, Int16, UInt64, UInt64);
329     scilab_fill_sub(Int16, ScalarBool, M_S, Int16, Bool, Int16);
330
331     //Scalar - Matrix
332     scilab_fill_sub(ScalarInt16, Double, S_M, Int16, Double, Int16);
333     scilab_fill_sub(ScalarInt16, Int8, S_M, Int16, Int8, Int16);
334     scilab_fill_sub(ScalarInt16, UInt8, S_M, Int16, UInt8, UInt16);
335     scilab_fill_sub(ScalarInt16, Int16, S_M, Int16, Int16, Int16);
336     scilab_fill_sub(ScalarInt16, UInt16, S_M, Int16, UInt16, UInt16);
337     scilab_fill_sub(ScalarInt16, Int32, S_M, Int16, Int32, Int32);
338     scilab_fill_sub(ScalarInt16, UInt32, S_M, Int16, UInt32, UInt32);
339     scilab_fill_sub(ScalarInt16, Int64, S_M, Int16, Int64, Int64);
340     scilab_fill_sub(ScalarInt16, UInt64, S_M, Int16, UInt64, UInt64);
341     scilab_fill_sub(ScalarInt16, Bool, S_M, Int16, Bool, Int16);
342     scilab_fill_sub(ScalarInt16, Empty, S_E, Int16, Double, Int16);
343
344     //Scalar - Scalar
345     scilab_fill_sub(ScalarInt16, ScalarDouble, S_S, Int16, Double, Int16);
346     scilab_fill_sub(ScalarInt16, ScalarInt8, S_S, Int16, Int8, Int16);
347     scilab_fill_sub(ScalarInt16, ScalarUInt8, S_S, Int16, UInt8, UInt16);
348     scilab_fill_sub(ScalarInt16, ScalarInt16, S_S, Int16, Int16, Int16);
349     scilab_fill_sub(ScalarInt16, ScalarUInt16, S_S, Int16, UInt16, UInt16);
350     scilab_fill_sub(ScalarInt16, ScalarInt32, S_S, Int16, Int32, Int32);
351     scilab_fill_sub(ScalarInt16, ScalarUInt32, S_S, Int16, UInt32, UInt32);
352     scilab_fill_sub(ScalarInt16, ScalarInt64, S_S, Int16, Int64, Int64);
353     scilab_fill_sub(ScalarInt16, ScalarUInt64, S_S, Int16, UInt64, UInt64);
354     scilab_fill_sub(ScalarInt16, ScalarBool, S_S, Int16, Bool, Int16);
355
356     //UInt16
357     //Matrix - Matrix
358     scilab_fill_sub(UInt16, Double, M_M, UInt16, Double, UInt16);
359     scilab_fill_sub(UInt16, Int8, M_M, UInt16, Int8, UInt16);
360     scilab_fill_sub(UInt16, UInt8, M_M, UInt16, UInt8, UInt16);
361     scilab_fill_sub(UInt16, Int16, M_M, UInt16, Int16, UInt16);
362     scilab_fill_sub(UInt16, UInt16, M_M, UInt16, UInt16, UInt16);
363     scilab_fill_sub(UInt16, Int32, M_M, UInt16, Int32, UInt32);
364     scilab_fill_sub(UInt16, UInt32, M_M, UInt16, UInt32, UInt32);
365     scilab_fill_sub(UInt16, Int64, M_M, UInt16, Int64, UInt64);
366     scilab_fill_sub(UInt16, UInt64, M_M, UInt16, UInt64, UInt64);
367     scilab_fill_sub(UInt16, Bool, M_M, UInt16, Bool, UInt16);
368     scilab_fill_sub(UInt16, Empty, M_E, UInt16, Double, UInt16);
369
370     //Matrix - Scalar
371     scilab_fill_sub(UInt16, ScalarDouble, M_S, UInt16, Double, UInt16);
372     scilab_fill_sub(UInt16, ScalarInt8, M_S, UInt16, Int8, UInt16);
373     scilab_fill_sub(UInt16, ScalarUInt8, M_S, UInt16, UInt8, UInt16);
374     scilab_fill_sub(UInt16, ScalarInt16, M_S, UInt16, Int16, UInt16);
375     scilab_fill_sub(UInt16, ScalarUInt16, M_S, UInt16, UInt16, UInt16);
376     scilab_fill_sub(UInt16, ScalarInt32, M_S, UInt16, Int32, UInt32);
377     scilab_fill_sub(UInt16, ScalarUInt32, M_S, UInt16, UInt32, UInt32);
378     scilab_fill_sub(UInt16, ScalarInt64, M_S, UInt16, Int64, UInt64);
379     scilab_fill_sub(UInt16, ScalarUInt64, M_S, UInt16, UInt64, UInt64);
380     scilab_fill_sub(UInt16, ScalarBool, M_S, UInt16, Bool, UInt16);
381
382     //Scalar - Matrix
383     scilab_fill_sub(ScalarUInt16, Double, S_M, UInt16, Double, UInt16);
384     scilab_fill_sub(ScalarUInt16, Int8, S_M, UInt16, Int8, UInt16);
385     scilab_fill_sub(ScalarUInt16, UInt8, S_M, UInt16, UInt8, UInt16);
386     scilab_fill_sub(ScalarUInt16, Int16, S_M, UInt16, Int16, UInt16);
387     scilab_fill_sub(ScalarUInt16, UInt16, S_M, UInt16, UInt16, UInt16);
388     scilab_fill_sub(ScalarUInt16, Int32, S_M, UInt16, Int32, UInt32);
389     scilab_fill_sub(ScalarUInt16, UInt32, S_M, UInt16, UInt32, UInt32);
390     scilab_fill_sub(ScalarUInt16, Int64, S_M, UInt16, Int64, UInt64);
391     scilab_fill_sub(ScalarUInt16, UInt64, S_M, UInt16, UInt64, UInt64);
392     scilab_fill_sub(ScalarUInt16, Bool, S_M, UInt16, Bool, UInt16);
393     scilab_fill_sub(ScalarUInt16, Empty, S_E, UInt16, Double, UInt16);
394
395     //Scalar - Scalar
396     scilab_fill_sub(ScalarUInt16, ScalarDouble, S_S, UInt16, Double, UInt16);
397     scilab_fill_sub(ScalarUInt16, ScalarInt8, S_S, UInt16, Int8, UInt16);
398     scilab_fill_sub(ScalarUInt16, ScalarUInt8, S_S, UInt16, UInt8, UInt16);
399     scilab_fill_sub(ScalarUInt16, ScalarInt16, S_S, UInt16, Int16, UInt16);
400     scilab_fill_sub(ScalarUInt16, ScalarUInt16, S_S, UInt16, UInt16, UInt16);
401     scilab_fill_sub(ScalarUInt16, ScalarInt32, S_S, UInt16, Int32, UInt32);
402     scilab_fill_sub(ScalarUInt16, ScalarUInt32, S_S, UInt16, UInt32, UInt32);
403     scilab_fill_sub(ScalarUInt16, ScalarInt64, S_S, UInt16, Int64, UInt64);
404     scilab_fill_sub(ScalarUInt16, ScalarUInt64, S_S, UInt16, UInt64, UInt64);
405     scilab_fill_sub(ScalarUInt16, ScalarBool, S_S, UInt16, Bool, UInt16);
406
407     //Int32
408     //Matrix - Matrix
409     scilab_fill_sub(Int32, Double, M_M, Int32, Double, Int32);
410     scilab_fill_sub(Int32, Int8, M_M, Int32, Int8, Int32);
411     scilab_fill_sub(Int32, UInt8, M_M, Int32, UInt8, UInt32);
412     scilab_fill_sub(Int32, Int16, M_M, Int32, Int16, Int32);
413     scilab_fill_sub(Int32, UInt16, M_M, Int32, UInt16, UInt32);
414     scilab_fill_sub(Int32, Int32, M_M, Int32, Int32, Int32);
415     scilab_fill_sub(Int32, UInt32, M_M, Int32, UInt32, UInt32);
416     scilab_fill_sub(Int32, Int64, M_M, Int32, Int64, Int64);
417     scilab_fill_sub(Int32, UInt64, M_M, Int32, UInt64, UInt64);
418     scilab_fill_sub(Int32, Bool, M_M, Int32, Bool, Int32);
419     scilab_fill_sub(Int32, Empty, M_E, Int32, Double, Int32);
420
421     //Matrix - Scalar
422     scilab_fill_sub(Int32, ScalarDouble, M_S, Int32, Double, Int32);
423     scilab_fill_sub(Int32, ScalarInt8, M_S, Int32, Int8, Int32);
424     scilab_fill_sub(Int32, ScalarUInt8, M_S, Int32, UInt8, UInt32);
425     scilab_fill_sub(Int32, ScalarInt16, M_S, Int32, Int16, Int32);
426     scilab_fill_sub(Int32, ScalarUInt16, M_S, Int32, UInt16, UInt32);
427     scilab_fill_sub(Int32, ScalarInt32, M_S, Int32, Int32, Int32);
428     scilab_fill_sub(Int32, ScalarUInt32, M_S, Int32, UInt32, UInt32);
429     scilab_fill_sub(Int32, ScalarInt64, M_S, Int32, Int64, Int64);
430     scilab_fill_sub(Int32, ScalarUInt64, M_S, Int32, UInt64, UInt64);
431     scilab_fill_sub(Int32, ScalarBool, M_S, Int32, Bool, Int32);
432
433     //Scalar - Matrix
434     scilab_fill_sub(ScalarInt32, Double, S_M, Int32, Double, Int32);
435     scilab_fill_sub(ScalarInt32, Int8, S_M, Int32, Int8, Int32);
436     scilab_fill_sub(ScalarInt32, UInt8, S_M, Int32, UInt8, UInt32);
437     scilab_fill_sub(ScalarInt32, Int16, S_M, Int32, Int16, Int32);
438     scilab_fill_sub(ScalarInt32, UInt16, S_M, Int32, UInt16, UInt32);
439     scilab_fill_sub(ScalarInt32, Int32, S_M, Int32, Int32, Int32);
440     scilab_fill_sub(ScalarInt32, UInt32, S_M, Int32, UInt32, UInt32);
441     scilab_fill_sub(ScalarInt32, Int64, S_M, Int32, Int64, Int64);
442     scilab_fill_sub(ScalarInt32, UInt64, S_M, Int32, UInt64, UInt64);
443     scilab_fill_sub(ScalarInt32, Bool, S_M, Int32, Bool, Int32);
444     scilab_fill_sub(ScalarInt32, Empty, S_E, Int32, Double, Int32);
445
446     //Scalar - Scalar
447     scilab_fill_sub(ScalarInt32, ScalarDouble, S_S, Int32, Double, Int32);
448     scilab_fill_sub(ScalarInt32, ScalarInt8, S_S, Int32, Int8, Int32);
449     scilab_fill_sub(ScalarInt32, ScalarUInt8, S_S, Int32, UInt8, UInt32);
450     scilab_fill_sub(ScalarInt32, ScalarInt16, S_S, Int32, Int16, Int32);
451     scilab_fill_sub(ScalarInt32, ScalarUInt16, S_S, Int32, UInt16, UInt32);
452     scilab_fill_sub(ScalarInt32, ScalarInt32, S_S, Int32, Int32, Int32);
453     scilab_fill_sub(ScalarInt32, ScalarUInt32, S_S, Int32, UInt32, UInt32);
454     scilab_fill_sub(ScalarInt32, ScalarInt64, S_S, Int32, Int64, Int64);
455     scilab_fill_sub(ScalarInt32, ScalarUInt64, S_S, Int32, UInt64, UInt64);
456     scilab_fill_sub(ScalarInt32, ScalarBool, S_S, Int32, Bool, Int32);
457
458     //UInt32
459     //Matrix - Matrix
460     scilab_fill_sub(UInt32, Double, M_M, UInt32, Double, UInt32);
461     scilab_fill_sub(UInt32, Int8, M_M, UInt32, Int8, UInt32);
462     scilab_fill_sub(UInt32, UInt8, M_M, UInt32, UInt8, UInt32);
463     scilab_fill_sub(UInt32, Int16, M_M, UInt32, Int16, UInt32);
464     scilab_fill_sub(UInt32, UInt16, M_M, UInt32, UInt16, UInt32);
465     scilab_fill_sub(UInt32, Int32, M_M, UInt32, Int32, UInt32);
466     scilab_fill_sub(UInt32, UInt32, M_M, UInt32, UInt32, UInt32);
467     scilab_fill_sub(UInt32, Int64, M_M, UInt32, Int64, UInt64);
468     scilab_fill_sub(UInt32, UInt64, M_M, UInt32, UInt64, UInt64);
469     scilab_fill_sub(UInt32, Bool, M_M, UInt32, Bool, UInt32);
470     scilab_fill_sub(UInt32, Empty, M_E, UInt32, Double, UInt32);
471
472     //Matrix - Scalar
473     scilab_fill_sub(UInt32, ScalarDouble, M_S, UInt32, Double, UInt32);
474     scilab_fill_sub(UInt32, ScalarInt8, M_S, UInt32, Int8, UInt32);
475     scilab_fill_sub(UInt32, ScalarUInt8, M_S, UInt32, UInt8, UInt32);
476     scilab_fill_sub(UInt32, ScalarInt16, M_S, UInt32, Int16, UInt32);
477     scilab_fill_sub(UInt32, ScalarUInt16, M_S, UInt32, UInt16, UInt32);
478     scilab_fill_sub(UInt32, ScalarInt32, M_S, UInt32, Int32, UInt32);
479     scilab_fill_sub(UInt32, ScalarUInt32, M_S, UInt32, UInt32, UInt32);
480     scilab_fill_sub(UInt32, ScalarInt64, M_S, UInt32, Int64, UInt64);
481     scilab_fill_sub(UInt32, ScalarUInt64, M_S, UInt32, UInt64, UInt64);
482     scilab_fill_sub(UInt32, ScalarBool, M_S, UInt32, Bool, UInt32);
483
484     //Scalar - Matrix
485     scilab_fill_sub(ScalarUInt32, Double, S_M, UInt32, Double, UInt32);
486     scilab_fill_sub(ScalarUInt32, Int8, S_M, UInt32, Int8, UInt32);
487     scilab_fill_sub(ScalarUInt32, UInt8, S_M, UInt32, UInt8, UInt32);
488     scilab_fill_sub(ScalarUInt32, Int16, S_M, UInt32, Int16, UInt32);
489     scilab_fill_sub(ScalarUInt32, UInt16, S_M, UInt32, UInt16, UInt32);
490     scilab_fill_sub(ScalarUInt32, Int32, S_M, UInt32, Int32, UInt32);
491     scilab_fill_sub(ScalarUInt32, UInt32, S_M, UInt32, UInt32, UInt32);
492     scilab_fill_sub(ScalarUInt32, Int64, S_M, UInt32, Int64, UInt64);
493     scilab_fill_sub(ScalarUInt32, UInt64, S_M, UInt32, UInt64, UInt64);
494     scilab_fill_sub(ScalarUInt32, Bool, S_M, UInt32, Bool, UInt32);
495     scilab_fill_sub(ScalarUInt32, Empty, S_E, UInt32, Double, UInt32);
496
497     //Scalar - Scalar
498     scilab_fill_sub(ScalarUInt32, ScalarDouble, S_S, UInt32, Double, UInt32);
499     scilab_fill_sub(ScalarUInt32, ScalarInt8, S_S, UInt32, Int8, UInt32);
500     scilab_fill_sub(ScalarUInt32, ScalarUInt8, S_S, UInt32, UInt8, UInt32);
501     scilab_fill_sub(ScalarUInt32, ScalarInt16, S_S, UInt32, Int16, UInt32);
502     scilab_fill_sub(ScalarUInt32, ScalarUInt16, S_S, UInt32, UInt16, UInt32);
503     scilab_fill_sub(ScalarUInt32, ScalarInt32, S_S, UInt32, Int32, UInt32);
504     scilab_fill_sub(ScalarUInt32, ScalarUInt32, S_S, UInt32, UInt32, UInt32);
505     scilab_fill_sub(ScalarUInt32, ScalarInt64, S_S, UInt32, Int64, UInt64);
506     scilab_fill_sub(ScalarUInt32, ScalarUInt64, S_S, UInt32, UInt64, UInt64);
507     scilab_fill_sub(ScalarUInt32, ScalarBool, S_S, UInt32, Bool, UInt32);
508
509     //Int64
510     //Matrix - Matrix
511     scilab_fill_sub(Int64, Double, M_M, Int64, Double, Int64);
512     scilab_fill_sub(Int64, Int8, M_M, Int64, Int8, Int64);
513     scilab_fill_sub(Int64, UInt8, M_M, Int64, UInt8, UInt64);
514     scilab_fill_sub(Int64, Int16, M_M, Int64, Int16, Int64);
515     scilab_fill_sub(Int64, UInt16, M_M, Int64, UInt16, UInt64);
516     scilab_fill_sub(Int64, Int32, M_M, Int64, Int32, Int64);
517     scilab_fill_sub(Int64, UInt32, M_M, Int64, UInt32, UInt64);
518     scilab_fill_sub(Int64, Int64, M_M, Int64, Int64, Int64);
519     scilab_fill_sub(Int64, UInt64, M_M, Int64, UInt64, UInt64);
520     scilab_fill_sub(Int64, Bool, M_M, Int64, Bool, Int64);
521     scilab_fill_sub(Int64, Empty, M_E, Int64, Double, Int64);
522
523     //Matrix - Scalar
524     scilab_fill_sub(Int64, ScalarDouble, M_S, Int64, Double, Int64);
525     scilab_fill_sub(Int64, ScalarInt8, M_S, Int64, Int8, Int64);
526     scilab_fill_sub(Int64, ScalarUInt8, M_S, Int64, UInt8, UInt64);
527     scilab_fill_sub(Int64, ScalarInt16, M_S, Int64, Int16, Int64);
528     scilab_fill_sub(Int64, ScalarUInt16, M_S, Int64, UInt16, UInt64);
529     scilab_fill_sub(Int64, ScalarInt32, M_S, Int64, Int32, Int64);
530     scilab_fill_sub(Int64, ScalarUInt32, M_S, Int64, UInt32, UInt64);
531     scilab_fill_sub(Int64, ScalarInt64, M_S, Int64, Int64, Int64);
532     scilab_fill_sub(Int64, ScalarUInt64, M_S, Int64, UInt64, UInt64);
533     scilab_fill_sub(Int64, ScalarBool, M_S, Int64, Bool, Int64);
534
535     //Scalar - Matrix
536     scilab_fill_sub(ScalarInt64, Double, S_M, Int64, Double, Int64);
537     scilab_fill_sub(ScalarInt64, Int8, S_M, Int64, Int8, Int64);
538     scilab_fill_sub(ScalarInt64, UInt8, S_M, Int64, UInt8, UInt64);
539     scilab_fill_sub(ScalarInt64, Int16, S_M, Int64, Int16, Int64);
540     scilab_fill_sub(ScalarInt64, UInt16, S_M, Int64, UInt16, UInt64);
541     scilab_fill_sub(ScalarInt64, Int32, S_M, Int64, Int32, Int64);
542     scilab_fill_sub(ScalarInt64, UInt32, S_M, Int64, UInt32, UInt64);
543     scilab_fill_sub(ScalarInt64, Int64, S_M, Int64, Int64, Int64);
544     scilab_fill_sub(ScalarInt64, UInt64, S_M, Int64, UInt64, UInt64);
545     scilab_fill_sub(ScalarInt64, Bool, S_M, Int64, Bool, Int64);
546     scilab_fill_sub(ScalarInt64, Empty, S_E, Int64, Double, Int64);
547
548     //Scalar - Scalar
549     scilab_fill_sub(ScalarInt64, ScalarDouble, S_S, Int64, Double, Int64);
550     scilab_fill_sub(ScalarInt64, ScalarInt8, S_S, Int64, Int8, Int64);
551     scilab_fill_sub(ScalarInt64, ScalarUInt8, S_S, Int64, UInt8, UInt64);
552     scilab_fill_sub(ScalarInt64, ScalarInt16, S_S, Int64, Int16, Int64);
553     scilab_fill_sub(ScalarInt64, ScalarUInt16, S_S, Int64, UInt16, UInt64);
554     scilab_fill_sub(ScalarInt64, ScalarInt32, S_S, Int64, Int32, Int64);
555     scilab_fill_sub(ScalarInt64, ScalarUInt32, S_S, Int64, UInt32, UInt64);
556     scilab_fill_sub(ScalarInt64, ScalarInt64, S_S, Int64, Int64, Int64);
557     scilab_fill_sub(ScalarInt64, ScalarUInt64, S_S, Int64, UInt64, UInt64);
558     scilab_fill_sub(ScalarInt64, ScalarBool, S_S, Int64, Bool, Int64);
559
560     //UInt64
561     //Matrix - Matrix
562     scilab_fill_sub(UInt64, Double, M_M, UInt64, Double, UInt64);
563     scilab_fill_sub(UInt64, Int8, M_M, UInt64, Int8, UInt64);
564     scilab_fill_sub(UInt64, UInt8, M_M, UInt64, UInt8, UInt64);
565     scilab_fill_sub(UInt64, Int16, M_M, UInt64, Int16, UInt64);
566     scilab_fill_sub(UInt64, UInt16, M_M, UInt64, UInt16, UInt64);
567     scilab_fill_sub(UInt64, Int32, M_M, UInt64, Int32, UInt64);
568     scilab_fill_sub(UInt64, UInt32, M_M, UInt64, UInt32, UInt64);
569     scilab_fill_sub(UInt64, Int64, M_M, UInt64, Int64, UInt64);
570     scilab_fill_sub(UInt64, UInt64, M_M, UInt64, UInt64, UInt64);
571     scilab_fill_sub(UInt64, Bool, M_M, UInt64, Bool, UInt64);
572     scilab_fill_sub(UInt64, Empty, M_E, UInt64, Double, UInt64);
573
574     //Matrix - Scalar
575     scilab_fill_sub(UInt64, ScalarDouble, M_S, UInt64, Double, UInt64);
576     scilab_fill_sub(UInt64, ScalarInt8, M_S, UInt64, Int8, UInt64);
577     scilab_fill_sub(UInt64, ScalarUInt8, M_S, UInt64, UInt8, UInt64);
578     scilab_fill_sub(UInt64, ScalarInt16, M_S, UInt64, Int16, UInt64);
579     scilab_fill_sub(UInt64, ScalarUInt16, M_S, UInt64, UInt16, UInt64);
580     scilab_fill_sub(UInt64, ScalarInt32, M_S, UInt64, Int32, UInt64);
581     scilab_fill_sub(UInt64, ScalarUInt32, M_S, UInt64, UInt32, UInt64);
582     scilab_fill_sub(UInt64, ScalarInt64, M_S, UInt64, Int64, UInt64);
583     scilab_fill_sub(UInt64, ScalarUInt64, M_S, UInt64, UInt64, UInt64);
584     scilab_fill_sub(UInt64, ScalarBool, M_S, UInt64, Bool, UInt64);
585
586     //Scalar - Matrix
587     scilab_fill_sub(ScalarUInt64, Double, S_M, UInt64, Double, UInt64);
588     scilab_fill_sub(ScalarUInt64, Int8, S_M, UInt64, Int8, UInt64);
589     scilab_fill_sub(ScalarUInt64, UInt8, S_M, UInt64, UInt8, UInt64);
590     scilab_fill_sub(ScalarUInt64, Int16, S_M, UInt64, Int16, UInt64);
591     scilab_fill_sub(ScalarUInt64, UInt16, S_M, UInt64, UInt16, UInt64);
592     scilab_fill_sub(ScalarUInt64, Int32, S_M, UInt64, Int32, UInt64);
593     scilab_fill_sub(ScalarUInt64, UInt32, S_M, UInt64, UInt32, UInt64);
594     scilab_fill_sub(ScalarUInt64, Int64, S_M, UInt64, Int64, UInt64);
595     scilab_fill_sub(ScalarUInt64, UInt64, S_M, UInt64, UInt64, UInt64);
596     scilab_fill_sub(ScalarUInt64, Bool, S_M, UInt64, Bool, UInt64);
597     scilab_fill_sub(ScalarUInt64, Empty, S_E, UInt64, Double, UInt64);
598
599     //Scalar - Scalar
600     scilab_fill_sub(ScalarUInt64, ScalarDouble, S_S, UInt64, Double, UInt64);
601     scilab_fill_sub(ScalarUInt64, ScalarInt8, S_S, UInt64, Int8, UInt64);
602     scilab_fill_sub(ScalarUInt64, ScalarUInt8, S_S, UInt64, UInt8, UInt64);
603     scilab_fill_sub(ScalarUInt64, ScalarInt16, S_S, UInt64, Int16, UInt64);
604     scilab_fill_sub(ScalarUInt64, ScalarUInt16, S_S, UInt64, UInt16, UInt64);
605     scilab_fill_sub(ScalarUInt64, ScalarInt32, S_S, UInt64, Int32, UInt64);
606     scilab_fill_sub(ScalarUInt64, ScalarUInt32, S_S, UInt64, UInt32, UInt64);
607     scilab_fill_sub(ScalarUInt64, ScalarInt64, S_S, UInt64, Int64, UInt64);
608     scilab_fill_sub(ScalarUInt64, ScalarUInt64, S_S, UInt64, UInt64, UInt64);
609     scilab_fill_sub(ScalarUInt64, ScalarBool, S_S, UInt64, Bool, UInt64);
610
611     //Bool
612     //Matrix - Matrix
613     scilab_fill_sub(Bool, Double, M_M, Bool, Double, Double);
614     scilab_fill_sub(Bool, Int8, M_M, Bool, Int8, Int8);
615     scilab_fill_sub(Bool, UInt8, M_M, Bool, UInt8, UInt8);
616     scilab_fill_sub(Bool, Int16, M_M, Bool, Int16, Int16);
617     scilab_fill_sub(Bool, UInt16, M_M, Bool, UInt16, UInt16);
618     scilab_fill_sub(Bool, Int32, M_M, Bool, Int32, Int32);
619     scilab_fill_sub(Bool, UInt32, M_M, Bool, UInt32, UInt32);
620     scilab_fill_sub(Bool, Int64, M_M, Bool, Int64, Int64);
621     scilab_fill_sub(Bool, UInt64, M_M, Bool, UInt64, UInt64);
622     scilab_fill_sub(Bool, Bool, M_M, Bool, Bool, Double);
623     scilab_fill_sub(Bool, Empty, M_E, Bool, Double, Double);
624
625     //Matrix - Scalar
626     scilab_fill_sub(Bool, ScalarDouble, M_S, Bool, Double, Double);
627     scilab_fill_sub(Bool, ScalarInt8, M_S, Bool, Int8, Int8);
628     scilab_fill_sub(Bool, ScalarUInt8, M_S, Bool, UInt8, UInt8);
629     scilab_fill_sub(Bool, ScalarInt16, M_S, Bool, Int16, Int16);
630     scilab_fill_sub(Bool, ScalarUInt16, M_S, Bool, UInt16, UInt16);
631     scilab_fill_sub(Bool, ScalarInt32, M_S, Bool, Int32, Int32);
632     scilab_fill_sub(Bool, ScalarUInt32, M_S, Bool, UInt32, UInt32);
633     scilab_fill_sub(Bool, ScalarInt64, M_S, Bool, Int64, Int64);
634     scilab_fill_sub(Bool, ScalarUInt64, M_S, Bool, UInt64, UInt64);
635     scilab_fill_sub(Bool, ScalarBool, M_S, Bool, Bool, Double);
636
637     //Scalar - Matrix
638     scilab_fill_sub(ScalarBool, Double, S_M, Bool, Double, Double);
639     scilab_fill_sub(ScalarBool, Int8, S_M, Bool, Int8, Int8);
640     scilab_fill_sub(ScalarBool, UInt8, S_M, Bool, UInt8, UInt8);
641     scilab_fill_sub(ScalarBool, Int16, S_M, Bool, Int16, Int16);
642     scilab_fill_sub(ScalarBool, UInt16, S_M, Bool, UInt16, UInt16);
643     scilab_fill_sub(ScalarBool, Int32, S_M, Bool, Int32, Int32);
644     scilab_fill_sub(ScalarBool, UInt32, S_M, Bool, UInt32, UInt32);
645     scilab_fill_sub(ScalarBool, Int64, S_M, Bool, Int64, Int64);
646     scilab_fill_sub(ScalarBool, UInt64, S_M, Bool, UInt64, UInt64);
647     scilab_fill_sub(ScalarBool, Bool, S_M, Bool, Bool, Double);
648     scilab_fill_sub(ScalarBool, Empty, S_E, Bool, Double, Double);
649
650     //Scalar - Scalar
651     scilab_fill_sub(ScalarBool, ScalarDouble, S_S, Bool, Double, Double);
652     scilab_fill_sub(ScalarBool, ScalarInt8, S_S, Bool, Int8, Int8);
653     scilab_fill_sub(ScalarBool, ScalarUInt8, S_S, Bool, UInt8, UInt8);
654     scilab_fill_sub(ScalarBool, ScalarInt16, S_S, Bool, Int16, Int16);
655     scilab_fill_sub(ScalarBool, ScalarUInt16, S_S, Bool, UInt16, UInt16);
656     scilab_fill_sub(ScalarBool, ScalarInt32, S_S, Bool, Int32, Int32);
657     scilab_fill_sub(ScalarBool, ScalarUInt32, S_S, Bool, UInt32, UInt32);
658     scilab_fill_sub(ScalarBool, ScalarInt64, S_S, Bool, Int64, Int64);
659     scilab_fill_sub(ScalarBool, ScalarUInt64, S_S, Bool, UInt64, UInt64);
660     scilab_fill_sub(ScalarBool, ScalarBool, S_S, Bool, Bool, Double);
661
662     //Identity
663     scilab_fill_sub(Identity, Double, I_M, Double, Double, Double);
664     scilab_fill_sub(Identity, DoubleComplex, I_MC, Double, Double, Double);
665     scilab_fill_sub(Identity, ScalarDouble, I_S, Double, Double, Double);
666     scilab_fill_sub(Identity, ScalarDoubleComplex, I_SC, Double, Double, Double);
667     scilab_fill_sub(Identity, Identity, I_I, Double, Double, Double);
668     scilab_fill_sub(Identity, IdentityComplex, I_IC, Double, Double, Double);
669     scilab_fill_sub(Identity, Empty, I_E, Double, Double, Double);
670
671     scilab_fill_sub(Identity, Polynom, I_M, Double, Polynom, Polynom);
672     scilab_fill_sub(Identity, PolynomComplex, I_MC, Double, Polynom, Polynom);
673     scilab_fill_sub(Identity, ScalarPolynom, I_M, Double, Polynom, Polynom);
674     scilab_fill_sub(Identity, ScalarPolynomComplex, I_MC, Double, Polynom, Polynom);
675     scilab_fill_sub(Identity, Sparse, M_M, Double, Sparse, Sparse);
676     scilab_fill_sub(Identity, SparseComplex, M_M, Double, Sparse, Sparse);
677
678     scilab_fill_sub(IdentityComplex, Double, IC_M, Double, Double, Double);
679     scilab_fill_sub(IdentityComplex, DoubleComplex, IC_MC, Double, Double, Double);
680     scilab_fill_sub(IdentityComplex, ScalarDouble, IC_S, Double, Double, Double);
681     scilab_fill_sub(IdentityComplex, ScalarDoubleComplex, IC_SC, Double, Double, Double);
682     scilab_fill_sub(IdentityComplex, Identity, IC_I, Double, Double, Double);
683     scilab_fill_sub(IdentityComplex, IdentityComplex, IC_IC, Double, Double, Double);
684     scilab_fill_sub(IdentityComplex, Empty, IC_E, Double, Double, Double);
685
686     scilab_fill_sub(IdentityComplex, Polynom, IC_M, Double, Polynom, Polynom);
687     scilab_fill_sub(IdentityComplex, PolynomComplex, IC_MC, Double, Polynom, Polynom);
688     scilab_fill_sub(IdentityComplex, ScalarPolynom, IC_M, Double, Polynom, Polynom);
689     scilab_fill_sub(IdentityComplex, ScalarPolynomComplex, IC_MC, Double, Polynom, Polynom);
690     scilab_fill_sub(IdentityComplex, Sparse, M_M, Double, Sparse, Sparse);
691     scilab_fill_sub(IdentityComplex, SparseComplex, M_M, Double, Sparse, Sparse);
692
693     scilab_fill_sub(Identity, ScalarInt8, I_S, Double, Int8, Int8);
694     scilab_fill_sub(Identity, ScalarUInt8, I_S, Double, UInt8, UInt8);
695     scilab_fill_sub(Identity, ScalarInt16, I_S, Double, Int16, Int16);
696     scilab_fill_sub(Identity, ScalarUInt16, I_S, Double, UInt16, UInt16);
697     scilab_fill_sub(Identity, ScalarInt32, I_S, Double, Int32, Int32);
698     scilab_fill_sub(Identity, ScalarUInt32, I_S, Double, UInt32, UInt32);
699     scilab_fill_sub(Identity, ScalarInt64, I_S, Double, Int64, Int64);
700     scilab_fill_sub(Identity, ScalarUInt64, I_S, Double, UInt64, UInt64);
701
702     scilab_fill_sub(Identity, Int8, I_M, Double, Int8, Int8);
703     scilab_fill_sub(Identity, UInt8, I_M, Double, UInt8, UInt8);
704     scilab_fill_sub(Identity, Int16, I_M, Double, Int16, Int16);
705     scilab_fill_sub(Identity, UInt16, I_M, Double, UInt16, UInt16);
706     scilab_fill_sub(Identity, Int32, I_M, Double, Int32, Int32);
707     scilab_fill_sub(Identity, UInt32, I_M, Double, UInt32, UInt32);
708     scilab_fill_sub(Identity, Int64, I_M, Double, Int64, Int64);
709     scilab_fill_sub(Identity, UInt64, I_M, Double, UInt64, UInt64);
710
711     //Polynom
712
713     //poly - poly
714     scilab_fill_sub(Polynom, Polynom, M_M, Polynom, Polynom, Polynom);
715     scilab_fill_sub(Polynom, PolynomComplex, M_M, Polynom, Polynom, Polynom);
716     scilab_fill_sub(PolynomComplex, Polynom, M_M, Polynom, Polynom, Polynom);
717     scilab_fill_sub(PolynomComplex, PolynomComplex, M_M, Polynom, Polynom, Polynom);
718
719     //poly - scalar poly
720     scilab_fill_sub(Polynom, ScalarPolynom, M_M, Polynom, Polynom, Polynom);
721     scilab_fill_sub(Polynom, ScalarPolynomComplex, M_M, Polynom, Polynom, Polynom);
722     scilab_fill_sub(PolynomComplex, ScalarPolynom, M_M, Polynom, Polynom, Polynom);
723     scilab_fill_sub(PolynomComplex, ScalarPolynomComplex, M_M, Polynom, Polynom, Polynom);
724
725     //poly - double
726     scilab_fill_sub(Polynom, Double, M_M, Polynom, Double, Polynom);
727     scilab_fill_sub(Polynom, DoubleComplex, M_M, Polynom, Double, Polynom);
728     scilab_fill_sub(PolynomComplex, Double, M_M, Polynom, Double, Polynom);
729     scilab_fill_sub(PolynomComplex, DoubleComplex, M_M, Polynom, Double, Polynom);
730
731     //poly - scalar double
732     scilab_fill_sub(Polynom, ScalarDouble, M_M, Polynom, Double, Polynom);
733     scilab_fill_sub(Polynom, ScalarDoubleComplex, M_M, Polynom, Double, Polynom);
734     scilab_fill_sub(PolynomComplex, ScalarDouble, M_M, Polynom, Double, Polynom);
735     scilab_fill_sub(PolynomComplex, ScalarDoubleComplex, M_M, Polynom, Double, Polynom);
736
737     //poly - []
738     scilab_fill_sub(Polynom, Empty, M_E, Polynom, Double, Polynom);
739     scilab_fill_sub(PolynomComplex, Empty, M_E, Polynom, Double, Polynom);
740
741     //poly - eye
742     scilab_fill_sub(Polynom, Identity, M_M, Polynom, Double, Polynom);
743     scilab_fill_sub(Polynom, IdentityComplex, M_M, Polynom, Double, Polynom);
744     scilab_fill_sub(PolynomComplex, Identity, M_M, Polynom, Double, Polynom);
745     scilab_fill_sub(PolynomComplex, IdentityComplex, M_M, Polynom, Double, Polynom);
746
747     //scalar poly - poly
748     scilab_fill_sub(ScalarPolynom, Polynom, M_M, Polynom, Polynom, Polynom);
749     scilab_fill_sub(ScalarPolynom, PolynomComplex, M_M, Polynom, Polynom, Polynom);
750     scilab_fill_sub(ScalarPolynomComplex, Polynom, M_M, Polynom, Polynom, Polynom);
751     scilab_fill_sub(ScalarPolynomComplex, PolynomComplex, M_M, Polynom, Polynom, Polynom);
752
753     //scalar poly - scalar poly
754     scilab_fill_sub(ScalarPolynom, ScalarPolynom, M_M, Polynom, Polynom, Polynom);
755     scilab_fill_sub(ScalarPolynom, ScalarPolynomComplex, M_M, Polynom, Polynom, Polynom);
756     scilab_fill_sub(ScalarPolynomComplex, ScalarPolynom, M_M, Polynom, Polynom, Polynom);
757     scilab_fill_sub(ScalarPolynomComplex, ScalarPolynomComplex, M_M, Polynom, Polynom, Polynom);
758
759     //scalar poly - double
760     scilab_fill_sub(ScalarPolynom, Double, M_M, Polynom, Double, Polynom);
761     scilab_fill_sub(ScalarPolynom, DoubleComplex, M_M, Polynom, Double, Polynom);
762     scilab_fill_sub(ScalarPolynomComplex, Double, M_M, Polynom, Double, Polynom);
763     scilab_fill_sub(ScalarPolynomComplex, DoubleComplex, M_M, Polynom, Double, Polynom);
764
765     //scalar poly - scalar double
766     scilab_fill_sub(ScalarPolynom, ScalarDouble, M_M, Polynom, Double, Polynom);
767     scilab_fill_sub(ScalarPolynom, ScalarDoubleComplex, M_M, Polynom, Double, Polynom);
768     scilab_fill_sub(ScalarPolynomComplex, ScalarDouble, M_M, Polynom, Double, Polynom);
769     scilab_fill_sub(ScalarPolynomComplex, ScalarDoubleComplex, M_M, Polynom, Double, Polynom);
770
771     //scalar poly - []
772     scilab_fill_sub(ScalarPolynom, Empty, M_E, Polynom, Double, Polynom);
773     scilab_fill_sub(ScalarPolynomComplex, Empty, M_E, Polynom, Double, Polynom);
774
775     //scalar poly - eye
776     scilab_fill_sub(ScalarPolynom, Identity, M_M, Polynom, Double, Polynom);
777     scilab_fill_sub(ScalarPolynom, IdentityComplex, M_M, Polynom, Double, Polynom);
778     scilab_fill_sub(ScalarPolynomComplex, Identity, M_M, Polynom, Double, Polynom);
779     scilab_fill_sub(ScalarPolynomComplex, IdentityComplex, M_M, Polynom, Double, Polynom);
780
781     //Sparse
782     scilab_fill_sub(Sparse, Sparse, M_M, Sparse, Sparse, Sparse);
783     scilab_fill_sub(Sparse, SparseComplex, M_M, Sparse, Sparse, Sparse);
784     scilab_fill_sub(Sparse, Double, M_M, Sparse, Double, Double);
785     scilab_fill_sub(Sparse, DoubleComplex, M_M, Sparse, Double, Double);
786     scilab_fill_sub(Sparse, ScalarDouble, M_M, Sparse, Double, Double);
787     scilab_fill_sub(Sparse, ScalarDoubleComplex, M_M, Sparse, Double, Double);
788
789     scilab_fill_sub(Sparse, Empty, M_E, Sparse, Double, Sparse);
790     scilab_fill_sub(Sparse, Identity, M_M, Sparse, Double, Sparse);
791     scilab_fill_sub(Sparse, IdentityComplex, M_M, Sparse, Double, Sparse);
792
793     scilab_fill_sub(SparseComplex, Sparse, M_M, Sparse, Sparse, Sparse);
794     scilab_fill_sub(SparseComplex, SparseComplex, M_M, Sparse, Sparse, Sparse);
795     scilab_fill_sub(SparseComplex, Double, M_M, Sparse, Double, Double);
796     scilab_fill_sub(SparseComplex, DoubleComplex, M_M, Sparse, Double, Double);
797     scilab_fill_sub(SparseComplex, ScalarDouble, M_M, Sparse, Double, Double);
798     scilab_fill_sub(SparseComplex, ScalarDoubleComplex, M_M, Sparse, Double, Double);
799
800     scilab_fill_sub(SparseComplex, Empty, M_E, Sparse, Double, Sparse);
801     scilab_fill_sub(SparseComplex, Identity, M_M, Sparse, Double, Sparse);
802     scilab_fill_sub(SparseComplex, IdentityComplex, M_M, Sparse, Double, Sparse);
803
804 #undef scilab_fill_sub
805 }
806
807 InternalType* GenericMinus(InternalType* _pLeftOperand, InternalType* _pRightOperand)
808 {
809     InternalType *pResult = NULL;
810
811     sub_function sub = pSubfunction[_pLeftOperand->getId()][_pRightOperand->getId()];
812     if (sub)
813     {
814         pResult = sub(_pLeftOperand, _pRightOperand);
815         if (pResult)
816         {
817             return pResult;
818         }
819     }
820
821     /*
822     ** Default case : Return NULL will Call Overloading.
823     */
824     return NULL;
825 }
826
827 template<class T, class U, class O>
828 InternalType* sub_M_M(T *_pL, U *_pR)
829 {
830     int iDimsL = _pL->getDims();
831     int iDimsR = _pR->getDims();
832
833     if (iDimsL != iDimsR)
834     {
835         return nullptr;
836     }
837
838     int* piDimsL = _pL->getDimsArray();
839     int* piDimsR = _pR->getDimsArray();
840
841     for (int i = 0 ; i < iDimsL ; ++i)
842     {
843         if (piDimsL[i] != piDimsR[i])
844         {
845             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
846         }
847     }
848
849     O* pOut = new O(iDimsL, piDimsL);
850
851     sub(_pL->get(), (size_t)_pL->getSize(), _pR->get(), pOut->get());
852     return pOut;
853 }
854
855 template<class T, class U, class O>
856 InternalType* sub_M_MC(T *_pL, U *_pR)
857 {
858     int iDimsL = _pL->getDims();
859     int iDimsR = _pR->getDims();
860
861     if (iDimsL != iDimsR)
862     {
863         return nullptr;
864     }
865
866     int* piDimsL = _pL->getDimsArray();
867     int* piDimsR = _pR->getDimsArray();
868
869     for (int i = 0 ; i < iDimsL ; ++i)
870     {
871         if (piDimsL[i] != piDimsR[i])
872         {
873             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
874         }
875     }
876
877     O* pOut = new O(iDimsL, piDimsL, true);
878
879     sub(_pL->get(), (size_t)_pL->getSize(), _pR->get(), _pR->getImg(), pOut->get(), pOut->getImg());
880     return pOut;
881 }
882
883 template<class T, class U, class O>
884 InternalType* sub_M_S(T *_pL, U *_pR)
885 {
886     O* pOut = new O(_pL->getDims(), _pL->getDimsArray());
887     sub(_pL->get(), (size_t)_pL->getSize(), _pR->get(0), pOut->get());
888     return pOut;
889 }
890
891 template<class T, class U, class O>
892 InternalType* sub_M_SC(T *_pL, U *_pR)
893 {
894     O* pOut = new O(_pL->getDims(), _pL->getDimsArray(), true);
895     sub(_pL->get(), (size_t)_pL->getSize(), _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
896     return pOut;
897 }
898
899 template<class T, class U, class O>
900 InternalType* sub_M_E(T *_pL, U * /*_pR*/)
901 {
902     if (ConfigVariable::getOldEmptyBehaviour())
903     {
904         Sciwarning(_("operation -: Warning adding a matrix with the empty matrix old behaviour.\n"));
905         return _pL;
906     }
907
908     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n"));
909     Double* pOut = Double::Empty();
910     sub();
911     return pOut;
912 }
913
914
915 template<class T, class U, class O>
916 InternalType* sub_MC_M(T *_pL, U *_pR)
917 {
918     int iDimsL = _pL->getDims();
919     int iDimsR = _pR->getDims();
920
921     if (iDimsL != iDimsR)
922     {
923         return nullptr;
924     }
925
926     int* piDimsL = _pL->getDimsArray();
927     int* piDimsR = _pR->getDimsArray();
928
929     for (int i = 0 ; i < iDimsL ; ++i)
930     {
931         if (piDimsL[i] != piDimsR[i])
932         {
933             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
934         }
935     }
936
937     O* pOut = new O(iDimsL, piDimsL, true);
938
939     sub(_pL->get(), _pL->getImg(), (size_t)_pL->getSize(), _pR->get(), pOut->get(), pOut->getImg());
940     return pOut;
941 }
942
943 template<class T, class U, class O>
944 InternalType* sub_MC_MC(T *_pL, U *_pR)
945 {
946     int iDimsL = _pL->getDims();
947     int iDimsR = _pR->getDims();
948
949     if (iDimsL != iDimsR)
950     {
951         return nullptr;
952     }
953
954     int* piDimsL = _pL->getDimsArray();
955     int* piDimsR = _pR->getDimsArray();
956
957     for (int i = 0 ; i < iDimsL ; ++i)
958     {
959         if (piDimsL[i] != piDimsR[i])
960         {
961             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
962         }
963     }
964
965     O* pOut = new O(iDimsL, piDimsL, true);
966
967     sub(_pL->get(), _pL->getImg(), (size_t)_pL->getSize(), _pR->get(), _pR->getImg(), pOut->get(), pOut->getImg());
968     return pOut;
969 }
970
971 template<class T, class U, class O>
972 InternalType* sub_MC_S(T *_pL, U *_pR)
973 {
974     O* pOut = new O(_pL->getDims(), _pL->getDimsArray(), true);
975     sub(_pL->get(), _pL->getImg(), (size_t)_pL->getSize(), _pR->get(0), pOut->get(), pOut->getImg());
976     return pOut;
977 }
978
979 template<class T, class U, class O>
980 InternalType* sub_MC_SC(T *_pL, U *_pR)
981 {
982     O* pOut = new O(_pL->getDims(), _pL->getDimsArray(), true);
983     sub(_pL->get(), _pL->getImg(), (size_t)_pL->getSize(), _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
984     return pOut;
985 }
986
987 template<class T, class U, class O>
988 InternalType* sub_MC_E(T *_pL, U * /*_pR*/)
989 {
990     if (ConfigVariable::getOldEmptyBehaviour())
991     {
992         Sciwarning(_("operation -: Warning adding a matrix with the empty matrix old behaviour.\n"));
993         return _pL;
994     }
995
996     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n"));
997     Double* pOut = Double::Empty();
998     sub();
999     return pOut;
1000 }
1001
1002
1003 template<class T, class U, class O>
1004 InternalType* sub_S_M(T *_pL, U *_pR)
1005 {
1006     O* pOut = new O(_pR->getDims(), _pR->getDimsArray());
1007     sub(_pL->get(0), (size_t)_pR->getSize(), _pR->get(), pOut->get());
1008     return pOut;
1009 }
1010
1011 template<class T, class U, class O>
1012 InternalType* sub_S_MC(T *_pL, U *_pR)
1013 {
1014     O* pOut = new O(_pR->getDims(), _pR->getDimsArray(), true);
1015     sub(_pL->get(0), (size_t)_pR->getSize(), _pR->get(), _pR->getImg(), pOut->get(), pOut->getImg());
1016     return pOut;
1017 }
1018
1019 template<class T, class U, class O>
1020 InternalType* sub_S_S(T *_pL, U *_pR)
1021 {
1022     O* pOut = new O(0);
1023     sub(_pL->get(0), _pR->get(0), pOut->get());
1024     return pOut;
1025 }
1026
1027 template<class T, class U, class O>
1028 InternalType* sub_S_SC(T *_pL, U *_pR)
1029 {
1030     O* pOut = new O(0.0, 0.0);
1031     sub(_pL->get(), (size_t)1, _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
1032     return pOut;
1033 }
1034
1035 template<class T, class U, class O>
1036 InternalType* sub_S_E(T *_pL, U * /*_pR*/)
1037 {
1038     if (ConfigVariable::getOldEmptyBehaviour())
1039     {
1040         Sciwarning(_("operation -: Warning adding a matrix with the empty matrix old behaviour.\n"));
1041         return _pL;
1042     }
1043
1044     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n"));
1045     Double* pOut = Double::Empty();
1046     sub();
1047     return pOut;
1048 }
1049
1050
1051 template<class T, class U, class O>
1052 InternalType* sub_SC_M(T *_pL, U *_pR)
1053 {
1054     O* pOut = new O(_pR->getDims(), _pR->getDimsArray(), true);
1055     sub(_pL->get(0), _pL->getImg(0), (size_t)_pR->getSize(), _pR->get(), pOut->get(), pOut->getImg());
1056     return pOut;
1057 }
1058
1059 template<class T, class U, class O>
1060 InternalType* sub_SC_MC(T *_pL, U *_pR)
1061 {
1062     O* pOut = new O(_pR->getDims(), _pR->getDimsArray(), true);
1063     sub(_pL->get(0), _pL->getImg(0), (size_t)_pR->getSize(), _pR->get(), _pR->getImg(), pOut->get(), pOut->getImg());
1064     return pOut;
1065 }
1066
1067 template<class T, class U, class O>
1068 InternalType* sub_SC_S(T *_pL, U *_pR)
1069 {
1070     O* pOut = new O(0.0, 0.0);
1071     sub(_pL->get(), _pL->getImg(), (size_t)1, _pR->get(0), pOut->get(), pOut->getImg());
1072     return pOut;
1073 }
1074
1075 template<class T, class U, class O>
1076 InternalType* sub_SC_SC(T *_pL, U *_pR)
1077 {
1078     O* pOut = new O(0.0, 0.0);
1079     sub(_pL->get(0), _pL->getImg(0), _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
1080     return pOut;
1081 }
1082
1083 template<class T, class U, class O>
1084 InternalType* sub_SC_E(T *_pL, U * /*_pR*/)
1085 {
1086     if (ConfigVariable::getOldEmptyBehaviour())
1087     {
1088         Sciwarning(_("operation -: Warning adding a matrix with the empty matrix old behaviour.\n"));
1089         return _pL;
1090     }
1091
1092     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n"));
1093     Double* pOut = Double::Empty();
1094     sub();
1095     return pOut;
1096 }
1097
1098
1099 template<class T, class U, class O>
1100 InternalType* sub_E_M(T * /*_pL*/, U *_pR)
1101 {
1102     if (ConfigVariable::getOldEmptyBehaviour())
1103     {
1104         Sciwarning(_("operation -: Warning adding a matrix with the empty matrix old behaviour.\n"));
1105         return opposite_M<U, O>(_pR);
1106     }
1107
1108     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n"));
1109     Double* pOut = Double::Empty();
1110     sub();
1111     return pOut;
1112 }
1113
1114 template<class T, class U, class O>
1115 InternalType* sub_E_MC(T * /*_pL*/, U *_pR)
1116 {
1117     if (ConfigVariable::getOldEmptyBehaviour())
1118     {
1119         Sciwarning(_("operation -: Warning adding a matrix with the empty matrix old behaviour.\n"));
1120         return opposite_MC<U, O>(_pR);
1121     }
1122
1123     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n"));
1124     Double* pOut = Double::Empty();
1125     sub();
1126     return pOut;
1127 }
1128
1129 template<class T, class U, class O>
1130 InternalType* sub_E_E(T * /*_pL*/, U * /*_pR*/)
1131 {
1132     Double* pOut = Double::Empty();
1133     sub();
1134     return pOut;
1135 }
1136
1137 template<class T, class U, class O>
1138 InternalType* sub_I_M(T *_pL, U *_pR)
1139 {
1140     int iDims = _pR->getDims();
1141     int* piDims = _pR->getDimsArray();
1142     O* pOut = (O*)opposite_M<U, O>(_pR);
1143     double dblLeft = _pL->get(0);
1144     int iLeadDims = piDims[0];
1145     int* piIndex = new int[iDims];
1146     piIndex[0] = 0;
1147
1148     //find smaller dims
1149     for (int i = 1 ; i < iDims ; ++i)
1150     {
1151         //init
1152         piIndex[i] = 0;
1153
1154         if (iLeadDims > piDims[i])
1155         {
1156             iLeadDims = piDims[i];
1157         }
1158     }
1159
1160     for (int i = 0 ; i < iLeadDims ; ++i)
1161     {
1162         for (int j = 0 ; j < iDims ; ++j)
1163         {
1164             piIndex[j] = i;
1165         }
1166
1167         int index = _pR->getIndex(piIndex);
1168         sub(dblLeft, _pR->get(index), pOut->get() + index);
1169     }
1170
1171     delete[] piIndex;
1172     return pOut;
1173 }
1174
1175 template<class T, class U, class O>
1176 InternalType* sub_I_MC(T *_pL, U *_pR)
1177 {
1178     int iDims = _pR->getDims();
1179     int* piDims = _pR->getDimsArray();
1180     O* pOut = (O*)opposite_MC<U, O>(_pR);
1181     double* pdblOut = pOut->get();
1182     double* pdblRight = _pR->get();
1183     double dblLeft = _pL->get(0);
1184     int iLeadDims = piDims[0];
1185     int* piIndex = new int[iDims];
1186     piIndex[0] = 0;
1187
1188     //find smaller dims
1189     for (int i = 1 ; i < iDims ; ++i)
1190     {
1191         //init
1192         piIndex[i] = 0;
1193
1194         if (iLeadDims > piDims[i])
1195         {
1196             iLeadDims = piDims[i];
1197         }
1198     }
1199
1200     for (int i = 0 ; i < iLeadDims ; ++i)
1201     {
1202         for (int j = 0 ; j < iDims ; ++j)
1203         {
1204             piIndex[j] = i;
1205         }
1206
1207         int index = _pR->getIndex(piIndex);
1208         sub(dblLeft, pdblRight[index], pdblOut + index);
1209     }
1210
1211     delete[] piIndex;
1212     return pOut;
1213 }
1214
1215 template<class T, class U, class O>
1216 InternalType* sub_IC_M(T *_pL, U *_pR)
1217 {
1218     int iDims = _pR->getDims();
1219     int* piDims = _pR->getDimsArray();
1220     O* pOut = (O*)opposite_M<U, O>(_pR);
1221     pOut->setComplex(true);
1222     double* pdblOutR = pOut->get();
1223     double* pdblOutI = pOut->getImg();
1224     double* pdblRight = _pR->get();
1225     double dblLeftR = _pL->get(0);
1226     double dblLeftI = _pL->getImg(0);
1227     int iLeadDims = piDims[0];
1228     int* piIndex = new int[iDims];
1229     piIndex[0] = 0;
1230
1231     //find smaller dims
1232     for (int i = 1 ; i < iDims ; ++i)
1233     {
1234         //init
1235         piIndex[i] = 0;
1236
1237         if (iLeadDims > piDims[i])
1238         {
1239             iLeadDims = piDims[i];
1240         }
1241     }
1242
1243     for (int i = 0 ; i < iLeadDims ; ++i)
1244     {
1245         for (int j = 0 ; j < iDims ; ++j)
1246         {
1247             piIndex[j] = i;
1248         }
1249
1250         int index = _pR->getIndex(piIndex);
1251         sub(&dblLeftR, &dblLeftI, (size_t)1, pdblRight[index], pdblOutR + index, pdblOutI + index);
1252     }
1253
1254     delete[] piIndex;
1255     return pOut;
1256 }
1257
1258 template<class T, class U, class O>
1259 InternalType* sub_IC_MC(T *_pL, U *_pR)
1260 {
1261     int iDims = _pR->getDims();
1262     int* piDims = _pR->getDimsArray();
1263     O* pOut = (O*)opposite_MC<U, O>(_pR);
1264     double dblLeftR = _pL->get(0);
1265     double dblLeftI = _pL->getImg(0);
1266     double* pdblOutR = pOut->get();
1267     double* pdblOutI = pOut->getImg();
1268     double* pdblRightR = _pR->get();
1269     double* pdblRightI = _pR->getImg();
1270     int iLeadDims = piDims[0];
1271     int* piIndex = new int[iDims];
1272     piIndex[0] = 0;
1273     //find smaller dims
1274     for (int i = 1 ; i < iDims ; ++i)
1275     {
1276         //init
1277         piIndex[i] = 0;
1278
1279         if (iLeadDims > piDims[i])
1280         {
1281             iLeadDims = piDims[i];
1282         }
1283     }
1284
1285     for (int i = 0 ; i < iLeadDims ; ++i)
1286     {
1287         for (int j = 0 ; j < iDims ; ++j)
1288         {
1289             piIndex[j] = i;
1290         }
1291
1292         int index = _pR->getIndex(piIndex);
1293         sub(dblLeftR, dblLeftI, pdblRightR[index], pdblRightI[index], pdblOutR + index, pdblOutI + index);
1294     }
1295
1296     delete[] piIndex;
1297     return pOut;
1298 }
1299
1300 template<class T, class U, class O>
1301 InternalType* sub_I_S(T *_pL, U *_pR)
1302 {
1303     return sub_S_S<T, U, O>(_pL, _pR);
1304 }
1305
1306 template<class T, class U, class O>
1307 InternalType* sub_IC_S(T *_pL, U *_pR)
1308 {
1309     return sub_SC_S<T, U, O>(_pL, _pR);
1310 }
1311
1312 template<class T, class U, class O>
1313 InternalType* sub_I_SC(T *_pL, U *_pR)
1314 {
1315     return sub_S_SC<T, U, O>(_pL, _pR);
1316 }
1317
1318 template<class T, class U, class O>
1319 InternalType* sub_IC_SC(T *_pL, U *_pR)
1320 {
1321     return sub_SC_SC<T, U, O>(_pL, _pR);
1322 }
1323
1324 template<class T, class U, class O> InternalType* sub_M_I(T *_pL, U *_pR)
1325 {
1326     int iDims = _pL->getDims();
1327     int* piDims = _pL->getDimsArray();
1328     O* pOut = (O*)_pL->clone();
1329     double* pdblOutR = pOut->get();
1330     double* pdblLeft = _pL->get();
1331     double dblRight = _pR->get(0);
1332     int iLeadDims = piDims[0];
1333     int* piIndex = new int[iDims];
1334     piIndex[0] = 0;
1335
1336     //find smaller dims
1337     for (int i = 1 ; i < iDims ; ++i)
1338     {
1339         //init
1340         piIndex[i] = 0;
1341
1342         if (iLeadDims > piDims[i])
1343         {
1344             iLeadDims = piDims[i];
1345         }
1346     }
1347
1348     for (int i = 0 ; i < iLeadDims ; ++i)
1349     {
1350         for (int j = 0 ; j < iDims ; ++j)
1351         {
1352             piIndex[j] = i;
1353         }
1354
1355         int index = _pL->getIndex(piIndex);
1356         sub(pdblLeft[index], (size_t)1, &dblRight, pdblOutR + index);
1357     }
1358
1359     delete[] piIndex;
1360     return pOut;
1361 }
1362
1363 template<class T, class U, class O> InternalType* sub_MC_I(T *_pL, U *_pR)
1364 {
1365     return sub_M_I<T, U, O>(_pL, _pR);
1366 }
1367
1368 template<class T, class U, class O> InternalType* sub_M_IC(T *_pL, U *_pR)
1369 {
1370     int iDims = _pL->getDims();
1371     int* piDims = _pL->getDimsArray();
1372     O* pOut = (O*)_pL->clone();
1373     pOut->setComplex(true);
1374     double* pdblOutR = pOut->get();
1375     double* pdblOutI = pOut->getImg();
1376     double* pdblLeft = _pL->get();
1377     double dblRightR = _pR->get(0);
1378     double dblRightI = _pR->getImg(0);
1379     int iLeadDims = piDims[0];
1380     int* piIndex = new int[iDims];
1381     piIndex[0] = 0;
1382
1383     //find smaller dims
1384     for (int i = 1 ; i < iDims ; ++i)
1385     {
1386         //init
1387         piIndex[i] = 0;
1388
1389         if (iLeadDims > piDims[i])
1390         {
1391             iLeadDims = piDims[i];
1392         }
1393     }
1394
1395     for (int i = 0 ; i < iLeadDims ; ++i)
1396     {
1397         for (int j = 0 ; j < iDims ; ++j)
1398         {
1399             piIndex[j] = i;
1400         }
1401
1402         int index = _pL->getIndex(piIndex);
1403         sub(pdblLeft[index], (size_t)1, &dblRightR, &dblRightI, pdblOutR + index, pdblOutI + index);
1404     }
1405
1406     delete[] piIndex;
1407     return pOut;
1408 }
1409
1410 template<class T, class U, class O> InternalType* sub_MC_IC(T *_pL, U *_pR)
1411 {
1412     int iDims = _pL->getDims();
1413     int* piDims = _pL->getDimsArray();
1414     O* pOut = (O*)_pL->clone();
1415     pOut->setComplex(true);
1416     double* pdblOutR = pOut->get();
1417     double* pdblOutI = pOut->getImg();
1418     double* pdblLeftR = _pL->get();
1419     double* pdblLeftI = _pL->getImg();
1420     double dblRightR = _pR->get(0);
1421     double dblRightI = _pR->getImg(0);
1422     int iLeadDims = piDims[0];
1423     int* piIndex = new int[iDims];
1424     piIndex[0] = 0;
1425
1426     //find smaller dims
1427     for (int i = 1 ; i < iDims ; ++i)
1428     {
1429         //init
1430         piIndex[i] = 0;
1431
1432         if (iLeadDims > piDims[i])
1433         {
1434             iLeadDims = piDims[i];
1435         }
1436     }
1437
1438     for (int i = 0 ; i < iLeadDims ; ++i)
1439     {
1440         for (int j = 0 ; j < iDims ; ++j)
1441         {
1442             piIndex[j] = i;
1443         }
1444
1445         int index = _pL->getIndex(piIndex);
1446         sub(pdblLeftR[index], pdblLeftI[index], (size_t)1, &dblRightR, &dblRightI, pdblOutR + index, pdblOutI + index);
1447     }
1448
1449     delete[] piIndex;
1450     return pOut;
1451 }
1452
1453 template<class T, class U, class O> InternalType* sub_S_I(T *_pL, U *_pR)
1454 {
1455     return sub_S_S<T, U, O>(_pL, _pR);
1456 }
1457
1458 template<class T, class U, class O> InternalType* sub_SC_I(T *_pL, U *_pR)
1459 {
1460     return sub_SC_SC<T, U, O>(_pL, _pR);
1461 }
1462
1463 template<class T, class U, class O> InternalType* sub_S_IC(T *_pL, U *_pR)
1464 {
1465     return sub_S_SC<T, U, O>(_pL, _pR);
1466 }
1467
1468 template<class T, class U, class O> InternalType* sub_SC_IC(T *_pL, U *_pR)
1469 {
1470     return sub_SC_SC<T, U, O>(_pL, _pR);
1471 }
1472
1473 template<class T, class U, class O> InternalType* sub_I_I(T *_pL, U *_pR)
1474 {
1475     O* pOut = types::Double::Identity(-1, -1);
1476     sub(_pL->get(0), _pR->get(0), pOut->get());
1477     return pOut;
1478 }
1479
1480 template<class T, class U, class O> InternalType* sub_I_IC(T *_pL, U *_pR)
1481 {
1482     O* pOut = types::Double::Identity(-1, -1);
1483     pOut->setComplex(true);
1484     sub(_pL->get(), (size_t)1, _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
1485     return pOut;
1486 }
1487
1488 template<class T, class U, class O> InternalType* sub_IC_I(T *_pL, U *_pR)
1489 {
1490     O* pOut = types::Double::Identity(-1, -1);
1491     pOut->setComplex(true);
1492     sub(_pL->get(0), _pL->getImg(0), _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
1493     return pOut;
1494 }
1495
1496 template<class T, class U, class O> InternalType* sub_IC_IC(T *_pL, U *_pR)
1497 {
1498     O* pOut = types::Double::Identity(-1, -1);
1499     pOut->setComplex(true);
1500     sub(_pL->get(0), _pL->getImg(0), _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
1501     return pOut;
1502 }
1503
1504 template<class T, class U, class O> types::InternalType* sub_I_E(T *_pL, U * /*_pR*/)
1505 {
1506     if (ConfigVariable::getOldEmptyBehaviour())
1507     {
1508         Sciwarning(_("operation -: Warning adding a matrix with the empty matrix old behaviour.\n"));
1509         return _pL;
1510     }
1511
1512     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n"));
1513     Double* pOut = Double::Empty();
1514     sub();
1515     return pOut;
1516 }
1517
1518 template<class T, class U, class O> types::InternalType* sub_IC_E(T *_pL, U * /*_pR*/)
1519 {
1520     if (ConfigVariable::getOldEmptyBehaviour())
1521     {
1522         Sciwarning(_("operation -: Warning adding a matrix with the empty matrix old behaviour.\n"));
1523         return _pL;
1524     }
1525
1526     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n"));
1527     Double* pOut = Double::Empty();
1528     sub();
1529     return pOut;
1530 }
1531
1532 template<class T, class U, class O> types::InternalType* sub_E_I(T * /*_pL*/, U *_pR)
1533 {
1534     if (ConfigVariable::getOldEmptyBehaviour())
1535     {
1536         Sciwarning(_("operation -: Warning adding a matrix with the empty matrix old behaviour.\n"));
1537         return opposite_I<U, O>(_pR);
1538     }
1539
1540     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n"));
1541     Double* pOut = Double::Empty();
1542     sub();
1543     return pOut;
1544 }
1545
1546 template<class T, class U, class O> types::InternalType* sub_E_IC(T * /*_pL*/, U *_pR)
1547 {
1548     if (ConfigVariable::getOldEmptyBehaviour())
1549     {
1550         Sciwarning(_("operation -: Warning adding a matrix with the empty matrix old behaviour.\n"));
1551         return opposite_IC<U, O>(_pR);
1552     }
1553
1554     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n"));
1555     Double* pOut = Double::Empty();
1556     sub();
1557     return pOut;
1558 }
1559
1560 template<> InternalType* sub_M_M<Polynom, Polynom, Polynom>(Polynom* _pL, Polynom* _pR)
1561 {
1562     Polynom* pOut = NULL;
1563
1564     //check varname
1565     if (_pL->getVariableName() != _pR->getVariableName())
1566     {
1567         //call overload
1568         return NULL;
1569     }
1570
1571     if (_pL->isScalar())
1572     {
1573         SinglePoly* p1Coef  = _pL->get(0);
1574         int iRank1          = p1Coef->getRank();
1575         int* pRank2         = new int[_pR->getSize()];
1576         int* pRankOut       = new int[_pR->getSize()];
1577
1578         _pR->getRank(pRank2);
1579         for (int i = 0 ; i < _pR->getSize() ; i++)
1580         {
1581             pRankOut[i] = std::max(iRank1, pRank2[i]);
1582         }
1583
1584         pOut = new Polynom(_pR->getVariableName(), _pR->getDims(), _pR->getDimsArray(), pRankOut);
1585         bool isOutComplex = _pL->isComplex() || _pR->isComplex();
1586
1587         //Result P1(0) - P2(i)
1588         double* p1R = p1Coef->get();
1589         if (isOutComplex)
1590         {
1591             double* p1I = p1Coef->getImg();
1592             for (int i = 0 ; i < _pR->getSize() ; i++)
1593             {
1594                 SinglePoly* p2Coef   = _pR->get(i);
1595                 double* p2R          = p2Coef->get();
1596                 double* p2I          = p2Coef->getImg();
1597
1598                 SinglePoly* pOutCoef = pOut->get(i);
1599                 pOutCoef->setComplex(isOutComplex);
1600                 double* pOutR        = pOutCoef->get();
1601                 double* pOutI        = pOutCoef->getImg();
1602
1603                 for (int j = 0 ; j < pRankOut[i] + 1 ; j++)
1604                 {
1605                     if (j > iRank1)
1606                     {
1607                         pOutR[j] = - p2R[j];
1608                         pOutI[j] = - (p2I ? p2I[j] : 0);
1609                     }
1610                     else if (j > pRank2[i])
1611                     {
1612                         pOutR[j] = p1R[j];
1613                         pOutI[j] = (p1I ? p1I[j] : 0);
1614                     }
1615                     else
1616                     {
1617                         pOutR[j] = p1R[j] - p2R[j];
1618                         pOutI[j] = (p1I ? p1I[j] : 0) - (p2I ? p2I[j] : 0);
1619                     }
1620                 }
1621             }
1622         }
1623         else
1624         {
1625             for (int i = 0 ; i < _pR->getSize() ; i++)
1626             {
1627                 SinglePoly* p2Coef   = _pR->get(i);
1628                 double* p2R          = p2Coef->get();
1629
1630                 SinglePoly* pOutCoef = pOut->get(i);
1631                 double* pOutR        = pOutCoef->get();
1632
1633                 for (int j = 0 ; j < pRankOut[i] + 1 ; j++)
1634                 {
1635                     if (j > iRank1)
1636                     {
1637                         pOutR[j] = - p2R[j];
1638                     }
1639                     else if (j > pRank2[i])
1640                     {
1641                         pOutR[j] = p1R[j];
1642                     }
1643                     else
1644                     {
1645                         pOutR[j] = p1R[j] - p2R[j];
1646                     }
1647                 }
1648             }
1649         }
1650
1651         delete[] pRankOut;
1652         delete[] pRank2;
1653
1654         if (pOut != NULL)
1655         {
1656             pOut->updateRank();
1657         }
1658
1659         return pOut;
1660     }
1661
1662     if (_pR->isScalar())
1663     {
1664         //size(p2) == 1
1665         int *pRank1     = new int[_pL->getSize()];
1666         int iRank2      = _pR->get(0)->getRank();
1667         int *pRankOut   = new int[_pL->getSize()];
1668
1669         _pL->getRank(pRank1);
1670         for (int i = 0 ; i < _pL->getSize() ; i++)
1671         {
1672             pRankOut[i] = std::max(pRank1[i], iRank2);
1673         }
1674
1675         pOut = new Polynom(_pL->getVariableName(), _pL->getDims(), _pL->getDimsArray(), pRankOut);
1676         bool isOutComplex = _pL->isComplex() || _pR->isComplex();
1677
1678         //Result P1(i) - P2(0)
1679         SinglePoly *p2Coef          = _pR->get(0);
1680         double *p2R                 = p2Coef->get();
1681
1682         if (isOutComplex)
1683         {
1684             double *p2I             = p2Coef->getImg();
1685             for (int i = 0 ; i < _pL->getSize() ; i++)
1686             {
1687                 SinglePoly *p1Coef      = _pL->get(i);
1688                 double *p1R             = p1Coef->get();
1689                 double *p1I             = p1Coef->getImg();
1690
1691                 SinglePoly *pOutCoef    = pOut->get(i);
1692                 pOutCoef->setComplex(isOutComplex);
1693                 double *pOutR           = pOutCoef->get();
1694                 double *pOutI           = pOutCoef->getImg();
1695
1696                 for (int j = 0 ; j < pRankOut[i] + 1 ; j++)
1697                 {
1698                     if (j > pRank1[j])
1699                     {
1700                         pOutR[j] = - p2R[j];
1701                         pOutI[j] = - (p2I ? p2I[j] : 0);
1702                     }
1703                     else if (j > iRank2)
1704                     {
1705                         pOutR[j] = p1R[j];
1706                         pOutI[j] = (p1I ? p1I[j] : 0);
1707                     }
1708                     else
1709                     {
1710                         pOutR[j] = p1R[j] - p2R[j];
1711                         pOutI[j] = (p1I ? p1I[j] : 0) - (p2I ? p2I[j] : 0);
1712                     }
1713                 }
1714             }
1715         }
1716         else
1717         {
1718             for (int i = 0 ; i < _pL->getSize() ; i++)
1719             {
1720                 SinglePoly *p1Coef      = _pL->get(i);
1721                 double *p1R             = p1Coef->get();
1722
1723                 SinglePoly *pOutCoef    = pOut->get(i);
1724                 double *pOutR           = pOutCoef->get();
1725
1726                 for (int j = 0 ; j < pRankOut[i] + 1 ; j++)
1727                 {
1728                     if (j > pRank1[j])
1729                     {
1730                         pOutR[j] = - p2R[j];
1731                     }
1732                     else if (j > iRank2)
1733                     {
1734                         pOutR[j] = p1R[j];
1735                     }
1736                     else
1737                     {
1738                         pOutR[j] = p1R[j] - p2R[j];
1739                     }
1740                 }
1741             }
1742         }
1743
1744         delete[] pRankOut;
1745         delete[] pRank1;
1746
1747         if (pOut != NULL)
1748         {
1749             pOut->updateRank();
1750         }
1751
1752         return pOut;
1753     }
1754
1755     //check dimension compatibilities
1756     int iDims1  = _pL->getDims();
1757     int iDims2  = _pR->getDims();
1758
1759     if (iDims1 != iDims2)
1760     {
1761         return nullptr;
1762     }
1763
1764     int* piDims1    = _pL->getDimsArray();
1765     int* piDims2    = _pR->getDimsArray();
1766
1767     for (int i = 0 ; i < iDims1 ; i++)
1768     {
1769         if (piDims1[i] != piDims2[i])
1770         {
1771             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
1772         }
1773     }
1774
1775     int *pRankOut   = new int[_pL->getSize()];
1776     int *pRank1     = new int[_pL->getSize()];
1777     int *pRank2     = new int[_pR->getSize()];
1778
1779     _pL->getRank(pRank1);
1780     _pR->getRank(pRank2);
1781     for (int i = 0 ; i < _pL->getSize() ; i++)
1782     {
1783         pRankOut[i] = std::max(pRank1[i], pRank2[i]);
1784     }
1785
1786     pOut = new Polynom(_pR->getVariableName(), iDims1, piDims1, pRankOut);
1787     bool isOutComplex =  _pL->isComplex() || _pR->isComplex();
1788
1789     //Result P1(i) - P2(i)
1790     for (int i = 0 ; i < _pL->getSize() ; i++)
1791     {
1792         SinglePoly* pOutCoef = pOut->get(i);
1793         pOutCoef->setComplex(isOutComplex);
1794
1795         double *p1R     = _pL->get(i)->get();
1796         double *p2R     = _pR->get(i)->get();
1797         double *p1I     = _pL->get(i)->getImg();
1798         double *p2I     = _pR->get(i)->getImg();
1799         double *pOutR   = pOutCoef->get();
1800         int iMin        = std::min(pRank1[i], pRank2[i]);
1801         int iMax        = std::max(pRank1[i], pRank2[i]);
1802
1803         for (int j = 0 ; j < iMin + 1 ; j++)
1804         {
1805             pOutR[j]    = p1R[j] - p2R[j];
1806         }
1807
1808         double *pTemp   = NULL;
1809         double *pTempI  = NULL;
1810         int iCoef       = 1;
1811         if (pRank1[i] > pRank2[i])
1812         {
1813             pTemp       = p1R;
1814             pTempI      = p1I;
1815             iCoef       = 1;
1816         }
1817         else
1818         {
1819             pTemp       = p2R;
1820             pTempI      = p2I;
1821             iCoef       = -1;
1822         }
1823
1824         for (int j = iMin + 1 ; j < iMax + 1 ; j++)
1825         {
1826             pOutR[j]    = pTemp[j] * iCoef;
1827         }
1828
1829         if (isOutComplex)
1830         {
1831             double *pOutI   = pOutCoef->getImg();
1832
1833             for (int j = 0 ; j < iMin + 1 ; j++)
1834             {
1835                 pOutI[j]    = (p1I == NULL ? 0 : p1I[j]) - (p2I == NULL ? 0 : p2I[j]);
1836             }
1837
1838             if (pTempI != NULL)
1839             {
1840                 for (int j = iMin + 1 ; j < iMax + 1 ; j++)
1841                 {
1842                     pOutI[j]  = pTempI[j] * iCoef;
1843                 }                
1844             }
1845         }
1846     }
1847
1848     delete[] pRankOut;
1849     delete[] pRank1;
1850     delete[] pRank2;
1851
1852     if (pOut != NULL)
1853     {
1854         pOut->updateRank();
1855     }
1856
1857     return pOut;
1858 }
1859
1860 // P - D
1861 template<> InternalType* sub_M_M<Polynom, Double, Polynom>(Polynom* _pL, Double* _pR)
1862 {
1863     Polynom* pOut = NULL;
1864
1865     bool bComplex1 = _pL->isComplex();
1866     bool bComplex2 = _pR->isComplex();
1867     bool isComplexOut = bComplex2 || bComplex1;
1868
1869     // P - []
1870     if (_pR->isEmpty())
1871     {
1872         return _pL;
1873     }
1874
1875     if (_pR->isIdentity())
1876     {
1877         double dblRightR = _pR->get(0);
1878         double dblRightI = 0;
1879         if (_pR->isComplex())
1880         {
1881             dblRightI = _pR->getImg(0);
1882         }
1883
1884         int iDims = _pL->getDims();
1885         int* piDims = _pL->getDimsArray();
1886         pOut = (Polynom*)_pL->clone();
1887         pOut->setComplex(isComplexOut);
1888         SinglePoly** pSPOut = pOut->get();
1889         SinglePoly** pSPLeft = _pL->get();
1890         int iLeadDims = piDims[0];
1891         int* piIndex = new int[iDims];
1892         piIndex[0] = 0;
1893
1894         //find smaller dims
1895         for (int i = 1 ; i < iDims ; ++i)
1896         {
1897             //init
1898             piIndex[i] = 0;
1899
1900             if (iLeadDims > piDims[i])
1901             {
1902                 iLeadDims = piDims[i];
1903             }
1904         }
1905
1906         if (isComplexOut)
1907         {
1908             for (int i = 0 ; i < iLeadDims ; ++i)
1909             {
1910                 for (int j = 0 ; j < iDims ; ++j)
1911                 {
1912                     piIndex[j] = i;
1913                 }
1914
1915                 int index = _pL->getIndex(piIndex);
1916                 sub(pSPLeft[index]->get(0), pSPLeft[index]->getImg(0), (size_t)1, &dblRightR, &dblRightI, pSPOut[index]->get(), pSPOut[index]->getImg());
1917             }
1918         }
1919         else
1920         {
1921             for (int i = 0 ; i < iLeadDims ; ++i)
1922             {
1923                 for (int j = 0 ; j < iDims ; ++j)
1924                 {
1925                     piIndex[j] = i;
1926                 }
1927
1928                 int index = _pL->getIndex(piIndex);
1929                 sub(pSPLeft[index]->get(0), (size_t)1, &dblRightR, pSPOut[index]->get());
1930             }
1931         }
1932
1933         delete[] piIndex;
1934         return pOut;
1935
1936     }
1937
1938     if (_pR->isScalar())
1939     {
1940         //subtract same value to all polynoms
1941         pOut = (Polynom*)_pL->clone();
1942         pOut->setComplex(isComplexOut);
1943
1944         SinglePoly** pSP = pOut->get();
1945         int iSize = pOut->getSize();
1946
1947         double dblR = _pR->get(0);
1948         if (isComplexOut)
1949         {
1950             double dblI = _pR->getImg(0);
1951             for (int i = 0 ; i < iSize ; ++i)
1952             {
1953                 pSP[i]->get()[0] -= dblR;
1954                 pSP[i]->getImg()[0] -= dblI;
1955             }
1956         }
1957         else
1958         {
1959             for (int i = 0 ; i < iSize ; ++i)
1960             {
1961                 pSP[i]->get()[0] -= dblR;
1962             }
1963         }
1964
1965         return pOut;
1966     }
1967
1968     if (_pL->isScalar())
1969     {
1970         //and _pR is not !
1971         int iDims = _pR->getDims();
1972         int* piDims = _pR->getDimsArray();
1973         int iSize = _pR->getSize();
1974
1975         pOut = new Polynom(_pL->getVariableName(), iDims, piDims);
1976
1977         SinglePoly* pSPL = _pL->get(0);
1978         if (_pR->isComplex())
1979         {
1980             double* pdblR = _pR->get();
1981             double* pdblI = _pR->getImg();
1982
1983             for (int i = 0 ; i < iSize ; ++i)
1984             {
1985                 SinglePoly* pSPOut = pSPL->clone();
1986                 //in case of original is not complex
1987                 pSPOut->setComplex(isComplexOut);
1988                 pSPOut->get()[0] -= pdblR[i];
1989                 pSPOut->getImg()[0] -= pdblI[i];
1990                 pOut->set(i, pSPOut);
1991                 delete pSPOut;
1992             }
1993         }
1994         else
1995         {
1996             double* pdblR = _pR->get();
1997
1998             for (int i = 0 ; i < iSize ; ++i)
1999             {
2000                 SinglePoly* pSPOut = pSPL->clone();
2001                 //update 0th rank value
2002                 pSPOut->get()[0] -= pdblR[i];
2003                 pOut->set(i, pSPOut);
2004                 delete pSPOut;
2005             }
2006         }
2007
2008         return pOut;
2009     }
2010
2011     //P - D
2012
2013     //check dimensions
2014     int iDims1 = _pR->getDims();
2015     int iDims2 = _pL->getDims();
2016
2017     if (iDims1 != iDims2)
2018     {
2019         return nullptr;
2020     }
2021
2022     int* piDims1 = _pR->getDimsArray();
2023     int* piDims2 = _pL->getDimsArray();
2024
2025     for (int i = 0 ; i < iDims1 ; i++)
2026     {
2027         if (piDims1[i] != piDims2[i])
2028         {
2029             wchar_t pMsg[bsiz];
2030             os_swprintf(pMsg, bsiz, _W("Error: operator %ls: Matrix dimensions must agree (op1 is %ls, op2 is %ls).\n").c_str(),  L"-", _pL->DimToString().c_str(), _pR->DimToString().c_str());
2031             throw ast::InternalError(pMsg);
2032         }
2033     }
2034
2035     double* pInDblR = _pR->get();
2036     pOut = (Polynom*)_pL->clone();
2037     if (bComplex1 && bComplex2)
2038     {
2039         double* pInDblI = _pR->getImg();
2040         for (int i = 0 ; i < pOut->getSize() ; i++)
2041         {
2042             SinglePoly *pSPOut   = pOut->get(i);
2043             double *pOutPolyR    = pSPOut->get();
2044             double *pOutPolyI    = pSPOut->getImg();
2045
2046             pOutPolyR[0] -= pInDblR[i];
2047             pOutPolyI[0] -= pInDblI[i];
2048         }
2049     }
2050     else if (bComplex2)
2051     {
2052         double* pInDblI = _pR->getImg();
2053         pOut->setComplex(true);
2054         for (int i = 0 ; i < pOut->getSize() ; i++)
2055         {
2056             SinglePoly *pSPOut   = pOut->get(i);
2057             double *pOutPolyR    = pSPOut->get();
2058             double *pOutPolyI    = pSPOut->getImg();
2059
2060             pOutPolyR[0] -= pInDblR[i];
2061             pOutPolyI[0] = -pInDblI[i];
2062         }
2063     }
2064     else
2065     {
2066         for (int i = 0 ; i < pOut->getSize() ; i++)
2067         {
2068             SinglePoly *pSPOut = pOut->get(i);
2069             double *pOutPolyR  = pSPOut->get();
2070
2071             pOutPolyR[0] -= pInDblR[i];
2072         }
2073     }
2074
2075     return pOut;
2076 }
2077
2078 template<> InternalType* sub_I_M<Double, Polynom, Polynom>(Double* _pL, Polynom* _pR)
2079 {
2080     Polynom* pOut = (Polynom*)opposite_M<Polynom, Polynom>(_pR);
2081     double dblLeft = _pL->get(0);
2082
2083     int iDims = _pR->getDims();
2084     int* piDims = _pR->getDimsArray();
2085     int iLeadDims = piDims[0];
2086     int* piIndex = new int[iDims];
2087     SinglePoly** pSP = _pR->get();
2088     SinglePoly** pSPOut = pOut->get();
2089     piIndex[0] = 0;
2090
2091     //find smaller dims
2092     for (int i = 1 ; i < iDims ; ++i)
2093     {
2094         //init
2095         piIndex[i] = 0;
2096
2097         if (iLeadDims > piDims[i])
2098         {
2099             iLeadDims = piDims[i];
2100         }
2101     }
2102
2103     for (int i = 0 ; i < iLeadDims ; ++i)
2104     {
2105         for (int j = 0 ; j < iDims ; ++j)
2106         {
2107             piIndex[j] = i;
2108         }
2109
2110         int index = _pR->getIndex(piIndex);
2111         sub(dblLeft, pSP[index]->get(0), pSPOut[index]->get());
2112     }
2113
2114     delete[] piIndex;
2115     return pOut;
2116 }
2117
2118 template<> InternalType* sub_I_MC<Double, Polynom, Polynom>(Double* _pL, Polynom* _pR)
2119 {
2120     Polynom* pOut = (Polynom*)opposite_MC<Polynom, Polynom>(_pR);
2121     double dblLeft = _pL->get(0);
2122
2123     int iDims = _pR->getDims();
2124     int* piDims = _pR->getDimsArray();
2125     int iLeadDims = piDims[0];
2126     int* piIndex = new int[iDims];
2127     SinglePoly** pSP = _pR->get();
2128     SinglePoly** pSPOut = pOut->get();
2129     piIndex[0] = 0;
2130
2131     for (int i = 0 ; i < iDims ; ++i)
2132     {
2133         //init
2134         piIndex[i] = 0;
2135
2136         if (iLeadDims > piDims[i])
2137         {
2138             iLeadDims = piDims[i];
2139         }
2140     }
2141
2142     for (int i = 0 ; i < iLeadDims ; ++i)
2143     {
2144         for (int j = 0 ; j < iDims ; ++j)
2145         {
2146             piIndex[j] = i;
2147         }
2148
2149         int index = _pR->getIndex(piIndex);
2150         sub(dblLeft, pSP[index]->get(0), pSPOut[index]->get());
2151     }
2152
2153     delete[] piIndex;
2154     return pOut;
2155 }
2156
2157 template<> InternalType* sub_IC_M<Double, Polynom, Polynom>(Double* _pL, Polynom* _pR)
2158 {
2159     Polynom* pOut = (Polynom*)opposite_M<Polynom, Polynom>(_pR);
2160     pOut->setComplex(true);
2161     double dblLeftR = _pL->get(0);
2162     double dblLeftI = _pL->getImg(0);
2163
2164     int iDims = _pR->getDims();
2165     int* piDims = _pR->getDimsArray();
2166     int iLeadDims = piDims[0];
2167     int* piIndex = new int[iDims];
2168     SinglePoly** pSP = _pR->get();
2169     SinglePoly** pSPOut = pOut->get();
2170     piIndex[0] = 0;
2171
2172     for (int i = 0 ; i < iDims ; ++i)
2173     {
2174         //init
2175         piIndex[i] = 0;
2176
2177         if (iLeadDims > piDims[i])
2178         {
2179             iLeadDims = piDims[i];
2180         }
2181     }
2182
2183     for (int i = 0 ; i < iLeadDims ; ++i)
2184     {
2185         for (int j = 0 ; j < iDims ; ++j)
2186         {
2187             piIndex[j] = i;
2188         }
2189
2190         int index = _pR->getIndex(piIndex);
2191         sub(&dblLeftR, &dblLeftI, (size_t)1, pSP[index]->get(0), pSPOut[index]->get(), pSPOut[index]->getImg());
2192     }
2193
2194     delete[] piIndex;
2195     return pOut;
2196 }
2197
2198 template<> InternalType* sub_IC_MC<Double, Polynom, Polynom>(Double* _pL, Polynom* _pR)
2199 {
2200     Polynom* pOut = (Polynom*)opposite_MC<Polynom, Polynom>(_pR);
2201     double dblLeftR = _pL->get(0);
2202     double dblLeftI = _pL->getImg(0);
2203
2204     int iDims = _pR->getDims();
2205     int* piDims = _pR->getDimsArray();
2206     int iLeadDims = piDims[0];
2207     int* piIndex = new int[iDims];
2208     SinglePoly** pSP = _pR->get();
2209     SinglePoly** pSPOut = pOut->get();
2210     piIndex[0] = 0;
2211
2212     for (int i = 0 ; i < iDims ; ++i)
2213     {
2214         //init
2215         piIndex[i] = 0;
2216
2217         if (iLeadDims > piDims[i])
2218         {
2219             iLeadDims = piDims[i];
2220         }
2221     }
2222
2223     for (int i = 0 ; i < iLeadDims ; ++i)
2224     {
2225         for (int j = 0 ; j < iDims ; ++j)
2226         {
2227             piIndex[j] = i;
2228         }
2229
2230         int index = _pR->getIndex(piIndex);
2231         sub(dblLeftR, dblLeftI, pSP[index]->get(0), pSP[index]->getImg(0), pSPOut[index]->get(), pSPOut[index]->getImg());
2232     }
2233
2234     delete[] piIndex;
2235     return pOut;
2236 }
2237
2238 template<> InternalType* sub_M_M<Double, Polynom, Polynom>(Double* _pL, Polynom* _pR)
2239 {
2240     Polynom* pOut = NULL;
2241     bool bComplex1 = _pL->isComplex();
2242     bool bComplex2 = _pR->isComplex();
2243
2244     double *pInDblR = _pL->getReal();
2245     double *pInDblI = _pL->getImg();
2246
2247     if (_pL->isEmpty())
2248     {
2249         return _pR;
2250     }
2251
2252     if (_pR->isScalar())
2253     {
2254         int *piRank = new int[_pL->getSize()];
2255         for (int i = 0 ; i < _pL->getSize() ; i++)
2256         {
2257             piRank[i] = _pR->get(0)->getRank();
2258         }
2259
2260         pOut = new Polynom(_pR->getVariableName(), _pL->getDims(), _pL->getDimsArray(), piRank);
2261         delete[] piRank;
2262         if (bComplex1 || bComplex2)
2263         {
2264             pOut->setComplex(true);
2265         }
2266
2267         for (int i = 0 ; i < pOut->getSize() ; i++)
2268         {
2269             SinglePoly *pInPoly  = _pR->get(0);
2270             SinglePoly *pOutPoly = pOut->get(i);
2271             double *pInPolyR     = pInPoly->get();
2272             double *pOutPolyR    = pOutPoly->get();
2273
2274             pOutPolyR[0] = pInDblR[i] - pInPolyR[0];
2275
2276             for (int j = 1 ; j < pInPoly->getSize() ; j++)
2277             {
2278                 pOutPolyR[j] = -pInPolyR[j];
2279             }
2280         }
2281
2282         if (pOut->isComplex())
2283         {
2284             for (int i = 0 ; i < pOut->getSize() ; i++)
2285             {
2286                 SinglePoly *pInPoly  = _pR->get(0);
2287                 SinglePoly *pOutPoly = pOut->get(i);
2288                 double *pInPolyI     = pInPoly->getImg();
2289                 double *pOutPolyI    = pOutPoly->getImg();
2290
2291                 pOutPolyI[0] = (pInDblI != NULL ? pInDblI[i] : 0) - (pInPolyI != NULL ? pInPolyI[0] : 0);
2292
2293                 for (int j = 1 ; j < pInPoly->getSize() ; j++)
2294                 {
2295                     pOutPolyI[j] = (pInPolyI != NULL ? -pInPolyI[j] : 0);
2296                 }
2297             }
2298         }
2299
2300         return pOut;
2301     }
2302
2303     if (_pL->isScalar())
2304     {
2305         if (bComplex2)
2306         {
2307             pOut = (Polynom*)opposite_MC<Polynom, Polynom>(_pR);
2308         }
2309         else
2310         {
2311             pOut = (Polynom*)opposite_M<Polynom, Polynom>(_pR);
2312         }
2313
2314         int iSize = pOut->getSize();
2315
2316         if (bComplex1 && bComplex2)
2317         {
2318             for (int i = 0 ; i < iSize ; i++)
2319             {
2320                 SinglePoly *pSPOut   = pOut->get(i);
2321                 double *pOutPolyR    = pSPOut->get();
2322                 double *pOutPolyI    = pSPOut->getImg();
2323
2324                 pOutPolyR[0] += pInDblR[0];
2325                 pOutPolyI[0] += pInDblI[0];
2326             }
2327         }
2328         else if (bComplex1)
2329         {
2330             pOut->setComplex(true);
2331             for (int i = 0 ; i < iSize ; i++)
2332             {
2333                 SinglePoly *pSPOut   = pOut->get(i);
2334                 double *pOutPolyR    = pSPOut->get();
2335                 double *pOutPolyI    = pSPOut->getImg();
2336
2337                 pOutPolyR[0] += pInDblR[0];
2338                 pOutPolyI[0] = pInDblI[0];
2339             }
2340         }
2341         else
2342         {
2343             for (int i = 0 ; i < iSize ; i++)
2344             {
2345                 SinglePoly *pSPOut = pOut->get(i);
2346                 double *pOutPolyR  = pSPOut->get();
2347
2348                 pOutPolyR[0] += pInDblR[0];
2349             }
2350         }
2351
2352         return pOut;
2353     }
2354
2355     int iDims1 = _pR->getDims();
2356     int iDims2 = _pL->getDims();
2357
2358     if (iDims1 != iDims2)
2359     {
2360         return nullptr;
2361     }
2362
2363     int* piDims1 = _pR->getDimsArray();
2364     int* piDims2 = _pL->getDimsArray();
2365
2366     for (int i = 0 ; i < iDims1 ; i++)
2367     {
2368         if (piDims1[i] != piDims2[i])
2369         {
2370             wchar_t pMsg[bsiz];
2371             os_swprintf(pMsg, bsiz, _W("Error: operator %ls: Matrix dimensions must agree (op1 is %ls, op2 is %ls).\n").c_str(),  L"+", _pL->DimToString().c_str(), _pR->DimToString().c_str());
2372             throw ast::InternalError(pMsg);
2373         }
2374     }
2375
2376     if (bComplex2)
2377     {
2378         pOut = (Polynom*)opposite_MC<Polynom, Polynom>(_pR);
2379     }
2380     else
2381     {
2382         pOut = (Polynom*)opposite_M<Polynom, Polynom>(_pR);
2383     }
2384
2385     if (bComplex1 && bComplex2)
2386     {
2387         for (int i = 0 ; i < pOut->getSize() ; i++)
2388         {
2389             SinglePoly *pSPOut   = pOut->get(i);
2390             double *pOutPolyR    = pSPOut->get();
2391             double *pOutPolyI    = pSPOut->getImg();
2392
2393             pOutPolyR[0] += pInDblR[i];
2394             pOutPolyI[0] += pInDblI[i];
2395         }
2396     }
2397     else if (bComplex1)
2398     {
2399         pOut->setComplex(true);
2400         for (int i = 0 ; i < pOut->getSize() ; i++)
2401         {
2402             SinglePoly *pSPOut   = pOut->get(i);
2403             double *pOutPolyR    = pSPOut->get();
2404             double *pOutPolyI    = pSPOut->getImg();
2405
2406             pOutPolyR[0] += pInDblR[i];
2407             pOutPolyI[0] = pInDblI[i];
2408         }
2409     }
2410     else
2411     {
2412         for (int i = 0 ; i < pOut->getSize() ; i++)
2413         {
2414             SinglePoly *pSPOut = pOut->get(i);
2415             double *pOutPolyR  = pSPOut->get();
2416
2417             pOutPolyR[0] += pInDblR[i];
2418         }
2419     }
2420
2421     return pOut;
2422 }
2423
2424 //sp - sp
2425 template<> InternalType* sub_M_M<Sparse, Sparse, Sparse>(Sparse* _pL, Sparse* _pR)
2426 {
2427     //check scalar hidden in a sparse ;)
2428     if (_pL->isScalar() || _pR->isScalar())
2429     {
2430         // scalar - sp  or  sp - scalar
2431         // call Overload
2432         return NULL;
2433     }
2434
2435     if (_pL->getRows() != _pR->getRows() || _pL->getCols() != _pR->getCols())
2436     {
2437         //dimensions not match
2438         throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
2439     }
2440
2441     types::Sparse* pSPOut = _pL->substract(*_pR);
2442     pSPOut->finalize();
2443     return pSPOut;
2444 }
2445
2446 //d - sp
2447 template<> InternalType* sub_M_M<Double, Sparse, Double>(Double* _pL, Sparse* _pR)
2448 {
2449     Double* pOut = NULL;
2450     int iOne = 1; //fortran
2451     bool bComplex1 = _pL->isComplex();
2452     bool bComplex2 = _pR->isComplex();
2453
2454     if (_pL->isScalar() && _pR->isScalar())
2455     {
2456         //d - sp
2457         pOut = (Double*)_pL->clone();
2458         pOut->setComplex(bComplex1 || bComplex2);
2459         if (bComplex2)
2460         {
2461             std::complex<double> dbl = _pR->getImg(0, 0);
2462             pOut->set(0, pOut->get(0) - dbl.real());
2463             pOut->setImg(0, pOut->getImg(0) - dbl.imag());
2464         }
2465         else
2466         {
2467             pOut->set(0, pOut->get(0) - _pR->get(0, 0));
2468         }
2469
2470         return pOut;
2471     }
2472
2473     if (_pL->isScalar())
2474     {
2475         //d - SP
2476         pOut = new Double(_pR->getRows(), _pR->getCols(), bComplex1 || bComplex2);
2477         int iSize = _pR->getSize();
2478         double dblVal = _pL->get(0);
2479         double dblValI = 0;
2480         C2F(dset)(&iSize, &dblVal, pOut->get(), &iOne);
2481         if (bComplex1)
2482         {
2483             //initialize imag part at 0
2484             dblValI = _pL->getImg(0);
2485             C2F(dset)(&iSize, &dblValI, pOut->getImg(), &iOne);
2486         }
2487         else if (bComplex2)
2488         {
2489             dblValI = 0;
2490             C2F(dset)(&iSize, &dblValI, pOut->getImg(), &iOne);
2491         }
2492
2493         int nonZeros = static_cast<int>(_pR->nonZeros());
2494         int* pRows = new int[nonZeros * 2];
2495         _pR->outputRowCol(pRows);
2496         int* pCols = pRows + nonZeros;
2497
2498         if (pOut->isComplex())
2499         {
2500             for (int i = 0 ; i < nonZeros ; i++)
2501             {
2502                 int iRow = static_cast<int>(pRows[i]) - 1;
2503                 int iCol = static_cast<int>(pCols[i]) - 1;
2504                 std::complex<double> dbl = _pR->getImg(iRow, iCol);
2505                 pOut->set(iRow, iCol, dblVal - dbl.real());
2506                 pOut->setImg(iRow, iCol, dblValI - dbl.imag());
2507             }
2508         }
2509         else
2510         {
2511             for (int i = 0 ; i < nonZeros ; i++)
2512             {
2513                 int iRow = static_cast<int>(pRows[i]) - 1;
2514                 int iCol = static_cast<int>(pCols[i]) - 1;
2515                 pOut->set(iRow, iCol, dblVal - _pR->get(iRow, iCol));
2516             }
2517         }
2518
2519         //clear
2520         delete[] pRows;
2521
2522         return pOut;
2523     }
2524
2525     if (_pR->isScalar())
2526     {
2527         //D - sp
2528         pOut = (Double*)_pL->clone();
2529         pOut->setComplex(bComplex1 || bComplex2);
2530
2531         if (pOut->isComplex())
2532         {
2533             double* pReal = pOut->get();
2534             double* pImg = pOut->getImg();
2535             int size = pOut->getSize();
2536             std::complex<double> dbl = _pR->getImg(0, 0);
2537             for (int i = 0 ; i < size ; i++)
2538             {
2539                 pReal[i] -= dbl.real();
2540                 pImg[i] -= dbl.imag();
2541             }
2542         }
2543         else
2544         {
2545             double* pReal = pOut->get();
2546             int size = pOut->getSize();
2547             double dblTmp = _pR->get(0);
2548             for (int i = 0 ; i < size ; i++)
2549             {
2550                 pReal[i] -= dblTmp;
2551             }
2552         }
2553
2554         return pOut;
2555     }
2556
2557     if (_pL->getDims() > 2)
2558     {
2559         return nullptr;
2560     }
2561
2562     if (_pL->getRows() == _pR->getRows() && _pL->getCols() == _pR->getCols())
2563     {
2564         //D - SP
2565         pOut = (Double*)_pL->clone();
2566         pOut->setComplex(bComplex1 || bComplex2);
2567
2568         int nonZeros = static_cast<int>(_pR->nonZeros());
2569         int* pRows = new int[nonZeros * 2];
2570         _pR->outputRowCol(pRows);
2571         int* pCols = pRows + nonZeros;
2572
2573         if (bComplex1 || bComplex2)
2574         {
2575             for (int i = 0 ; i < nonZeros ; i++)
2576             {
2577                 int iRow = static_cast<int>(pRows[i]) - 1;
2578                 int iCol = static_cast<int>(pCols[i]) - 1;
2579                 std::complex<double> dbl = _pR->getImg(iRow, iCol);
2580                 pOut->set(iRow, iCol, pOut->get(iRow, iCol) - dbl.real());
2581                 pOut->setImg(iRow, iCol, pOut->getImg(iRow, iCol) - dbl.imag());
2582             }
2583         }
2584         else
2585         {
2586             for (int i = 0 ; i < nonZeros ; i++)
2587             {
2588                 int iRow = static_cast<int>(pRows[i]) - 1;
2589                 int iCol = static_cast<int>(pCols[i]) - 1;
2590                 pOut->set(iRow, iCol, pOut->get(iRow, iCol) - _pR->get(iRow, iCol));
2591             }
2592         }
2593
2594         //clear
2595         delete[] pRows;
2596         return pOut;
2597     }
2598     else
2599     {
2600         throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
2601     }
2602 }
2603
2604 //sp - d
2605 template<> InternalType* sub_M_M<Sparse, Double, Double>(Sparse* _pL, Double* _pR)
2606 {
2607     Double* pOut = NULL;
2608     int iOne = 1; //fortran
2609     bool bComplex1 = _pL->isComplex();
2610     bool bComplex2 = _pR->isComplex();
2611     bool bComplexOut = bComplex1 || bComplex2;
2612
2613     if (_pL->isScalar() && _pR->isScalar())
2614     {
2615         //sp - d
2616         pOut = (Double*)_pR->clone();
2617         pOut->setComplex(bComplexOut);
2618         if (bComplex1)
2619         {
2620             std::complex<double> dbl = _pL->getImg(0, 0);
2621             pOut->set(0, dbl.real() - pOut->get(0));
2622             pOut->setImg(0, dbl.imag() - pOut->getImg(0));
2623         }
2624         else if (bComplex2)
2625         {
2626             pOut->set(0, _pL->get(0, 0) - pOut->get(0));
2627             pOut->setImg(0, - pOut->getImg(0));
2628         }
2629         else
2630         {
2631             pOut->set(0, _pL->get(0, 0) - pOut->get(0));
2632         }
2633
2634         return pOut;
2635     }
2636
2637     if (_pR->isScalar())
2638     {
2639         //SP - d
2640         pOut = new Double(_pL->getRows(), _pL->getCols(), bComplexOut);
2641         int iSize = _pL->getSize();
2642         double dblVal = -_pR->get(0);
2643         double dblValI = 0;
2644         C2F(dset)(&iSize, &dblVal, pOut->get(), &iOne);
2645         if (bComplex2)
2646         {
2647             dblValI = -_pR->getImg(0);
2648             C2F(dset)(&iSize, &dblValI, pOut->getImg(), &iOne);
2649         }
2650         else if (bComplex1)
2651         {
2652             //initialize imag part at 0
2653             dblValI = 0;
2654             C2F(dset)(&iSize, &dblValI, pOut->getImg(), &iOne);
2655         }
2656
2657         int nonZeros = static_cast<int>(_pL->nonZeros());
2658         int* pRows = new int[nonZeros * 2];
2659         _pL->outputRowCol(pRows);
2660         int* pCols = pRows + nonZeros;
2661
2662         if (bComplex1)
2663         {
2664             for (int i = 0 ; i < nonZeros ; i++)
2665             {
2666                 int iRow = static_cast<int>(pRows[i]) - 1;
2667                 int iCol = static_cast<int>(pCols[i]) - 1;
2668                 std::complex<double> dbl = _pL->getImg(iRow, iCol);
2669                 pOut->set(iRow, iCol, dbl.real() - (-dblVal));
2670                 pOut->setImg(iRow, iCol, dbl.imag() - (-dblValI));
2671             }
2672         }
2673         else
2674         {
2675             for (int i = 0 ; i < nonZeros ; i++)
2676             {
2677                 int iRow = static_cast<int>(pRows[i]) - 1;
2678                 int iCol = static_cast<int>(pCols[i]) - 1;
2679                 pOut->set(iRow, iCol, _pL->get(iRow, iCol) - (-dblVal));
2680             }
2681         }
2682
2683         //clear
2684         delete[] pRows;
2685
2686         return pOut;
2687     }
2688
2689     if (_pL->isScalar())
2690     {
2691         //sp - D
2692         pOut = (Double*)_pR->clone();
2693         pOut->setComplex(bComplexOut);
2694
2695         if (bComplex1)
2696         {
2697             double* pReal = pOut->get();
2698             double* pImg = pOut->getImg();
2699             int size = pOut->getSize();
2700             for (int i = 0 ; i < size ; i++)
2701             {
2702                 std::complex<double> dbl = _pL->getImg(0, 0);
2703                 pReal[i] = dbl.real() - pReal[i];
2704                 pImg[i] = dbl.imag() - pImg[i];
2705             }
2706         }
2707         else
2708         {
2709             double* pReal = pOut->get();
2710             int size = pOut->getSize();
2711             for (int i = 0 ; i < size ; i++)
2712             {
2713                 pReal[i] = _pL->get(0, 0) - pReal[i];
2714             }
2715         }
2716
2717         return pOut;
2718     }
2719
2720
2721     if (_pR->getDims() > 2)
2722     {
2723         return nullptr;
2724     }
2725
2726     if (_pL->getRows() != _pR->getRows() || _pL->getCols() != _pR->getCols())
2727     {
2728         throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
2729     }
2730
2731     //SP - D
2732     pOut = new types::Double(_pL->getRows(), _pL->getCols(), _pL->isComplex());
2733     _pL->fill(*pOut);
2734     int iSize = pOut->getSize();
2735     pOut->setComplex(bComplexOut);
2736     double* pdblOutR = pOut->get();
2737     double* pdblOutI = pOut->getImg(); //can be null for non complex output
2738     double* pdblRR = _pR->get();
2739     double* pdblRI = _pR->getImg(); //can be null for non complex output
2740
2741     if (bComplex1)
2742     {
2743         if (bComplex2)
2744         {
2745             sub(pdblOutR, pdblOutI, iSize, pdblRR, pdblRI, pdblOutR, pdblOutI);
2746         }
2747         else
2748         {
2749             sub(pdblOutR, pdblOutI, iSize, pdblRR, pdblOutR, pdblOutI);
2750         }
2751     }
2752     else
2753     {
2754         if (bComplex2)
2755         {
2756             sub(pdblOutR, iSize, pdblRR, pdblRI, pdblOutR, pdblOutI);
2757         }
2758         else
2759         {
2760             sub(pdblOutR, iSize, pdblRR, pdblOutR);
2761         }
2762     }
2763
2764     return pOut;
2765 }
2766
2767 //Identity - sp
2768 template<> InternalType* sub_M_M<Double, Sparse, Sparse>(Double* _pL, Sparse* _pR)
2769 {
2770     Sparse* pOut = NULL;
2771     if (_pL->isIdentity())
2772     {
2773         //convert to _pL
2774         Sparse* pS = new Sparse(_pR->getRows(), _pR->getCols(), _pL->isComplex());
2775         int size = std::min(_pR->getRows(), _pR->getCols());
2776         double dblLeftR = _pL->get(0);
2777
2778
2779         if (_pL->isComplex())
2780         {
2781             std::complex<double> complexLeft(dblLeftR, _pL->getImg(0));
2782             for (int i = 0 ; i < size ; i++)
2783             {
2784                 pS->set(i, i, complexLeft);
2785             }
2786         }
2787         else
2788         {
2789             for (int i = 0 ; i < size ; i++)
2790             {
2791                 pS->set(i, i, dblLeftR);
2792             }
2793         }
2794         pS->finalize();
2795
2796
2797         pOut = pS->substract(*_pR);
2798         delete pS;
2799         return pOut;
2800     }
2801     else
2802     {
2803         // Call overload if the matrix is not identity
2804         return NULL;
2805     }
2806 }
2807
2808 //sp - Identity
2809 template<> InternalType* sub_M_M<Sparse, Double, Sparse>(Sparse* _pL, Double* _pR)
2810 {
2811     Sparse* pOut = NULL;
2812     if (_pR->isIdentity())
2813     {
2814         //convert to _pL
2815         Sparse* pS = new Sparse(_pL->getRows(), _pL->getCols(), _pR->isComplex());
2816         int size = std::min(_pL->getRows(), _pL->getCols());
2817         double dblRightR = _pR->get(0);
2818
2819         if (_pR->isComplex())
2820         {
2821             std::complex<double> complexRight(dblRightR, _pR->getImg(0));
2822             for (int i = 0 ; i < size ; i++)
2823             {
2824                 pS->set(i, i, complexRight, false);
2825             }
2826         }
2827         else
2828         {
2829             for (int i = 0 ; i < size ; i++)
2830             {
2831                 pS->set(i, i, dblRightR, false);
2832             }
2833         }
2834         pS->finalize();
2835
2836
2837         pOut = _pL->substract(*pS);
2838         delete pS;
2839         return pOut;
2840     }
2841     else
2842     {
2843         // Call overload if the matrix is not identity
2844         return NULL;
2845     }
2846 }