* Bug 15599 fixed: now degree of zero polynomial is -Inf
[scilab.git] / scilab / modules / ast / includes / operations / types_dotdivide.hxx
1 /*
2 *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  *  Copyright (C) 2014 - Scilab Enterprises - Antoine ELIAS
4  *  Copyright (C) 2014 - Scilab Enterprises - Sylvain GENIN
5 *
6  * Copyright (C) 2012 - 2016 - Scilab Enterprises
7  *
8  * This file is hereby licensed under the terms of the GNU GPL v2.0,
9  * pursuant to article 5.3.4 of the CeCILL v.2.1.
10  * This file was originally licensed under the terms of the CeCILL v2.1,
11  * and continues to be available under such terms.
12  * For more information, see the COPYING file which you should have received
13  * along with this program.
14 *
15 */
16 #ifndef __TYPES_DOTDIVIDE_HXX__
17 #define __TYPES_DOTDIVIDE_HXX__
18
19 #include "generic_operations.hxx"
20 #include "configvariable.hxx"
21 #include "double.hxx"
22 #include "polynom.hxx"
23 #include "sparse.hxx"
24
25 extern"C"
26 {
27 #include "abs.h"
28 }
29
30 void fillDotDivFunction();
31
32 //define arrays on operation functions
33 typedef types::InternalType*(*dotdiv_function)(types::InternalType*, types::InternalType*);
34
35 #define DECLARE_DOTDIV_PROTO(x) template<class T, class U, class O> types::InternalType* x(T *_pL, U *_pR)
36
37 DECLARE_DOTDIV_PROTO(dotdiv_M_M);
38 DECLARE_DOTDIV_PROTO(dotdiv_M_MC);
39 DECLARE_DOTDIV_PROTO(dotdiv_M_S);
40 DECLARE_DOTDIV_PROTO(dotdiv_M_SC);
41 DECLARE_DOTDIV_PROTO(dotdiv_M_E);
42 DECLARE_DOTDIV_PROTO(dotdiv_M_I);
43 DECLARE_DOTDIV_PROTO(dotdiv_M_IC);
44
45 DECLARE_DOTDIV_PROTO(dotdiv_MC_M);
46 DECLARE_DOTDIV_PROTO(dotdiv_MC_MC);
47 DECLARE_DOTDIV_PROTO(dotdiv_MC_S);
48 DECLARE_DOTDIV_PROTO(dotdiv_MC_SC);
49 DECLARE_DOTDIV_PROTO(dotdiv_MC_I);
50 DECLARE_DOTDIV_PROTO(dotdiv_MC_IC);
51
52 DECLARE_DOTDIV_PROTO(dotdiv_S_M);
53 DECLARE_DOTDIV_PROTO(dotdiv_S_MC);
54 DECLARE_DOTDIV_PROTO(dotdiv_S_S);
55 DECLARE_DOTDIV_PROTO(dotdiv_S_SC);
56 DECLARE_DOTDIV_PROTO(dotdiv_S_I);
57 DECLARE_DOTDIV_PROTO(dotdiv_S_IC);
58
59 DECLARE_DOTDIV_PROTO(dotdiv_SC_M);
60 DECLARE_DOTDIV_PROTO(dotdiv_SC_MC);
61 DECLARE_DOTDIV_PROTO(dotdiv_SC_S);
62 DECLARE_DOTDIV_PROTO(dotdiv_SC_SC);
63 DECLARE_DOTDIV_PROTO(dotdiv_SC_I);
64 DECLARE_DOTDIV_PROTO(dotdiv_SC_IC);
65
66 //[]
67 DECLARE_DOTDIV_PROTO(dotdiv_E_M);
68
69 //eye
70 DECLARE_DOTDIV_PROTO(dotdiv_I_M);
71 DECLARE_DOTDIV_PROTO(dotdiv_I_MC);
72 DECLARE_DOTDIV_PROTO(dotdiv_I_S);
73 DECLARE_DOTDIV_PROTO(dotdiv_I_SC);
74 DECLARE_DOTDIV_PROTO(dotdiv_I_I);
75 DECLARE_DOTDIV_PROTO(dotdiv_I_IC);
76
77 DECLARE_DOTDIV_PROTO(dotdiv_IC_M);
78 DECLARE_DOTDIV_PROTO(dotdiv_IC_MC);
79 DECLARE_DOTDIV_PROTO(dotdiv_IC_S);
80 DECLARE_DOTDIV_PROTO(dotdiv_IC_SC);
81 DECLARE_DOTDIV_PROTO(dotdiv_IC_I);
82 DECLARE_DOTDIV_PROTO(dotdiv_IC_IC);
83
84 #undef DECLARE_DOTDIV_PROTO
85
86 template<> types::InternalType* dotdiv_M_M<types::Sparse, types::Sparse, types::Sparse>(types::Sparse* _pL, types::Sparse* _pR);
87 template<> types::InternalType* dotdiv_M_M<types::Double, types::Sparse, types::Double>(types::Double* _pL, types::Sparse* _pR);
88 template<> types::InternalType* dotdiv_M_M<types::Sparse, types::Double, types::Double>(types::Sparse* _pL, types::Double* _pR);
89 template<> types::InternalType* dotdiv_M_M<types::Double, types::Sparse, types::Sparse>(types::Double* _pL, types::Sparse* _pR);
90 template<> types::InternalType* dotdiv_M_M<types::Sparse, types::Double, types::Sparse>(types::Sparse* _pL, types::Double* _pR);
91
92 template<> types::InternalType* dotdiv_M_M<types::Polynom, types::Polynom, types::Polynom>(types::Polynom* _pL, types::Polynom* _pR);
93 template<> types::InternalType* dotdiv_M_M<types::Polynom, types::Double, types::Polynom>(types::Polynom* _pL, types::Double* _pR);
94 template<> types::InternalType* dotdiv_M_M<types::Double, types::Polynom, types::Polynom>(types::Double* _pL, types::Polynom* _pR);
95
96 //x1 ./ x1
97 template<typename T, typename U, typename O> inline static void dotdiv(T l, U r, O* o)
98 {
99     if ((O)r == 0)
100     {
101         ConfigVariable::setDivideByZero(true);
102     }
103     *o = (O)l / (O)r;
104 }
105
106 //x1 ./ x1c
107 template<typename T, typename U, typename O> inline static void dotdiv(T l, size_t size, U r, U rc, O* o, O* oc)
108 {
109     dotdiv((O)l * (O)r, (O)rc * (O)rc + (O)r * (O)r, o);
110     dotdiv((O)l * -(O)rc, (O)rc * (O)rc + (O)r * (O)r, oc);
111     //*o  = ((O)l * (O)r) / ((O)rc * (O)rc + (O)r * (O)r) ;
112     //*oc = ((O)l * -(O)rc) / ((O)rc * (O)rc + (O)r * (O)r) ;
113 }
114
115 template<> inline void dotdiv<double, double, double>(double l, size_t size, double r, double rc, double* o, double* oc)
116 {
117     if (rc == 0)
118     {
119         dotdiv<double, double, double>(l, r, o);
120         *oc = 0;
121     }
122     else if (r == 0)
123     {
124         *o = 0;
125         dotdiv<double, double, double>(-l, rc, oc);
126     }
127     else
128     {
129         double dblAbsSum    = dabss(r) + dabss(rc);
130         double dblReal1Sum  = l  / dblAbsSum;
131         double dblReal2Sum  = r  / dblAbsSum;
132         double dblImg2Sum   = rc / dblAbsSum;
133         double dblSum       = pow(dblReal2Sum, 2) + pow(dblImg2Sum, 2);
134         *o                  = (dblReal1Sum * dblReal2Sum) / dblSum;
135         *oc                 = (-dblReal1Sum * dblImg2Sum) / dblSum;
136     }
137 }
138
139
140
141 //x1c ./ x1
142 template<typename T, typename U, typename O> inline static void dotdiv(T l, T lc, size_t size, U r, O* o, O* oc)
143 {
144     dotdiv<T, U, O>(l, r, o);
145     dotdiv<T, U, O>(lc, r, oc);
146     //*o = (O)l / (O)r;
147     //*oc = (O)lc / (O)r;
148 }
149
150 //x1c ./ x1c
151 template<typename T, typename U, typename O> inline static void dotdiv(T l, T lc, size_t size, U r, U rc, O* o, O* oc)
152 {
153     dotdiv((O)l * (O)r + (O)lc * (O)rc, (O)rc * (O)rc + (O)r * (O)r, o);
154     dotdiv((O)r * (O)lc - (O)rc * (O)l, (O)rc * (O)rc + (O)r * (O)r, oc);
155     //*o = ((O)l * (O)r + (O)lc * (O)rc) / ((O)rc * (O)rc + (O)r * (O)r);
156     //*oc = ((O)r * (O)lc - (O)rc * (O)l) / ((O)rc * (O)rc + (O)r * (O)r) ;
157 }
158
159 template<> inline void dotdiv<double, double, double>(double l, double lc, size_t size, double r, double rc, double* o, double* oc)
160 {
161     if (rc == 0)
162     {
163         if (r == 0)
164         {
165             //got NaN + i NaN
166             dotdiv<double, double, double>(rc, r, o);
167             *oc = *o;
168         }
169         else
170         {
171             *o  = l  / r;
172             *oc = lc / r;
173         }
174     }
175     else if (r == 0)
176     {
177         *o  = lc / rc;
178         *oc = -l / rc;
179     }
180     else
181     {
182         if (dabss(r) >= dabss(rc))
183         {
184             double oRatio = rc / r;
185             double oSum = r + oRatio * rc;
186             *o = (l + lc * oRatio) / oSum;
187             *oc = (lc - l * oRatio) / oSum;
188         }
189         else
190         {
191             double oRatio = r / rc;
192             double oSum = rc + oRatio * r;
193             *o = (l * oRatio + lc) / oSum;
194             *oc = (lc * oRatio - l) / oSum;
195         }
196     }
197 }
198
199
200 //x ./ x
201 template<typename T, typename U, typename O> inline static void dotdiv(T* l, size_t size, U* r, O* o)
202 {
203     for (size_t i = 0; i < size ; ++i)
204     {
205         //dotdiv(T l, U r, O* o)
206         dotdiv<T, U, O>(l[i], r[i], &o[i]);
207         //o[i] = (O)l[i] / (O)r[i];
208     }
209 }
210
211 //xC ./ x
212 template<typename T, typename U, typename O> inline static void dotdiv(T* l, T* lc, size_t size, U* r, O* o, O* oc)
213 {
214     for (size_t i = 0; i < size ; ++i)
215     {
216         //dotdiv(T l, T lc, size_t size, U r, O* o, O* oc)
217         dotdiv<T, U, O>(l[i], lc[i], (size_t)1, r[i], &o[i], &oc[i]);
218         //o[i] = (O)l[i] / (O)r[i];
219         //oc[i] = (O)lc[i] / (O)r[i];
220     }
221 }
222
223 //x ./ xC
224 template<typename T, typename U, typename O> inline static void dotdiv(T* l, size_t size, U* r, U* rc, O* o, O* oc)
225 {
226     for (size_t i = 0; i < size ; ++i)
227     {
228         //dotdiv(T l, size_t size, U r, U rc, O* o, O* oc)
229         dotdiv<T, U, O>(l[i], 1, r[i], rc[i], &o[i], &oc[i]);
230         //o[i] = ((O)l[i] * (O)r[i]) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]) ;
231         //oc[i] = ((O)l[i] * -(O)rc[i]) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]) ;
232     }
233 }
234
235 //xC ./ xC
236 template<typename T, typename U, typename O> inline static void dotdiv(T* l, T* lc, size_t size, U* r, U* rc, O* o, O* oc)
237 {
238     for (size_t i = 0; i < size ; ++i)
239     {
240         //dotdiv(T l, T lc, size_t size, U r, U rc, O* o, O* oc)
241         dotdiv<T, U, O>(l[i], lc[i], 1, r[i], rc[i], &o[i], &oc[i]);
242         //o[i] =  ((O)l[i] * (O)r[i] + (O)lc[i] * (O)rc[i] ) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]) ;
243         //oc[i] = ((O)r[i] * (O)lc[i] - (O)rc[i] * (O)l[i] ) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]) ;
244     }
245 }
246
247 //x ./ x1
248 template<typename T, typename U, typename O> inline static void dotdiv(T* l, size_t size, U r, O* o)
249 {
250     for (size_t i = 0; i < size ; ++i)
251     {
252         //dotdiv(T l, U r, O* o)
253         dotdiv<T, U, O>(l[i], r, &o[i]);
254         //o[i] = (O)l[i] / (O)r;
255     }
256 }
257
258 //x1 ./ x
259 template<typename T, typename U, typename O> inline static void dotdiv(T l, size_t size, U* r, O* o)
260 {
261     for (size_t i = 0; i < size ; ++i)
262     {
263         //dotdiv(T l, U r, O* o)
264         dotdiv<T, U, O>(l, r[i], &o[i]);
265         //o[i] = (O)l / (O)r[i];
266     }
267 }
268
269
270 //x ./ x1c
271 template<typename T, typename U, typename O> inline static void dotdiv(T* l, size_t size, U r, U rc, O* o, O* oc)
272 {
273     //O denum = ((O)rc * (O)rc + (O)r * (O)r);
274     for (size_t i = 0; i < size ; ++i)
275     {
276         //dotdiv(T l, size_t size, U r, U rc, O* o, O* oc)
277         dotdiv<T, U, O >(l[i], (size_t)1, r, rc, &o[i], &oc[i]);
278         //o[i] = ((O)l[i] * (O)r) / denum;
279         //oc[i] = ((O)l[i] * -(O)rc) / denum;
280     }
281 }
282
283 //x1 ./ xc
284 template<typename T, typename U, typename O> inline static void dotdiv(T l, size_t size, U* r, U* rc, O* o, O* oc)
285 {
286     for (size_t i = 0; i < size ; ++i)
287     {
288         //dotdiv(T l, size_t size, U r, U rc, O* o, O* oc)
289         dotdiv<T, U, O>(l, (size_t)1, r[i], rc[i], &o[i], &oc[i]);
290         //o[i] = ((O)l * (O)r[i]) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]) ;
291         //oc[i] = ((O)l * -(O)rc[i]) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]) ;
292     }
293 }
294
295 //xC ./ x1
296 template<typename T, typename U, typename O> inline static void dotdiv(T* l, T* lc, size_t size, U r, O* o, O* oc)
297 {
298     for (size_t i = 0; i < size ; ++i)
299     {
300         //dotdiv(T l, T lc, size_t size, U r, O* o, O* oc)
301         dotdiv<T, U, O>(l[i], lc[i], (size_t)1, r, &o[i], &oc[i]);
302         //o[i] = (O)l[i] / (O)r;
303         //oc[i] = (O)lc[i] / (O)r;
304     }
305 }
306
307 //x1C ./ x
308 template<typename T, typename U, typename O> inline static void dotdiv(T l, T lc, size_t size, U* r, O* o, O* oc)
309 {
310     for (size_t i = 0; i < size ; ++i)
311     {
312         //dotdiv(T l, T lc, size_t size, U r, O* o, O* oc)
313         dotdiv<T, U, O>(l, lc, (size_t)1, r[i], &o[i], &oc[i]);
314         //o[i] = (O)l / (O)r[i];
315         //oc[i] = (O)lc / (O)r[i];
316     }
317 }
318
319
320 //xC ./ x1c
321 template<typename T, typename U, typename O> inline static void dotdiv(T* l, T* lc, size_t size, U r, U rc, O* o, O* oc)
322 {
323     //    O denum = ((O)rc * (O)rc + (O)r * (O)r);
324     for (size_t i = 0; i < size ; ++i)
325     {
326         //dotdiv(T l, T lc, size_t size, U r, O* o, O* oc)
327         dotdiv<T, U, O>(l[i], lc[i], (size_t)1, r, rc, &o[i], &oc[i]);
328         //o[i] = ((O)l[i] * (O)r + (O)lc[i] * (O)rc) / denum;
329         //oc[i] = ((O)r * (O)lc[i] - (O)rc * (O)l[i] ) / denum ;
330     }
331 }
332
333 //x1C ./ xc
334 template<typename T, typename U, typename O> inline static void dotdiv(T l, T lc, size_t size, U* r, U* rc, O* o, O* oc)
335 {
336     for (size_t i = 0; i < size ; ++i)
337     {
338         //dotdiv(T l, T lc, size_t size, U r, O* o, O* oc)
339         dotdiv<T, U, O>(l, lc, (size_t)1, r[i], rc[i], &o[i], &oc[i]);
340         //o[i] = ((O)l * (O)r[i] + (O)lc * (O)rc[i]) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]);
341         //oc[i] = ((O)r[i] * (O)lc - (O)rc[i] * (O)l ) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]) ;
342     }
343 }
344
345
346 #endif /* !__TYPES_DOTDIVIDE_HXX__ */