    /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
    /*
     * This file is part of the LibreOffice project.
     *
     * This Source Code Form is subject to the terms of the Mozilla Public
     * License, v. 2.0. If a copy of the MPL was not distributed with this
     * file, You can obtain one at http://mozilla.org/MPL/2.0/.
     *
     * This file incorporates work covered by the following license notice:
     *
     *   Licensed to the Apache Software Foundation (ASF) under one or more
     *   contributor license agreements. See the NOTICE file distributed
     *   with this work for additional information regarding copyright
     *   ownership. The ASF licenses this file to you under the Apache
     *   License, Version 2.0 (the "License"); you may not use this file
     *   except in compliance with the License. You may obtain a copy of
     *   the License at http://www.apache.org/licenses/LICENSE-2.0 .
     */

    #include <tools/solar.h>
    #include <stdlib.h>
    #include <string.h>

    #include <interpre.hxx>
    #include <global.hxx>
    #include <compiler.hxx>
    #include <formulacell.hxx>
    #include <document.hxx>
    #include <dociter.hxx>
    #include <matrixoperators.hxx>
    #include <scmatrix.hxx>

    #include <math.h>
    #include <memory>
    #include <vector>
    #include <algorithm>
    #include <comphelper/random.hxx>

    using ::std::vector;
    using namespace formula;

    /// Two columns of data should be sortable with GetSortArray() and QuickSort()
    // This is an arbitrary limit.
    #define MAX_COUNT_DOUBLE_FOR_SORT (MAXROWCOUNT * 2)

    const double ScInterpreter::fMaxGammaArgument = 171.624376956302;  // found experimental
    const double fMachEps = ::std::numeric_limits<double>::epsilon();

    class ScDistFunc
    {
    public:
        virtual double GetValue(double x) const = 0;

    protected:
        ~ScDistFunc() {}
    };

    //  iteration for inverse distributions

    //template< class T > double lcl_IterateInverse( const T& rFunction, double x0, double x1, bool& rConvError )

    /** u*w<0.0 fails for values near zero */
    static inline bool lcl_HasChangeOfSign( double u, double w )
    {
        return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0);
    }

    static double lcl_IterateInverse( const ScDistFunc& rFunction, double fAx, double fBx, bool& rConvError )
    {
        rConvError = false;
        const double fYEps = 1.0E-307;
        const double fXEps = ::std::numeric_limits<double>::epsilon();

        OSL_ENSURE(fAx<fBx, "IterateInverse: wrong interval");

        //  find enclosing interval

        double fAy = rFunction.GetValue(fAx);
        double fBy = rFunction.GetValue(fBx);
        double fTemp;
        unsigned short nCount;
        for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy); nCount++)
        {
            if (fabs(fAy) <= fabs(fBy))
            {
                fTemp = fAx;
                fAx += 2.0 * (fAx - fBx);
                if (fAx < 0.0)
                    fAx = 0.0;
                fBx = fTemp;
                fBy = fAy;
                fAy = rFunction.GetValue(fAx);
            }
            else
            {
                fTemp = fBx;
                fBx += 2.0 * (fBx - fAx);
                fAx = fTemp;
                fAy = fBy;
                fBy = rFunction.GetValue(fBx);
            }
        }

        if (fAy == 0.0)
            return fAx;
        if (fBy == 0.0)
            return fBx;
        if (!lcl_HasChangeOfSign( fAy, fBy))
        {
            rConvError = true;
            return 0.0;
        }
        // inverse quadric interpolation with additional brackets
        // set three points
        double fPx = fAx;
        double fPy = fAy;
        double fQx = fBx;
        double fQy = fBy;
        double fRx = fAx;
        double fRy = fAy;
        double fSx = 0.5 * (fAx + fBx); // potential next point
        bool bHasToInterpolate = true;
        nCount = 0;
        while ( nCount < 500 && fabs(fRy) > fYEps &&
                (fBx-fAx) > ::std::max( fabs(fAx), fabs(fBx)) * fXEps )
        {
            if (bHasToInterpolate)
            {
                if (fPy!=fQy && fQy!=fRy && fRy!=fPy)
                {
                    fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy)
                        + fRx * fQy * fPy / (fQy-fRy) / (fPy-fRy)
                        + fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy);
                    bHasToInterpolate = (fAx < fSx) && (fSx < fBx); // inside the brackets?
                }
                else
                    bHasToInterpolate = false;
            }
            if(!bHasToInterpolate)
            {
                fSx = 0.5 * (fAx + fBx);
                // reset points
                fQx = fBx; fQy = fBy;
                bHasToInterpolate = true;
            }
            // shift points for next interpolation
            fPx = fQx; fQx = fRx; fRx = fSx;
            fPy = fQy; fQy = fRy; fRy = rFunction.GetValue(fSx);
            // update brackets
            if (lcl_HasChangeOfSign( fAy, fRy))
            {
                fBx = fRx; fBy = fRy;
            }
            else
            {
                fAx = fRx; fAy = fRy;
            }
            // if last iteration brought to small advance, then do bisection next
            // time, for safety
            bHasToInterpolate = bHasToInterpolate && (fabs(fRy) * 2.0 <= fabs(fQy));
            ++nCount;
        }
        return fRx;
    }

    // General functions

    void ScInterpreter::ScNoName()
    {
        PushError(FormulaError::NoName);
    }

    void ScInterpreter::ScBadName()
    {
        short nParamCount = GetByte();
        while (nParamCount-- > 0)
        {
            PopError();
        }
        PushError( FormulaError::NoName);
    }

    double ScInterpreter::phi(double x)
    {
        return  0.39894228040143268 * exp(-(x * x) / 2.0);
    }

    double ScInterpreter::integralPhi(double x)
    { // Using gauss(x)+0.5 has severe cancellation errors for x<-4
        return 0.5 * ::rtl::math::erfc(-x * 0.7071067811865475); // * 1/sqrt(2)
    }

    double ScInterpreter::taylor(const double* pPolynom, sal_uInt16 nMax, double x)
    {
        double nVal = pPolynom[nMax];
        for (short i = nMax-1; i >= 0; i--)
        {
            nVal = pPolynom[i] + (nVal * x);
        }
        return nVal;
    }

    double ScInterpreter::gauss(double x)
    {

        double xAbs = fabs(x);
        sal_uInt16 xShort = static_cast<sal_uInt16>(::rtl::math::approxFloor(xAbs));
        double nVal = 0.0;
        if (xShort == 0)
        {
            static const double t0[] =
            { 0.39894228040143268, -0.06649038006690545,  0.00997355701003582,
             -0.00118732821548045,  0.00011543468761616, -0.00000944465625950,
              0.00000066596935163, -0.00000004122667415,  0.00000000227352982,
              0.00000000011301172,  0.00000000000511243, -0.00000000000021218 };
            nVal = taylor(t0, 11, (xAbs * xAbs)) * xAbs;
        }
        else if ((xShort >= 1) && (xShort <= 2))
        {
            static const double t2[] =
            { 0.47724986805182079,  0.05399096651318805, -0.05399096651318805,
              0.02699548325659403, -0.00449924720943234, -0.00224962360471617,
              0.00134977416282970, -0.00011783742691370, -0.00011515930357476,
              0.00003704737285544,  0.00000282690796889, -0.00000354513195524,
              0.00000037669563126,  0.00000019202407921, -0.00000005226908590,
             -0.00000000491799345,  0.00000000366377919, -0.00000000015981997,
             -0.00000000017381238,  0.00000000002624031,  0.00000000000560919,
             -0.00000000000172127, -0.00000000000008634,  0.00000000000007894 };
            nVal = taylor(t2, 23, (xAbs - 2.0));
        }
        else if ((xShort >= 3) && (xShort <= 4))
        {
            static const double t4[] =
           { 0.49996832875816688,  0.00013383022576489, -0.00026766045152977,
             0.00033457556441221, -0.00028996548915725,  0.00018178605666397,
            -0.00008252863922168,  0.00002551802519049, -0.00000391665839292,
            -0.00000074018205222,  0.00000064422023359, -0.00000017370155340,
             0.00000000909595465,  0.00000000944943118, -0.00000000329957075,
             0.00000000029492075,  0.00000000011874477, -0.00000000004420396,
             0.00000000000361422,  0.00000000000143638, -0.00000000000045848 };
            nVal = taylor(t4, 20, (xAbs - 4.0));
        }
        else
        {
            static const double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 };
            nVal = 0.5 + phi(xAbs) * taylor(asympt, 4, 1.0 / (xAbs * xAbs)) / xAbs;
        }
        if (x < 0.0)
            return -nVal;
        else
            return nVal;
    }

    //  #i26836# new gaussinv implementation by Martin Eitzenberger <m.eitzenberger@unix.net>

    double ScInterpreter::gaussinv(double x)
    {
        double q,t,z;

        q=x-0.5;

        if(fabs(q)<=.425)
        {
            t=0.180625-q*q;

            z=
            q*
            (
                (
                    (
                        (
                            (
                                (
                                    (
                                        t*2509.0809287301226727+33430.575583588128105
                                    )
                                    *t+67265.770927008700853
                                )
                                *t+45921.953931549871457
                            )
                            *t+13731.693765509461125
                        )
                        *t+1971.5909503065514427
                    )
                    *t+133.14166789178437745
                )
                *t+3.387132872796366608
            )
            /
            (
                (
                    (
                        (
                            (
                                (
                                    (
                                        t*5226.495278852854561+28729.085735721942674
                                    )
                                    *t+39307.89580009271061
                                )
                                *t+21213.794301586595867
                            )
                            *t+5394.1960214247511077
                        )
                        *t+687.1870074920579083
                    )
                    *t+42.313330701600911252
                )
                *t+1.0
            );

        }
        else
        {
            if(q>0) t=1-x;
            else        t=x;

            t=sqrt(-log(t));

            if(t<=5.0)
            {
                t+=-1.6;

                z=
                (
                    (
                        (
                            (
                                (
                                    (
                                        (
                                            t*7.7454501427834140764e-4+0.0227238449892691845833
                                        )
                                        *t+0.24178072517745061177
                                    )
                                    *t+1.27045825245236838258
                                )
                                *t+3.64784832476320460504
                            )
                            *t+5.7694972214606914055
                        )
                        *t+4.6303378461565452959
                    )
                    *t+1.42343711074968357734
                )
                /
                (
                    (
                        (
                            (
                                (
                                    (
                                        (
                                            t*1.05075007164441684324e-9+5.475938084995344946e-4
                                        )
                                        *t+0.0151986665636164571966
                                    )
                                    *t+0.14810397642748007459
                                )
                                *t+0.68976733498510000455
                            )
                            *t+1.6763848301838038494
                        )
                        *t+2.05319162663775882187
                    )
                    *t+1.0
                );

            }
            else
            {
                t+=-5.0;

                z=
                (
                    (
                        (
                            (
                                (
                                    (
                                        (
                                            t*2.01033439929228813265e-7+2.71155556874348757815e-5
                                        )
                                        *t+0.0012426609473880784386
                                    )
                                    *t+0.026532189526576123093
                                )
                                *t+0.29656057182850489123
                            )
                            *t+1.7848265399172913358
                        )
                        *t+5.4637849111641143699
                    )
                    *t+6.6579046435011037772
                )
                /
                (
                    (
                        (
                            (
                                (
                                    (
                                        (
                                            t*2.04426310338993978564e-15+1.4215117583164458887e-7
                                        )
                                        *t+1.8463183175100546818e-5
                                    )
                                    *t+7.868691311456132591e-4
                                )
                                *t+0.0148753612908506148525
                            )
                            *t+0.13692988092273580531
                        )
                        *t+0.59983220655588793769
                    )
                    *t+1.0
                );

            }

            if(q<0.0) z=-z;
        }

        return z;
    }

    double ScInterpreter::Fakultaet(double x)
    {
        x = ::rtl::math::approxFloor(x);
        if (x < 0.0)
            return 0.0;
        else if (x == 0.0)
            return 1.0;
        else if (x <= 170.0)
        {
            double fTemp = x;
            while (fTemp > 2.0)
            {
                fTemp--;
                x *= fTemp;
            }
        }
        else
            SetError(FormulaError::NoValue);
        return x;
    }

    double ScInterpreter::BinomKoeff(double n, double k)
    {
        // this method has been duplicated as BinomialCoefficient()
        // in scaddins/source/analysis/analysishelper.cxx

        double nVal = 0.0;
        k = ::rtl::math::approxFloor(k);
        if (n < k)
            nVal = 0.0;
        else if (k == 0.0)
            nVal = 1.0;
        else
        {
            nVal = n/k;
            n--;
            k--;
            while (k > 0.0)
            {
                nVal *= n/k;
                k--;
                n--;
            }

        }
        return nVal;
    }

    // The algorithm is based on lanczos13m53 in lanczos.hpp
    // in math library from http://www.boost.org
    /** you must ensure fZ>0
        Uses a variant of the Lanczos sum with a rational function. */
    static double lcl_getLanczosSum(double fZ)
    {
        static const double fNum[13] ={
            23531376880.41075968857200767445163675473,
            42919803642.64909876895789904700198885093,
            35711959237.35566804944018545154716670596,
            17921034426.03720969991975575445893111267,
            6039542586.35202800506429164430729792107,
            1439720407.311721673663223072794912393972,
            248874557.8620541565114603864132294232163,
            31426415.58540019438061423162831820536287,
            2876370.628935372441225409051620849613599,
            186056.2653952234950402949897160456992822,
            8071.672002365816210638002902272250613822,
            210.8242777515793458725097339207133627117,
            2.506628274631000270164908177133837338626
            };
        static const double fDenom[13] = {
            0,
            39916800,
            120543840,
            150917976,
            105258076,
            45995730,
            13339535,
            2637558,
            357423,
            32670,
            1925,
            66,
            1
            };
        // Horner scheme
        double fSumNum;
        double fSumDenom;
        int nI;
        if (fZ<=1.0)
        {
            fSumNum = fNum[12];
            fSumDenom = fDenom[12];
            for (nI = 11; nI >= 0; --nI)
            {
                fSumNum *= fZ;
                fSumNum += fNum[nI];
                fSumDenom *= fZ;
                fSumDenom += fDenom[nI];
            }
        }
        else
        // Cancel down with fZ^12; Horner scheme with reverse coefficients
        {
            double fZInv = 1/fZ;
            fSumNum = fNum[0];
            fSumDenom = fDenom[0];
            for (nI = 1; nI <=12; ++nI)
            {
                fSumNum *= fZInv;
                fSumNum += fNum[nI];
                fSumDenom *= fZInv;
                fSumDenom += fDenom[nI];
            }
        }
        return fSumNum/fSumDenom;
    }

    // The algorithm is based on tgamma in gamma.hpp
    // in math library from http://www.boost.org
    /** You must ensure fZ>0; fZ>171.624376956302 will overflow. */
    static double lcl_GetGammaHelper(double fZ)
    {
        double fGamma = lcl_getLanczosSum(fZ);
        const double fg = 6.024680040776729583740234375;
        double fZgHelp = fZ + fg - 0.5;
        // avoid intermediate overflow
        double fHalfpower = pow( fZgHelp, fZ / 2 - 0.25);
        fGamma *= fHalfpower;
        fGamma /= exp(fZgHelp);
        fGamma *= fHalfpower;
        if (fZ <= 20.0 && fZ == ::rtl::math::approxFloor(fZ))
            fGamma = ::rtl::math::round(fGamma);
        return fGamma;
    }

    // The algorithm is based on tgamma in gamma.hpp
    // in math library from http://www.boost.org
    /** You must ensure fZ>0 */
    static double lcl_GetLogGammaHelper(double fZ)
    {
        const double fg = 6.024680040776729583740234375;
        double fZgHelp = fZ + fg - 0.5;
        return log( lcl_getLanczosSum(fZ)) + (fZ-0.5) * log(fZgHelp) - fZgHelp;
    }

    /** You must ensure non integer arguments for fZ<1 */
    double ScInterpreter::GetGamma(double fZ)
    {
        const double fLogPi = log(F_PI);
        const double fLogDblMax = log( ::std::numeric_limits<double>::max());

        if (fZ > fMaxGammaArgument)
        {
            SetError(FormulaError::IllegalFPOperation);
            return HUGE_VAL;
        }

        if (fZ >= 1.0)
            return lcl_GetGammaHelper(fZ);

        if (fZ >= 0.5)  // shift to x>=1 using Gamma(x)=Gamma(x+1)/x
            return lcl_GetGammaHelper(fZ+1) / fZ;

        if (fZ >= -0.5) // shift to x>=1, might overflow
        {
            double fLogTest = lcl_GetLogGammaHelper(fZ+2) - rtl::math::log1p(fZ) - log( fabs(fZ));
            if (fLogTest >= fLogDblMax)
            {
                SetError( FormulaError::IllegalFPOperation);
                return HUGE_VAL;
            }
            return lcl_GetGammaHelper(fZ+2) / (fZ+1) / fZ;
        }
        // fZ<-0.5
        // Use Euler's reflection formula: gamma(x)= pi/ ( gamma(1-x)*sin(pi*x) )
        double fLogDivisor = lcl_GetLogGammaHelper(1-fZ) + log( fabs( ::rtl::math::sin( F_PI*fZ)));
        if (fLogDivisor - fLogPi >= fLogDblMax)     // underflow
            return 0.0;

        if (fLogDivisor<0.0)
            if (fLogPi - fLogDivisor > fLogDblMax)  // overflow
            {
                SetError(FormulaError::IllegalFPOperation);
                return HUGE_VAL;
            }

        return exp( fLogPi - fLogDivisor) * ((::rtl::math::sin( F_PI*fZ) < 0.0) ? -1.0 : 1.0);
    }

    /** You must ensure fZ>0 */
    double ScInterpreter::GetLogGamma(double fZ)
    {
        if (fZ >= fMaxGammaArgument)
            return lcl_GetLogGammaHelper(fZ);
        if (fZ >= 1.0)
            return log(lcl_GetGammaHelper(fZ));
        if (fZ >= 0.5)
            return log( lcl_GetGammaHelper(fZ+1) / fZ);
        return lcl_GetLogGammaHelper(fZ+2) - rtl::math::log1p(fZ) - log(fZ);
    }

    double ScInterpreter::GetFDist(double x, double fF1, double fF2)
    {
        double arg = fF2/(fF2+fF1*x);
        double alpha = fF2/2.0;
        double beta = fF1/2.0;
        return GetBetaDist(arg, alpha, beta);
    }

    double ScInterpreter::GetTDist( double T, double fDF, int nType )
    {
        switch ( nType )
        {
            case 1 : // 1-tailed T-distribution
                return 0.5 * GetBetaDist( fDF / ( fDF + T * T ), fDF / 2.0, 0.5 );
            case 2 : // 2-tailed T-distribution
                return GetBetaDist( fDF / ( fDF + T * T ), fDF / 2.0, 0.5);
            case 3 : // left-tailed T-distribution (probability density function)
                return pow( 1 + ( T * T / fDF ), -( fDF + 1 ) / 2 ) / ( sqrt( fDF ) * GetBeta( 0.5, fDF / 2.0 ) );
            case 4 : // left-tailed T-distribution (cumulative distribution function)
                double X = fDF / ( T * T + fDF );
                double R = 0.5 * GetBetaDist( X, 0.5 * fDF, 0.5 );
                return ( T < 0 ? R : 1 - R );
        }
        SetError( FormulaError::IllegalArgument );
        return HUGE_VAL;
    }

    // for LEGACY.CHIDIST, returns right tail, fDF=degrees of freedom
    /** You must ensure fDF>0.0 */
    double ScInterpreter::GetChiDist(double fX, double fDF)
    {
        if (fX <= 0.0)
            return 1.0; // see ODFF
        else
            return GetUpRegIGamma( fDF/2.0, fX/2.0);
    }

    // ready for ODF 1.2
    // for ODF CHISQDIST; cumulative distribution function, fDF=degrees of freedom
    // returns left tail
    /** You must ensure fDF>0.0 */
    double ScInterpreter::GetChiSqDistCDF(double fX, double fDF)
    {
        if (fX <= 0.0)
            return 0.0; // see ODFF
        else
            return GetLowRegIGamma( fDF/2.0, fX/2.0);
    }

    double ScInterpreter::GetChiSqDistPDF(double fX, double fDF)
    {
        // you must ensure fDF is positive integer
        double fValue;
        if (fX <= 0.0)
            return 0.0; // see ODFF
        if (fDF*fX > 1391000.0)
        {
            // intermediate invalid values, use log
            fValue = exp((0.5*fDF - 1) * log(fX*0.5) - 0.5 * fX - log(2.0) - GetLogGamma(0.5*fDF));
        }
        else // fDF is small in most cases, we can iterate
        {
            double fCount;
            if (fmod(fDF,2.0)<0.5)
            {
                // even
                fValue = 0.5;
                fCount = 2.0;
            }
            else
            {
                fValue = 1/sqrt(fX*2*F_PI);
                fCount = 1.0;
            }
            while ( fCount < fDF)
            {
                fValue *= (fX / fCount);
                fCount += 2.0;
            }
            if (fX>=1425.0) // underflow in e^(-x/2)
                fValue = exp(log(fValue)-fX/2);
            else
                fValue *= exp(-fX/2);
        }
        return fValue;
    }

    void ScInterpreter::ScChiSqDist()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
            return;
        bool bCumulative;
        if (nParamCount == 3)
            bCumulative = GetBool();
        else
            bCumulative = true;
        double fDF = ::rtl::math::approxFloor(GetDouble());
        if (fDF < 1.0)
            PushIllegalArgument();
        else
        {
            double fX = GetDouble();
            if (bCumulative)
                PushDouble(GetChiSqDistCDF(fX,fDF));
            else
                PushDouble(GetChiSqDistPDF(fX,fDF));
        }
    }

    void ScInterpreter::ScChiSqDist_MS()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, 3, 3 ) )
            return;
        bool bCumulative = GetBool();
        double fDF = ::rtl::math::approxFloor( GetDouble() );
        if ( fDF < 1.0 || fDF > 1E10 )
            PushIllegalArgument();
        else
        {
            double fX = GetDouble();
            if ( fX < 0 )
                PushIllegalArgument();
            else
            {
                if ( bCumulative )
                    PushDouble( GetChiSqDistCDF( fX, fDF ) );
                else
                    PushDouble( GetChiSqDistPDF( fX, fDF ) );
            }
        }
    }

    void ScInterpreter::ScGamma()
    {
        double x = GetDouble();
        if (x <= 0.0 && x == ::rtl::math::approxFloor(x))
            PushIllegalArgument();
        else
        {
            double fResult = GetGamma(x);
            if (nGlobalError != FormulaError::NONE)
            {
                PushError( nGlobalError);
                return;
            }
            PushDouble(fResult);
        }
    }

    void ScInterpreter::ScLogGamma()
    {
        double x = GetDouble();
        if (x > 0.0)    // constraint from ODFF
            PushDouble( GetLogGamma(x));
        else
            PushIllegalArgument();
    }

    double ScInterpreter::GetBeta(double fAlpha, double fBeta)
    {
        double fA;
        double fB;
        if (fAlpha > fBeta)
        {
            fA = fAlpha; fB = fBeta;
        }
        else
        {
            fA = fBeta; fB = fAlpha;
        }
        if (fA+fB < fMaxGammaArgument) // simple case
            return GetGamma(fA)/GetGamma(fA+fB)*GetGamma(fB);
        // need logarithm
        // GetLogGamma is not accurate enough, back to Lanczos for all three
        // GetGamma and arrange factors newly.
        const double fg = 6.024680040776729583740234375; //see GetGamma
        double fgm = fg - 0.5;
        double fLanczos = lcl_getLanczosSum(fA);
        fLanczos /= lcl_getLanczosSum(fA+fB);
        fLanczos *= lcl_getLanczosSum(fB);
        double fABgm = fA+fB+fgm;
        fLanczos *= sqrt((fABgm/(fA+fgm))/(fB+fgm));
        double fTempA = fB/(fA+fgm); // (fA+fgm)/fABgm = 1 / ( 1 + fB/(fA+fgm))
        double fTempB = fA/(fB+fgm);
        double fResult = exp(-fA * ::rtl::math::log1p(fTempA)
                                -fB * ::rtl::math::log1p(fTempB)-fgm);
        fResult *= fLanczos;
        return fResult;
    }

    // Same as GetBeta but with logarithm
    double ScInterpreter::GetLogBeta(double fAlpha, double fBeta)
    {
        double fA;
        double fB;
        if (fAlpha > fBeta)
        {
            fA = fAlpha; fB = fBeta;
        }
        else
        {
            fA = fBeta; fB = fAlpha;
        }
        const double fg = 6.024680040776729583740234375; //see GetGamma
        double fgm = fg - 0.5;
        double fLanczos = lcl_getLanczosSum(fA);
        fLanczos /= lcl_getLanczosSum(fA+fB);
        fLanczos *= lcl_getLanczosSum(fB);
        double fLogLanczos = log(fLanczos);
        double fABgm = fA+fB+fgm;
        fLogLanczos += 0.5*(log(fABgm)-log(fA+fgm)-log(fB+fgm));
        double fTempA = fB/(fA+fgm); // (fA+fgm)/fABgm = 1 / ( 1 + fB/(fA+fgm))
        double fTempB = fA/(fB+fgm);
        double fResult = -fA * ::rtl::math::log1p(fTempA)
                            -fB * ::rtl::math::log1p(fTempB)-fgm;
        fResult += fLogLanczos;
        return fResult;
    }

    // beta distribution probability density function
    double ScInterpreter::GetBetaDistPDF(double fX, double fA, double fB)
    {
        // special cases
        if (fA == 1.0) // result b*(1-x)^(b-1)
        {
            if (fB == 1.0)
                return 1.0;
            if (fB == 2.0)
                return -2.0*fX + 2.0;
            if (fX == 1.0 && fB < 1.0)
            {
                SetError(FormulaError::IllegalArgument);
                return HUGE_VAL;
            }
            if (fX <= 0.01)
                return fB + fB * ::rtl::math::expm1((fB-1.0) * ::rtl::math::log1p(-fX));
            else
                return fB * pow(0.5-fX+0.5,fB-1.0);
        }
        if (fB == 1.0) // result a*x^(a-1)
        {
            if (fA == 2.0)
                return fA * fX;
            if (fX == 0.0 && fA < 1.0)
            {
                SetError(FormulaError::IllegalArgument);
                return HUGE_VAL;
            }
            return fA * pow(fX,fA-1);
        }
        if (fX <= 0.0)
        {
            if (fA < 1.0 && fX == 0.0)
            {
                SetError(FormulaError::IllegalArgument);
                return HUGE_VAL;
            }
            else
                return 0.0;
        }
        if (fX >= 1.0)
        {
            if (fB < 1.0 && fX == 1.0)
            {
                SetError(FormulaError::IllegalArgument);
                return HUGE_VAL;
            }
            else
                return 0.0;
        }

        // normal cases; result x^(a-1)*(1-x)^(b-1)/Beta(a,b)
        const double fLogDblMax = log( ::std::numeric_limits<double>::max());
        const double fLogDblMin = log( ::std::numeric_limits<double>::min());
        double fLogY = (fX < 0.1) ? ::rtl::math::log1p(-fX) : log(0.5-fX+0.5);
        double fLogX = log(fX);
        double fAm1LogX = (fA-1.0) * fLogX;
        double fBm1LogY = (fB-1.0) * fLogY;
        double fLogBeta = GetLogBeta(fA,fB);
        // check whether parts over- or underflow
        if (   fAm1LogX < fLogDblMax  && fAm1LogX > fLogDblMin
            && fBm1LogY < fLogDblMax  && fBm1LogY > fLogDblMin
            && fLogBeta < fLogDblMax  && fLogBeta > fLogDblMin
            && fAm1LogX + fBm1LogY < fLogDblMax && fAm1LogX + fBm1LogY > fLogDblMin)
            return pow(fX,fA-1.0) * pow(0.5-fX+0.5,fB-1.0) / GetBeta(fA,fB);
        else // need logarithm;
            // might overflow as a whole, but seldom, not worth to pre-detect it
            return exp( fAm1LogX + fBm1LogY - fLogBeta);
    }

    /*
                    x^a * (1-x)^b
        I_x(a,b) = ----------------  * result of ContFrac
                    a * Beta(a,b)
    */
    static double lcl_GetBetaHelperContFrac(double fX, double fA, double fB)
    {   // like old version
        double a1, b1, a2, b2, fnorm, cfnew, cf;
        a1 = 1.0; b1 = 1.0;
        b2 = 1.0 - (fA+fB)/(fA+1.0)*fX;
        if (b2 == 0.0)
        {
            a2 = 0.0;
            fnorm = 1.0;
            cf = 1.0;
        }
        else
        {
            a2 = 1.0;
            fnorm = 1.0/b2;
            cf = a2*fnorm;
        }
        cfnew = 1.0;
        double rm = 1.0;

        const double fMaxIter = 50000.0;
        // loop security, normal cases converge in less than 100 iterations.
        // FIXME: You will get so much iterations for fX near mean,
        // I do not know a better algorithm.
        bool bfinished = false;
        do
        {
            const double apl2m = fA + 2.0*rm;
            const double d2m = rm*(fB-rm)*fX/((apl2m-1.0)*apl2m);
            const double d2m1 = -(fA+rm)*(fA+fB+rm)*fX/(apl2m*(apl2m+1.0));
            a1 = (a2+d2m*a1)*fnorm;
            b1 = (b2+d2m*b1)*fnorm;
            a2 = a1 + d2m1*a2*fnorm;
            b2 = b1 + d2m1*b2*fnorm;
            if (b2 != 0.0)
            {
                fnorm = 1.0/b2;
                cfnew = a2*fnorm;
                bfinished = (fabs(cf-cfnew) < fabs(cf)*fMachEps);
            }
            cf = cfnew;
            rm += 1.0;
        }
        while (rm < fMaxIter && !bfinished);
        return cf;
    }

    // cumulative distribution function, normalized
    double ScInterpreter::GetBetaDist(double fXin, double fAlpha, double fBeta)
    {
    // special cases
        if (fXin <= 0.0)  // values are valid, see spec
            return 0.0;
        if (fXin >= 1.0)  // values are valid, see spec
            return 1.0;
        if (fBeta == 1.0)
            return pow(fXin, fAlpha);
        if (fAlpha == 1.0)
        //            1.0 - pow(1.0-fX,fBeta) is not accurate enough
            return -::rtl::math::expm1(fBeta * ::rtl::math::log1p(-fXin));
        //FIXME: need special algorithm for fX near fP for large fA,fB
        double fResult;
        // I use always continued fraction, power series are neither
        // faster nor more accurate.
        double fY = (0.5-fXin)+0.5;
        double flnY = ::rtl::math::log1p(-fXin);
        double fX = fXin;
        double flnX = log(fXin);
        double fA = fAlpha;
        double fB = fBeta;
        bool bReflect = fXin > fAlpha/(fAlpha+fBeta);
        if (bReflect)
        {
            fA = fBeta;
            fB = fAlpha;
            fX = fY;
            fY = fXin;
            flnX = flnY;
            flnY = log(fXin);
        }
        fResult = lcl_GetBetaHelperContFrac(fX,fA,fB);
        fResult = fResult/fA;
        double fP = fA/(fA+fB);
        double fQ = fB/(fA+fB);
        double fTemp;
        if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97) //found experimental
            fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY;
        else
            fTemp = exp(fA*flnX + fB*flnY - GetLogBeta(fA,fB));
        fResult *= fTemp;
        if (bReflect)
            fResult = 0.5 - fResult + 0.5;
        if (fResult > 1.0) // ensure valid range
            fResult = 1.0;
        if (fResult < 0.0)
            fResult = 0.0;
        return fResult;
    }

    void ScInterpreter::ScBetaDist()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, 3, 6 ) ) // expanded, see #i91547#
            return;
        double fLowerBound, fUpperBound;
        double alpha, beta, x;
        bool bIsCumulative;
        if (nParamCount == 6)
            bIsCumulative = GetBool();
        else
            bIsCumulative = true;
        if (nParamCount >= 5)
            fUpperBound = GetDouble();
        else
            fUpperBound = 1.0;
        if (nParamCount >= 4)
            fLowerBound = GetDouble();
        else
            fLowerBound = 0.0;
        beta = GetDouble();
        alpha = GetDouble();
        x = GetDouble();
        double fScale = fUpperBound - fLowerBound;
        if (fScale <= 0.0 || alpha <= 0.0 || beta <= 0.0)
        {
            PushIllegalArgument();
            return;
        }
        if (bIsCumulative) // cumulative distribution function
        {
            // special cases
            if (x < fLowerBound)
            {
                PushDouble(0.0); return; //see spec
            }
            if (x > fUpperBound)
            {
                PushDouble(1.0); return; //see spec
            }
            // normal cases
            x = (x-fLowerBound)/fScale;  // convert to standard form
            PushDouble(GetBetaDist(x, alpha, beta));
            return;
        }
        else // probability density function
        {
            if (x < fLowerBound || x > fUpperBound)
            {
                PushDouble(0.0);
                return;
            }
            x = (x-fLowerBound)/fScale;
            PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale);
            return;
        }
    }

    /**
      Microsoft version has parameters in different order
      Also, upper and lowerbound are optional and have default values
      and different constraints apply.
      Basically, function is identical with ScInterpreter::ScBetaDist()
    */
    void ScInterpreter::ScBetaDist_MS()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, 4, 6 ) )
            return;
        double fLowerBound, fUpperBound;
        double alpha, beta, x;
        bool bIsCumulative;
        if (nParamCount == 6)
            fUpperBound = GetDouble();
        else
            fUpperBound = 1.0;
        if (nParamCount >= 5)
            fLowerBound = GetDouble();
        else
            fLowerBound = 0.0;
        bIsCumulative = GetBool();
        beta = GetDouble();
        alpha = GetDouble();
        x = GetDouble();
        if (alpha <= 0.0 || beta <= 0.0 || x < fLowerBound || x > fUpperBound)
        {
            PushIllegalArgument();
            return;
        }
        double fScale = fUpperBound - fLowerBound;
        if (bIsCumulative) // cumulative distribution function
        {
            x = (x-fLowerBound)/fScale;  // convert to standard form
            PushDouble(GetBetaDist(x, alpha, beta));
            return;
        }
        else // probability density function
        {
            x = (x-fLowerBound)/fScale;
            PushDouble(GetBetaDistPDF(x, alpha, beta)/fScale);
            return;
        }
    }

    void ScInterpreter::ScPhi()
    {
        PushDouble(phi(GetDouble()));
    }

    void ScInterpreter::ScGauss()
    {
        PushDouble(gauss(GetDouble()));
    }

    void ScInterpreter::ScFisher()
    {
        double fVal = GetDouble();
        if (fabs(fVal) >= 1.0)
            PushIllegalArgument();
        else
            PushDouble( ::rtl::math::atanh( fVal));
    }

    void ScInterpreter::ScFisherInv()
    {
        PushDouble( tanh( GetDouble()));
    }

    void ScInterpreter::ScFact()
    {
        double nVal = GetDouble();
        if (nVal < 0.0)
            PushIllegalArgument();
        else
            PushDouble(Fakultaet(nVal));
    }

    void ScInterpreter::ScCombin()
    {
        if ( MustHaveParamCount( GetByte(), 2 ) )
        {
            double k = ::rtl::math::approxFloor(GetDouble());
            double n = ::rtl::math::approxFloor(GetDouble());
            if (k < 0.0 || n < 0.0 || k > n)
                PushIllegalArgument();
            else
                PushDouble(BinomKoeff(n, k));
        }
    }

    void ScInterpreter::ScCombinA()
    {
        if ( MustHaveParamCount( GetByte(), 2 ) )
        {
            double k = ::rtl::math::approxFloor(GetDouble());
            double n = ::rtl::math::approxFloor(GetDouble());
            if (k < 0.0 || n < 0.0 || k > n)
                PushIllegalArgument();
            else
                PushDouble(BinomKoeff(n + k - 1, k));
        }
    }

    void ScInterpreter::ScPermut()
    {
        if ( MustHaveParamCount( GetByte(), 2 ) )
        {
            double k = ::rtl::math::approxFloor(GetDouble());
            double n = ::rtl::math::approxFloor(GetDouble());
            if (n < 0.0 || k < 0.0 || k > n)
                PushIllegalArgument();
            else if (k == 0.0)
                PushInt(1);     // (n! / (n - 0)!) == 1
            else
            {
                double nVal = n;
                for (sal_uLong i = static_cast<sal_uLong>(k)-1; i >= 1; i--)
                    nVal *= n-static_cast<double>(i);
                PushDouble(nVal);
            }
        }
    }

    void ScInterpreter::ScPermutationA()
    {
        if ( MustHaveParamCount( GetByte(), 2 ) )
        {
            double k = ::rtl::math::approxFloor(GetDouble());
            double n = ::rtl::math::approxFloor(GetDouble());
            if (n < 0.0 || k < 0.0)
                PushIllegalArgument();
            else
                PushDouble(pow(n,k));
        }
    }

    double ScInterpreter::GetBinomDistPMF(double x, double n, double p)
    // used in ScB and ScBinomDist
    // preconditions: 0.0 <= x <= n, 0.0 < p < 1.0;  x,n integral although double
    {
        double q = (0.5 - p) + 0.5;
        double fFactor = pow(q, n);
        if (fFactor <=::std::numeric_limits<double>::min())
        {
            fFactor = pow(p, n);
            if (fFactor <= ::std::numeric_limits<double>::min())
                return GetBetaDistPDF(p, x+1.0, n-x+1.0)/(n+1.0);
            else
            {
                sal_uInt32 max = static_cast<sal_uInt32>(n - x);
                for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
                    fFactor *= (n-i)/(i+1)*q/p;
                return fFactor;
            }
        }
        else
        {
            sal_uInt32 max = static_cast<sal_uInt32>(x);
            for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
                fFactor *= (n-i)/(i+1)*p/q;
            return fFactor;
        }
    }

    double lcl_GetBinomDistRange(double n, double xs,double xe,
                double fFactor /* q^n */, double p, double q)
    //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
        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);
    }

    void ScInterpreter::ScB()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
            return ;
        if (nParamCount == 3)   // mass function
        {
            double x = ::rtl::math::approxFloor(GetDouble());
            double p = GetDouble();
            double n = ::rtl::math::approxFloor(GetDouble());
            if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
                PushIllegalArgument();
            else if (p == 0.0)
                PushDouble( (x == 0.0) ? 1.0 : 0.0 );
            else if ( p == 1.0)
                PushDouble( (x == n) ? 1.0 : 0.0);
            else
                PushDouble(GetBinomDistPMF(x,n,p));
        }
        else
        {   // nParamCount == 4
            double xe = ::rtl::math::approxFloor(GetDouble());
            double xs = ::rtl::math::approxFloor(GetDouble());
            double p = GetDouble();
            double n = ::rtl::math::approxFloor(GetDouble());
            double q = (0.5 - p) + 0.5;
            bool bIsValidX = ( 0.0 <= xs && xs <= xe && xe <= n);
            if ( bIsValidX && 0.0 < p && p < 1.0)
            {
                if (xs == xe)       // mass function
                    PushDouble(GetBinomDistPMF(xs,n,p));
                else
                {
                    double fFactor = pow(q, n);
                    if (fFactor > ::std::numeric_limits<double>::min())
                        PushDouble(lcl_GetBinomDistRange(n,xs,xe,fFactor,p,q));
                    else
                    {
                        fFactor = pow(p, n);
                        if (fFactor > ::std::numeric_limits<double>::min())
                        {
                            // sum from j=xs to xe {(n choose j) * p^j * q^(n-j)}
                            // = sum from i = n-xe to n-xs { (n choose i) * q^i * p^(n-i)}
                            PushDouble(lcl_GetBinomDistRange(n,n-xe,n-xs,fFactor,q,p));
                        }
                        else
                            PushDouble(GetBetaDist(q,n-xe,xe+1.0)-GetBetaDist(q,n-xs+1,xs) );
                    }
                }
            }
            else
            {
                if ( bIsValidX ) // not(0<p<1)
                {
                    if ( p == 0.0 )
                        PushDouble( (xs == 0.0) ? 1.0 : 0.0 );
                    else if ( p == 1.0 )
                        PushDouble( (xe == n) ? 1.0 : 0.0 );
                    else
                        PushIllegalArgument();
                }
                else
                    PushIllegalArgument();
            }
        }
    }

    void ScInterpreter::ScBinomDist()
    {
        if ( MustHaveParamCount( GetByte(), 4 ) )
        {
            bool bIsCum   = GetBool();     // false=mass function; true=cumulative
            double p      = GetDouble();
            double n      = ::rtl::math::approxFloor(GetDouble());
            double x      = ::rtl::math::approxFloor(GetDouble());
            double q = (0.5 - p) + 0.5;           // get one bit more for p near 1.0
            if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0)
            {
                PushIllegalArgument();
                return;
            }
            if ( p == 0.0)
            {
                PushDouble( (x==0.0 || bIsCum) ? 1.0 : 0.0 );
                return;
            }
            if ( p == 1.0)
            {
                PushDouble( (x==n) ? 1.0 : 0.0);
                return;
            }
            if (!bIsCum)
                PushDouble( GetBinomDistPMF(x,n,p));
            else
            {
                if (x == n)
                    PushDouble(1.0);
                else
                {
                    double fFactor = pow(q, n);
                    if (x == 0.0)
                        PushDouble(fFactor);
                    else if (fFactor <= ::std::numeric_limits<double>::min())
                    {
                        fFactor = pow(p, n);
                        if (fFactor <= ::std::numeric_limits<double>::min())
                            PushDouble(GetBetaDist(q,n-x,x+1.0));
                        else
                        {
                            if (fFactor > fMachEps)
                            {
                                double fSum = 1.0 - fFactor;
                                sal_uInt32 max = static_cast<sal_uInt32> (n - x) - 1;
                                for (sal_uInt32 i = 0; i < max && fFactor > 0.0; i++)
                                {
                                    fFactor *= (n-i)/(i+1)*q/p;
                                    fSum -= fFactor;
                                }
                                PushDouble( (fSum < 0.0) ? 0.0 : fSum );
                            }
                            else
                                PushDouble(lcl_GetBinomDistRange(n,n-x,n,fFactor,q,p));
                        }
                    }
                    else
                        PushDouble( lcl_GetBinomDistRange(n,0.0,x,fFactor,p,q)) ;
                }
            }
        }
    }

    void ScInterpreter::ScCritBinom()
    {
        if ( MustHaveParamCount( GetByte(), 3 ) )
        {
            double alpha  = GetDouble();
            double p      = GetDouble();
            double n      = ::rtl::math::approxFloor(GetDouble());
            if (n < 0.0 || alpha < 0.0 || alpha > 1.0 || p < 0.0 || p > 1.0)
                PushIllegalArgument();
            else if ( alpha == 0.0 )
                PushDouble( 0.0 );
            else if ( alpha == 1.0 )
                PushDouble( p == 0 ? 0.0 : n );
            else
            {
                double fFactor;
                double q = (0.5 - p) + 0.5;           // get one bit more for p near 1.0
                if ( q > p )                          // work from the side where the cumulative curve is
                {
                    // work from 0 upwards
                    fFactor = pow(q,n);
                    if (fFactor > ::std::numeric_limits<double>::min())
                    {
                        double fSum = fFactor;
                        sal_uInt32 max = static_cast<sal_uInt32> (n), i;
                        for (i = 0; i < max && fSum < alpha; i++)
                        {
                            fFactor *= (n-i)/(i+1)*p/q;
                            fSum += fFactor;
                        }
                        PushDouble(i);
                    }
                    else
                    {
                        // accumulate BinomDist until accumulated BinomDist reaches alpha
                        double 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();
                                return;
                            }
                        }
                        PushDouble( i - 1 );
                    }
                }
                else
                {
                    // work from n backwards
                    fFactor = pow(p, n);
                    if (fFactor > ::std::numeric_limits<double>::min())
                    {
                        double fSum = 1.0 - fFactor;
                        sal_uInt32 max = static_cast<sal_uInt32> (n), i;
                        for (i = 0; i < max && fSum >= alpha; i++)
                        {
                            fFactor *= (n-i)/(i+1)*q/p;
                            fSum -= fFactor;
                        }
                        PushDouble(n-i);
                    }
                    else
                    {
                        // accumulate BinomDist until accumulated BinomDist reaches alpha
                        double 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();
                                return;
                            }
                        }
                        PushDouble( n - i + 1 );
                    }
                }
            }
        }
    }

    void ScInterpreter::ScNegBinomDist()
    {
        if ( MustHaveParamCount( GetByte(), 3 ) )
        {
            double p = GetDouble();                            // probability
            double s = ::rtl::math::approxFloor(GetDouble());  // No of successes
            double f = ::rtl::math::approxFloor(GetDouble());  // No of failures
            if ((f + s) <= 1.0 || p < 0.0 || p > 1.0)
                PushIllegalArgument();
            else
            {
                double q = 1.0 - p;
                double fFactor = pow(p,s);
                for (double i = 0.0; i < f; i++)
                    fFactor *= (i+s)/(i+1.0)*q;
                PushDouble(fFactor);
            }
        }
    }

    void ScInterpreter::ScNegBinomDist_MS()
    {
        if ( MustHaveParamCount( GetByte(), 4 ) )
        {
            bool bCumulative = GetBool();
            double p = GetDouble();                            // probability
            double s = ::rtl::math::approxFloor(GetDouble());  // No of successes
            double f = ::rtl::math::approxFloor(GetDouble());  // No of failures
            if ( s < 1.0 || f < 0.0 || p < 0.0 || p > 1.0 )
                PushIllegalArgument();
            else
            {
                double q = 1.0 - p;
                if ( bCumulative )
                    PushDouble( 1.0 - GetBetaDist( q, f + 1, s ) );
                else
                {
                    double fFactor = pow( p, s );
                    for ( double i = 0.0; i < f; i++ )
                        fFactor *= ( i + s ) / ( i + 1.0 ) * q;
                    PushDouble( fFactor );
                }
            }
        }
    }

    void ScInterpreter::ScNormDist( int nMinParamCount )
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) )
            return;
        bool bCumulative = nParamCount != 4 || GetBool();
        double sigma = GetDouble();                 // standard deviation
        double mue = GetDouble();                   // mean
        double x = GetDouble();                     // x
        if (sigma <= 0.0)
        {
            PushIllegalArgument();
            return;
        }
        if (bCumulative)
            PushDouble(integralPhi((x-mue)/sigma));
        else
            PushDouble(phi((x-mue)/sigma)/sigma);
    }

    void ScInterpreter::ScLogNormDist( int nMinParamCount ) //expanded, see #i100119# and fdo72158
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) )
            return;
        bool bCumulative = nParamCount != 4 || GetBool(); // cumulative
        double sigma = nParamCount >= 3 ? GetDouble() : 1.0; // standard deviation
        double mue = nParamCount >= 2 ? GetDouble() : 0.0;   // mean
        double x = GetDouble();                              // x
        if (sigma <= 0.0)
        {
            PushIllegalArgument();
            return;
        }
        if (bCumulative)
        { // cumulative
            if (x <= 0.0)
                PushDouble(0.0);
            else
                PushDouble(integralPhi((log(x)-mue)/sigma));
        }
        else
        { // density
            if (x <= 0.0)
                PushIllegalArgument();
            else
                PushDouble(phi((log(x)-mue)/sigma)/sigma/x);
        }
    }

    void ScInterpreter::ScStdNormDist()
    {
        PushDouble(integralPhi(GetDouble()));
    }

    void ScInterpreter::ScStdNormDist_MS()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, 2 ) )
            return;
        bool bCumulative = GetBool();                        // cumulative
        double x = GetDouble();                              // x

        if ( bCumulative )
            PushDouble( integralPhi( x ) );
        else
            PushDouble( exp( - pow( x, 2 ) / 2 ) / sqrt( 2 * F_PI ) );
    }

    void ScInterpreter::ScExpDist()
    {
        if ( MustHaveParamCount( GetByte(), 3 ) )
        {
            double kum    = GetDouble();                    // 0 or 1
            double lambda = GetDouble();                    // lambda
            double x      = GetDouble();                    // x
            if (lambda <= 0.0)
                PushIllegalArgument();
            else if (kum == 0.0)                        // density
            {
                if (x >= 0.0)
                    PushDouble(lambda * exp(-lambda*x));
                else
                    PushInt(0);
            }
            else                                        // distribution
            {
                if (x > 0.0)
                    PushDouble(1.0 - exp(-lambda*x));
                else
                    PushInt(0);
            }
        }
    }

    void ScInterpreter::ScTDist()
    {
        if ( !MustHaveParamCount( GetByte(), 3 ) )
            return;
        double fFlag = ::rtl::math::approxFloor(GetDouble());
        double fDF   = ::rtl::math::approxFloor(GetDouble());
        double T     = GetDouble();
        if (fDF < 1.0 || T < 0.0 || (fFlag != 1.0 && fFlag != 2.0) )
        {
            PushIllegalArgument();
            return;
        }
        PushDouble( GetTDist( T, fDF, static_cast<int>(fFlag) ) );
    }

    void ScInterpreter::ScTDist_T( int nTails )
    {
        if ( !MustHaveParamCount( GetByte(), 2 ) )
            return;
        double fDF = ::rtl::math::approxFloor( GetDouble() );
        double fT  = GetDouble();
        if ( fDF < 1.0 || ( nTails == 2 && fT < 0.0 ) )
        {
            PushIllegalArgument();
            return;
        }
        double fRes = GetTDist( fT, fDF, nTails );
        if ( nTails == 1 && fT < 0.0 )
            PushDouble( 1.0 - fRes ); // tdf#105937, right tail, negative X
        else
            PushDouble( fRes );
    }

    void ScInterpreter::ScTDist_MS()
    {
        if ( !MustHaveParamCount( GetByte(), 3 ) )
            return;
        bool   bCumulative = GetBool();
        double fDF = ::rtl::math::approxFloor( GetDouble() );
        double T   = GetDouble();
        if ( fDF < 1.0 )
        {
            PushIllegalArgument();
            return;
        }
        PushDouble( GetTDist( T, fDF, ( bCumulative ? 4 : 3 ) ) );
    }

    void ScInterpreter::ScFDist()
    {
        if ( !MustHaveParamCount( GetByte(), 3 ) )
            return;
        double fF2 = ::rtl::math::approxFloor(GetDouble());
        double fF1 = ::rtl::math::approxFloor(GetDouble());
        double fF  = GetDouble();
        if (fF < 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10)
        {
            PushIllegalArgument();
            return;
        }
        PushDouble(GetFDist(fF, fF1, fF2));
    }

    void ScInterpreter::ScFDist_LT()
    {
        int nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
            return;
        bool bCum;
        if ( nParamCount == 3 )
            bCum = true;
        else if ( IsMissing() )
        {
            bCum = true;
            Pop();
        }
        else
            bCum = GetBool();
        double fF2 = ::rtl::math::approxFloor( GetDouble() );
        double fF1 = ::rtl::math::approxFloor( GetDouble() );
        double fF  = GetDouble();
        if ( fF < 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 )
        {
            PushIllegalArgument();
            return;
        }
        if ( bCum )
        {
            // left tail cumulative distribution
            PushDouble( 1.0 - GetFDist( fF, fF1, fF2 ) );
        }
        else
        {
            // probability density function
            PushDouble( pow( fF1 / fF2, fF1 / 2 ) * pow( fF, ( fF1 / 2 ) - 1 ) /
                        ( pow( ( 1 + ( fF * fF1 / fF2 ) ), ( fF1 + fF2 ) / 2 ) *
                          GetBeta( fF1 / 2, fF2 / 2 ) ) );
        }
    }

    void ScInterpreter::ScChiDist( bool bODFF )
    {
        double fResult;
        if ( !MustHaveParamCount( GetByte(), 2 ) )
            return;
        double fDF  = ::rtl::math::approxFloor(GetDouble());
        double fChi = GetDouble();
        if ( fDF < 1.0 // x<=0 returns 1, see ODFF1.2 6.18.11
           || ( !bODFF && fChi < 0 ) ) // Excel does not accept negative fChi
        {
            PushIllegalArgument();
            return;
        }
        fResult = GetChiDist( fChi, fDF);
        if (nGlobalError != FormulaError::NONE)
        {
            PushError( nGlobalError);
            return;
        }
        PushDouble(fResult);
    }

    void ScInterpreter::ScWeibull()
    {
        if ( MustHaveParamCount( GetByte(), 4 ) )
        {
            double kum   = GetDouble();                 // 0 or 1
            double beta  = GetDouble();                 // beta
            double alpha = GetDouble();                 // alpha
            double x     = GetDouble();                 // x
            if (alpha <= 0.0 || beta <= 0.0 || x < 0.0)
                PushIllegalArgument();
            else if (kum == 0.0)                        // Density
                PushDouble(alpha/pow(beta,alpha)*pow(x,alpha-1.0)*
                           exp(-pow(x/beta,alpha)));
            else                                        // Distribution
                PushDouble(1.0 - exp(-pow(x/beta,alpha)));
        }
    }

    void ScInterpreter::ScPoissonDist( bool bODFF )
    {
        sal_uInt8 nParamCount = GetByte();
        if ( MustHaveParamCount( nParamCount, ( bODFF ? 2 : 3 ), 3 ) )
        {
            bool bCumulative = nParamCount != 3 || GetBool();         // default cumulative
            double lambda    = GetDouble();                           // Mean
            double x         = ::rtl::math::approxFloor(GetDouble()); // discrete distribution
            if (lambda <= 0.0 || x < 0.0)
                PushIllegalArgument();
            else if (!bCumulative)                            // Probability mass function
            {
                if (lambda >712.0)    // underflow in exp(-lambda)
                {   // accuracy 11 Digits
                    PushDouble( exp(x*log(lambda)-lambda-GetLogGamma(x+1.0)));
                }
                else
                {
                    double fPoissonVar = 1.0;
                    for ( double f = 0.0; f < x; ++f )
                        fPoissonVar *= lambda / ( f + 1.0 );
                    PushDouble( fPoissonVar * exp( -lambda ) );
                }
            }
            else                                // Cumulative distribution function
            {
                if (lambda > 712.0)  // underflow in exp(-lambda)
                {   // accuracy 12 Digits
                    PushDouble(GetUpRegIGamma(x+1.0,lambda));
                }
                else
                {
                    if (x >= 936.0) // result is always indistinguishable from 1
                        PushDouble (1.0);
                    else
                    {
                        double fSummand = exp(-lambda);
                        double 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);
                    }
                }
            }
        }
    }

    /** Local function used in the calculation of the hypergeometric distribution.
     */
    static void lcl_PutFactorialElements( ::std::vector< double >& cn, double fLower, double fUpper, double fBase )
    {
        for ( double i = fLower; i <= fUpper; ++i )
        {
            double fVal = fBase - i;
            if ( fVal > 1.0 )
                cn.push_back( fVal );
        }
    }

    /** Calculates a value of the hypergeometric distribution.

        @see #i47296#

        This function has an extra argument bCumulative,
        which only calculates the non-cumulative distribution and
        which is optional in Calc and mandatory with Excel's HYPGEOM.DIST()

        @see fdo#71722
        @see tdf#102948, make Calc function ODFF1.2-compliant
        @see tdf#117041, implement note at bottom of ODFF1.2 par.6.18.37
     */
    void ScInterpreter::ScHypGeomDist( int nMinParamCount )
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, nMinParamCount, 5 ) )
            return;

        bool bCumulative = ( nParamCount == 5 && GetBool() );
        double N = ::rtl::math::approxFloor(GetDouble());
        double M = ::rtl::math::approxFloor(GetDouble());
        double n = ::rtl::math::approxFloor(GetDouble());
        double x = ::rtl::math::approxFloor(GetDouble());

        if ( (x < 0.0) || (n < x) || (N < n) || (N < M) || (M < 0.0) )
        {
            PushIllegalArgument();
            return;
        }

        double fVal = 0.0;

        for ( int i = ( bCumulative ? 0 : x ); i <= x && nGlobalError == FormulaError::NONE; i++ )
        {
            if ( (n - i <= N - M) && (i <= M) )
                fVal +=  GetHypGeomDist( i, n, M, N );
        }

        PushDouble( fVal );
    }

    /** Calculates a value of the hypergeometric distribution.

        The algorithm is designed to avoid unnecessary multiplications and division
        by expanding all factorial elements (9 of them).  It is done by excluding
        those ranges that overlap in the numerator and the denominator.  This allows
        for a fast calculation for large values which would otherwise cause an overflow
        in the intermediate values.

        @see #i47296#
     */
    double ScInterpreter::GetHypGeomDist( double x, double n, double M, double N )
    {
        const size_t nMaxArraySize = 500000; // arbitrary max array size

        std::vector<double> cnNumer, cnDenom;

        size_t nEstContainerSize = static_cast<size_t>( x + ::std::min( n, M ) );
        size_t nMaxSize = ::std::min( cnNumer.max_size(), nMaxArraySize );
        if ( nEstContainerSize > nMaxSize )
        {
            PushNoValue();
            return 0;
        }
        cnNumer.reserve( nEstContainerSize + 10 );
        cnDenom.reserve( nEstContainerSize + 10 );

        // Trim coefficient C first
        double fCNumVarUpper = N - n - M + x - 1.0;
        double fCDenomVarLower = 1.0;
        if ( N - n - M + x >= M - x + 1.0 )
        {
            fCNumVarUpper = M - x - 1.0;
            fCDenomVarLower = N - n - 2.0*(M - x) + 1.0;
        }

        double fCNumLower = N - n - fCNumVarUpper;
        double fCDenomUpper = N - n - M + x + 1.0 - fCDenomVarLower;

        double fDNumVarLower = n - M;

        if ( n >= M + 1.0 )
        {
            if ( N - M < n + 1.0 )
            {
                // Case 1

                if ( N - n < n + 1.0 )
                {
                    // no overlap
                    lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
                    lcl_PutFactorialElements( cnDenom, 0.0, N - n - 1.0, N );
                }
                else
                {
                    // overlap
                    OSL_ENSURE( fCNumLower < n + 1.0, "ScHypGeomDist: wrong assertion" );
                    lcl_PutFactorialElements( cnNumer, N - 2.0*n, fCNumVarUpper, N - n );
                    lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
                }

                OSL_ENSURE( fCDenomUpper <= N - M, "ScHypGeomDist: wrong assertion" );

                if ( fCDenomUpper < n - x + 1.0 )
                    // no overlap
                    lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 );
                else
                {
                    // overlap
                    lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 );

                    fCDenomUpper = n - x;
                    fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
                }
            }
            else
            {
                // Case 2

                if ( n > M - 1.0 )
                {
                    // no overlap
                    lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
                    lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N );
                }
                else
                {
                    lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n );
                    lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
                }

                OSL_ENSURE( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" );

                if ( fCDenomUpper < n - x + 1.0 )
                    // no overlap
                    lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - n + x, N - M + 1.0 );
                else
                {
                    lcl_PutFactorialElements( cnNumer, N - M - n + 1.0, N - M - fCDenomUpper, N - M + 1.0 );
                    fCDenomUpper = n - x;
                    fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
                }
            }

            OSL_ENSURE( fCDenomUpper <= M, "ScHypGeomDist: wrong assertion" );
        }
        else
        {
            if ( N - M < M + 1.0 )
            {
                // Case 3

                if ( N - n < M + 1.0 )
                {
                    // No overlap
                    lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
                    lcl_PutFactorialElements( cnDenom, 0.0, N - M - 1.0, N );
                }
                else
                {
                    lcl_PutFactorialElements( cnNumer, N - n - M, fCNumVarUpper, N - n );
                    lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
                }

                if ( n - x + 1.0 > fCDenomUpper )
                    // No overlap
                    lcl_PutFactorialElements( cnNumer, 1.0, N - M - n + x, N - M + 1.0 );
                else
                {
                    // Overlap
                    lcl_PutFactorialElements( cnNumer, 1.0, N - M - fCDenomUpper, N - M + 1.0 );

                    fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
                    fCDenomUpper = n - x;
                }
            }
            else
            {
                // Case 4

                OSL_ENSURE( M >= n - x, "ScHypGeomDist: wrong assertion" );
                OSL_ENSURE( M - x <= N - M + 1.0, "ScHypGeomDist: wrong assertion" );

                if ( N - n < N - M + 1.0 )
                {
                    // No overlap
                    lcl_PutFactorialElements( cnNumer, 0.0, fCNumVarUpper, N - n );
                    lcl_PutFactorialElements( cnDenom, 0.0, M - 1.0, N );
                }
                else
                {
                    // Overlap
                    OSL_ENSURE( fCNumLower <= N - M + 1.0, "ScHypGeomDist: wrong assertion" );
                    lcl_PutFactorialElements( cnNumer, M - n, fCNumVarUpper, N - n );
                    lcl_PutFactorialElements( cnDenom, 0.0, n - 1.0, N );
                }

                if ( n - x + 1.0 > fCDenomUpper )
                    // No overlap
                    lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - n + x, N - M + 1.0 );
                else if ( M >= fCDenomUpper )
                {
                    lcl_PutFactorialElements( cnNumer, N - 2.0*M + 1.0, N - M - fCDenomUpper, N - M + 1.0 );

                    fCDenomUpper = n - x;
                    fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
                }
                else
                {
                    OSL_ENSURE( M <= fCDenomUpper, "ScHypGeomDist: wrong assertion" );
                    lcl_PutFactorialElements( cnDenom, fCDenomVarLower, N - n - 2.0*M + x,
                            N - n - M + x + 1.0 );

                    fCDenomUpper = n - x;
                    fCDenomVarLower = N - M - 2.0*(n - x) + 1.0;
                }
            }

            OSL_ENSURE( fCDenomUpper <= n, "ScHypGeomDist: wrong assertion" );

            fDNumVarLower = 0.0;
        }

        double nDNumVarUpper   = fCDenomUpper < x + 1.0 ? n - x - 1.0     : n - fCDenomUpper - 1.0;
        double nDDenomVarLower = fCDenomUpper < x + 1.0 ? fCDenomVarLower : N - n - M + 1.0;
        lcl_PutFactorialElements( cnNumer, fDNumVarLower, nDNumVarUpper, n );
        lcl_PutFactorialElements( cnDenom, nDDenomVarLower, N - n - M + x, N - n - M + x + 1.0 );

        ::std::sort( cnNumer.begin(), cnNumer.end() );
        ::std::sort( cnDenom.begin(), cnDenom.end() );
        auto it1 = cnNumer.rbegin(), it1End = cnNumer.rend();
        auto it2 = cnDenom.rbegin(), it2End = cnDenom.rend();

        double fFactor = 1.0;
        for ( ; it1 != it1End || it2 != it2End; )
        {
            double fEnum = 1.0, fDenom = 1.0;
            if ( it1 != it1End )
                fEnum  = *it1++;
            if ( it2 != it2End )
                fDenom = *it2++;
            fFactor *= fEnum / fDenom;
        }

        return fFactor;
    }

    void ScInterpreter::ScGammaDist( bool bODFF )
    {
        sal_uInt8 nMinParamCount = ( bODFF ? 3 : 4 );
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, nMinParamCount, 4 ) )
            return;
        bool bCumulative;
        if (nParamCount == 4)
            bCumulative = GetBool();
        else
            bCumulative = true;
        double fBeta = GetDouble();                 // scale
        double fAlpha = GetDouble();                // shape
        double fX = GetDouble();                    // x
        if ((!bODFF && fX < 0) || fAlpha <= 0.0 || fBeta <= 0.0)
            PushIllegalArgument();
        else
        {
            if (bCumulative)                        // distribution
                PushDouble( GetGammaDist( fX, fAlpha, fBeta));
            else                                    // density
                PushDouble( GetGammaDistPDF( fX, fAlpha, fBeta));
        }
    }

    void ScInterpreter::ScNormInv()
    {
        if ( MustHaveParamCount( GetByte(), 3 ) )
        {
            double sigma = GetDouble();
            double mue   = GetDouble();
            double x     = GetDouble();
            if (sigma <= 0.0 || x < 0.0 || x > 1.0)
                PushIllegalArgument();
            else if (x == 0.0 || x == 1.0)
                PushNoValue();
            else
                PushDouble(gaussinv(x)*sigma + mue);
        }
    }

    void ScInterpreter::ScSNormInv()
    {
        double x = GetDouble();
        if (x < 0.0 || x > 1.0)
            PushIllegalArgument();
        else if (x == 0.0 || x == 1.0)
            PushNoValue();
        else
            PushDouble(gaussinv(x));
    }

    void ScInterpreter::ScLogNormInv()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( MustHaveParamCount( nParamCount, 1, 3 ) )
        {
            double fSigma = ( nParamCount == 3 ? GetDouble() : 1.0 );  // Stddev
            double fMue = ( nParamCount >= 2 ? GetDouble() : 0.0 );    // Mean
            double fP = GetDouble();                                   // p
            if ( fSigma <= 0.0 || fP <= 0.0 || fP >= 1.0 )
                PushIllegalArgument();
            else
                PushDouble( exp( fMue + fSigma * gaussinv( fP ) ) );
        }
    }

    class ScGammaDistFunction : public ScDistFunc
    {
        ScInterpreter&  rInt;
        double          fp, fAlpha, fBeta;

    public:
                ScGammaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
                    rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}

        virtual ~ScGammaDistFunction() {}

        double  GetValue( double x ) const override  { return fp - rInt.GetGammaDist(x, fAlpha, fBeta); }
    };

    void ScInterpreter::ScGammaInv()
    {
        if ( !MustHaveParamCount( GetByte(), 3 ) )
            return;
        double fBeta  = GetDouble();
        double fAlpha = GetDouble();
        double fP = GetDouble();
        if (fAlpha <= 0.0 || fBeta <= 0.0 || fP < 0.0 || fP >= 1.0 )
        {
            PushIllegalArgument();
            return;
        }
        if (fP == 0.0)
            PushInt(0);
        else
        {
            bool bConvError;
            ScGammaDistFunction aFunc( *this, fP, fAlpha, fBeta );
            double fStart = fAlpha * fBeta;
            double fVal = lcl_IterateInverse( aFunc, fStart*0.5, fStart, bConvError );
            if (bConvError)
                SetError(FormulaError::NoConvergence);
            PushDouble(fVal);
        }
    }

    class ScBetaDistFunction : public ScDistFunc
    {
        ScInterpreter&  rInt;
        double          fp, fAlpha, fBeta;

    public:
                ScBetaDistFunction( ScInterpreter& rI, double fpVal, double fAlphaVal, double fBetaVal ) :
                    rInt(rI), fp(fpVal), fAlpha(fAlphaVal), fBeta(fBetaVal) {}

        virtual ~ScBetaDistFunction() {}

        double  GetValue( double x ) const override  { return fp - rInt.GetBetaDist(x, fAlpha, fBeta); }
    };

    void ScInterpreter::ScBetaInv()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, 3, 5 ) )
            return;
        double fP, fA, fB, fAlpha, fBeta;
        if (nParamCount == 5)
            fB = GetDouble();
        else
            fB = 1.0;
        if (nParamCount >= 4)
            fA = GetDouble();
        else
            fA = 0.0;
        fBeta  = GetDouble();
        fAlpha = GetDouble();
        fP     = GetDouble();
        if (fP < 0.0 || fP > 1.0 || fA >= fB || fAlpha <= 0.0 || fBeta <= 0.0)
        {
            PushIllegalArgument();
            return;
        }

        bool bConvError;
        ScBetaDistFunction aFunc( *this, fP, fAlpha, fBeta );
        // 0..1 as range for iteration so it isn't extended beyond the valid range
        double fVal = lcl_IterateInverse( aFunc, 0.0, 1.0, bConvError );
        if (bConvError)
            PushError( FormulaError::NoConvergence);
        else
            PushDouble(fA + fVal*(fB-fA));                  // scale to (A,B)
    }

    // Note: T, F, and Chi are
    // monotonically decreasing,
    // therefore 1-Dist as function

    class ScTDistFunction : public ScDistFunc
    {
        ScInterpreter&  rInt;
        double          fp, fDF;
        int             nT;

    public:
                ScTDistFunction( ScInterpreter& rI, double fpVal, double fDFVal, int nType ) :
                    rInt( rI ), fp( fpVal ), fDF( fDFVal ), nT( nType ) {}

        virtual ~ScTDistFunction() {}

        double  GetValue( double x ) const override  { return fp - rInt.GetTDist( x, fDF, nT ); }
    };

    void ScInterpreter::ScTInv( int nType )
    {
        if ( !MustHaveParamCount( GetByte(), 2 ) )
            return;
        double fDF  = ::rtl::math::approxFloor(GetDouble());
        double fP = GetDouble();
        if (fDF < 1.0 || fP <= 0.0 || fP > 1.0 )
        {
            PushIllegalArgument();
            return;
        }
        if ( nType == 4 ) // left-tailed cumulative t-distribution
        {
            if ( fP == 1.0 )
                PushIllegalArgument();
            else if ( fP < 0.5 )
                PushDouble( -GetTInv( 1 - fP, fDF, nType ) );
            else
                PushDouble( GetTInv( fP, fDF, nType ) );
        }
        else
            PushDouble( GetTInv( fP, fDF, nType ) );
    };

    double ScInterpreter::GetTInv( double fAlpha, double fSize, int nType )
    {
        bool bConvError;
        ScTDistFunction aFunc( *this, fAlpha, fSize, nType );
        double fVal = lcl_IterateInverse( aFunc, fSize * 0.5, fSize, bConvError );
        if (bConvError)
            SetError(FormulaError::NoConvergence);
        return fVal;
    }

    class ScFDistFunction : public ScDistFunc
    {
        ScInterpreter&  rInt;
        double          fp, fF1, fF2;

    public:
                ScFDistFunction( ScInterpreter& rI, double fpVal, double fF1Val, double fF2Val ) :
                    rInt(rI), fp(fpVal), fF1(fF1Val), fF2(fF2Val) {}

        virtual ~ScFDistFunction() {}

        double  GetValue( double x ) const override  { return fp - rInt.GetFDist(x, fF1, fF2); }
    };

    void ScInterpreter::ScFInv()
    {
        if ( !MustHaveParamCount( GetByte(), 3 ) )
            return;
        double fF2 = ::rtl::math::approxFloor(GetDouble());
        double fF1 = ::rtl::math::approxFloor(GetDouble());
        double fP  = GetDouble();
        if (fP <= 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 || fP > 1.0)
        {
            PushIllegalArgument();
            return;
        }

        bool bConvError;
        ScFDistFunction aFunc( *this, fP, fF1, fF2 );
        double fVal = lcl_IterateInverse( aFunc, fF1*0.5, fF1, bConvError );
        if (bConvError)
            SetError(FormulaError::NoConvergence);
        PushDouble(fVal);
    }

    void ScInterpreter::ScFInv_LT()
    {
        if ( !MustHaveParamCount( GetByte(), 3 ) )
            return;
        double fF2 = ::rtl::math::approxFloor(GetDouble());
        double fF1 = ::rtl::math::approxFloor(GetDouble());
        double fP  = GetDouble();
        if (fP <= 0.0 || fF1 < 1.0 || fF2 < 1.0 || fF1 >= 1.0E10 || fF2 >= 1.0E10 || fP > 1.0)
        {
            PushIllegalArgument();
            return;
        }

        bool bConvError;
        ScFDistFunction aFunc( *this, ( 1.0 - fP ), fF1, fF2 );
        double fVal = lcl_IterateInverse( aFunc, fF1*0.5, fF1, bConvError );
        if (bConvError)
            SetError(FormulaError::NoConvergence);
        PushDouble(fVal);
    }

    class ScChiDistFunction : public ScDistFunc
    {
        ScInterpreter&  rInt;
        double          fp, fDF;

    public:
                ScChiDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
                    rInt(rI), fp(fpVal), fDF(fDFVal) {}

        virtual ~ScChiDistFunction() {}

        double  GetValue( double x ) const override  { return fp - rInt.GetChiDist(x, fDF); }
    };

    void ScInterpreter::ScChiInv()
    {
        if ( !MustHaveParamCount( GetByte(), 2 ) )
            return;
        double fDF  = ::rtl::math::approxFloor(GetDouble());
        double fP = GetDouble();
        if (fDF < 1.0 || fP <= 0.0 || fP > 1.0 )
        {
            PushIllegalArgument();
            return;
        }

        bool bConvError;
        ScChiDistFunction aFunc( *this, fP, fDF );
        double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
        if (bConvError)
            SetError(FormulaError::NoConvergence);
        PushDouble(fVal);
    }

    /***********************************************/
    class ScChiSqDistFunction : public ScDistFunc
    {
        ScInterpreter&  rInt;
        double          fp, fDF;

    public:
                ScChiSqDistFunction( ScInterpreter& rI, double fpVal, double fDFVal ) :
                    rInt(rI), fp(fpVal), fDF(fDFVal) {}

        virtual ~ScChiSqDistFunction() {}

        double  GetValue( double x ) const override  { return fp - rInt.GetChiSqDistCDF(x, fDF); }
    };

    void ScInterpreter::ScChiSqInv()
    {
        if ( !MustHaveParamCount( GetByte(), 2 ) )
            return;
        double fDF  = ::rtl::math::approxFloor(GetDouble());
        double fP = GetDouble();
        if (fDF < 1.0 || fP < 0.0 || fP >= 1.0 )
        {
            PushIllegalArgument();
            return;
        }

        bool bConvError;
        ScChiSqDistFunction aFunc( *this, fP, fDF );
        double fVal = lcl_IterateInverse( aFunc, fDF*0.5, fDF, bConvError );
        if (bConvError)
            SetError(FormulaError::NoConvergence);
        PushDouble(fVal);
    }

    void ScInterpreter::ScConfidence()
    {
        if ( MustHaveParamCount( GetByte(), 3 ) )
        {
            double n     = ::rtl::math::approxFloor(GetDouble());
            double sigma = GetDouble();
            double alpha = GetDouble();
            if (sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1.0)
                PushIllegalArgument();
            else
                PushDouble( gaussinv(1.0-alpha/2.0) * sigma/sqrt(n) );
        }
    }

    void ScInterpreter::ScConfidenceT()
    {
        if ( MustHaveParamCount( GetByte(), 3 ) )
        {
            double n     = ::rtl::math::approxFloor(GetDouble());
            double sigma = GetDouble();
            double alpha = GetDouble();
            if (sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1.0)
                PushIllegalArgument();
            else if (n == 1.0) // for interoperability with Excel
                PushError(FormulaError::DivisionByZero);
            else
                PushDouble( sigma * GetTInv( alpha, n - 1, 2 ) / sqrt( n ) );
        }
    }

    void ScInterpreter::ScZTest()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
            return;
        double sigma = 0.0, x;
        if (nParamCount == 3)
        {
            sigma = GetDouble();
            if (sigma <= 0.0)
            {
                PushIllegalArgument();
                return;
            }
        }
        x = GetDouble();

        double fSum    = 0.0;
        double fSumSqr = 0.0;
        double fVal;
        double rValCount = 0.0;
        switch (GetStackType())
        {
            case svDouble :
            {
                fVal = GetDouble();
                fSum    += fVal;
                fSumSqr += fVal*fVal;
                rValCount++;
            }
            break;
            case svSingleRef :
            {
                ScAddress aAdr;
                PopSingleRef( aAdr );
                ScRefCellValue aCell(*pDok, aAdr);
                if (aCell.hasNumeric())
                {
                    fVal = GetCellValue(aAdr, aCell);
                    fSum += fVal;
                    fSumSqr += fVal*fVal;
                    rValCount++;
                }
            }
            break;
            case svRefList :
            case svDoubleRef :
            {
                short nParam = 1;
                size_t nRefInList = 0;
                while (nParam-- > 0)
                {
                    ScRange aRange;
                    FormulaError nErr = FormulaError::NONE;
                    PopDoubleRef( aRange, nParam, nRefInList);
                    ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
                    if (aValIter.GetFirst(fVal, nErr))
                    {
                        fSum += fVal;
                        fSumSqr += fVal*fVal;
                        rValCount++;
                        while ((nErr == FormulaError::NONE) && aValIter.GetNext(fVal, nErr))
                        {
                            fSum += fVal;
                            fSumSqr += fVal*fVal;
                            rValCount++;
                        }
                        SetError(nErr);
                    }
                }
            }
            break;
            case svMatrix :
            case svExternalSingleRef:
            case svExternalDoubleRef:
            {
                ScMatrixRef pMat = GetMatrix();
                if (pMat)
                {
                    SCSIZE nCount = pMat->GetElementCount();
                    if (pMat->IsNumeric())
                    {
                        for ( SCSIZE i = 0; i < nCount; i++ )
                        {
                            fVal= pMat->GetDouble(i);
                            fSum += fVal;
                            fSumSqr += fVal * fVal;
                            rValCount++;
                        }
                    }
                    else
                    {
                        for (SCSIZE i = 0; i < nCount; i++)
                            if (!pMat->IsStringOrEmpty(i))
                            {
                                fVal= pMat->GetDouble(i);
                                fSum += fVal;
                                fSumSqr += fVal * fVal;
                                rValCount++;
                            }
                    }
                }
            }
            break;
            default : SetError(FormulaError::IllegalParameter); break;
        }
        if (rValCount <= 1.0)
            PushError( FormulaError::DivisionByZero);
        else
        {
            double mue = fSum/rValCount;

            if (nParamCount != 3)
            {
                sigma = (fSumSqr - fSum*fSum/rValCount)/(rValCount-1.0);
                if (sigma == 0.0)
                {
                    PushError(FormulaError::DivisionByZero);
                    return;
                }
                PushDouble(0.5 - gauss((mue-x)/sqrt(sigma/rValCount)));
            }
            else
                PushDouble(0.5 - gauss((mue-x)*sqrt(rValCount)/sigma));
        }
    }

    bool ScInterpreter::CalculateTest(bool _bTemplin
                                      ,const SCSIZE nC1, const SCSIZE nC2,const SCSIZE nR1,const SCSIZE nR2
                                      ,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 fVal;
        SCSIZE i,j;
        for (i = 0; i < nC1; i++)
            for (j = 0; j < nR1; j++)
            {
                if (!pMat1->IsStringOrEmpty(i,j))
                {
                    fVal = pMat1->GetDouble(i,j);
                    fSum1    += fVal;
                    fSumSqr1 += fVal * fVal;
                    fCount1++;
                }
            }
        for (i = 0; i < nC2; i++)
            for (j = 0; j < nR2; j++)
            {
                if (!pMat2->IsStringOrEmpty(i,j))
                {
                    fVal = pMat2->GetDouble(i,j);
                    fSum2    += fVal;
                    fSumSqr2 += fVal * fVal;
                    fCount2++;
                }
            }
        if (fCount1 < 2.0 || fCount2 < 2.0)
        {
            PushNoValue();
            return false;
        } // 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;
            if (fS1 + fS2 == 0.0)
            {
                PushNoValue();
                return false;
            }
            fT = fabs(fSum1/fCount1 - fSum2/fCount2)/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
            fF = 1.0/(c*c/(fCount1-1.0)+(1.0-c)*(1.0-c)/(fCount2-1.0));
        }
        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 = fabs( fSum1/fCount1 - fSum2/fCount2 ) /
                 sqrt( (fCount1-1.0)*fS1 + (fCount2-1.0)*fS2 ) *
                 sqrt( fCount1*fCount2*(fCount1+fCount2-2)/(fCount1+fCount2) );
            fF = fCount1 + fCount2 - 2;
        }
        return true;
    }
    void ScInterpreter::ScTTest()
    {
        if ( !MustHaveParamCount( GetByte(), 4 ) )
            return;
        double fTyp   = ::rtl::math::approxFloor(GetDouble());
        double fTails = ::rtl::math::approxFloor(GetDouble());
        if (fTails != 1.0 && fTails != 2.0)
        {
            PushIllegalArgument();
            return;
        }

        ScMatrixRef pMat2 = GetMatrix();
        ScMatrixRef pMat1 = GetMatrix();
        if (!pMat1 || !pMat2)
        {
            PushIllegalParameter();
            return;
        }
        double fT, fF;
        SCSIZE nC1, nC2;
        SCSIZE nR1, nR2;
        SCSIZE i, j;
        pMat1->GetDimensions(nC1, nR1);
        pMat2->GetDimensions(nC2, nR2);
        if (fTyp == 1.0)
        {
            if (nC1 != nC2 || nR1 != nR2)
            {
                PushIllegalArgument();
                return;
            }
            double fCount   = 0.0;
            double fSum1    = 0.0;
            double fSum2    = 0.0;
            double fSumSqrD = 0.0;
            double fVal1, fVal2;
            for (i = 0; i < nC1; i++)
                for (j = 0; j < nR1; j++)
                {
                    if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
                    {
                        fVal1 = pMat1->GetDouble(i,j);
                        fVal2 = pMat2->GetDouble(i,j);
                        fSum1    += fVal1;
                        fSum2    += fVal2;
                        fSumSqrD += (fVal1 - fVal2)*(fVal1 - fVal2);
                        fCount++;
                    }
                }
            if (fCount < 1.0)
            {
                PushNoValue();
                return;
            }
            double fSumD = fSum1 - fSum2;
            double fDivider = (fCount*fSumSqrD - fSumD*fSumD);
            if ( fDivider == 0.0 )
            {
                PushError(FormulaError::DivisionByZero);
                return;
            }
            fT = fabs(fSumD) * sqrt((fCount-1.0) / fDivider);
            fF = fCount - 1.0;
        }
        else if (fTyp == 2.0)
        {
            if (!CalculateTest(false,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF))
                return;     // error was pushed
        }
        else if (fTyp == 3.0)
        {
            if (!CalculateTest(true,nC1, nC2,nR1, nR2,pMat1,pMat2,fT,fF))
                return;     // error was pushed
        }
        else
        {
            PushIllegalArgument();
            return;
        }
        PushDouble( GetTDist( fT, fF, static_cast<int>(fTails) ) );
    }

    void ScInterpreter::ScFTest()
    {
        if ( !MustHaveParamCount( GetByte(), 2 ) )
            return;
        ScMatrixRef pMat2 = GetMatrix();
        ScMatrixRef pMat1 = GetMatrix();
        if (!pMat1 || !pMat2)
        {
            PushIllegalParameter();
            return;
        }
        SCSIZE nC1, nC2;
        SCSIZE nR1, nR2;
        pMat1->GetDimensions(nC1, nR1);
        pMat2->GetDimensions(nC2, nR2);
        double fCount1  = 0.0;
        double fCount2  = 0.0;
        double fSum1    = 0.0;
        double fSumSqr1 = 0.0;
        double fSum2    = 0.0;
        double fSumSqr2 = 0.0;

        std::vector<std::unique_ptr<sc::op::Op>> aOp;
        aOp.emplace_back(new sc::op::Op(0.0, [](double& rAccum, double fVal){rAccum += fVal;}));
        aOp.emplace_back(new sc::op::Op(0.0, [](double& rAccum, double fVal){rAccum += fVal * fVal;}));

        auto aVal1 = pMat1->Collect(aOp);
        fSum1 = aVal1[0].mfFirst + aVal1[0].mfRest;
        fSumSqr1 = aVal1[1].mfFirst + aVal1[1].mfRest;
        fCount1 = aVal1[2].mnCount;

        auto aVal2 = pMat2->Collect(aOp);
        fSum2 = aVal2[0].mfFirst + aVal2[0].mfRest;
        fSumSqr2 = aVal2[1].mfFirst + aVal2[1].mfRest;
        fCount2 = aVal2[2].mnCount;

        if (fCount1 < 2.0 || fCount2 < 2.0)
        {
            PushNoValue();
            return;
        }
        double fS1 = (fSumSqr1-fSum1*fSum1/fCount1)/(fCount1-1.0);
        double fS2 = (fSumSqr2-fSum2*fSum2/fCount2)/(fCount2-1.0);
        if (fS1 == 0.0 || fS2 == 0.0)
        {
            PushNoValue();
            return;
        }
        double fF, fF1, fF2;
        if (fS1 > fS2)
        {
            fF = fS1/fS2;
            fF1 = fCount1-1.0;
            fF2 = fCount2-1.0;
        }
        else
        {
            fF = fS2/fS1;
            fF1 = fCount2-1.0;
            fF2 = fCount1-1.0;
        }
        double fFcdf = GetFDist(fF, fF1, fF2);
        PushDouble(2.0*std::min(fFcdf, 1.0 - fFcdf));
    }

    void ScInterpreter::ScChiTest()
    {
        if ( !MustHaveParamCount( GetByte(), 2 ) )
            return;
        ScMatrixRef pMat2 = GetMatrix();
        ScMatrixRef pMat1 = GetMatrix();
        if (!pMat1 || !pMat2)
        {
            PushIllegalParameter();
            return;
        }
        SCSIZE nC1, nC2;
        SCSIZE nR1, nR2;
        pMat1->GetDimensions(nC1, nR1);
        pMat2->GetDimensions(nC2, nR2);
        if (nR1 != nR2 || nC1 != nC2)
        {
            PushIllegalArgument();
            return;
        }
        double fChi = 0.0;
        bool bEmpty = true;
        for (SCSIZE i = 0; i < nC1; i++)
        {
            for (SCSIZE j = 0; j < nR1; j++)
            {
                if (!(pMat1->IsEmpty(i,j) || pMat2->IsEmpty(i,j)))
                {
                    bEmpty = false;
                    if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
                    {
                        double fValX = pMat1->GetDouble(i,j);
                        double fValE = pMat2->GetDouble(i,j);
                        if ( fValE == 0.0 )
                        {
                            PushError(FormulaError::DivisionByZero);
                            return;
                        }
                        // These fTemp values guard against a failure when compiled
                        // with optimization (using g++ 4.8.2 on tinderbox 71-TDF),
                        // where ((fValX - fValE) * (fValX - fValE)) with
                        // fValE==1e+308 should had produced Infinity but did
                        // not, instead the result of divide() then was 1e+308.
                        volatile double fTemp1 = (fValX - fValE) * (fValX - fValE);
                        double fTemp2 = fTemp1;
                        fChi += sc::divide( fTemp2, fValE);
                    }
                    else
                    {
                        PushIllegalArgument();
                        return;
                    }
                }
            }
        }
        if ( bEmpty )
        {
            // not in ODFF1.2, but for interoperability with Excel
            PushIllegalArgument();
            return;
        }
        double fDF;
        if (nC1 == 1 || nR1 == 1)
        {
            fDF = static_cast<double>(nC1*nR1 - 1);
            if (fDF == 0.0)
            {
                PushNoValue();
                return;
            }
        }
        else
            fDF = static_cast<double>(nC1-1)*static_cast<double>(nR1-1);
        PushDouble(GetChiDist(fChi, fDF));
    }

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

        if (fCount == 0.0)
        {
            PushError( FormulaError::DivisionByZero);
            return;
        }

        double fMean = fSum / fCount;

        for (double v : values)
            vSum += (v - fMean) * (v - fMean);

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

        if (fStdDev == 0.0)
        {
            PushError( FormulaError::DivisionByZero);
            return;
        }

        for (double v : values)
        {
            double dx = (v - fMean) / fStdDev;
            xpower4 = 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);
    }

    void ScInterpreter::ScHarMean()
    {
        short nParamCount = GetByte();
        double nVal = 0.0;
        double nValCount = 0.0;
        ScAddress aAdr;
        ScRange aRange;
        size_t nRefInList = 0;
        while ((nGlobalError == FormulaError::NONE) && (nParamCount-- > 0))
        {
            switch (GetStackType())
            {
                case svDouble    :
                {
                    double x = GetDouble();
                    if (x > 0.0)
                    {
                        nVal += 1.0/x;
                        nValCount++;
                    }
                    else
                        SetError( FormulaError::IllegalArgument);
                    break;
                }
                case svSingleRef :
                {
                    PopSingleRef( aAdr );
                    ScRefCellValue aCell(*pDok, aAdr);
                    if (aCell.hasNumeric())
                    {
                        double x = GetCellValue(aAdr, aCell);
                        if (x > 0.0)
                        {
                            nVal += 1.0/x;
                            nValCount++;
                        }
                        else
                            SetError( FormulaError::IllegalArgument);
                    }
                    break;
                }
                case svDoubleRef :
                case svRefList :
                {
                    FormulaError nErr = FormulaError::NONE;
                    PopDoubleRef( aRange, nParamCount, nRefInList);
                    double nCellVal;
                    ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
                    if (aValIter.GetFirst(nCellVal, nErr))
                    {
                        if (nCellVal > 0.0)
                        {
                            nVal += 1.0/nCellVal;
                            nValCount++;
                        }
                        else
                            SetError( FormulaError::IllegalArgument);
                        SetError(nErr);
                        while ((nErr == FormulaError::NONE) && aValIter.GetNext(nCellVal, nErr))
                        {
                            if (nCellVal > 0.0)
                            {
                                nVal += 1.0/nCellVal;
                                nValCount++;
                            }
                            else
                                SetError( FormulaError::IllegalArgument);
                        }
                        SetError(nErr);
                    }
                }
                break;
                case svMatrix :
                case svExternalSingleRef:
                case svExternalDoubleRef:
                {
                    ScMatrixRef pMat = GetMatrix();
                    if (pMat)
                    {
                        SCSIZE nCount = pMat->GetElementCount();
                        if (pMat->IsNumeric())
                        {
                            for (SCSIZE nElem = 0; nElem < nCount; nElem++)
                            {
                                double x = pMat->GetDouble(nElem);
                                if (x > 0.0)
                                {
                                    nVal += 1.0/x;
                                    nValCount++;
                                }
                                else
                                    SetError( FormulaError::IllegalArgument);
                            }
                        }
                        else
                        {
                            for (SCSIZE nElem = 0; nElem < nCount; nElem++)
                                if (!pMat->IsStringOrEmpty(nElem))
                                {
                                    double x = pMat->GetDouble(nElem);
                                    if (x > 0.0)
                                    {
                                        nVal += 1.0/x;
                                        nValCount++;
                                    }
                                    else
                                        SetError( FormulaError::IllegalArgument);
                                }
                        }
                    }
                }
                break;
                default : SetError(FormulaError::IllegalParameter); break;
            }
        }
        if (nGlobalError == FormulaError::NONE)
            PushDouble( nValCount / nVal );
        else
            PushError( nGlobalError);
    }

    void ScInterpreter::ScGeoMean()
    {
        short nParamCount = GetByte();
        double nVal = 0.0;
        double nValCount = 0.0;
        ScAddress aAdr;
        ScRange aRange;

        size_t nRefInList = 0;
        while ((nGlobalError == FormulaError::NONE) && (nParamCount-- > 0))
        {
            switch (GetStackType())
            {
                case svDouble    :
                {
                    double x = GetDouble();
                    if (x > 0.0)
                    {
                        nVal += log(x);
                        nValCount++;
                    }
                    else if ( x == 0.0 )
                    {
                        // value of 0 means that function result will be 0
                        while ( nParamCount-- > 0 )
                            PopError();
                        PushDouble( 0.0 );
                        return;
                    }
                    else
                        SetError( FormulaError::IllegalArgument);
                    break;
                }
                case svSingleRef :
                {
                    PopSingleRef( aAdr );
                    ScRefCellValue aCell(*pDok, aAdr);
                    if (aCell.hasNumeric())
                    {
                        double x = GetCellValue(aAdr, aCell);
                        if (x > 0.0)
                        {
                            nVal += log(x);
                            nValCount++;
                        }
                        else if ( x == 0.0 )
                        {
                            // value of 0 means that function result will be 0
                            while ( nParamCount-- > 0 )
                                PopError();
                            PushDouble( 0.0 );
                            return;
                        }
                        else
                            SetError( FormulaError::IllegalArgument);
                    }
                    break;
                }
                case svDoubleRef :
                case svRefList :
                {
                    FormulaError nErr = FormulaError::NONE;
                    PopDoubleRef( aRange, nParamCount, nRefInList);
                    double nCellVal;
                    ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
                    if (aValIter.GetFirst(nCellVal, nErr))
                    {
                        if (nCellVal > 0.0)
                        {
                            nVal += log(nCellVal);
                            nValCount++;
                        }
                        else if ( nCellVal == 0.0 )
                        {
                            // value of 0 means that function result will be 0
                            while ( nParamCount-- > 0 )
                                PopError();
                            PushDouble( 0.0 );
                            return;
                        }
                        else
                            SetError( FormulaError::IllegalArgument);
                        SetError(nErr);
                        while ((nErr == FormulaError::NONE) && aValIter.GetNext(nCellVal, nErr))
                        {
                            if (nCellVal > 0.0)
                            {
                                nVal += log(nCellVal);
                                nValCount++;
                            }
                            else if ( nCellVal == 0.0 )
                            {
                                // value of 0 means that function result will be 0
                                while ( nParamCount-- > 0 )
                                    PopError();
                                PushDouble( 0.0 );
                                return;
                            }
                            else
                                SetError( FormulaError::IllegalArgument);
                        }
                        SetError(nErr);
                    }
                }
                break;
                case svMatrix :
                case svExternalSingleRef:
                case svExternalDoubleRef:
                {
                    ScMatrixRef pMat = GetMatrix();
                    if (pMat)
                    {
                        SCSIZE nCount = pMat->GetElementCount();
                        if (pMat->IsNumeric())
                        {
                            for (SCSIZE ui = 0; ui < nCount; ui++)
                            {
                                double x = pMat->GetDouble(ui);
                                if (x > 0.0)
                                {
                                    nVal += log(x);
                                    nValCount++;
                                }
                                else if ( x == 0.0 )
                                {
                                    // value of 0 means that function result will be 0
                                    while ( nParamCount-- > 0 )
                                        PopError();
                                    PushDouble( 0.0 );
                                    return;
                                }
                                else
                                    SetError( FormulaError::IllegalArgument);
                            }
                        }
                        else
                        {
                            for (SCSIZE ui = 0; ui < nCount; ui++)
                            {
                                if (!pMat->IsStringOrEmpty(ui))
                                {
                                    double x = pMat->GetDouble(ui);
                                    if (x > 0.0)
                                    {
                                        nVal += log(x);
                                        nValCount++;
                                    }
                                    else if ( x == 0.0 )
                                    {
                                        // value of 0 means that function result will be 0
                                        while ( nParamCount-- > 0 )
                                            PopError();
                                        PushDouble( 0.0 );
                                        return;
                                    }
                                    else
                                        SetError( FormulaError::IllegalArgument);
                                }
                            }
                        }
                    }
                }
                break;
                default : SetError(FormulaError::IllegalParameter); break;
            }
        }
        if (nGlobalError == FormulaError::NONE)
            PushDouble(exp(nVal / nValCount));
        else
            PushError( nGlobalError);
    }

    void ScInterpreter::ScStandard()
    {
        if ( MustHaveParamCount( GetByte(), 3 ) )
        {
            double sigma = GetDouble();
            double mue   = GetDouble();
            double x     = GetDouble();
            if (sigma < 0.0)
                PushError( FormulaError::IllegalArgument);
            else if (sigma == 0.0)
                PushError( FormulaError::DivisionByZero);
            else
                PushDouble((x-mue)/sigma);
        }
    }
    bool ScInterpreter::CalculateSkew(double& fSum,double& fCount,double& vSum,std::vector<double>& values)
    {
        short nParamCount = GetByte();
        if ( !MustHaveParamCountMin( nParamCount, 1 )  )
            return false;

        fSum   = 0.0;
        fCount = 0.0;
        vSum   = 0.0;
        double fVal = 0.0;
        ScAddress aAdr;
        ScRange aRange;
        size_t nRefInList = 0;
        while (nParamCount-- > 0)
        {
            switch (GetStackType())
            {
                case svDouble :
                {
                    fVal = GetDouble();
                    fSum += fVal;
                    values.push_back(fVal);
                    fCount++;
                }
                    break;
                case svSingleRef :
                {
                    PopSingleRef( aAdr );
                    ScRefCellValue aCell(*pDok, aAdr);
                    if (aCell.hasNumeric())
                    {
                        fVal = GetCellValue(aAdr, aCell);
                        fSum += fVal;
                        values.push_back(fVal);
                        fCount++;
                    }
                }
                break;
                case svDoubleRef :
                case svRefList :
                {
                    PopDoubleRef( aRange, nParamCount, nRefInList);
                    FormulaError nErr = FormulaError::NONE;
                    ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
                    if (aValIter.GetFirst(fVal, nErr))
                    {
                        fSum += fVal;
                        values.push_back(fVal);
                        fCount++;
                        SetError(nErr);
                        while ((nErr == FormulaError::NONE) && aValIter.GetNext(fVal, nErr))
                        {
                            fSum += fVal;
                            values.push_back(fVal);
                            fCount++;
                        }
                        SetError(nErr);
                    }
                }
                break;
                case svMatrix :
                case svExternalSingleRef:
                case svExternalDoubleRef:
                {
                    ScMatrixRef pMat = GetMatrix();
                    if (pMat)
                    {
                        SCSIZE nCount = pMat->GetElementCount();
                        if (pMat->IsNumeric())
                        {
                            for (SCSIZE nElem = 0; nElem < nCount; nElem++)
                            {
                                fVal = pMat->GetDouble(nElem);
                                fSum += fVal;
                                values.push_back(fVal);
                                fCount++;
                            }
                        }
                        else
                        {
                            for (SCSIZE nElem = 0; nElem < nCount; nElem++)
                                if (!pMat->IsStringOrEmpty(nElem))
                                {
                                    fVal = pMat->GetDouble(nElem);
                                    fSum += fVal;
                                    values.push_back(fVal);
                                    fCount++;
                                }
                        }
                    }
                }
                break;
                default :
                    SetError(FormulaError::IllegalParameter);
                break;
            }
        }

        if (nGlobalError != FormulaError::NONE)
        {
            PushError( nGlobalError);
            return false;
        } // if (nGlobalError != FormulaError::NONE)
        return true;
    }

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

        double fMean = fSum / fCount;

        for (double v : values)
            vSum += (v - fMean) * (v - fMean);

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

        if (fStdDev == 0)
        {
            PushIllegalArgument();
            return;
        }

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

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

    void ScInterpreter::ScSkew()
    {
        CalculateSkewOrSkewp( false );
    }

    void ScInterpreter::ScSkewp()
    {
        CalculateSkewOrSkewp( true );
    }

    double ScInterpreter::GetMedian( vector<double> & rArray )
    {
        size_t nSize = rArray.size();
        if (nSize == 0 || nGlobalError != FormulaError::NONE)
        {
            SetError( FormulaError::NoValue);
            return 0.0;
        }

        // Upper median.
        size_t nMid = nSize / 2;
        vector<double>::iterator iMid = rArray.begin() + nMid;
        ::std::nth_element( rArray.begin(), iMid, rArray.end());
        if (nSize & 1)
            return *iMid;   // Lower and upper median are equal.
        else
        {
            double fUp = *iMid;
            // Lower median.
            iMid = ::std::max_element( rArray.begin(), rArray.begin() + nMid);
            return (fUp + *iMid) / 2;
        }
    }

    void ScInterpreter::ScMedian()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCountMin( nParamCount, 1 )  )
            return;
        vector<double> aArray;
        GetNumberSequenceArray( nParamCount, aArray, false );
        PushDouble( GetMedian( aArray));
    }

    double ScInterpreter::GetPercentile( vector<double> & rArray, double fPercentile )
    {
        size_t nSize = rArray.size();
        if (nSize == 1)
            return rArray[0];
        else
        {
            size_t nIndex = static_cast<size_t>(::rtl::math::approxFloor( fPercentile * (nSize-1)));
            double fDiff = fPercentile * (nSize-1) - ::rtl::math::approxFloor( fPercentile * (nSize-1));
            OSL_ENSURE(nIndex < nSize, "GetPercentile: wrong index(1)");
            vector<double>::iterator iter = rArray.begin() + nIndex;
            ::std::nth_element( rArray.begin(), iter, rArray.end());
            if (fDiff == 0.0)
                return *iter;
            else
            {
                OSL_ENSURE(nIndex < nSize-1, "GetPercentile: wrong index(2)");
                double fVal = *iter;
                iter = ::std::min_element( rArray.begin() + nIndex + 1, rArray.end());
                return fVal + fDiff * (*iter - fVal);
            }
        }
    }

    double ScInterpreter::GetPercentileExclusive( vector<double> & rArray, double fPercentile )
    {
        size_t nSize1 = rArray.size() + 1;
        if ( rArray.empty() || nSize1 == 1 || nGlobalError != FormulaError::NONE)
        {
            SetError( FormulaError::NoValue );
            return 0.0;
        }
        if ( fPercentile * nSize1 < 1.0 || fPercentile * nSize1 > static_cast<double>( nSize1 - 1 ) )
        {
            SetError( FormulaError::IllegalParameter );
            return 0.0;
        }

        size_t nIndex = static_cast<size_t>(::rtl::math::approxFloor( fPercentile * nSize1 - 1 ));
        double fDiff = fPercentile *  nSize1 - 1 - ::rtl::math::approxFloor( fPercentile * nSize1 - 1 );
        OSL_ENSURE(nIndex < ( nSize1 - 1 ), "GetPercentile: wrong index(1)");
        vector<double>::iterator iter = rArray.begin() + nIndex;
        ::std::nth_element( rArray.begin(), iter, rArray.end());
        if (fDiff == 0.0)
            return *iter;
        else
        {
            OSL_ENSURE(nIndex < nSize1, "GetPercentile: wrong index(2)");
            double fVal = *iter;
            iter = ::std::min_element( rArray.begin() + nIndex + 1, rArray.end());
            return fVal + fDiff * (*iter - fVal);
        }
    }

    void ScInterpreter::ScPercentile( bool bInclusive )
    {
        if ( !MustHaveParamCount( GetByte(), 2 ) )
            return;
        double alpha = GetDouble();
        if ( bInclusive ? ( alpha < 0.0 || alpha > 1.0 ) : ( alpha <= 0.0 || alpha >= 1.0 ) )
        {
            PushIllegalArgument();
            return;
        }
        vector<double> aArray;
        GetNumberSequenceArray( 1, aArray, false );
        if ( aArray.empty() || nGlobalError != FormulaError::NONE )
        {
            SetError( FormulaError::NoValue );
            return;
        }
        if ( bInclusive )
            PushDouble( GetPercentile( aArray, alpha ));
        else
            PushDouble( GetPercentileExclusive( aArray, alpha ));
    }

    void ScInterpreter::ScQuartile( bool bInclusive )
    {
        if ( !MustHaveParamCount( GetByte(), 2 ) )
            return;
        double fFlag = ::rtl::math::approxFloor(GetDouble());
        if ( bInclusive ? ( fFlag < 0.0 || fFlag > 4.0 ) : ( fFlag <= 0.0 || fFlag >= 4.0 ) )
        {
            PushIllegalArgument();
            return;
        }
        vector<double> aArray;
        GetNumberSequenceArray( 1, aArray, false );
        if ( aArray.empty() || nGlobalError != FormulaError::NONE )
        {
            SetError( FormulaError::NoValue );
            return;
        }
        if ( bInclusive )
            PushDouble( fFlag == 2.0 ? GetMedian( aArray ) : GetPercentile( aArray, 0.25 * fFlag ) );
        else
            PushDouble( fFlag == 2.0 ? GetMedian( aArray ) : GetPercentileExclusive( aArray, 0.25 * fFlag ) );
    }

    void ScInterpreter::ScModalValue()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCountMin( nParamCount, 1 ) )
            return;
        vector<double> aSortArray;
        GetSortArray( nParamCount, aSortArray, nullptr, false, false );
        SCSIZE nSize = aSortArray.size();
        if (nSize == 0 || nGlobalError != FormulaError::NONE)
            PushNoValue();
        else
        {
            SCSIZE nMaxIndex = 0, nMax = 1, nCount = 1;
            double nOldVal = aSortArray[0];
            SCSIZE i;
            for ( i = 1; i < nSize; i++)
            {
                if (aSortArray[i] == nOldVal)
                    nCount++;
                else
                {
                    nOldVal = aSortArray[i];
                    if (nCount > nMax)
                    {
                        nMax = nCount;
                        nMaxIndex = i-1;
                    }
                    nCount = 1;
                }
            }
            if (nCount > nMax)
            {
                nMax = nCount;
                nMaxIndex = i-1;
            }
            if (nMax == 1 && nCount == 1)
                PushNoValue();
            else if (nMax == 1)
                PushDouble(nOldVal);
            else
                PushDouble(aSortArray[nMaxIndex]);
        }
    }

    void ScInterpreter::ScModalValue_MS( bool bSingle )
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCountMin( nParamCount, 1 ) )
            return;
        vector<double> aArray;
        GetNumberSequenceArray( nParamCount, aArray, false );
        vector< double > aSortArray( aArray );
        QuickSort( aSortArray, nullptr );
        SCSIZE nSize = aSortArray.size();
        if ( nSize == 0 || nGlobalError != FormulaError::NONE )
            PushNoValue();
        else
        {
            SCSIZE nMax = 1, nCount = 1;
            double nOldVal = aSortArray[ 0 ];
            vector< double > aResultArray( 1 );
            SCSIZE i;
            for ( i = 1; i < nSize; i++ )
            {
                if ( aSortArray[ i ] == nOldVal )
                    nCount++;
                else
                {
                    if ( nCount >= nMax && nCount > 1 )
                    {
                        if ( nCount > nMax )
                        {
                            nMax = nCount;
                            if ( aResultArray.size() != 1 )
                                vector< double >( 1 ).swap( aResultArray );
                            aResultArray[ 0 ] = nOldVal;
                        }
                        else
                            aResultArray.emplace_back( nOldVal );
                    }
                    nOldVal = aSortArray[ i ];
                    nCount = 1;
                }
            }
            if ( nCount >= nMax && nCount > 1 )
            {
                if ( nCount > nMax )
                    vector< double >().swap( aResultArray );
                aResultArray.emplace_back( nOldVal );
            }
            if ( nMax == 1 && nCount == 1 )
                PushNoValue();
            else if ( nMax == 1 )
                PushDouble( nOldVal ); // there is only 1 result, no reordering needed
            else
            {
                // sort resultArray according to ordering of aArray
                vector< vector< double > > aOrder;
                aOrder.resize( aResultArray.size(), vector< double >( 2 ) );
                for ( i = 0; i < aResultArray.size(); i++ )
                {
                    for ( SCSIZE j = 0; j < nSize; j++ )
                    {
                        if ( aArray[ j ] == aResultArray[ i ] )
                        {
                            aOrder[ i ][ 0 ] = aResultArray[ i ];
                            aOrder[ i ][ 1 ] = j;
                            break;
                        }
                    }
                }
                sort( aOrder.begin(), aOrder.end(), []( const std::vector< double >& lhs,
                                                        const std::vector< double >& rhs )
                                                        { return lhs[ 1 ] < rhs[ 1 ]; } );

                if ( bSingle )
                    PushDouble( aOrder[ 0 ][ 0 ] );
                else
                {
                    // put result in correct order in aResultArray
                    for ( i = 0; i < aResultArray.size(); i++ )
                        aResultArray[ i ] = aOrder[ i ][ 0 ];
                    ScMatrixRef pResMatrix = GetNewMat( 1, aResultArray.size(), true );
                    pResMatrix->PutDoubleVector( aResultArray, 0, 0 );
                    PushMatrix( pResMatrix );
                }
            }
        }
    }

    void ScInterpreter::CalculateSmallLarge(bool bSmall)
    {
        if ( !MustHaveParamCount( GetByte(), 2 )  )
            return;
        double f = ::rtl::math::approxFloor(GetDouble());
        if (f < 1.0)
        {
            PushIllegalArgument();
            return;
        }
        SCSIZE k = static_cast<SCSIZE>(f);
        vector<double> aSortArray;
        /* TODO: using nth_element() is best for one single value, but LARGE/SMALL
         * actually are defined to return an array of values if an array of
         * positions was passed, in which case, depending on the number of values,
         * we may or will need a real sorted array again, see #i32345. */
        GetNumberSequenceArray(1, aSortArray, false );
        SCSIZE nSize = aSortArray.size();
        if (nSize == 0 || nGlobalError != FormulaError::NONE || nSize < k)
            PushNoValue();
        else
        {
            // TODO: the sorted case for array: PushDouble( aSortArray[ bSmall ? k-1 : nSize-k ] );
            vector<double>::iterator iPos = aSortArray.begin() + (bSmall ? k-1 : nSize-k);
            ::std::nth_element( aSortArray.begin(), iPos, aSortArray.end());
            PushDouble( *iPos);
        }
    }

    void ScInterpreter::ScLarge()
    {
        CalculateSmallLarge(false);
    }

    void ScInterpreter::ScSmall()
    {
        CalculateSmallLarge(true);
    }

    void ScInterpreter::ScPercentrank( bool bInclusive )
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
            return;
        double fSignificance = ( nParamCount == 3 ? ::rtl::math::approxFloor( GetDouble() ) : 3.0 );
        if ( fSignificance < 1.0 )
        {
            PushIllegalArgument();
            return;
        }
        double fNum = GetDouble();
        vector<double> aSortArray;
        GetSortArray( 1, aSortArray, nullptr, false, false );
        SCSIZE nSize = aSortArray.size();
        if ( nSize == 0 || nGlobalError != FormulaError::NONE )
            PushNoValue();
        else
        {
            if ( fNum < aSortArray[ 0 ] || fNum > aSortArray[ nSize - 1 ] )
                PushNoValue();
            else
            {
                double fRes;
                if ( nSize == 1 )
                    fRes = 1.0;            // fNum == aSortArray[ 0 ], see test above
                else
                    fRes = GetPercentrank( aSortArray, fNum, bInclusive );
                if ( fRes != 0.0 )
                {
                    double fExp = ::rtl::math::approxFloor( log10( fRes ) ) + 1.0 - fSignificance;
                    fRes = ::rtl::math::round( fRes * pow( 10, -fExp ) ) / pow( 10, -fExp );
                }
                PushDouble( fRes );
            }
        }
    }

    double ScInterpreter::GetPercentrank( ::std::vector<double> & rArray, double fVal, bool bInclusive )
    {
        SCSIZE nSize = rArray.size();
        double fRes;
        if ( fVal == rArray[ 0 ] )
        {
            if ( bInclusive )
                fRes = 0.0;
            else
                fRes = 1.0 / static_cast<double>( nSize + 1 );
        }
        else
        {
            SCSIZE nOldCount = 0;
            double fOldVal = rArray[ 0 ];
            SCSIZE i;
            for ( i = 1; i < nSize && rArray[ i ] < fVal; i++ )
            {
                if ( rArray[ i ] != fOldVal )
                {
                    nOldCount = i;
                    fOldVal = rArray[ i ];
                }
            }
            if ( rArray[ i ] != fOldVal )
                nOldCount = i;
            if ( fVal == rArray[ i ] )
            {
                if ( bInclusive )
                    fRes = div( nOldCount, nSize - 1 );
                else
                    fRes = static_cast<double>( i + 1 ) / static_cast<double>( nSize + 1 );
            }
            else
            {
                //  nOldCount is the count of smaller entries
                //  fVal is between rArray[ nOldCount - 1 ] and rArray[ nOldCount ]
                //  use linear interpolation to find a position between the entries
                if ( nOldCount == 0 )
                {
                    OSL_FAIL( "should not happen" );
                    fRes = 0.0;
                }
                else
                {
                    double fFract = ( fVal - rArray[ nOldCount - 1 ] ) /
                        ( rArray[ nOldCount ] - rArray[ nOldCount - 1 ] );
                    if ( bInclusive )
                        fRes = div( static_cast<double>( nOldCount - 1 ) + fFract, nSize - 1 );
                    else
                        fRes = ( static_cast<double>(nOldCount) + fFract ) / static_cast<double>( nSize + 1 );
                }
            }
        }
        return fRes;
    }

    void ScInterpreter::ScTrimMean()
    {
        if ( !MustHaveParamCount( GetByte(), 2 ) )
            return;
        double alpha = GetDouble();
        if (alpha < 0.0 || alpha >= 1.0)
        {
            PushIllegalArgument();
            return;
        }
        vector<double> aSortArray;
        GetSortArray( 1, aSortArray, nullptr, false, false );
        SCSIZE nSize = aSortArray.size();
        if (nSize == 0 || nGlobalError != FormulaError::NONE)
            PushNoValue();
        else
        {
            sal_uLong nIndex = static_cast<sal_uLong>(::rtl::math::approxFloor(alpha*static_cast<double>(nSize)));
            if (nIndex % 2 != 0)
                nIndex--;
            nIndex /= 2;
            OSL_ENSURE(nIndex < nSize, "ScTrimMean: wrong index");
            double fSum = 0.0;
            for (SCSIZE i = nIndex; i < nSize-nIndex; i++)
                fSum += aSortArray[i];
            PushDouble(fSum/static_cast<double>(nSize-2*nIndex));
        }
    }

    void ScInterpreter::GetNumberSequenceArray( sal_uInt8 nParamCount, vector<double>& rArray, bool bConvertTextInArray )
    {
        ScAddress aAdr;
        ScRange aRange;
        short nParam = nParamCount;
        size_t nRefInList = 0;
        ReverseStack( nParamCount );
        while (nParam-- > 0)
        {
            const StackVar eStackType = GetStackType();
            switch (eStackType)
            {
                case svDouble :
                    rArray.push_back( PopDouble());
                break;
                case svSingleRef :
                {
                    PopSingleRef( aAdr );
                    ScRefCellValue aCell(*pDok, aAdr);
                    if (aCell.hasNumeric())
                        rArray.push_back(GetCellValue(aAdr, aCell));
                }
                break;
                case svDoubleRef :
                case svRefList :
                {
                    PopDoubleRef( aRange, nParam, nRefInList);
                    if (nGlobalError != FormulaError::NONE)
                        break;

                    aRange.PutInOrder();
                    SCSIZE nCellCount = aRange.aEnd.Col() - aRange.aStart.Col() + 1;
                    nCellCount *= aRange.aEnd.Row() - aRange.aStart.Row() + 1;
                    rArray.reserve( rArray.size() + nCellCount);

                    FormulaError nErr = FormulaError::NONE;
                    double fCellVal;
                    ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
                    if (aValIter.GetFirst( fCellVal, nErr))
                    {
                        rArray.push_back( fCellVal);
                        SetError(nErr);
                        while ((nErr == FormulaError::NONE) && aValIter.GetNext( fCellVal, nErr))
                            rArray.push_back( fCellVal);
                        SetError(nErr);
                    }
                }
                break;
                case svMatrix :
                case svExternalSingleRef:
                case svExternalDoubleRef:
                {
                    ScMatrixRef pMat = GetMatrix();
                    if (!pMat)
                        break;

                    SCSIZE nCount = pMat->GetElementCount();
                    rArray.reserve( rArray.size() + nCount);
                    if (pMat->IsNumeric())
                    {
                        for (SCSIZE i = 0; i < nCount; ++i)
                            rArray.push_back( pMat->GetDouble(i));
                    }
                    else if (bConvertTextInArray && eStackType == svMatrix)
                    {
                        for (SCSIZE i = 0; i < nCount; ++i)
                        {
                            if ( pMat->IsValue( i ) )
                                rArray.push_back( pMat->GetDouble(i));
                            else
                            {
                                // tdf#88547 try to convert string to (date)value
                                OUString aStr = pMat->GetString( i ).getString();
                                if ( aStr.getLength() > 0 )
                                {
                                    FormulaError nErr = nGlobalError;
                                    nGlobalError = FormulaError::NONE;
                                    double fVal = ConvertStringToValue( aStr );
                                    if ( nGlobalError == FormulaError::NONE )
                                    {
                                        rArray.push_back( fVal );
                                        nGlobalError = nErr;
                                    }
                                    else
                                    {
                                        rArray.push_back( CreateDoubleError( FormulaError::NoValue));
                                        // Propagate previous error if any, else
                                        // the current #VALUE! error.
                                        if (nErr != FormulaError::NONE)
                                            nGlobalError = nErr;
                                        else
                                            nGlobalError = FormulaError::NoValue;
                                    }
                                }
                            }
                        }
                    }
                    else
                    {
                        for (SCSIZE i = 0; i < nCount; ++i)
                        {
                            if ( pMat->IsValue( i ) )
                                rArray.push_back( pMat->GetDouble(i));
                        }
                    }
                }
                break;
                default :
                    PopError();
                    SetError( FormulaError::IllegalParameter);
                break;
            }
            if (nGlobalError != FormulaError::NONE)
                break;  // while
        }
        // nParam > 0 in case of error, clean stack environment and obtain earlier
        // error if there was one.
        while (nParam-- > 0)
            PopError();
    }

    void ScInterpreter::GetSortArray( sal_uInt8 nParamCount, vector<double>& rSortArray, vector<long>* pIndexOrder, bool bConvertTextInArray, bool bAllowEmptyArray )
    {
        GetNumberSequenceArray( nParamCount, rSortArray, bConvertTextInArray );
        if (rSortArray.size() > MAX_COUNT_DOUBLE_FOR_SORT)
            SetError( FormulaError::MatrixSize);
        else if ( rSortArray.empty() )
        {
            if ( bAllowEmptyArray )
                return;
            SetError( FormulaError::NoValue);
        }

        if (nGlobalError == FormulaError::NONE)
            QuickSort( rSortArray, pIndexOrder);
    }

    static void lcl_QuickSort( long nLo, long nHi, vector<double>& rSortArray, vector<long>* pIndexOrder )
    {
        // If pIndexOrder is not NULL, we assume rSortArray.size() == pIndexOrder->size().

        using ::std::swap;

        if (nHi - nLo == 1)
        {
            if (rSortArray[nLo] > rSortArray[nHi])
            {
                swap(rSortArray[nLo],  rSortArray[nHi]);
                if (pIndexOrder)
                    swap(pIndexOrder->at(nLo), pIndexOrder->at(nHi));
            }
            return;
        }

        long ni = nLo;
        long nj = nHi;
        do
        {
            double fLo = rSortArray[nLo];
            while (ni <= nHi && rSortArray[ni] < fLo) ni++;
            while (nj >= nLo && fLo < rSortArray[nj]) nj--;
            if (ni <= nj)
            {
                if (ni != nj)
                {
                    swap(rSortArray[ni],  rSortArray[nj]);
                    if (pIndexOrder)
                        swap(pIndexOrder->at(ni), pIndexOrder->at(nj));
                }

                ++ni;
                --nj;
            }
        }
        while (ni < nj);

        if ((nj - nLo) < (nHi - ni))
        {
            if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
            if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
        }
        else
        {
            if (ni < nHi) lcl_QuickSort(ni, nHi, rSortArray, pIndexOrder);
            if (nLo < nj) lcl_QuickSort(nLo, nj, rSortArray, pIndexOrder);
        }
    }

    void ScInterpreter::QuickSort( vector<double>& rSortArray, vector<long>* pIndexOrder )
    {
        long n = static_cast<long>(rSortArray.size());

        if (pIndexOrder)
        {
            pIndexOrder->clear();
            pIndexOrder->reserve(n);
            for (long i = 0; i < n; ++i)
                pIndexOrder->push_back(i);
        }

        if (n < 2)
            return;

        size_t nValCount = rSortArray.size();
        for (size_t i = 0; (i + 4) <= nValCount-1; i += 4)
        {
            size_t nInd = comphelper::rng::uniform_size_distribution(0, nValCount-2);
            ::std::swap( rSortArray[i], rSortArray[nInd]);
            if (pIndexOrder)
                ::std::swap( pIndexOrder->at(i), pIndexOrder->at(nInd));
        }

        lcl_QuickSort(0, n-1, rSortArray, pIndexOrder);
    }

    void ScInterpreter::ScRank( bool bAverage )
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, 2, 3 ) )
            return;
        bool bAscending;
        if ( nParamCount == 3 )
            bAscending = GetBool();
        else
            bAscending = false;

        vector<double> aSortArray;
        GetSortArray( 1, aSortArray, nullptr, false, false );
        double fVal = GetDouble();
        SCSIZE nSize = aSortArray.size();
        if ( nSize == 0 || nGlobalError != FormulaError::NONE )
            PushNoValue();
        else
        {
            if ( fVal < aSortArray[ 0 ] || fVal > aSortArray[ nSize - 1 ] )
                PushNoValue();
            else
            {
                double fLastPos = 0;
                double fFirstPos = -1.0;
                bool bFinished = false;
                SCSIZE i;
                for ( i = 0; i < nSize && !bFinished && nGlobalError == FormulaError::NONE; i++ )
                {
                    if ( aSortArray[ i ] == fVal )
                    {
                        if ( fFirstPos < 0 )
                            fFirstPos = i + 1.0;
                    }
                    else
                    {
                        if ( aSortArray[ i ] > fVal )
                        {
                            fLastPos = i;
                            bFinished = true;
                        }
                    }
                }
                if ( !bFinished )
                    fLastPos = i;
                if ( !bAverage )
                {
                    if ( bAscending )
                        PushDouble( fFirstPos );
                    else
                        PushDouble( nSize + 1.0 - fLastPos );
                }
                else
                {
                    if ( bAscending )
                        PushDouble( ( fFirstPos + fLastPos ) / 2.0 );
                    else
                        PushDouble( nSize + 1.0 - ( fFirstPos + fLastPos ) / 2.0 );
                }
            }
        }
    }

    void ScInterpreter::ScAveDev()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCountMin( nParamCount, 1 ) )
            return;
        sal_uInt16 SaveSP = sp;
        double nMiddle = 0.0;
        double rVal = 0.0;
        double rValCount = 0.0;
        ScAddress aAdr;
        ScRange aRange;
        short nParam = nParamCount;
        size_t nRefInList = 0;
        while (nParam-- > 0)
        {
            switch (GetStackType())
            {
                case svDouble :
                    rVal += GetDouble();
                    rValCount++;
                    break;
                case svSingleRef :
                {
                    PopSingleRef( aAdr );
                    ScRefCellValue aCell(*pDok, aAdr);
                    if (aCell.hasNumeric())
                    {
                        rVal += GetCellValue(aAdr, aCell);
                        rValCount++;
                    }
                }
                break;
                case svDoubleRef :
                case svRefList :
                {
                    FormulaError nErr = FormulaError::NONE;
                    double nCellVal;
                    PopDoubleRef( aRange, nParam, nRefInList);
                    ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
                    if (aValIter.GetFirst(nCellVal, nErr))
                    {
                        rVal += nCellVal;
                        rValCount++;
                        SetError(nErr);
                        while ((nErr == FormulaError::NONE) && aValIter.GetNext(nCellVal, nErr))
                        {
                            rVal += nCellVal;
                            rValCount++;
                        }
                        SetError(nErr);
                    }
                }
                break;
                case svMatrix :
                case svExternalSingleRef:
                case svExternalDoubleRef:
                {
                    ScMatrixRef pMat = GetMatrix();
                    if (pMat)
                    {
                        SCSIZE nCount = pMat->GetElementCount();
                        if (pMat->IsNumeric())
                        {
                            for (SCSIZE nElem = 0; nElem < nCount; nElem++)
                            {
                                rVal += pMat->GetDouble(nElem);
                                rValCount++;
                            }
                        }
                        else
                        {
                            for (SCSIZE nElem = 0; nElem < nCount; nElem++)
                                if (!pMat->IsStringOrEmpty(nElem))
                                {
                                    rVal += pMat->GetDouble(nElem);
                                    rValCount++;
                                }
                        }
                    }
                }
                break;
                default :
                    SetError(FormulaError::IllegalParameter);
                break;
            }
        }
        if (nGlobalError != FormulaError::NONE)
        {
            PushError( nGlobalError);
            return;
        }
        nMiddle = rVal / rValCount;
        sp = SaveSP;
        rVal = 0.0;
        nParam = nParamCount;
        nRefInList = 0;
        while (nParam-- > 0)
        {
            switch (GetStackType())
            {
                case svDouble :
                    rVal += fabs(GetDouble() - nMiddle);
                    break;
                case svSingleRef :
                {
                    PopSingleRef( aAdr );
                    ScRefCellValue aCell(*pDok, aAdr);
                    if (aCell.hasNumeric())
                        rVal += fabs(GetCellValue(aAdr, aCell) - nMiddle);
                }
                break;
                case svDoubleRef :
                case svRefList :
                {
                    FormulaError nErr = FormulaError::NONE;
                    double nCellVal;
                    PopDoubleRef( aRange, nParam, nRefInList);
                    ScValueIterator aValIter( pDok, aRange, mnSubTotalFlags );
                    if (aValIter.GetFirst(nCellVal, nErr))
                    {
                        rVal += (fabs(nCellVal - nMiddle));
                        while (aValIter.GetNext(nCellVal, nErr))
                             rVal += fabs(nCellVal - nMiddle);
                    }
                }
                break;
                case svMatrix :
                case svExternalSingleRef:
                case svExternalDoubleRef:
                {
                    ScMatrixRef pMat = GetMatrix();
                    if (pMat)
                    {
                        SCSIZE nCount = pMat->GetElementCount();
                        if (pMat->IsNumeric())
                        {
                            for (SCSIZE nElem = 0; nElem < nCount; nElem++)
                            {
                                rVal += fabs(pMat->GetDouble(nElem) - nMiddle);
                            }
                        }
                        else
                        {
                            for (SCSIZE nElem = 0; nElem < nCount; nElem++)
                            {
                                if (!pMat->IsStringOrEmpty(nElem))
                                    rVal += fabs(pMat->GetDouble(nElem) - nMiddle);
                            }
                        }
                    }
                }
                break;
                default : SetError(FormulaError::IllegalParameter); break;
            }
        }
        PushDouble(rVal / rValCount);
    }

    void ScInterpreter::ScDevSq()
    {
        auto VarResult = []( double fVal, size_t /*nValCount*/ )
        {
            return fVal;
        };
        GetStVarParams( false /*bTextAsZero*/, VarResult);
    }

    void ScInterpreter::ScProbability()
    {
        sal_uInt8 nParamCount = GetByte();
        if ( !MustHaveParamCount( nParamCount, 3, 4 ) )
            return;
        double fUp, fLo;
        fUp = GetDouble();
        if (nParamCount == 4)
            fLo = GetDouble();
        else
            fLo = fUp;
        if (fLo > fUp)
        {
            double fTemp = fLo;
            fLo = fUp;
            fUp = fTemp;
        }
        ScMatrixRef pMatP = GetMatrix();
        ScMatrixRef pMatW = GetMatrix();
        if (!pMatP || !pMatW)
            PushIllegalParameter();
        else
        {
            SCSIZE nC1, nC2;
            SCSIZE nR1, nR2;
            pMatP->GetDimensions(nC1, nR1);
            pMatW->GetDimensions(nC2, nR2);
            if (nC1 != nC2 || nR1 != nR2 || nC1 == 0 || nR1 == 0 ||
                nC2 == 0 || nR2 == 0)
                PushNA();
            else
            {
                double fSum = 0.0;
                double fRes = 0.0;
                bool bStop = false;
                double fP, fW;
                for ( SCSIZE i = 0; i < nC1 && !bStop; i++ )
                {
                    for (SCSIZE j = 0; j < nR1 && !bStop; ++j )
                    {
                        if (pMatP->IsValue(i,j) && pMatW->IsValue(i,j))
                        {
                            fP = pMatP->GetDouble(i,j);
                            fW = pMatW->GetDouble(i,j);
                            if (fP < 0.0 || fP > 1.0)
                                bStop = true;
                            else
                            {
                                fSum += fP;
                                if (fW >= fLo && fW <= fUp)
                                    fRes += fP;
                            }
                        }
                        else
                            SetError( FormulaError::IllegalArgument);
                    }
                }
                if (bStop || fabs(fSum -1.0) > 1.0E-7)
                    PushNoValue();
                else
                    PushDouble(fRes);
            }
        }
    }

    void ScInterpreter::ScCorrel()
    {
        // This is identical to ScPearson()
        ScPearson();
    }

    void ScInterpreter::ScCovarianceP()
    {
        CalculatePearsonCovar( false, false, false );
    }

    void ScInterpreter::ScCovarianceS()
    {
        CalculatePearsonCovar( false, false, true );
    }

    void ScInterpreter::ScPearson()
    {
        CalculatePearsonCovar( true, false, false );
    }

    void ScInterpreter::CalculatePearsonCovar( bool _bPearson, bool _bStexy, bool _bSample )
    {
        if ( !MustHaveParamCount( GetByte(), 2 ) )
            return;
        ScMatrixRef pMat1 = GetMatrix();
        ScMatrixRef pMat2 = GetMatrix();
        if (!pMat1 || !pMat2)
        {
            PushIllegalParameter();
            return;
        }
        SCSIZE nC1, nC2;
        SCSIZE nR1, nR2;
        pMat1->GetDimensions(nC1, nR1);
        pMat2->GetDimensions(nC2, nR2);
        if (nR1 != nR2 || nC1 != nC2)
        {
            PushIllegalArgument();
            return;
        }
        /* #i78250#
         * (sum((X-MeanX)(Y-MeanY)))/N equals (SumXY)/N-MeanX*MeanY mathematically,
         * but the latter produces wrong results if the absolute values are high,
         * for example above 10^8
         */
        double fCount           = 0.0;
        double fSumX            = 0.0;
        double fSumY            = 0.0;

        for (SCSIZE i = 0; i < nC1; i++)
        {
            for (SCSIZE j = 0; j < nR1; j++)
            {
                if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
                {
                    double fValX = pMat1->GetDouble(i,j);
                    double fValY = pMat2->GetDouble(i,j);
                    fSumX += fValX;
                    fSumY += fValY;
                    fCount++;
                }
            }
        }
        if (fCount < (_bStexy ? 3.0 : (_bSample ? 2.0 : 1.0)))
            PushNoValue();
        else
        {
            double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
            double fSumSqrDeltaX    = 0.0; // sum of (ValX-MeanX)^2
            double fSumSqrDeltaY    = 0.0; // sum of (ValY-MeanY)^2
            const double fMeanX = fSumX / fCount;
            const double fMeanY = fSumY / fCount;
            for (SCSIZE i = 0; i < nC1; i++)
            {
                for (SCSIZE j = 0; j < nR1; j++)
                {
                    if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
                    {
                        const double fValX = pMat1->GetDouble(i,j);
                        const double fValY = pMat2->GetDouble(i,j);
                        fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
                        if ( _bPearson )
                        {
                            fSumSqrDeltaX    += (fValX - fMeanX) * (fValX - fMeanX);
                            fSumSqrDeltaY    += (fValY - fMeanY) * (fValY - fMeanY);
                        }
                    }
                }
            }
            if ( _bPearson )
            {
                if (fSumSqrDeltaX == 0.0 || ( !_bStexy && fSumSqrDeltaY == 0.0) )
                    PushError( FormulaError::DivisionByZero);
                else if ( _bStexy )
                    PushDouble( sqrt( (fSumSqrDeltaY - fSumDeltaXDeltaY *
                                fSumDeltaXDeltaY / fSumSqrDeltaX) / (fCount-2)));
                else
                    PushDouble( fSumDeltaXDeltaY / sqrt( fSumSqrDeltaX * fSumSqrDeltaY));
            }
            else
            {
                if ( _bSample )
                    PushDouble( fSumDeltaXDeltaY / ( fCount - 1 ) );
                else
                    PushDouble( fSumDeltaXDeltaY / fCount);
            }
        }
    }

    void ScInterpreter::ScRSQ()
    {
        // Same as ScPearson()*ScPearson()
        ScPearson();
        if (nGlobalError == FormulaError::NONE)
        {
            switch (GetStackType())
            {
                case svDouble:
                    {
                        double fVal = PopDouble();
                        PushDouble( fVal * fVal);
                    }
                    break;
                default:
                    PopError();
                    PushNoValue();
            }
        }
    }

    void ScInterpreter::ScSTEYX()
    {
        CalculatePearsonCovar( true, true, false );
    }
    void ScInterpreter::CalculateSlopeIntercept(bool bSlope)
    {
        if ( !MustHaveParamCount( GetByte(), 2 ) )
            return;
        ScMatrixRef pMat1 = GetMatrix();
        ScMatrixRef pMat2 = GetMatrix();
        if (!pMat1 || !pMat2)
        {
            PushIllegalParameter();
            return;
        }
        SCSIZE nC1, nC2;
        SCSIZE nR1, nR2;
        pMat1->GetDimensions(nC1, nR1);
        pMat2->GetDimensions(nC2, nR2);
        if (nR1 != nR2 || nC1 != nC2)
        {
            PushIllegalArgument();
            return;
        }
        // #i78250# numerical stability improved
        double fCount           = 0.0;
        double fSumX            = 0.0;
        double fSumY            = 0.0;

        for (SCSIZE i = 0; i < nC1; i++)
        {
            for (SCSIZE j = 0; j < nR1; j++)
            {
                if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
                {
                    double fValX = pMat1->GetDouble(i,j);
                    double fValY = pMat2->GetDouble(i,j);
                    fSumX += fValX;
                    fSumY += fValY;
                    fCount++;
                }
            }
        }
        if (fCount < 1.0)
            PushNoValue();
        else
        {
            double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
            double fSumSqrDeltaX    = 0.0; // sum of (ValX-MeanX)^2
            double fMeanX = fSumX / fCount;
            double fMeanY = fSumY / fCount;
            for (SCSIZE i = 0; i < nC1; i++)
            {
                for (SCSIZE j = 0; j < nR1; j++)
                {
                    if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
                    {
                        double fValX = pMat1->GetDouble(i,j);
                        double fValY = pMat2->GetDouble(i,j);
                        fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
                        fSumSqrDeltaX    += (fValX - fMeanX) * (fValX - fMeanX);
                    }
                }
            }
            if (fSumSqrDeltaX == 0.0)
                PushError( FormulaError::DivisionByZero);
            else
            {
                if ( bSlope )
                    PushDouble( fSumDeltaXDeltaY / fSumSqrDeltaX);
                else
                    PushDouble( fMeanY - fSumDeltaXDeltaY / fSumSqrDeltaX * fMeanX);
            }
        }
    }

    void ScInterpreter::ScSlope()
    {
        CalculateSlopeIntercept(true);
    }

    void ScInterpreter::ScIntercept()
    {
        CalculateSlopeIntercept(false);
    }

    void ScInterpreter::ScForecast()
    {
        if ( !MustHaveParamCount( GetByte(), 3 ) )
            return;
        ScMatrixRef pMat1 = GetMatrix();
        ScMatrixRef pMat2 = GetMatrix();
        if (!pMat1 || !pMat2)
        {
            PushIllegalParameter();
            return;
        }
        SCSIZE nC1, nC2;
        SCSIZE nR1, nR2;
        pMat1->GetDimensions(nC1, nR1);
        pMat2->GetDimensions(nC2, nR2);
        if (nR1 != nR2 || nC1 != nC2)
        {
            PushIllegalArgument();
            return;
        }
        double fVal = GetDouble();
        // #i78250# numerical stability improved
        double fCount           = 0.0;
        double fSumX            = 0.0;
        double fSumY            = 0.0;

        for (SCSIZE i = 0; i < nC1; i++)
        {
            for (SCSIZE j = 0; j < nR1; j++)
            {
                if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
                {
                    double fValX = pMat1->GetDouble(i,j);
                    double fValY = pMat2->GetDouble(i,j);
                    fSumX += fValX;
                    fSumY += fValY;
                    fCount++;
                }
            }
        }
        if (fCount < 1.0)
            PushNoValue();
        else
        {
            double fSumDeltaXDeltaY = 0.0; // sum of (ValX-MeanX)*(ValY-MeanY)
            double fSumSqrDeltaX    = 0.0; // sum of (ValX-MeanX)^2
            double fMeanX = fSumX / fCount;
            double fMeanY = fSumY / fCount;
            for (SCSIZE i = 0; i < nC1; i++)
            {
                for (SCSIZE j = 0; j < nR1; j++)
                {
                    if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
                    {
                        double fValX = pMat1->GetDouble(i,j);
                        double fValY = pMat2->GetDouble(i,j);
                        fSumDeltaXDeltaY += (fValX - fMeanX) * (fValY - fMeanY);
                        fSumSqrDeltaX    += (fValX - fMeanX) * (fValX - fMeanX);
                    }
                }
            }
            if (fSumSqrDeltaX == 0.0)
                PushError( FormulaError::DivisionByZero);
            else
                PushDouble( fMeanY + fSumDeltaXDeltaY / fSumSqrDeltaX * (fVal - fMeanX));
        }
    }

    /* vim:set shiftwidth=4 softtabstop=4 expandtab: */

↑ V1023 A pointer without owner is added to the 'aOp' container by the 'emplace_back' method. A memory leak will occur in case of an exception.

↑ V560 A part of conditional expression is always true: (xShort >= 1).

↑ V1023 A pointer without owner is added to the 'aOp' container by the 'emplace_back' method. A memory leak will occur in case of an exception.

↑ V560 A part of conditional expression is always true: nGlobalError == FormulaError::NONE.

↑ V560 A part of conditional expression is always true: (xShort >= 3).