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