Quelle  interpr3.cxx   Sprache: C

/* -*- 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 <interpre.hxx>
#include <global.hxx>
#include <document.hxx>
#include <dociter.hxx>
#include <matrixoperators.hxx>
#include <scmatrix.hxx>
#include <columniterator.hxx>
#include <unotools/collatorwrapper.hxx>

#include <cassert>
#include <cmath>
#include <memory>
#include <set>
#include <vector>
#include <algorithm>
#include <comphelper/random.hxx>
#include <o3tl/float_int_conversion.hxx>
#include <osl/diagnose.h>

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

/// Two columns of data should be sortable with GetSortArray() and QuickSort()
// This is an arbitrary limit.
static size_t MAX_COUNT_DOUBLE_FOR_SORT(const ScSheetLimits& rSheetLimits)
{
    return rSheetLimits.GetMaxRowCount() * 2;
}

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

namespace {

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 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

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

    fAx = fkAx.get();
    fBx = fkBx.get();
    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 && std::abs(fRy) > fYEps &&
            (fBx-fAx) > ::std::max( std::abs(fAx), std::abs(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 too small advance, then do bisection next
        // time, for safety
        bHasToInterpolate = bHasToInterpolate && (std::abs(fRy) * 2.0 <= std::abs(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 * std::erfc(-x * M_SQRT1_2);
}

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

double ScInterpreter::gauss(double x)
{

    double xAbs = std::abs(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 <= 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 <= 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(std::abs(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::BinomCoeff(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(M_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) - std::log1p(fZ) - log( std::abs(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( std::abs( ::rtl::math::sin( M_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( M_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) - std::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*M_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 * std::log1p(fTempA)
                            -fB * std::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 * std::log1p(fTempA)
                        -fB * std::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 * std::expm1((fB-1.0) * std::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) ? std::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 = (std::abs(cf-cfnew) < std::abs(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 -std::expm1(fBeta * std::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 = std::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 (std::abs(fVal) >= 1.0)
        PushIllegalArgument();
    else
        PushDouble(::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(BinomCoeff(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(BinomCoeff(n + k - 1, k));
    }
}

void ScInterpreter::ScPermut()
{
    if ( !MustHaveParamCount( GetByte(), 2 ) )
        return;

    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;
    }
}

static 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;
    // 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;
    KahanSum fSum = fFactor; // Summand xs
    sal_uInt32 nXe = static_cast<sal_uInt32>(xe);
    for (i = nXs+1; i <= nXe && fFactor > 0.0; i++)
    {
        fFactor *= (n-i+1)/i * p/q;
        fSum += fFactor;
    }
    return std::min(fSum.get(), 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 ) )
        return;

    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 ) )
        return;

    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())
            {
                KahanSum 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
                KahanSum fSum = 0.0;
                sal_uInt32 max = static_cast<sal_uInt32> (n), i;
                for (i = 0; i < max && fSum < alpha; i++)
                {
                    const double x = GetBetaDistPDF( p, ( i + 1 ), ( n - i + 1 ) )/( n + 1 );
                    if ( nGlobalError == FormulaError::NONE )
                        fSum += x;
                    else
                    {
                        PushNoValue();
                        return;
                    }
                }
                assert(i > 0 && "coverity 2023.12.2");
                PushDouble( i - 1 );
            }
        }
        else
        {
            // work from n backwards
            fFactor = pow(p, n);
            if (fFactor > ::std::numeric_limits<double>::min())
            {
                KahanSum 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
                KahanSum fSum = 0.0;
                sal_uInt32 max = static_cast<sal_uInt32> (n), i;
                alpha = 1 - alpha;
                for (i = 0; i < max && fSum < alpha; i++)
                {
                    const double x = GetBetaDistPDF( q, ( i + 1 ), ( n - i + 1 ) )/( n + 1 );
                    if ( nGlobalError == FormulaError::NONE )
                        fSum += x;
                    else
                    {
                        PushNoValue();
                        return;
                    }
                }
                PushDouble( n - i + 1 );
            }
        }
    }
}

void ScInterpreter::ScNegBinomDist()
{
    if ( !MustHaveParamCount( GetByte(), 3 ) )
        return;

    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 ) )
        return;

    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 * M_PI ) );
}

void ScInterpreter::ScExpDist()
{
    if ( !MustHaveParamCount( GetByte(), 3 ) )
        return;

    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 ) )
        return;

    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 ) )
        return;

    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 = std::exp(-lambda);
                KahanSum fSum = fSummand;
                int nEnd = sal::static_int_cast<int>( x );
                for (int i = 1; i <= nEnd; i++)
                {
                    fSummand = (fSummand * lambda)/static_cast<double>(i);
                    fSum += fSummand;
                }
                PushDouble(fSum.get());
            }
        }
    }
}

/** 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;
    }

    KahanSum 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.get() );
}

/** 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();
--> --------------------

--> maximum size reached

--> --------------------