dotdiv fixed
[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 <limits>
20
21 #include "generic_operations.hxx"
22 #include "configvariable.hxx"
23 #include "double.hxx"
24 #include "polynom.hxx"
25 #include "sparse.hxx"
26
27 extern"C"
28 {
29 #include "abs.h"
30 }
31
32 void fillDotDivFunction();
33
34 //define arrays on operation functions
35 typedef types::InternalType*(*dotdiv_function)(types::InternalType*, types::InternalType*);
36
37 #define DECLARE_DOTDIV_PROTO(x) template<class T, class U, class O> types::InternalType* x(T *_pL, U *_pR)
38
39 DECLARE_DOTDIV_PROTO(dotdiv_M_M);
40 DECLARE_DOTDIV_PROTO(dotdiv_M_MC);
41 DECLARE_DOTDIV_PROTO(dotdiv_M_S);
42 DECLARE_DOTDIV_PROTO(dotdiv_M_SC);
43 DECLARE_DOTDIV_PROTO(dotdiv_M_E);
44 DECLARE_DOTDIV_PROTO(dotdiv_M_I);
45 DECLARE_DOTDIV_PROTO(dotdiv_M_IC);
46
47 DECLARE_DOTDIV_PROTO(dotdiv_MC_M);
48 DECLARE_DOTDIV_PROTO(dotdiv_MC_MC);
49 DECLARE_DOTDIV_PROTO(dotdiv_MC_S);
50 DECLARE_DOTDIV_PROTO(dotdiv_MC_SC);
51 DECLARE_DOTDIV_PROTO(dotdiv_MC_I);
52 DECLARE_DOTDIV_PROTO(dotdiv_MC_IC);
53
54 DECLARE_DOTDIV_PROTO(dotdiv_S_M);
55 DECLARE_DOTDIV_PROTO(dotdiv_S_MC);
56 DECLARE_DOTDIV_PROTO(dotdiv_S_S);
57 DECLARE_DOTDIV_PROTO(dotdiv_S_SC);
58 DECLARE_DOTDIV_PROTO(dotdiv_S_I);
59 DECLARE_DOTDIV_PROTO(dotdiv_S_IC);
60
61 DECLARE_DOTDIV_PROTO(dotdiv_SC_M);
62 DECLARE_DOTDIV_PROTO(dotdiv_SC_MC);
63 DECLARE_DOTDIV_PROTO(dotdiv_SC_S);
64 DECLARE_DOTDIV_PROTO(dotdiv_SC_SC);
65 DECLARE_DOTDIV_PROTO(dotdiv_SC_I);
66 DECLARE_DOTDIV_PROTO(dotdiv_SC_IC);
67
68 //[]
69 DECLARE_DOTDIV_PROTO(dotdiv_E_M);
70
71 //eye
72 DECLARE_DOTDIV_PROTO(dotdiv_I_M);
73 DECLARE_DOTDIV_PROTO(dotdiv_I_MC);
74 DECLARE_DOTDIV_PROTO(dotdiv_I_S);
75 DECLARE_DOTDIV_PROTO(dotdiv_I_SC);
76 DECLARE_DOTDIV_PROTO(dotdiv_I_I);
77 DECLARE_DOTDIV_PROTO(dotdiv_I_IC);
78
79 DECLARE_DOTDIV_PROTO(dotdiv_IC_M);
80 DECLARE_DOTDIV_PROTO(dotdiv_IC_MC);
81 DECLARE_DOTDIV_PROTO(dotdiv_IC_S);
82 DECLARE_DOTDIV_PROTO(dotdiv_IC_SC);
83 DECLARE_DOTDIV_PROTO(dotdiv_IC_I);
84 DECLARE_DOTDIV_PROTO(dotdiv_IC_IC);
85
86 #undef DECLARE_DOTDIV_PROTO
87
88 template<> types::InternalType* dotdiv_M_M<types::Sparse, types::Sparse, types::Sparse>(types::Sparse* _pL, types::Sparse* _pR);
89 template<> types::InternalType* dotdiv_M_M<types::Double, types::Sparse, types::Double>(types::Double* _pL, types::Sparse* _pR);
90 template<> types::InternalType* dotdiv_M_M<types::Sparse, types::Double, types::Double>(types::Sparse* _pL, types::Double* _pR);
91 template<> types::InternalType* dotdiv_M_M<types::Double, types::Sparse, types::Sparse>(types::Double* _pL, types::Sparse* _pR);
92 template<> types::InternalType* dotdiv_M_M<types::Sparse, types::Double, types::Sparse>(types::Sparse* _pL, types::Double* _pR);
93
94 template<> types::InternalType* dotdiv_M_M<types::Polynom, types::Polynom, types::Polynom>(types::Polynom* _pL, types::Polynom* _pR);
95 template<> types::InternalType* dotdiv_M_M<types::Polynom, types::Double, types::Polynom>(types::Polynom* _pL, types::Double* _pR);
96 template<> types::InternalType* dotdiv_M_M<types::Double, types::Polynom, types::Polynom>(types::Double* _pL, types::Polynom* _pR);
97
98 //x1 ./ x1
99 template<typename T, typename U> inline static void dotdiv(T l, U r, double* o)
100 {
101     if (r == 0)
102     {
103         ConfigVariable::setDivideByZero(true);
104     }
105     *o = (double)l / (double)r;
106 }
107
108 template<typename T, typename U, typename O> inline static void dotdiv(T l, U r, O* o)
109 {
110     if ((O)r == 0)
111     {
112         ConfigVariable::setDivideByZero(true);
113         double tmp = (double)l / (double)r;
114         if (std::isnan(tmp))
115         {
116             *o = 0;
117         }
118         else if (std::isinf(tmp))
119         {
120             tmp < 0 ? *o = std::numeric_limits<O>::min() : *o = std::numeric_limits<O>::max();
121         }
122         // no else because division by zero result to nan or inf
123     }
124     else
125     {
126         *o = (O)l / (O)r;
127     }
128 }
129
130 //x1 ./ x1c
131 template<typename T, typename U, typename O> inline static void dotdiv(T l, size_t size, U r, U rc, O* o, O* oc)
132 {
133     dotdiv((O)l * (O)r, (O)rc * (O)rc + (O)r * (O)r, o);
134     dotdiv((O)l * -(O)rc, (O)rc * (O)rc + (O)r * (O)r, oc);
135     //*o  = ((O)l * (O)r) / ((O)rc * (O)rc + (O)r * (O)r) ;
136     //*oc = ((O)l * -(O)rc) / ((O)rc * (O)rc + (O)r * (O)r) ;
137 }
138
139 template<> inline void dotdiv<double, double, double>(double l, size_t size, double r, double rc, double* o, double* oc)
140 {
141     if (rc == 0)
142     {
143         dotdiv<double, double, double>(l, r, o);
144         *oc = 0;
145     }
146     else if (r == 0)
147     {
148         *o = 0;
149         dotdiv<double, double, double>(-l, rc, oc);
150     }
151     else
152     {
153         double dblAbsSum    = dabss(r) + dabss(rc);
154         double dblReal1Sum  = l  / dblAbsSum;
155         double dblReal2Sum  = r  / dblAbsSum;
156         double dblImg2Sum   = rc / dblAbsSum;
157         double dblSum       = pow(dblReal2Sum, 2) + pow(dblImg2Sum, 2);
158         *o                  = (dblReal1Sum * dblReal2Sum) / dblSum;
159         *oc                 = (-dblReal1Sum * dblImg2Sum) / dblSum;
160     }
161 }
162
163
164
165 //x1c ./ x1
166 template<typename T, typename U, typename O> inline static void dotdiv(T l, T lc, size_t size, U r, O* o, O* oc)
167 {
168     dotdiv(l, r, o);
169     dotdiv(lc, r, oc);
170     //*o = (O)l / (O)r;
171     //*oc = (O)lc / (O)r;
172 }
173
174 //x1c ./ x1c
175 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)
176 {
177     dotdiv((O)l * (O)r + (O)lc * (O)rc, (O)rc * (O)rc + (O)r * (O)r, o);
178     dotdiv((O)r * (O)lc - (O)rc * (O)l, (O)rc * (O)rc + (O)r * (O)r, oc);
179     //*o = ((O)l * (O)r + (O)lc * (O)rc) / ((O)rc * (O)rc + (O)r * (O)r);
180     //*oc = ((O)r * (O)lc - (O)rc * (O)l) / ((O)rc * (O)rc + (O)r * (O)r) ;
181 }
182
183 template<> inline void dotdiv<double, double, double>(double l, double lc, size_t size, double r, double rc, double* o, double* oc)
184 {
185     if (rc == 0)
186     {
187         if (r == 0)
188         {
189             //got NaN + i NaN
190             dotdiv<double, double, double>(rc, r, o);
191             *oc = *o;
192         }
193         else
194         {
195             *o  = l  / r;
196             *oc = lc / r;
197         }
198     }
199     else if (r == 0)
200     {
201         *o  = lc / rc;
202         *oc = -l / rc;
203     }
204     else
205     {
206         if (dabss(r) >= dabss(rc))
207         {
208             double oRatio = rc / r;
209             double oSum = r + oRatio * rc;
210             *o = (l + lc * oRatio) / oSum;
211             *oc = (lc - l * oRatio) / oSum;
212         }
213         else
214         {
215             double oRatio = r / rc;
216             double oSum = rc + oRatio * r;
217             *o = (l * oRatio + lc) / oSum;
218             *oc = (lc * oRatio - l) / oSum;
219         }
220     }
221 }
222
223
224 //x ./ x
225 template<typename T, typename U, typename O> inline static void dotdiv(T* l, size_t size, U* r, O* o)
226 {
227     for (size_t i = 0; i < size ; ++i)
228     {
229         //dotdiv(T l, U r, O* o)
230         dotdiv(l[i], r[i], &o[i]);
231         //o[i] = (O)l[i] / (O)r[i];
232     }
233 }
234
235 //xC ./ x
236 template<typename T, typename U, typename O> inline static void dotdiv(T* l, T* lc, size_t size, U* r, 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, O* o, O* oc)
241         dotdiv(l[i], lc[i], (size_t)1, r[i], &o[i], &oc[i]);
242         //o[i] = (O)l[i] / (O)r[i];
243         //oc[i] = (O)lc[i] / (O)r[i];
244     }
245 }
246
247 //x ./ xC
248 template<typename T, typename U, typename O> inline static void dotdiv(T* l, size_t size, U* r, U* rc, O* o, O* oc)
249 {
250     for (size_t i = 0; i < size ; ++i)
251     {
252         //dotdiv(T l, size_t size, U r, U rc, O* o, O* oc)
253         dotdiv(l[i], 1, r[i], rc[i], &o[i], &oc[i]);
254         //o[i] = ((O)l[i] * (O)r[i]) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]) ;
255         //oc[i] = ((O)l[i] * -(O)rc[i]) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]) ;
256     }
257 }
258
259 //xC ./ xC
260 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)
261 {
262     for (size_t i = 0; i < size ; ++i)
263     {
264         //dotdiv(T l, T lc, size_t size, U r, U rc, O* o, O* oc)
265         dotdiv(l[i], lc[i], 1, r[i], rc[i], &o[i], &oc[i]);
266         //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]) ;
267         //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]) ;
268     }
269 }
270
271 //x ./ x1
272 template<typename T, typename U, typename O> inline static void dotdiv(T* l, size_t size, U r, O* o)
273 {
274     for (size_t i = 0; i < size ; ++i)
275     {
276         //dotdiv(T l, U r, O* o)
277         dotdiv(l[i], r, &o[i]);
278         //o[i] = (O)l[i] / (O)r;
279     }
280 }
281
282 //x1 ./ x
283 template<typename T, typename U, typename O> inline static void dotdiv(T l, size_t size, U* r, O* o)
284 {
285     for (size_t i = 0; i < size ; ++i)
286     {
287         //dotdiv(T l, U r, O* o)
288         dotdiv(l, r[i], &o[i]);
289         //o[i] = (O)l / (O)r[i];
290     }
291 }
292
293
294 //x ./ x1c
295 template<typename T, typename U, typename O> inline static void dotdiv(T* l, size_t size, U r, U rc, O* o, O* oc)
296 {
297     //O denum = ((O)rc * (O)rc + (O)r * (O)r);
298     for (size_t i = 0; i < size ; ++i)
299     {
300         //dotdiv(T l, size_t size, U r, U rc, O* o, O* oc)
301         dotdiv(l[i], (size_t)1, r, rc, &o[i], &oc[i]);
302         //o[i] = ((O)l[i] * (O)r) / denum;
303         //oc[i] = ((O)l[i] * -(O)rc) / denum;
304     }
305 }
306
307 //x1 ./ xc
308 template<typename T, typename U, typename O> inline static void dotdiv(T l, size_t size, U* r, U* rc, O* o, O* oc)
309 {
310     for (size_t i = 0; i < size ; ++i)
311     {
312         //dotdiv(T l, size_t size, U r, U rc, O* o, O* oc)
313         dotdiv(l, (size_t)1, r[i], rc[i], &o[i], &oc[i]);
314         //o[i] = ((O)l * (O)r[i]) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]) ;
315         //oc[i] = ((O)l * -(O)rc[i]) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]) ;
316     }
317 }
318
319 //xC ./ x1
320 template<typename T, typename U, typename O> inline static void dotdiv(T* l, T* lc, size_t size, U r, O* o, O* oc)
321 {
322     for (size_t i = 0; i < size ; ++i)
323     {
324         //dotdiv(T l, T lc, size_t size, U r, O* o, O* oc)
325         dotdiv(l[i], lc[i], (size_t)1, r, &o[i], &oc[i]);
326         //o[i] = (O)l[i] / (O)r;
327         //oc[i] = (O)lc[i] / (O)r;
328     }
329 }
330
331 //x1C ./ x
332 template<typename T, typename U, typename O> inline static void dotdiv(T l, T lc, size_t size, U* r, O* o, O* oc)
333 {
334     for (size_t i = 0; i < size ; ++i)
335     {
336         //dotdiv(T l, T lc, size_t size, U r, O* o, O* oc)
337         dotdiv(l, lc, (size_t)1, r[i], &o[i], &oc[i]);
338         //o[i] = (O)l / (O)r[i];
339         //oc[i] = (O)lc / (O)r[i];
340     }
341 }
342
343
344 //xC ./ x1c
345 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)
346 {
347     //    O denum = ((O)rc * (O)rc + (O)r * (O)r);
348     for (size_t i = 0; i < size ; ++i)
349     {
350         //dotdiv(T l, T lc, size_t size, U r, O* o, O* oc)
351         dotdiv(l[i], lc[i], (size_t)1, r, rc, &o[i], &oc[i]);
352         //o[i] = ((O)l[i] * (O)r + (O)lc[i] * (O)rc) / denum;
353         //oc[i] = ((O)r * (O)lc[i] - (O)rc * (O)l[i] ) / denum ;
354     }
355 }
356
357 //x1C ./ xc
358 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)
359 {
360     for (size_t i = 0; i < size ; ++i)
361     {
362         //dotdiv(T l, T lc, size_t size, U r, O* o, O* oc)
363         dotdiv(l, lc, (size_t)1, r[i], rc[i], &o[i], &oc[i]);
364         //o[i] = ((O)l * (O)r[i] + (O)lc * (O)rc[i]) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]);
365         //oc[i] = ((O)r[i] * (O)lc - (O)rc[i] * (O)l ) / ((O)rc[i] * (O)rc[i] + (O)r[i] * (O)r[i]) ;
366     }
367 }
368
369
370 #endif /* !__TYPES_DOTDIVIDE_HXX__ */