Polynomials multiplication & addition fixed by copying Scilab 5.X algorithm. 48/16048/5
Vincent COUVERT [Thu, 26 Feb 2015 13:32:30 +0000 (14:32 +0100)]
To test: test_run("polynomials", "bug_12679");

Revert "fixe : polynomials : bug_12679"
This reverts commit 9fa5a0bde7255d26dbddbabc6a50478a64fec825.

Change-Id: I5a4fffd58000e62994b357720bf4e354c977629f

scilab/modules/ast/src/c/operations/matrix_addition.c
scilab/modules/ast/src/c/operations/matrix_multiplication.c
scilab/modules/polynomials/tests/nonreg_tests/bug_12679.dia.ref
scilab/modules/polynomials/tests/nonreg_tests/bug_12679.tst

index 1a59d41..f6acea5 100644 (file)
 */
 
 #include <stdio.h>
+#include <string.h>
 #include "matrix_addition.h"
 #include "core_math.h"
 #include "operations_tools.h"
-#include <string.h>
+#include "elem_common.h"
 
 /* SinglePoly + SinglePoly*/
 int iAddScilabPolynomToScilabPolynom(double* _pCoef1R, int _iRank1, double* _pCoef2R, int _iRank2, double* _pCoefOutR, int _iRanOut)
 {
     int iRankMin                       = Min(_iRank1, _iRank2);
     int iRankMax                       = Max(_iRank1, _iRank2);
-    int iRank                                  = 0;
+    int iRank                          = 0;
+    double dblSum                      = 0.0;
     double* pCoefMaxR  = (_iRank1 > _iRank2) ? _pCoef1R : _pCoef2R;
 
     for (iRank = 0; iRank < iRankMin ; iRank++)
     {
-        _pCoefOutR[iRank] = _pCoef1R[iRank] + _pCoef2R[iRank];
+        dblSum = _pCoef1R[iRank] + _pCoef2R[iRank];
+        if (fabs(dblSum) > Max(fabs(_pCoef1R[iRank]), fabs(_pCoef2R[iRank])) * 2 * getRelativeMachinePrecision())
+        {
+            _pCoefOutR[iRank] = dblSum;
+        }
+        else
+        {
+            _pCoefOutR[iRank] = 0.0;
+        }
     }
 
     for (iRank = iRankMin ; iRank < iRankMax ; iRank++)
index a19b85c..ce45a9d 100644 (file)
@@ -156,13 +156,25 @@ int iMultiScilabPolynomByScilabPolynom(
 {
     int i1             = 0;
     int i2             = 0;
+    double dblMult = 0.0;
+    double dblAdd = 0.0;
 
     memset(_pdblRealOut, 0x00, _iRankOut * sizeof(double));
+
     for (i1 = 0 ; i1 < _iRank1 ; i1++)
     {
         for (i2 = 0 ; i2 < _iRank2 ; i2++)
         {
-            _pdblRealOut[i1 + i2] += _pdblReal1[i1] * _pdblReal2[i2];
+            dblMult = _pdblReal1[i1] * _pdblReal2[i2];
+            dblAdd = _pdblRealOut[i1 + i2] + dblMult;
+            if (fabs(dblAdd) > 2 * getRelativeMachinePrecision() * Max(fabs(_pdblRealOut[i1 + i2]), fabs(dblMult)))
+            {
+                _pdblRealOut[i1 + i2] = dblAdd;
+            }
+            else
+            {
+                _pdblRealOut[i1 + i2] = 0.0;
+            }
         }
     }
     return 0;
index a9c37f0..1a295ae 100644 (file)
@@ -20,8 +20,7 @@
 s = poly(0, "s");
 p = [s, s*(s+1)^2, 2*s^2+s^3];
 [pgcd, u] = gcd(p);
-pout = clean(p*u);
-assert_checkequal(pout, [0 0 s]);
+assert_checkequal(p*u, [0 0 s]);
 // Complex polynomials should yield an error
 // Normal behavior, with integers
 V = [2^2*3^5 2^3*3^2 2^2*3^4*5];
index 8254bdf..88f98dd 100644 (file)
@@ -22,8 +22,7 @@
 s = poly(0, "s");
 p = [s, s*(s+1)^2, 2*s^2+s^3];
 [pgcd, u] = gcd(p);
-pout = clean(p*u);
-assert_checkequal(pout, [0 0 s]);
+assert_checkequal(p*u, [0 0 s]);
 // Complex polynomials should yield an error
 
 // Normal behavior, with integers