tdf#137679 Use Kahan summation for interpr5.cxx

Change-Id: Iefd09f2af4c77a464e47367c75b9e7e303f2c2a9
Reviewed-on: https://gerrit.libreoffice.org/c/core/+/115128
Tested-by: Jenkins
Reviewed-by: Mike Kaganski <mike.kaganski@collabora.com>

1 files changed, 25 insertions, 25 deletions

diff --git a/sc/source/core/tool/interpr5.cxx b/sc/source/core/tool/interpr5.cxx
index 4f10bb8961d2..f7dd527780ea 100644
--- a/sc/source/core/tool/interpr5.cxx
+++ b/sc/source/core/tool/interpr5.cxx
@@ -90,15 +90,14 @@ struct MatrixPow
 void lcl_MFastMult(const ScMatrixRef& pA, const ScMatrixRef& pB, const ScMatrixRef& pR,
                    SCSIZE n, SCSIZE m, SCSIZE l)
 {
-    double sum;
     for (SCSIZE row = 0; row < n; row++)
     {
         for (SCSIZE col = 0; col < l; col++)
         {   // result element(col, row) =sum[ (row of A) * (column of B)]
-            sum = 0.0;
+            KahanSum fSum = 0.0;
             for (SCSIZE k = 0; k < m; k++)
-                sum += pA->GetDouble(k,row) * pB->GetDouble(col,k);
-            pR->PutDouble(sum, col, row);
+                fSum += pA->GetDouble(k,row) * pB->GetDouble(col,k);
+            pR->PutDouble(fSum.get(), col, row);
         }
     }
 }
@@ -886,7 +885,7 @@ static void lcl_LUP_solve( const ScMatrix* mLU, const SCSIZE n,
     // Define y=Ux and solve for y in Ly=Pb using forward substitution.
     for (SCSIZE i=0; i < n; ++i)
     {
-        double fSum = B[P[i]];
+        KahanSum fSum = B[P[i]];
         // Matrix inversion comes with a lot of zeros in the B vectors, we
         // don't have to do all the computing with results multiplied by zero.
         // Until then, simply lookout for the position of the first nonzero
@@ -894,19 +893,19 @@ static void lcl_LUP_solve( const ScMatrix* mLU, const SCSIZE n,
         if (nFirst != SCSIZE_MAX)
         {
             for (SCSIZE j = nFirst; j < i; ++j)
-                fSum -= mLU->GetDouble( j, i) * X[j];   // X[j] === y[j]
+                fSum -= mLU->GetDouble( j, i) * X[j];         // X[j] === y[j]
         }
-        else if (fSum)
+        else if (fSum != 0)
             nFirst = i;
-        X[i] = fSum;                                    // X[i] === y[i]
+        X[i] = fSum.get();                                    // X[i] === y[i]
     }
     // Solve for x in Ux=y using back substitution.
     for (SCSIZE i = n; i--; )
     {
-        double fSum = X[i];                             // X[i] === y[i]
+        KahanSum fSum = X[i];                                 // X[i] === y[i]
         for (SCSIZE j = i+1; j < n; ++j)
-            fSum -= mLU->GetDouble( j, i) * X[j];       // X[j] === x[j]
-        X[i] = fSum / mLU->GetDouble( i, i);            // X[i] === x[i]
+            fSum -= mLU->GetDouble( j, i) * X[j];             // X[j] === x[j]
+        X[i] = fSum.get() / mLU->GetDouble( i, i);            // X[i] === x[i]
     }
 #ifdef DEBUG_SC_LUP_DECOMPOSITION
     fprintf( stderr, "\n%s\n", "lcl_LUP_solve():");
@@ -1098,17 +1097,17 @@ void ScInterpreter::ScMatMult()
                 pRMat = GetNewMat(nC2, nR1, /*bEmpty*/true);
                 if (pRMat)
                 {
-                    double sum;
+                    KahanSum fSum;
                     for (SCSIZE i = 0; i < nR1; i++)
                     {
                         for (SCSIZE j = 0; j < nC2; j++)
                         {
-                            sum = 0.0;
+                            fSum = 0.0;
                             for (SCSIZE k = 0; k < nC1; k++)
                             {
-                                sum += pMat1->GetDouble(k,i)*pMat2->GetDouble(j,k);
+                                fSum += pMat1->GetDouble(k,i)*pMat2->GetDouble(j,k);
                             }
-                            pRMat->PutDouble(sum, j, i);
+                            pRMat->PutDouble(fSum.get(), j, i);
                         }
                     }
                     PushMatrix(pRMat);
@@ -1804,7 +1803,8 @@ void ScInterpreter::CalculateSumX2MY2SumX2DY2(bool _bSumX2DY2)
         PushNoValue();
         return;
     }
-    double fVal, fSum = 0.0;
+    double fVal;
+    KahanSum fSum = 0.0;
     for (i = 0; i < nC1; i++)
         for (j = 0; j < nR1; j++)
             if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
@@ -1817,7 +1817,7 @@ void ScInterpreter::CalculateSumX2MY2SumX2DY2(bool _bSumX2DY2)
                 else
                     fSum -= fVal * fVal;
             }
-    PushDouble(fSum);
+    PushDouble(fSum.get());
 }
 
 void ScInterpreter::ScSumX2DY2()
