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