commit 4283fb9d4a6152643364bfe1f98ee1f36aabbb78
author dante <dante19031999@gmail.com> Wed Apr 28 19:37:21 2021 +0200
committer Mike Kaganski <mike.kaganski@collabora.com> Fri Apr 30 23:29:59 2021 +0200
tree 7be85479adc740afa41ac360903bb51b86987d46
parent 6a113a4f14808ac7f4bbdb4a5baff9383541d49a

tdf#137679 Use kahan summation for ScInterpreter::CalculateSkew

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

sc/source/core/inc/interpre.hxx

diff --git a/sc/source/core/inc/interpre.hxx b/sc/source/core/inc/interpre.hxx
index bdde25a..c7b9379 100644
--- a/sc/source/core/inc/interpre.hxx
+++ b/sc/source/core/inc/interpre.hxx
@@ -835,7 +835,7 @@ private:
    void ScSumX2DY2();
    void ScSumXMY2();
    void ScGrowth();
    bool CalculateSkew(double& fSum,double& fCount,double& vSum,std::vector<double>& values);
    bool CalculateSkew(KahanSum& fSum, double& fCount, std::vector<double>& values);
    void CalculateSkewOrSkewp( bool bSkewp );
    void CalculateSlopeIntercept(bool bSlope);
    void CalculateSmallLarge(bool bSmall);

sc/source/core/tool/interpr3.cxx

diff --git a/sc/source/core/tool/interpr3.cxx b/sc/source/core/tool/interpr3.cxx
index 7adee32..d41f4f9 100644
--- a/sc/source/core/tool/interpr3.cxx
+++ b/sc/source/core/tool/interpr3.cxx
@@ -2872,9 +2872,10 @@ void ScInterpreter::ScChiTest()

void ScInterpreter::ScKurt()
{
    double fSum,fCount,vSum;
    KahanSum fSum;
    double fCount;
    std::vector<double> values;
    if ( !CalculateSkew(fSum,fCount,vSum,values) )
    if ( !CalculateSkew(fSum, fCount, values) )
        return;

    // ODF 1.2 constraints: # of numbers >= 4
@@ -2885,31 +2886,30 @@ void ScInterpreter::ScKurt()
        return;
    }

    double fMean = fSum / fCount;

    KahanSum vSum;
    double fMean = fSum.get() / fCount;
    for (double v : values)
        vSum += (v - fMean) * (v - fMean);

    double fStdDev = sqrt(vSum / (fCount - 1.0));
    double xpower4 = 0.0;

    double fStdDev = sqrt(vSum.get() / (fCount - 1.0));
    if (fStdDev == 0.0)
    {
        PushError( FormulaError::DivisionByZero);
        return;
    }

    KahanSum xpower4 = 0.0;
    for (double v : values)
    {
        double dx = (v - fMean) / fStdDev;
        xpower4 = xpower4 + (dx * dx * dx * dx);
        xpower4 += (dx * dx) * (dx * dx);
    }

    double k_d = (fCount - 2.0) * (fCount - 3.0);
    double k_l = fCount * (fCount + 1.0) / ((fCount - 1.0) * k_d);
    double k_t = 3.0 * (fCount - 1.0) * (fCount - 1.0) / k_d;

    PushDouble(xpower4 * k_l - k_t);
    PushDouble(xpower4.get() * k_l - k_t);
}

void ScInterpreter::ScHarMean()
@@ -3220,7 +3220,7 @@ void ScInterpreter::ScStandard()
            PushDouble((x-mue)/sigma);
    }
}
bool ScInterpreter::CalculateSkew(double& fSum,double& fCount,double& vSum,std::vector<double>& values)
bool ScInterpreter::CalculateSkew(KahanSum& fSum, double& fCount, std::vector<double>& values)
{
    short nParamCount = GetByte();
    if ( !MustHaveParamCountMin( nParamCount, 1 )  )
@@ -3228,7 +3228,6 @@ bool ScInterpreter::CalculateSkew(double& fSum,double& fCount,double& vSum,std::

    fSum   = 0.0;
    fCount = 0.0;
    vSum   = 0.0;
    double fVal = 0.0;
    ScAddress aAdr;
    ScRange aRange;
@@ -3328,9 +3327,10 @@ bool ScInterpreter::CalculateSkew(double& fSum,double& fCount,double& vSum,std::

void ScInterpreter::CalculateSkewOrSkewp( bool bSkewp )
{
    double fSum, fCount, vSum;
    KahanSum fSum;
    double fCount;
    std::vector<double> values;
    if (!CalculateSkew( fSum, fCount, vSum, values))
    if (!CalculateSkew( fSum, fCount, values))
        return;
     // SKEW/SKEWP's constraints: they require at least three numbers
    if (fCount < 3.0)
@@ -3340,30 +3340,29 @@ void ScInterpreter::CalculateSkewOrSkewp( bool bSkewp )
        return;
    }

    double fMean = fSum / fCount;

    KahanSum vSum;
    double fMean = fSum.get() / fCount;
    for (double v : values)
        vSum += (v - fMean) * (v - fMean);

    double fStdDev = sqrt( vSum / (bSkewp ? fCount : (fCount - 1.0)));
    double xcube = 0.0;

    double fStdDev = sqrt( vSum.get() / (bSkewp ? fCount : (fCount - 1.0)));
    if (fStdDev == 0)
    {
        PushIllegalArgument();
        return;
    }

    KahanSum xcube = 0.0;
    for (double v : values)
    {
        double dx = (v - fMean) / fStdDev;
        xcube = xcube + (dx * dx * dx);
        xcube += dx * dx * dx;
    }

    if (bSkewp)
        PushDouble( xcube / fCount );
        PushDouble( xcube.get() / fCount );
    else
        PushDouble( ((xcube * fCount) / (fCount - 1.0)) / (fCount - 2.0) );
        PushDouble( ((xcube.get() * fCount) / (fCount - 1.0)) / (fCount - 2.0) );
}

void ScInterpreter::ScSkew()