From: Cedric Delamarre Date: Wed, 1 Apr 2015 15:57:44 +0000 (+0200) Subject: dot divide fixed about precision error X-Git-Tag: 6.0.0-alpha-1~391 X-Git-Url: http://gitweb.scilab.org/?p=scilab.git;a=commitdiff_plain;h=8b043eff9972ed827e6fcd574acc11c705f24fb7 dot divide fixed about precision error r = (0.5857864376269046324808 - 1.41421356237309492343*%i) ./ (3.4142135623730953675192 + 1.41421356237309492343*%i); real(r) test_run signal_processing bug_9851 Change-Id: I22864f09cecc8596e4f141e7130494474153be65 --- diff --git a/scilab/modules/ast/includes/operations/types_dotdivide.hxx b/scilab/modules/ast/includes/operations/types_dotdivide.hxx index d8d0f2e..976e6e9 100644 --- a/scilab/modules/ast/includes/operations/types_dotdivide.hxx +++ b/scilab/modules/ast/includes/operations/types_dotdivide.hxx @@ -19,6 +19,11 @@ #include "polynom.hxx" #include "sparse.hxx" +extern"C" +{ +#include "abs.h" +} + void fillDotDivFunction(); //define arrays on operation functions @@ -104,6 +109,32 @@ template inline static void dotdiv(T l, size //*oc = ((O)l * -(O)rc) / ((O)rc * (O)rc + (O)r * (O)r) ; } +template<> inline void dotdiv(double l, size_t size, double r, double rc, double* o, double* oc) +{ + if (rc == 0) + { + dotdiv(l, r, o); + *oc = 0; + } + else if (r == 0) + { + *o = 0; + dotdiv(-l, rc, oc); + } + else + { + double dblAbsSum = dabss(r) + dabss(rc); + double dblReal1Sum = l / dblAbsSum; + double dblReal2Sum = r / dblAbsSum; + double dblImg2Sum = rc / dblAbsSum; + double dblSum = pow(dblReal2Sum, 2) + pow(dblImg2Sum, 2); + *o = (dblReal1Sum * dblReal2Sum) / dblSum; + *oc = (-dblReal1Sum * dblImg2Sum) / dblSum; + } +} + + + //x1c ./ x1 template inline static void dotdiv(T l, T lc, size_t size, U r, O* o, O* oc) { @@ -122,6 +153,46 @@ template inline static void dotdiv(T l, T lc //*oc = ((O)r * (O)lc - (O)rc * (O)l) / ((O)rc * (O)rc + (O)r * (O)r) ; } +template<> inline void dotdiv(double l, double lc, size_t size, double r, double rc, double* o, double* oc) +{ + if (rc == 0) + { + if (r == 0) + { + //got NaN + i NaN + dotdiv(rc, r, o); + *oc = *o; + } + else + { + *o = l / r; + *oc = lc / r; + } + } + else if (r == 0) + { + *o = lc / rc; + *oc = -l / rc; + } + else + { + if (dabss(r) >= dabss(rc)) + { + double oRatio = rc / r; + double oSum = r + oRatio * rc; + *o = (l + lc * oRatio) / oSum; + *oc = (lc - l * oRatio) / oSum; + } + else + { + double oRatio = r / rc; + double oSum = rc + oRatio * r; + *o = (l * oRatio + lc) / oSum; + *oc = (lc * oRatio - l) / oSum; + } + } +} + //x ./ x template inline static void dotdiv(T* l, size_t size, U* r, O* o)