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@collabora.com>
diff --git a/sc/inc/kahan.hxx b/sc/inc/kahan.hxx
index 23a29f6e..dffd74d 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 648f8b8..aaffa4a 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 67f803a..f72361a 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 c339a68..bc4265b 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()
