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 <dante19031...@gmail.com>
AuthorDate: Sat May 1 14:32:14 2021 +0200
Commit: Mike Kaganski <mike.kagan...@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.kagan...@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()
_______________________________________________
Libreoffice-commits mailing list
libreoffice-comm...@lists.freedesktop.org
https://lists.freedesktop.org/mailman/listinfo/libreoffice-commits
