[Libreoffice-commits] core.git: sc/inc sc/qa sc/source

dante (via logerrit) logerrit at kemper.freedesktop.org
Sat May 8 07:01:06 UTC 2021


 sc/inc/kahan.hxx                                           |   16 ++
 sc/qa/unit/data/functions/statistical/fods/chisq.test.fods |   16 +-
 sc/qa/unit/data/functions/statistical/fods/chitest.fods    |   16 +-
 sc/source/core/tool/interpr3.cxx                           |   80 ++++++-------
 4 files changed, 72 insertions(+), 56 deletions(-)

New commits:
commit 5545d8f515a2f7af45f0d2ab3ee40dfbde7fa20f
Author:     dante <dante19031999 at gmail.com>
AuthorDate: Sat May 1 14:32:14 2021 +0200
Commit:     Mike Kaganski <mike.kaganski at collabora.com>
CommitDate: Sat May 8 09:00:31 2021 +0200

    tdf#137679 Use Kahan summation for interpr3.cxx (2/2)
    
    Change-Id: I64113449f70d300feace6663c670db3448f43acf
    Reviewed-on: https://gerrit.libreoffice.org/c/core/+/114970
    Tested-by: Jenkins
    Reviewed-by: Mike Kaganski <mike.kaganski at collabora.com>

diff --git a/sc/inc/kahan.hxx b/sc/inc/kahan.hxx
index 23a29f6ee13f..dffd74d13f0f 100644
--- a/sc/inc/kahan.hxx
+++ b/sc/inc/kahan.hxx
@@ -133,6 +133,15 @@ public:
         m_fError /= fDivides;
     }
 
+    inline KahanSum operator*(const KahanSum& fTimes) const
+    {
+        KahanSum fSum(m_fSum * fTimes.m_fSum);
+        fSum += m_fSum * fTimes.m_fError;
+        fSum += m_fError * fTimes.m_fSum;
+        fSum += m_fError * fTimes.m_fError;
+        return fSum;
+    }
+
     constexpr KahanSum operator*(double fTimes) const
     {
         KahanSum fSum(*this);
@@ -140,6 +149,13 @@ public:
         return fSum;
     }
 