@@ -2701,7 +2701,7 @@ void ScInterpreter::CalculateRGPRKP(bool _bRKP)
                     // = RMSE * sqrt( Xmean * (R' R)^(-1) * Xmean' + 1/N)
                     // (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as
                     // a whole matrix, but iterate over unit vectors.
-                    double fSigmaIntercept = 0.0;
+                    KahanSum aSigmaIntercept = 0.0;
                     double fPart; // for Xmean * single column of (R' R)^(-1)
                     for (SCSIZE col = 0; col < K; col++)
                     {
@@ -2719,13 +2719,13 @@ void ScInterpreter::CalculateRGPRKP(bool _bRKP)
                         if (bConstant)
                         {
                             fPart = lcl_GetSumProduct(pMeans, pMatZ, K);
-                            fSigmaIntercept += fPart * pMeans->GetDouble(col);
+                            aSigmaIntercept += fPart * pMeans->GetDouble(col);
                         }
                     }
                     if (bConstant)
                     {
-                        fSigmaIntercept = fRMSE
-                                          * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N));
+                        double fSigmaIntercept = fRMSE
+                                          * sqrt( (aSigmaIntercept + 1.0 / static_cast<double>(N) ).get() );
                         pResMat->PutDouble(fSigmaIntercept, K, 1);
                     }
                     else
@@ -2859,7 +2859,7 @@ void ScInterpreter::CalculateRGPRKP(bool _bRKP)
                     // (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as
                     // a whole matrix, but iterate over unit vectors.
                     // (R' R) ^(-1) is symmetric
-                    double fSigmaIntercept = 0.0;
+                    KahanSum aSigmaIntercept = 0.0;
                     double fPart; // for Xmean * single col of (R' R)^(-1)
                     for (SCSIZE row = 0; row < K; row++)
                     {
@@ -2876,13 +2876,13 @@ void ScInterpreter::CalculateRGPRKP(bool _bRKP)
                         if (bConstant)
                         {
                             fPart = lcl_GetSumProduct(pMeans, pMatZ, K);
-                            fSigmaIntercept += fPart * pMeans->GetDouble(row);
+                            aSigmaIntercept += fPart * pMeans->GetDouble(row);
                         }
                     }
                     if (bConstant)
                     {
-                        fSigmaIntercept = fRMSE
-                                          * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N));
+                        double fSigmaIntercept = fRMSE
+                                          * sqrt( (aSigmaIntercept + 1.0 / static_cast<double>(N) ).get() );
                         pResMat->PutDouble(fSigmaIntercept, K, 1);
                     }
                     else