+    inline KahanSum operator/(const KahanSum& fDivides) const
+    {
+        KahanSum fSum(m_fSum / fDivides.get());
+        fSum += m_fError / fDivides.get();
+        return fSum;
+    }
+
     constexpr KahanSum operator/(double fTimes) const
     {
         KahanSum fSum(*this);
diff --git a/sc/qa/unit/data/functions/statistical/fods/chisq.test.fods b/sc/qa/unit/data/functions/statistical/fods/chisq.test.fods
index 648f8b864679..aaffa4a02a5e 100644
--- a/sc/qa/unit/data/functions/statistical/fods/chisq.test.fods
+++ b/sc/qa/unit/data/functions/statistical/fods/chisq.test.fods
@@ -3797,11 +3797,11 @@
      <table:table-cell table:style-name="ce46"/>
     </table:table-row>
     <table:table-row table:style-name="ro6">
-     <table:table-cell table:style-name="ce19" table:formula="of:=COM.MICROSOFT.CHISQ.TEST([.F1:.F11];[.G1:.G11])" office:value-type="float" office:value="1.16440189336067E-029" calcext:value-type="float">
-      <text:p>1.16440189336067E-29</text:p>
+     <table:table-cell table:style-name="ce19" table:formula="of:=COM.MICROSOFT.CHISQ.TEST([.F1:.F11];[.G1:.G11])" office:value-type="float" office:value="1.16440189336065E-029" calcext:value-type="float">
+      <text:p>1.16440189336065E-29</text:p>
      </table:table-cell>
-     <table:table-cell table:style-name="ce19" office:value-type="float" office:value="1.16440189336067E-029" calcext:value-type="float">
-      <text:p>1.16440189336067E-29</text:p>
+     <table:table-cell table:style-name="ce19" office:value-type="float" office:value="1.16440189336065E-029" calcext:value-type="float">
+      <text:p>1.16440189336065E-29</text:p>
      </table:table-cell>
      <table:table-cell table:style-name="ce54" table:formula="of:=ORG.LIBREOFFICE.ROUNDSIG([.A2];15)=ORG.LIBREOFFICE.ROUNDSIG([.B2];15)" office:value-type="boolean" office:boolean-value="true" calcext:value-type="boolean">
       <text:p>TRUE</text:p>
@@ -3860,11 +3860,11 @@
      <table:table-cell table:style-name="ce46"/>
     </table:table-row>
     <table:table-row table:style-name="ro6">
-     <table:table-cell table:style-name="ce66" table:formula="of:=COM.MICROSOFT.CHISQ.TEST([.F1:.F11];[.H1:.H11])" office:value-type="float" office:value="0.000000000903490480352966" calcext:value-type="float">
-      <text:p>9.03490480352966E-010</text:p>
+     <table:table-cell table:style-name="ce66" table:formula="of:=COM.MICROSOFT.CHISQ.TEST([.F1:.F11];[.H1:.H11])" office:value-type="float" office:value="0.000000000903490480352969" calcext:value-type="float">
+      <text:p>9.03490480352969E-010</text:p>
      </table:table-cell>
-     <table:table-cell table:style-name="ce19" office:value-type="float" office:value="0.000000000903490480352966" calcext:value-type="float">
-      <text:p>9.03490480352966E-10</text:p>
+     <table:table-cell table:style-name="ce19" office:value-type="float" office:value="0.000000000903490480352969" calcext:value-type="float">
+      <text:p>9.03490480352969E-10</text:p>
      </table:table-cell>
      <table:table-cell table:style-name="ce54" table:formula="of:=ORG.LIBREOFFICE.ROUNDSIG([.A3];15)=ORG.LIBREOFFICE.ROUNDSIG([.B3];15)" office:value-type="boolean" office:boolean-value="true" calcext:value-type="boolean">
       <text:p>TRUE</text:p>
diff --git a/sc/qa/unit/data/functions/statistical/fods/chitest.fods b/sc/qa/unit/data/functions/statistical/fods/chitest.fods
index 67f803af72ee..f72361ae35f8 100644
--- a/sc/qa/unit/data/functions/statistical/fods/chitest.fods
+++ b/sc/qa/unit/data/functions/statistical/fods/chitest.fods
@@ -3797,11 +3797,11 @@
      <table:table-cell table:style-name="ce46"/>
     </table:table-row>
     <table:table-row table:style-name="ro6">
-     <table:table-cell table:style-name="ce19" table:formula="of:=LEGACY.CHITEST([.F1:.F11];[.G1:.G11])" office:value-type="float" office:value="1.16440189336067E-029" calcext:value-type="float">
-      <text:p>1.16440189336067E-29</text:p>
+     <table:table-cell table:style-name="ce19" table:formula="of:=LEGACY.CHITEST([.F1:.F11];[.G1:.G11])" office:value-type="float" office:value="1.16440189336065E-029" calcext:value-type="float">
+      <text:p>1.16440189336065E-29</text:p>
      </table:table-cell>
-     <table:table-cell table:style-name="ce19" office:value-type="float" office:value="1.16440189336067E-029" calcext:value-type="float">
-      <text:p>1.16440189336067E-29</text:p>
+     <table:table-cell table:style-name="ce19" office:value-type="float" office:value="1.16440189336065E-029" calcext:value-type="float">
+      <text:p>1.16440189336065E-29</text:p>
      </table:table-cell>
      <table:table-cell table:style-name="ce54" table:formula="of:=ORG.LIBREOFFICE.ROUNDSIG([.A2];15)=ORG.LIBREOFFICE.ROUNDSIG([.B2];15)" office:value-type="boolean" office:boolean-value="true" calcext:value-type="boolean">
       <text:p>TRUE</text:p>
@@ -3860,11 +3860,11 @@
      <table:table-cell table:style-name="ce46"/>
     </table:table-row>
     <table:table-row table:style-name="ro6">
-     <table:table-cell table:style-name="ce66" table:formula="of:=LEGACY.CHITEST([.F1:.F11];[.H1:.H11])" office:value-type="float" office:value="0.000000000903490480352966" calcext:value-type="float">
-      <text:p>9.03490480352966E-010</text:p>
+     <table:table-cell table:style-name="ce66" table:formula="of:=LEGACY.CHITEST([.F1:.F11];[.H1:.H11])" office:value-type="float" office:value="0.000000000903490480352969" calcext:value-type="float">
+      <text:p>9.03490480352969E-010</text:p>
      </table:table-cell>
-     <table:table-cell table:style-name="ce19" office:value-type="float" office:value="0.000000000903490480352966" calcext:value-type="float">
-      <text:p>9.03490480352966E-10</text:p>
+     <table:table-cell table:style-name="ce19" office:value-type="float" office:value="0.000000000903490480352969" calcext:value-type="float">
+      <text:p>9.03490480352969E-10</text:p>
      </table:table-cell>
      <table:table-cell table:style-name="ce54" table:formula="of:=ORG.LIBREOFFICE.ROUNDSIG([.A3];15)=ORG.LIBREOFFICE.ROUNDSIG([.B3];15)" office:value-type="boolean" office:boolean-value="true" calcext:value-type="boolean">
       <text:p>TRUE</text:p>
diff --git a/sc/source/core/tool/interpr3.cxx b/sc/source/core/tool/interpr3.cxx
index c339a68dd80f..bc4265be0b67 100644
--- a/sc/source/core/tool/interpr3.cxx
+++ b/sc/source/core/tool/interpr3.cxx
@@ -1263,19 +1263,18 @@ static double lcl_GetBinomDistRange(double n, double xs,double xe,
 //preconditions: 0.0 <= xs < xe <= n; xs,xe,n integral although double
 {
     sal_uInt32 i;
-    double fSum;
     // skip summands index 0 to xs-1, start sum with index xs
     sal_uInt32 nXs = static_cast<sal_uInt32>( xs );
     for (i = 1; i <= nXs && fFactor > 0.0; i++)
         fFactor *= (n-i+1)/i * p/q;
-    fSum = fFactor; // Summand xs
+    KahanSum fSum = fFactor; // Summand xs
     sal_uInt32 nXe = static_cast<sal_uInt32>(xe);
     for (i = nXs+1; i <= nXe && fFactor > 0.0; i++)
     {
         fFactor *= (n-i+1)/i * p/q;
         fSum += fFactor;
     }
-    return std::min(fSum,1.0);
+    return std::min(fSum.get(), 1.0);
 }
 
 void ScInterpreter::ScB()
@@ -1433,7 +1432,7 @@ void ScInterpreter::ScCritBinom()
             fFactor = pow(q,n);
             if (fFactor > ::std::numeric_limits<double>::min())
             {
-                double fSum = fFactor;
+                KahanSum fSum = fFactor;
                 sal_uInt32 max = static_cast<sal_uInt32> (n), i;
                 for (i = 0; i < max && fSum < alpha; i++)
                 {
@@ -1445,15 +1444,13 @@ void ScInterpreter::ScCritBinom()
             else
             {
                 // accumulate BinomDist until accumulated BinomDist reaches alpha
-                double fSum = 0.0;
+                KahanSum fSum = 0.0;
                 sal_uInt32 max = static_cast<sal_uInt32> (n), i;
                 for (i = 0; i < max && fSum < alpha; i++)
                 {
                     const double x = GetBetaDistPDF( p, ( i + 1 ), ( n - i + 1 ) )/( n + 1 );
                     if ( nGlobalError == FormulaError::NONE )
-                    {
                         fSum += x;
-                    }
                     else
                     {
                         PushNoValue();
@@ -1469,7 +1466,7 @@ void ScInterpreter::ScCritBinom()
             fFactor = pow(p, n);
             if (fFactor > ::std::numeric_limits<double>::min())
             {
-                double fSum = 1.0 - fFactor;
+                KahanSum fSum = 1.0 - fFactor;
                 sal_uInt32 max = static_cast<sal_uInt32> (n), i;
                 for (i = 0; i < max && fSum >= alpha; i++)
                 {
@@ -1481,16 +1478,14 @@ void ScInterpreter::ScCritBinom()
             else
             {
                 // accumulate BinomDist until accumulated BinomDist reaches alpha
-                double fSum = 0.0;
+                KahanSum fSum = 0.0;
                 sal_uInt32 max = static_cast<sal_uInt32> (n), i;
                 alpha = 1 - alpha;
                 for (i = 0; i < max && fSum < alpha; i++)
                 {
                     const double x = GetBetaDistPDF( q, ( i + 1 ), ( n - i + 1 ) )/( n + 1 );
                     if ( nGlobalError == FormulaError::NONE )
-                    {
                         fSum += x;
-                    }
                     else
                     {
                         PushNoValue();
@@ -1821,15 +1816,15 @@ void ScInterpreter::ScPoissonDist( bool bODFF )
                 PushDouble (1.0);
             else
             {
-                double fSummand = exp(-lambda);
-                double fSum = fSummand;
+                double fSummand = std::exp(-lambda);
+                KahanSum fSum = fSummand;
                 int nEnd = sal::static_int_cast<int>( x );
                 for (int i = 1; i <= nEnd; i++)
                 {
                     fSummand = (fSummand * lambda)/static_cast<double>(i);
                     fSum += fSummand;
                 }
-                PushDouble(fSum);
+                PushDouble(fSum.get());
             }
         }
     }
@@ -1877,7 +1872,7 @@ void ScInterpreter::ScHypGeomDist( int nMinParamCount )
         return;
     }
 
-    double fVal = 0.0;
+    KahanSum fVal = 0.0;
 
     for ( int i = ( bCumulative ? 0 : x ); i <= x && nGlobalError == FormulaError::NONE; i++ )
     {
@@ -1885,7 +1880,7 @@ void ScInterpreter::ScHypGeomDist( int nMinParamCount )
             fVal +=  GetHypGeomDist( i, n, M, N );
     }
 
-    PushDouble( fVal );
+    PushDouble( fVal.get() );
 }
 
 /** Calculates a value of the hypergeometric distribution.
@@ -2473,8 +2468,8 @@ void ScInterpreter::ScZTest()
     }
     x = GetDouble();
 
-    double fSum    = 0.0;
-    double fSumSqr = 0.0;
+    KahanSum fSum    = 0.0;
+    KahanSum fSumSqr = 0.0;
     double fVal;
     double rValCount = 0.0;
     switch (GetStackType())
@@ -2566,11 +2561,11 @@ void ScInterpreter::ScZTest()
         PushError( FormulaError::DivisionByZero);
     else
     {
-        double mue = fSum/rValCount;
+        double mue = fSum.get()/rValCount;
 
         if (nParamCount != 3)
         {
-            sigma = (fSumSqr - fSum*fSum/rValCount)/(rValCount-1.0);
+            sigma = (fSumSqr - fSum*fSum/rValCount).get()/(rValCount-1.0);
             if (sigma == 0.0)
             {
                 PushError(FormulaError::DivisionByZero);
@@ -2588,12 +2583,12 @@ bool ScInterpreter::CalculateTest(bool _bTemplin
                                   ,const ScMatrixRef& pMat1,const ScMatrixRef& pMat2
                                   ,double& fT,double& fF)
 {
-    double fCount1  = 0.0;
-    double fCount2  = 0.0;
-    double fSum1    = 0.0;
-    double fSumSqr1 = 0.0;
-    double fSum2    = 0.0;
-    double fSumSqr2 = 0.0;
+    double fCount1    = 0.0;
+    double fCount2    = 0.0;
+    KahanSum fSum1    = 0.0;
+    KahanSum fSumSqr1 = 0.0;
+    KahanSum fSum2    = 0.0;
+    KahanSum fSumSqr2 = 0.0;
     double fVal;
     SCSIZE i,j;
     for (i = 0; i < nC1; i++)
@@ -2625,14 +2620,14 @@ bool ScInterpreter::CalculateTest(bool _bTemplin
     } // if (fCount1 < 2.0 || fCount2 < 2.0)
     if ( _bTemplin )
     {
-        double fS1 = (fSumSqr1-fSum1*fSum1/fCount1)/(fCount1-1.0)/fCount1;
-        double fS2 = (fSumSqr2-fSum2*fSum2/fCount2)/(fCount2-1.0)/fCount2;
+        double fS1 = (fSumSqr1-fSum1*fSum1/fCount1).get() / (fCount1-1.0) / fCount1;
+        double fS2 = (fSumSqr2-fSum2*fSum2/fCount2).get() / (fCount2-1.0) / fCount2;
         if (fS1 + fS2 == 0.0)
         {
             PushNoValue();
             return false;
         }
-        fT = std::abs(fSum1/fCount1 - fSum2/fCount2)/sqrt(fS1+fS2);
+        fT = std::abs(( fSum1/fCount1 - fSum2/fCount2 ).get())/sqrt(fS1+fS2);
         double c = fS1/(fS1+fS2);
     //  GetTDist is calculated via GetBetaDist and also works with non-integral
     // degrees of freedom. The result matches Excel
@@ -2641,9 +2636,9 @@ bool ScInterpreter::CalculateTest(bool _bTemplin
     else
     {
         //  according to Bronstein-Semendjajew
-        double fS1 = (fSumSqr1 - fSum1*fSum1/fCount1) / (fCount1 - 1.0);    // Variance
-        double fS2 = (fSumSqr2 - fSum2*fSum2/fCount2) / (fCount2 - 1.0);
-        fT = std::abs( fSum1/fCount1 - fSum2/fCount2 ) /
+        double fS1 = (fSumSqr1 - fSum1*fSum1/fCount1).get() / (fCount1 - 1.0);    // Variance
+        double fS2 = (fSumSqr2 - fSum2*fSum2/fCount2).get() / (fCount2 - 1.0);
+        fT = std::abs( fSum1.get()/fCount1 - fSum2.get()/fCount2 ) /
              sqrt( (fCount1-1.0)*fS1 + (fCount2-1.0)*fS2 ) *
              sqrt( fCount1*fCount2*(fCount1+fCount2-2)/(fCount1+fCount2) );
         fF = fCount1 + fCount2 - 2;
@@ -2683,9 +2678,9 @@ void ScInterpreter::ScTTest()
             return;
         }
         double fCount   = 0.0;
-        double fSum1    = 0.0;
-        double fSum2    = 0.0;
-        double fSumSqrD = 0.0;
+        KahanSum fSum1    = 0.0;
+        KahanSum fSum2    = 0.0;
+        KahanSum fSumSqrD = 0.0;
         double fVal1, fVal2;
         for (i = 0; i < nC1; i++)
             for (j = 0; j < nR1; j++)
@@ -2705,14 +2700,14 @@ void ScInterpreter::ScTTest()
             PushNoValue();
             return;
         }
-        double fSumD = fSum1 - fSum2;
-        double fDivider = fCount*fSumSqrD - fSumD*fSumD;
+        KahanSum fSumD = fSum1 - fSum2;
+        double fDivider = ( fSumSqrD*fCount - fSumD*fSumD ).get();
         if ( fDivider == 0.0 )
         {
             PushError(FormulaError::DivisionByZero);
             return;
         }
-        fT = std::abs(fSumD) * sqrt((fCount-1.0) / fDivider);
+        fT = std::abs(fSumD.get()) * sqrt((fCount-1.0) / fDivider);
         fF = fCount - 1.0;
     }
     else if (fTyp == 2.0)
@@ -2818,7 +2813,7 @@ void ScInterpreter::ScChiTest()
         PushIllegalArgument();
         return;
     }
-    double fChi = 0.0;
+    KahanSum fChi = 0.0;
     bool bEmpty = true;
     for (SCSIZE i = 0; i < nC1; i++)
     {
@@ -2843,6 +2838,11 @@ void ScInterpreter::ScChiTest()
                     // not, instead the result of divide() then was 1e+308.
                     volatile double fTemp1 = (fValX - fValE) * (fValX - fValE);
                     double fTemp2 = fTemp1;
+                    if (std::isinf(fTemp2))
+                    {
+                        PushError(FormulaError::NoConvergence);
+                        return;
+                    }
                     fChi += sc::divide( fTemp2, fValE);
                 }
                 else
@@ -2871,7 +2871,7 @@ void ScInterpreter::ScChiTest()
     }
     else
         fDF = static_cast<double>(nC1-1)*static_cast<double>(nR1-1);
-    PushDouble(GetChiDist(fChi, fDF));
+    PushDouble(GetChiDist(fChi.get(), fDF));
 }
 
 void ScInterpreter::ScKurt()


More information about the Libreoffice-commits mailing list